OmniSciDB  06b3bd477c
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
ExtensionFunctions.hpp
Go to the documentation of this file.
1 #include "../Shared/funcannotations.h"
2 #ifndef __CUDACC__
3 #include <cstdint>
4 #endif
5 #include <cmath>
6 #include <cstdlib>
7 
8 #include "OmniSciTypes.h"
9 
10 /* Example extension functions:
11  *
12  * EXTENSION_NOINLINE
13  * int32_t diff(const int32_t x, const int32_t y) {
14  * return x - y;
15  *
16  * Arrays map to a pair of pointer and element count. The pointer type must be
17  * consistent with the element type of the array column. Only array of numbers
18  * are supported for the moment. For example, an ARRAY_AT function which works
19  * for arrays of INTEGER and SMALLINT can be implemented as follows:
20  *
21  * EXTENSION_NOINLINE
22  * int64_t array_at(const int32_t* arr, const size_t elem_count, const size_t idx) {
23  * return idx < elem_count ? arr[idx] : -1;
24  * }
25  *
26  * EXTENSION_NOINLINE
27  * int64_t array_at__(const int16_t* arr, const size_t elem_count, const size_t idx) {
28  * return idx < elem_count ? arr[idx] : -1;
29  * }
30  *
31  * Note that the return type cannot be specialized and must be the same across
32  * all specializations. We can remove this constraint in the future if necessary.
33  */
34 
36 double Acos(const double x) {
37  return acos(x);
38 }
39 
41 double Asin(const double x) {
42  return asin(x);
43 }
44 
46 double Atan(const double x) {
47  return atan(x);
48 }
49 
51 double Atan2(const double y, const double x) {
52  return atan2(y, x);
53 }
54 
56 double Ceil(double x) {
57  return ceil(x);
58 }
59 
61 float Ceil__(float x) {
62  return ceil(x);
63 }
64 
66 int16_t Ceil__1(int16_t x) {
67  return x;
68 }
69 
71 int32_t Ceil__2(int32_t x) {
72  return x;
73 }
74 
76 int64_t Ceil__3(int64_t x) {
77  return x;
78 }
79 
81 double Cos(const double x) {
82  return cos(x);
83 }
84 
86 double Cot(const double x) {
87  return 1 / tan(x);
88 }
89 
91 double degrees(double x) {
92  return x * (180.0 / M_PI);
93 }
94 
96 double Exp(double x) {
97  return exp(x);
98 }
99 
101 double Floor(double x) {
102  return floor(x);
103 }
104 
106 float Floor__(float x) {
107  return floor(x);
108 }
109 
111 int16_t Floor__1(int16_t x) {
112  return x;
113 }
114 
116 int32_t Floor__2(int32_t x) {
117  return x;
118 }
119 
121 int64_t Floor__3(int64_t x) {
122  return x;
123 }
124 
126 double ln(const double x) {
127  return log(x);
128 }
129 
131 double ln__(const float x) {
132  return logf(x);
133 }
134 
136 double Log(const double x) {
137  return log(x);
138 }
139 
141 double Log__(const float x) {
142  return logf(x);
143 }
144 
146 double Log10(const double x) {
147  return log10(x);
148 }
149 
151 double Log10__(const float x) {
152  return log10f(x);
153 }
154 
156 double pi() {
157  return M_PI;
158 }
159 
161 double power(const double x, const double y) {
162  return pow(x, y);
163 }
164 
166 double radians(const double x) {
167  return x * (M_PI / 180.0);
168 }
169 
171 double Round(const double x, const int32_t y) {
172  if (y == 0) {
173  return round(x) + 0.0;
174  }
175 
176  double exp = pow(10, y);
177 #if defined(__powerpc__) && !defined(__CUDACC__)
178  int32_t yy = y - 1;
179  exp = 10 * powf((float)10L, yy);
180 #endif
181  return (round(x * exp) / exp) + 0.0;
182 }
183 
185 float Round__(const float x, const int32_t y) {
186  if (y == 0) {
187  return roundf(x) + 0.0f;
188  }
189 
190  float exp = powf((float)10L, y);
191 #if defined(__powerpc__) && !defined(__CUDACC__)
192  int32_t yy = y - 1;
193  exp = 10 * powf((float)10L, yy);
194 #endif
195  return roundf(x * exp) / exp + 0.0f;
196 }
197 
199 int16_t Round__1(const int16_t x, const int32_t y) {
200  if (y >= 0) {
201  return x;
202  }
203 
204  int32_t p = pow((float)10L, std::abs(y));
205  int32_t p_half = p >> 1;
206 
207  int64_t temp = x;
208 #if defined(__powerpc__) && !defined(__CUDACC__)
209  int16_t xx = x;
210  xx += 1;
211  temp = xx;
212  temp -= 1;
213 #endif
214  temp = temp >= 0 ? temp + p_half : temp - p_half;
215  temp = temp / p;
216  return temp * p;
217 }
218 
220 int32_t Round__2(const int32_t x, const int32_t y) {
221  if (y >= 0) {
222  return x;
223  }
224 
225  int32_t p = pow((float)10L, std::abs(y));
226  int32_t p_half = p >> 1;
227 
228  int64_t temp = x;
229 #if defined(__powerpc__) && !defined(__CUDACC__)
230  int32_t xx = x;
231  xx += 1;
232  temp = xx;
233  temp -= 1;
234 #endif
235  temp = temp >= 0 ? temp + p_half : temp - p_half;
236  temp = temp / p;
237  return temp * p;
238 }
239 
241 int64_t Round__3(const int64_t x, const int32_t y) {
242  if (y >= 0) {
243  return x;
244  }
245 
246  int64_t p = pow((double)10L, std::abs(y));
247  int64_t p_half = p >> 1;
248 
249  int64_t temp = x;
250  temp = temp >= 0 ? temp + p_half : temp - p_half;
251  temp = temp / p;
252  return temp * p;
253 }
254 
256 int64_t Round__4(const int64_t x, const int32_t y0, const int32_t scale) {
257  int32_t y = y0 - scale;
258 
259  if (y >= 0) {
260  return x;
261  }
262 
263  int64_t p = pow((double)10L, std::abs(y));
264  int64_t p_half = p >> 1;
265 
266  int64_t temp = x;
267  temp = temp >= 0 ? temp + p_half : temp - p_half;
268  temp = temp / p;
269  return temp * p;
270 }
271 
273 double Round2_to_digit(const double x, const int32_t y) {
274  double exp = pow(10, y);
275  return round(x * exp) / exp;
276 }
277 
279 double round_to_digit(const double x, const int32_t y) {
280  double exp = pow(10, y);
281  return round(x * exp) / exp;
282 }
283 
285 double Sin(const double x) {
286  return sin(x);
287 }
288 
290 double Tan(const double x) {
291  return tan(x);
292 }
293 
295 double Tan__(const float x) {
296  return tanf(x);
297 }
298 
300 double Truncate(const double x, const int32_t y) {
301  double p = pow((double)10L, y);
302  int64_t temp = x * p;
303  return temp / p;
304 }
305 
307 float Truncate__(const float x, const int32_t y) {
308  float p = powf((float)10L, y);
309  int64_t temp = x * p;
310  return temp / p;
311 }
312 
314 int16_t Truncate__1(const int16_t x, const int32_t y) {
315  if (y >= 0) {
316  return x;
317  }
318  int32_t p = pow((float)10L, std::abs(y));
319  int64_t temp = x / p;
320 #if defined(__powerpc__) && !defined(__CUDACC__)
321  int16_t xx = x;
322  xx += 1;
323  temp = xx;
324  temp -= 1;
325  temp /= p;
326 #endif
327  return temp * p;
328 }
329 
331 int32_t Truncate__2(const int32_t x, const int32_t y) {
332  if (y >= 0) {
333  return x;
334  }
335  int32_t p = pow((float)10L, std::abs(y));
336  int64_t temp = x / p;
337  return temp * p;
338 }
339 
341 int64_t Truncate__3(const int64_t x, const int32_t y) {
342  if (y >= 0) {
343  return x;
344  }
345  int64_t p = pow((double)10L, std::abs(y));
346  int64_t temp = x / p;
347  return temp * p;
348 }
349 
351 bool isNan(const double x) {
352  return std::isnan(x);
353 }
354 
356 bool isNan__(const float x) {
357  return std::isnan(x);
358 }
359 
361 double conv_4326_900913_x(const double x) {
362  return x * 111319.490778;
363 }
364 
366 double conv_4326_900913_y(const double y) {
367  return 6378136.99911 * log(tan(.00872664626 * y + .785398163397));
368 }
369 
392 double distance_in_meters(const double fromlon,
393  const double fromlat,
394  const double tolon,
395  const double tolat) {
396  double latitudeArc = (fromlat - tolat) * 0.017453292519943295769236907684886;
397  double longitudeArc = (fromlon - tolon) * 0.017453292519943295769236907684886;
398  double latitudeH = sin(latitudeArc * 0.5);
399  latitudeH *= latitudeH;
400  double lontitudeH = sin(longitudeArc * 0.5);
401  lontitudeH *= lontitudeH;
402  double tmp = cos(fromlat * 0.017453292519943295769236907684886) *
403  cos(tolat * 0.017453292519943295769236907684886);
404  return 6372797.560856 * (2.0 * asin(sqrt(latitudeH + tmp * lontitudeH)));
405 }
406 
408 double distance_in_meters__(const float fromlon,
409  const float fromlat,
410  const float tolon,
411  const float tolat) {
412  float latitudeArc = (fromlat - tolat) * 0.017453292519943295769236907684886;
413  float longitudeArc = (fromlon - tolon) * 0.017453292519943295769236907684886;
414  float latitudeH = sinf(latitudeArc * 0.5);
415  latitudeH *= latitudeH;
416  float lontitudeH = sinf(longitudeArc * 0.5);
417  lontitudeH *= lontitudeH;
418  float tmp = cosf(fromlat * 0.017453292519943295769236907684886) *
419  cosf(tolat * 0.017453292519943295769236907684886);
420  return 6372797.560856 * (2.0 * asinf(sqrtf(latitudeH + tmp * lontitudeH)));
421 }
422 
424 double approx_distance_in_meters(const float fromlon,
425  const float fromlat,
426  const float tolon,
427  const float tolat) {
428 #ifdef __CUDACC__
429  float latitudeArc = (fromlat - tolat) * 0.017453292519943295769236907684886;
430  float longitudeArc = (fromlon - tolon) * 0.017453292519943295769236907684886;
431  float latitudeH = __sinf(latitudeArc * 0.5);
432  latitudeH *= latitudeH;
433  float lontitudeH = __sinf(longitudeArc * 0.5);
434  lontitudeH *= lontitudeH;
435  float tmp = __cosf(fromlat * 0.017453292519943295769236907684886) *
436  __cosf(tolat * 0.017453292519943295769236907684886);
437  return 6372797.560856 * (2.0 * asinf(__fsqrt_rd(latitudeH + tmp * lontitudeH)));
438 #else
439  return distance_in_meters__(fromlon, fromlat, tolon, tolat);
440 #endif
441 }
442 
444 float rect_pixel_bin(const double val,
445  const double min,
446  const double max,
447  const int32_t numbins,
448  const int32_t dimensionsize) {
450  float numbinsf = float(numbins);
451  return float(int32_t(float((val - min) / (max - min)) * numbinsf)) *
452  float(dimensionsize) / numbinsf;
453 }
454 
456 float rect_pixel_bin_x(const double valx,
457  const double minx,
458  const double maxx,
459  const double rectwidth,
460  const double offsetx,
461  const int32_t imgwidth) {
462  const float imgwidthf = float(imgwidth);
463  const float rectwidthf = float(rectwidth);
464  double min = minx;
465  float offset = offsetx;
466  if (offset != 0) {
467  offset = fmodf(offset, rectwidthf);
468  if (offset > 0) {
469  offset -= rectwidthf;
470  }
471  min += offset * (maxx - minx) / imgwidthf;
472  }
473  return float(int32_t(float((valx - min) / (maxx - min)) * (imgwidthf - offset) /
474  rectwidthf)) *
475  rectwidthf +
476  offset + rectwidthf / 2.0f;
477 }
478 
480 float rect_pixel_bin_y(const double valy,
481  const double miny,
482  const double maxy,
483  const double rectheight,
484  const double offsety,
485  const int32_t imgheight) {
486  const float imgheightf = float(imgheight);
487  const float rectheightf = rectheight;
488  double min = miny;
489  float offset = offsety;
490  if (offset != 0) {
491  offset = fmodf(offset, rectheightf);
492  if (offset > 0) {
493  offset -= rectheightf;
494  }
495  min += offset * (maxy - miny) / imgheightf;
496  }
497  return float(int32_t(float((valy - min) / (maxy - min)) * (imgheightf - offset) /
498  rectheightf)) *
499  rectheightf +
500  offset + rectheightf / 2.0f;
501 }
502 
504 int32_t rect_pixel_bin_packed(const double valx,
505  const double minx,
506  const double maxx,
507  const double valy,
508  const double miny,
509  const double maxy,
510  const double rectwidth,
511  const double rectheight,
512  const double offsetx,
513  const double offsety,
514  const int32_t imgwidth,
515  const int32_t imgheight) {
516  const float imgwidthf = float(imgwidth);
517  const float rectwidthf = float(rectwidth);
518  double min = minx;
519  float offset = offsetx;
520  if (offset != 0) {
521  offset = fmodf(offset, rectwidthf);
522  if (offset > 0) {
523  offset -= rectwidthf;
524  }
525  min += offset * (maxx - minx) / imgwidthf;
526  }
527  float rx = float(int32_t(float((valx - min) / (maxx - min)) * (imgwidthf - offset) /
528  rectwidthf)) *
529  rectwidthf +
530  offset + rectwidthf / 2.0f;
531 
532  const float imgheightf = float(imgheight);
533  const float rectheightf = rectheight;
534  min = miny;
535  offset = offsety;
536  if (offset != 0) {
537  offset = fmodf(offset, rectheightf);
538  if (offset > 0) {
539  offset -= rectheightf;
540  }
541  min += offset * (maxy - miny) / imgheightf;
542  }
543  float ry = float(int32_t(float((valy - min) / (maxy - min)) * (imgheightf - offset) /
544  rectheightf)) *
545  rectheightf +
546  offset + rectheightf / 2.0f;
547 
548  // and pack as two 14.2 fixed-point values into 32bits
549  int32_t ux = static_cast<int32_t>(rx * 4.0f);
550  int32_t uy = static_cast<int32_t>(ry * 4.0f);
551  return (ux & 0x7FFF) | ((uy & 0x7FFF) << 16);
552 }
553 
555 float reg_hex_horiz_pixel_bin_x(const double valx,
556  const double minx,
557  const double maxx,
558  const double valy,
559  const double miny,
560  const double maxy,
561  const double hexwidth,
562  const double hexheight,
563  const double offsetx,
564  const double offsety,
565  const int32_t imgwidth,
566  const int32_t imgheight) {
567  const float sqrt3 = 1.7320508075688772;
568  const float imgwidthf = float(imgwidth);
569  const float imgheightf = float(imgheight);
570  const float hexwidthf = float(hexwidth);
571  const float hexheightf = float(hexheight);
572 
573  // expand the bounds of the data according
574  // to the input offsets. This is done because
575  // we also expand the image size according to the
576  // offsets because this algorithm layers the hexagon
577  // bins starting at the bottom left corner
578  double xmin = minx;
579  float xoffset = offsetx;
580  if (xoffset != 0) {
581  xoffset = fmodf(xoffset, hexwidthf);
582  if (xoffset > 0) {
583  xoffset -= hexwidthf;
584  }
585  xmin += xoffset * (maxx - xmin) / imgwidthf;
586  }
587 
588  double ymin = miny;
589  float yoffset = offsety;
590  if (yoffset != 0) {
591  yoffset = fmodf(yoffset, 1.5f * hexheightf);
592  if (yoffset > 0) {
593  yoffset -= 1.5f * hexheightf;
594  }
595  ymin += yoffset * (maxy - ymin) / imgheightf;
596  }
597 
598  // get the pixel position of the point
599  // assumes a linear scale here
600  // Rounds to the nearest pixel.
601  const float pix_x =
602  roundf((imgwidthf - xoffset) * float((valx - xmin) / (maxx - xmin)));
603  const float pix_y =
604  roundf((imgheightf - yoffset) * float((valy - ymin) / (maxy - ymin)));
605 
606  // Now convert the pixel position into a
607  // cube-coordinate system representation
608  const float hexsize = hexheightf / 2.0f;
609  const float cube_x = ((pix_x / sqrt3) - (pix_y / 3.0f)) / hexsize;
610  const float cube_z = (pix_y * 2.0f / 3.0f) / hexsize;
611  const float cube_y = -cube_x - cube_z;
612 
613  // need to round the cube coordinates above
614  float rx = round(cube_x);
615  float ry = round(cube_y);
616  float rz = round(cube_z);
617  const float x_diff = fabs(rx - cube_x);
618  const float y_diff = fabs(ry - cube_y);
619  const float z_diff = fabs(rz - cube_z);
620  if (x_diff > y_diff && x_diff > z_diff) {
621  rx = -ry - rz;
622  } else if (y_diff > z_diff) {
623  ry = -rx - rz;
624  } else {
625  rz = -rx - ry;
626  }
627 
628  // now convert the cube/hex coord to a pixel location
629  return hexsize * sqrt3 * (rx + rz / 2.0f) + xoffset;
630 }
631 
633 float reg_hex_horiz_pixel_bin_y(const double valx,
634  const double minx,
635  const double maxx,
636  const double valy,
637  const double miny,
638  const double maxy,
639  const double hexwidth,
640  const double hexheight,
641  const double offsetx,
642  const double offsety,
643  const int32_t imgwidth,
644  const int32_t imgheight) {
645  const float sqrt3 = 1.7320508075688772;
646  const float imgwidthf = float(imgwidth);
647  const float imgheightf = float(imgheight);
648  const float hexwidthf = float(hexwidth);
649  const float hexheightf = float(hexheight);
650 
651  // expand the bounds of the data according
652  // to the input offsets. This is done because
653  // we also expand the image size according to the
654  // offsets because this algorithm layers the hexagon
655  // bins starting at the bottom left corner
656  double xmin = minx;
657  float xoffset = offsetx;
658  if (xoffset != 0) {
659  xoffset = fmodf(xoffset, hexwidthf);
660  if (xoffset > 0) {
661  xoffset -= hexwidthf;
662  }
663  xmin += xoffset * (maxx - xmin) / imgwidthf;
664  }
665 
666  double ymin = miny;
667  float yoffset = offsety;
668  if (yoffset != 0) {
669  yoffset = fmodf(yoffset, 1.5f * hexheightf);
670  if (yoffset > 0) {
671  yoffset -= 1.5f * hexheightf;
672  }
673  ymin += yoffset * (maxy - ymin) / imgheightf;
674  }
675 
676  // get the pixel position of the point
677  // assumes a linear scale here
678  // Rounds to the nearest pixel.
679  const float pix_x =
680  roundf((imgwidthf - xoffset) * float((valx - xmin) / (maxx - xmin)));
681  const float pix_y =
682  roundf((imgheightf - yoffset) * float((valy - ymin) / (maxy - ymin)));
683 
684  // Now convert the pixel position into a
685  // cube-coordinate system representation
686  const float hexsize = hexheightf / 2.0f;
687  const float cube_x = ((pix_x / sqrt3) - (pix_y / 3.0f)) / hexsize;
688  const float cube_z = (pix_y * 2.0f / 3.0f) / hexsize;
689  const float cube_y = -cube_x - cube_z;
690 
691  // need to round the cube coordinates above
692  float rx = round(cube_x);
693  float ry = round(cube_y);
694  float rz = round(cube_z);
695  const float x_diff = fabs(rx - cube_x);
696  const float y_diff = fabs(ry - cube_y);
697  const float z_diff = fabs(rz - cube_z);
698  if ((x_diff <= y_diff || x_diff <= z_diff) && y_diff <= z_diff) {
699  rz = -rx - ry;
700  }
701 
702  // now convert the cube/hex coord to a pixel location
703  return hexsize * 3.0f / 2.0f * rz + yoffset;
704 }
705 
707 int32_t reg_hex_horiz_pixel_bin_packed(const double valx,
708  const double minx,
709  const double maxx,
710  const double valy,
711  const double miny,
712  const double maxy,
713  const double hexwidth,
714  const double hexheight,
715  const double offsetx,
716  const double offsety,
717  const int32_t imgwidth,
718  const int32_t imgheight) {
719  const float sqrt3 = 1.7320508075688772;
720  const float imgwidthf = float(imgwidth);
721  const float imgheightf = float(imgheight);
722  const float hexwidthf = float(hexwidth);
723  const float hexheightf = float(hexheight);
724 
725  // expand the bounds of the data according
726  // to the input offsets. This is done because
727  // we also expand the image size according to the
728  // offsets because this algorithm layers the hexagon
729  // bins starting at the bottom left corner
730  double xmin = minx;
731  float xoffset = offsetx;
732  if (xoffset != 0) {
733  xoffset = fmodf(xoffset, hexwidthf);
734  if (xoffset > 0) {
735  xoffset -= hexwidthf;
736  }
737  xmin += xoffset * (maxx - xmin) / imgwidthf;
738  }
739 
740  double ymin = miny;
741  float yoffset = offsety;
742  if (yoffset != 0) {
743  yoffset = fmodf(yoffset, 1.5f * hexheightf);
744  if (yoffset > 0) {
745  yoffset -= 1.5f * hexheightf;
746  }
747  ymin += yoffset * (maxy - ymin) / imgheightf;
748  }
749 
750  // get the pixel position of the point
751  // assumes a linear scale here
752  // Rounds to the nearest pixel.
753  const float pix_x =
754  roundf((imgwidthf - xoffset) * float((valx - xmin) / (maxx - xmin)));
755  const float pix_y =
756  roundf((imgheightf - yoffset) * float((valy - ymin) / (maxy - ymin)));
757 
758  // Now convert the pixel position into a
759  // cube-coordinate system representation
760  const float hexsize = hexheightf / 2.0f;
761  const float cube_x = ((pix_x / sqrt3) - (pix_y / 3.0f)) / hexsize;
762  const float cube_z = (pix_y * 2.0f / 3.0f) / hexsize;
763  const float cube_y = -cube_x - cube_z;
764 
765  // need to round the cube coordinates above
766  float rx = round(cube_x);
767  float ry = round(cube_y);
768  float rz = round(cube_z);
769  const float x_diff = fabs(rx - cube_x);
770  const float y_diff = fabs(ry - cube_y);
771  const float z_diff = fabs(rz - cube_z);
772  if (x_diff > y_diff && x_diff > z_diff) {
773  rx = -ry - rz;
774  } else if (y_diff <= z_diff) {
775  rz = -rx - ry;
776  }
777 
778  // now convert the cube/hex coord to pixel locations
779  float hx = hexsize * sqrt3 * (rx + rz / 2.0f) + xoffset;
780  float hy = hexsize * 3.0f / 2.0f * rz + yoffset;
781 
782  // and pack as two 14.2 fixed-point values into 32bits
783  int32_t ux = static_cast<int32_t>(hx * 4.0f);
784  int32_t uy = static_cast<int32_t>(hy * 4.0f);
785  return (ux & 0x7FFF) | ((uy & 0x7FFF) << 16);
786 }
787 
789 float reg_hex_vert_pixel_bin_x(const double valx,
790  const double minx,
791  const double maxx,
792  const double valy,
793  const double miny,
794  const double maxy,
795  const double hexwidth,
796  const double hexheight,
797  const double offsetx,
798  const double offsety,
799  const int32_t imgwidth,
800  const int32_t imgheight) {
801  const float sqrt3 = 1.7320508075688772;
802  const float imgwidthf = float(imgwidth);
803  const float imgheightf = float(imgheight);
804  const float hexwidthf = float(hexwidth);
805  const float hexheightf = float(hexheight);
806 
807  // expand the bounds of the data according
808  // to the input offsets. This is done because
809  // we also expand the image size according to the
810  // offsets because this algorithm layers the hexagon
811  // bins starting at the bottom left corner
812  double xmin = minx;
813  float xoffset = offsetx;
814  if (xoffset != 0) {
815  xoffset = fmodf(xoffset, 1.5f * hexwidthf);
816  if (xoffset > 0) {
817  xoffset -= 1.5f * hexwidthf;
818  }
819  xmin += xoffset * (maxx - xmin) / imgwidthf;
820  }
821 
822  double ymin = miny;
823  float yoffset = offsety;
824  if (yoffset != 0) {
825  yoffset = fmodf(yoffset, hexheightf);
826  if (yoffset > 0) {
827  yoffset -= hexheightf;
828  }
829  ymin += yoffset * (maxy - ymin) / imgheightf;
830  }
831 
832  // get the pixel position of the point
833  // assumes a linear scale here
834  // Rounds to the nearest pixel.
835  const float pix_x =
836  roundf((imgwidthf - xoffset) * float((valx - xmin) / (maxx - xmin)));
837  const float pix_y =
838  roundf((imgheightf - yoffset) * float((valy - ymin) / (maxy - ymin)));
839 
840  // Now convert the pixel position into a
841  // cube-coordinate system representation
842  const float hexsize = hexwidthf / 2.0f;
843  const float cube_x = (pix_x * 2.0f / 3.0f) / hexsize;
844  const float cube_z = ((pix_y / sqrt3) - (pix_x / 3.0f)) / hexsize;
845  const float cube_y = -cube_x - cube_z;
846 
847  // need to round the cube coordinates above
848  float rx = round(cube_x);
849  float ry = round(cube_y);
850  float rz = round(cube_z);
851  const float x_diff = fabs(rx - cube_x);
852  const float y_diff = fabs(ry - cube_y);
853  const float z_diff = fabs(rz - cube_z);
854  if (x_diff > y_diff && x_diff > z_diff) {
855  rx = -ry - rz;
856  }
857 
858  // now convert the cube/hex coord to a pixel location
859  return hexsize * 3.0f / 2.0f * rx + xoffset;
860 }
861 
863 float reg_hex_vert_pixel_bin_y(const double valx,
864  const double minx,
865  const double maxx,
866  const double valy,
867  const double miny,
868  const double maxy,
869  const double hexwidth,
870  const double hexheight,
871  const double offsetx,
872  const double offsety,
873  const int32_t imgwidth,
874  const int32_t imgheight) {
875  const float sqrt3 = 1.7320508075688772;
876  const float imgwidthf = float(imgwidth);
877  const float imgheightf = float(imgheight);
878  const float hexwidthf = float(hexwidth);
879  const float hexheightf = float(hexheight);
880 
881  // expand the bounds of the data according
882  // to the input offsets. This is done because
883  // we also expand the image size according to the
884  // offsets because this algorithm layers the hexagon
885  // bins starting at the bottom left corner
886  float xmin = minx;
887  float xoffset = offsetx;
888  if (xoffset != 0) {
889  xoffset = fmodf(xoffset, 1.5f * hexwidthf);
890  if (xoffset > 0) {
891  xoffset -= 1.5f * hexwidthf;
892  }
893  xmin += xoffset * (maxx - xmin) / imgwidthf;
894  }
895 
896  float ymin = miny;
897  float yoffset = offsety;
898  if (yoffset != 0) {
899  yoffset = fmodf(yoffset, hexheightf);
900  if (yoffset > 0) {
901  yoffset -= hexheightf;
902  }
903  ymin += yoffset * (maxy - ymin) / imgheightf;
904  }
905 
906  // get the pixel position of the point
907  // assumes a linear scale here
908  // Rounds to the nearest pixel.
909  const float pix_x = roundf((imgwidthf - xoffset) * (valx - xmin) / (maxx - xmin));
910  const float pix_y = roundf((imgheightf - yoffset) * (valy - ymin) / (maxy - ymin));
911 
912  // Now convert the pixel position into a
913  // cube-coordinate system representation
914  const float hexsize = hexwidthf / 2.0f;
915  const float cube_x = (pix_x * 2.0f / 3.0f) / hexsize;
916  const float cube_z = ((pix_y / sqrt3) - (pix_x / 3.0f)) / hexsize;
917  const float cube_y = -cube_x - cube_z;
918 
919  // need to round the cube coordinates above
920  float rx = round(cube_x);
921  float ry = round(cube_y);
922  float rz = round(cube_z);
923  const float x_diff = fabs(rx - cube_x);
924  const float y_diff = fabs(ry - cube_y);
925  const float z_diff = fabs(rz - cube_z);
926  if (x_diff > y_diff && x_diff > z_diff) {
927  rx = -ry - rz;
928  } else if (y_diff > z_diff) {
929  ry = -rx - rz;
930  } else {
931  rz = -rx - ry;
932  }
933 
934  // now convert the cube/hex coord to a pixel location
935  return hexsize * sqrt3 * (rz + rx / 2.0f) + yoffset;
936 }
937 
939 int32_t reg_hex_vert_pixel_bin_packed(const double valx,
940  const double minx,
941  const double maxx,
942  const double valy,
943  const double miny,
944  const double maxy,
945  const double hexwidth,
946  const double hexheight,
947  const double offsetx,
948  const double offsety,
949  const int32_t imgwidth,
950  const int32_t imgheight) {
951  const float sqrt3 = 1.7320508075688772;
952  const float imgwidthf = float(imgwidth);
953  const float imgheightf = float(imgheight);
954  const float hexwidthf = float(hexwidth);
955  const float hexheightf = float(hexheight);
956 
957  // expand the bounds of the data according
958  // to the input offsets. This is done because
959  // we also expand the image size according to the
960  // offsets because this algorithm layers the hexagon
961  // bins starting at the bottom left corner
962  double xmin = minx;
963  float xoffset = offsetx;
964  if (xoffset != 0) {
965  xoffset = fmodf(xoffset, 1.5f * hexwidthf);
966  if (xoffset > 0) {
967  xoffset -= 1.5f * hexwidthf;
968  }
969  xmin += xoffset * (maxx - xmin) / imgwidthf;
970  }
971 
972  double ymin = miny;
973  float yoffset = offsety;
974  if (yoffset != 0) {
975  yoffset = fmodf(yoffset, hexheightf);
976  if (yoffset > 0) {
977  yoffset -= hexheightf;
978  }
979  ymin += yoffset * (maxy - ymin) / imgheightf;
980  }
981 
982  // get the pixel position of the point
983  // assumes a linear scale here
984  // Rounds to the nearest pixel.
985  const float pix_x = roundf((imgwidthf - xoffset) * float(valx - xmin) / (maxx - xmin));
986  const float pix_y = roundf((imgheightf - yoffset) * float(valy - ymin) / (maxy - ymin));
987 
988  // Now convert the pixel position into a
989  // cube-coordinate system representation
990  const float hexsize = hexwidthf / 2.0f;
991  const float cube_x = (pix_x * 2.0f / 3.0f) / hexsize;
992  const float cube_z = ((pix_y / sqrt3) - (pix_x / 3.0f)) / hexsize;
993  const float cube_y = -cube_x - cube_z;
994 
995  // need to round the cube coordinates above
996  float rx = round(cube_x);
997  float ry = round(cube_y);
998  float rz = round(cube_z);
999  const float x_diff = fabs(rx - cube_x);
1000  const float y_diff = fabs(ry - cube_y);
1001  const float z_diff = fabs(rz - cube_z);
1002  if (x_diff > y_diff && x_diff > z_diff) {
1003  rx = -ry - rz;
1004  } else if (y_diff <= z_diff) {
1005  rz = -rx - ry;
1006  }
1007 
1008  // now convert the cube/hex coord to a pixel location
1009  float hx = hexsize * 3.0f / 2.0f * rx + xoffset;
1010  float hy = hexsize * sqrt3 * (rz + rx / 2.0f) + yoffset;
1011 
1012  // and pack as two 14.2 fixed-point values into 32bits
1013  int32_t ux = static_cast<int32_t>(hx * 4.0f);
1014  int32_t uy = static_cast<int32_t>(hy * 4.0f);
1015  return (ux & 0x7FFF) | ((uy & 0x7FFF) << 16);
1016 }
1017 
1019 double convert_meters_to_merc_pixel_width(const double meters,
1020  const double lon,
1021  const double lat,
1022  const double min_lon,
1023  const double max_lon,
1024  const int32_t img_width,
1025  const double min_width) {
1026  const double const1 = 0.017453292519943295769236907684886;
1027  const double const2 = 6372797.560856;
1028  double t1 = sinf(meters / (2.0 * const2));
1029  double t2 = cosf(const1 * lat);
1030  const double newlon = lon - (2.0 * asinf(t1 / t2)) / const1;
1031  t1 = conv_4326_900913_x(lon);
1032  t2 = conv_4326_900913_x(newlon);
1033  const double min_merc_x = conv_4326_900913_x(min_lon);
1034  const double max_merc_x = conv_4326_900913_x(max_lon);
1035  const double merc_diff = max_merc_x - min_merc_x;
1036  t1 = ((t1 - min_merc_x) / merc_diff) * static_cast<double>(img_width);
1037  t2 = ((t2 - min_merc_x) / merc_diff) * static_cast<double>(img_width);
1038 
1039  // TODO(croot): need to account for edge cases, such as getting close to the poles.
1040  const double sz = fabs(t1 - t2);
1041  return (sz < min_width ? min_width : sz);
1042 }
1043 
1045 double convert_meters_to_merc_pixel_height(const double meters,
1046  const double lon,
1047  const double lat,
1048  const double min_lat,
1049  const double max_lat,
1050  const int32_t img_height,
1051  const double min_height) {
1052  const double const1 = 0.017453292519943295769236907684886;
1053  const double const2 = 6372797.560856;
1054  const double latdiff = meters / (const1 * const2);
1055  const double newlat =
1056  (lat < 0) ? lat + latdiff : lat - latdiff; // assumes a lat range of [-90, 90]
1057  double t1 = conv_4326_900913_y(lat);
1058  double t2 = conv_4326_900913_y(newlat);
1059  const double min_merc_y = conv_4326_900913_y(min_lat);
1060  const double max_merc_y = conv_4326_900913_y(max_lat);
1061  const double merc_diff = max_merc_y - min_merc_y;
1062  t1 = ((t1 - min_merc_y) / merc_diff) * static_cast<double>(img_height);
1063  t2 = ((t2 - min_merc_y) / merc_diff) * static_cast<double>(img_height);
1064 
1065  // TODO(croot): need to account for edge cases, such as getting close to the poles.
1066  const double sz = fabs(t1 - t2);
1067  return (sz < min_height ? min_height : sz);
1068 }
1069 
1071  const double lat,
1072  const double min_lon,
1073  const double max_lon,
1074  const double min_lat,
1075  const double max_lat) {
1076  return !(lon < min_lon || lon > max_lon || lat < min_lat || lat > max_lat);
1077 }
1078 
1080  const double lat,
1081  const double meters,
1082  const double min_lon,
1083  const double max_lon,
1084  const double min_lat,
1085  const double max_lat) {
1086  const double const1 = 0.017453292519943295769236907684886;
1087  const double const2 = 6372797.560856;
1088  const double latdiff = meters / (const1 * const2);
1089  const double t1 = sinf(meters / (2.0 * const2));
1090  const double t2 = cosf(const1 * lat);
1091  const double londiff = (2.0 * asinf(t1 / t2)) / const1;
1092  return !(lon + londiff < min_lon || lon - londiff > max_lon ||
1093  lat + latdiff < min_lat || lat - latdiff > max_lat);
1094 }
1095 
1096 #include "ExtensionFunctionsArray.hpp"
1097 #include "ExtensionFunctionsGeo.hpp"
EXTENSION_NOINLINE double ln(const double x)
EXTENSION_NOINLINE int32_t reg_hex_horiz_pixel_bin_packed(const double valx, const double minx, const double maxx, const double valy, const double miny, const double maxy, const double hexwidth, const double hexheight, const double offsetx, const double offsety, const int32_t imgwidth, const int32_t imgheight)
EXTENSION_NOINLINE float reg_hex_horiz_pixel_bin_x(const double valx, const double minx, const double maxx, const double valy, const double miny, const double maxy, const double hexwidth, const double hexheight, const double offsetx, const double offsety, const int32_t imgwidth, const int32_t imgheight)
EXTENSION_NOINLINE double Tan__(const float x)
EXTENSION_INLINE bool is_point_in_merc_view(const double lon, const double lat, const double min_lon, const double max_lon, const double min_lat, const double max_lat)
EXTENSION_NOINLINE double Sin(const double x)
EXTENSION_NOINLINE int32_t Truncate__2(const int32_t x, const int32_t y)
EXTENSION_NOINLINE double power(const double x, const double y)
EXTENSION_NOINLINE double Cot(const double x)
EXTENSION_NOINLINE int32_t rect_pixel_bin_packed(const double valx, const double minx, const double maxx, const double valy, const double miny, const double maxy, const double rectwidth, const double rectheight, const double offsetx, const double offsety, const int32_t imgwidth, const int32_t imgheight)
EXTENSION_NOINLINE int64_t Round__3(const int64_t x, const int32_t y)
EXTENSION_NOINLINE double distance_in_meters__(const float fromlon, const float fromlat, const float tolon, const float tolat)
EXTENSION_NOINLINE int16_t Round__1(const int16_t x, const int32_t y)
EXTENSION_NOINLINE double Asin(const double x)
EXTENSION_NOINLINE double Log10(const double x)
EXTENSION_NOINLINE double convert_meters_to_merc_pixel_width(const double meters, const double lon, const double lat, const double min_lon, const double max_lon, const int32_t img_width, const double min_width)
EXTENSION_NOINLINE float Floor__(float x)
EXTENSION_NOINLINE int64_t Floor__3(int64_t x)
EXTENSION_NOINLINE double Ceil(double x)
#define EXTENSION_NOINLINE
Definition: OmniSciTypes.h:25
EXTENSION_NOINLINE int16_t Ceil__1(int16_t x)
EXTENSION_NOINLINE int64_t Truncate__3(const int64_t x, const int32_t y)
EXTENSION_NOINLINE int16_t Floor__1(int16_t x)
EXTENSION_NOINLINE double Exp(double x)
EXTENSION_NOINLINE bool isNan__(const float x)
EXTENSION_NOINLINE double radians(const double x)
EXTENSION_NOINLINE double Cos(const double x)
EXTENSION_NOINLINE int64_t Round__4(const int64_t x, const int32_t y0, const int32_t scale)
EXTENSION_NOINLINE int32_t Ceil__2(int32_t x)
EXTENSION_NOINLINE float rect_pixel_bin_y(const double valy, const double miny, const double maxy, const double rectheight, const double offsety, const int32_t imgheight)
EXTENSION_NOINLINE int32_t Round__2(const int32_t x, const int32_t y)
EXTENSION_NOINLINE bool isNan(const double x)
EXTENSION_NOINLINE double Floor(double x)
EXTENSION_NOINLINE double Round2_to_digit(const double x, const int32_t y)
EXTENSION_NOINLINE double degrees(double x)
EXTENSION_NOINLINE double round_to_digit(const double x, const int32_t y)
EXTENSION_NOINLINE double convert_meters_to_merc_pixel_height(const double meters, const double lon, const double lat, const double min_lat, const double max_lat, const int32_t img_height, const double min_height)
EXTENSION_NOINLINE double Log__(const float x)
EXTENSION_NOINLINE float reg_hex_horiz_pixel_bin_y(const double valx, const double minx, const double maxx, const double valy, const double miny, const double maxy, const double hexwidth, const double hexheight, const double offsetx, const double offsety, const int32_t imgwidth, const int32_t imgheight)
EXTENSION_NOINLINE double ln__(const float x)
EXTENSION_NOINLINE float rect_pixel_bin(const double val, const double min, const double max, const int32_t numbins, const int32_t dimensionsize)
EXTENSION_NOINLINE float reg_hex_vert_pixel_bin_x(const double valx, const double minx, const double maxx, const double valy, const double miny, const double maxy, const double hexwidth, const double hexheight, const double offsetx, const double offsety, const int32_t imgwidth, const int32_t imgheight)
EXTENSION_NOINLINE int64_t Ceil__3(int64_t x)
EXTENSION_NOINLINE float Round__(const float x, const int32_t y)
EXTENSION_NOINLINE float Truncate__(const float x, const int32_t y)
EXTENSION_NOINLINE double distance_in_meters(const double fromlon, const double fromlat, const double tolon, const double tolat)
Computes the distance, in meters, between two WGS-84 positions.
EXTENSION_NOINLINE double Round(const double x, const int32_t y)
EXTENSION_NOINLINE float Ceil__(float x)
EXTENSION_NOINLINE float reg_hex_vert_pixel_bin_y(const double valx, const double minx, const double maxx, const double valy, const double miny, const double maxy, const double hexwidth, const double hexheight, const double offsetx, const double offsety, const int32_t imgwidth, const int32_t imgheight)
EXTENSION_NOINLINE int32_t Floor__2(int32_t x)
EXTENSION_NOINLINE double Truncate(const double x, const int32_t y)
EXTENSION_NOINLINE double Log(const double x)
#define EXTENSION_INLINE
Definition: OmniSciTypes.h:24
EXTENSION_NOINLINE double approx_distance_in_meters(const float fromlon, const float fromlat, const float tolon, const float tolat)
EXTENSION_NOINLINE float rect_pixel_bin_x(const double valx, const double minx, const double maxx, const double rectwidth, const double offsetx, const int32_t imgwidth)
EXTENSION_NOINLINE double pi()
EXTENSION_NOINLINE int32_t reg_hex_vert_pixel_bin_packed(const double valx, const double minx, const double maxx, const double valy, const double miny, const double maxy, const double hexwidth, const double hexheight, const double offsetx, const double offsety, const int32_t imgwidth, const int32_t imgheight)
EXTENSION_NOINLINE double Tan(const double x)
EXTENSION_NOINLINE double Atan2(const double y, const double x)
EXTENSION_NOINLINE double conv_4326_900913_y(const double y)
EXTENSION_NOINLINE double Atan(const double x)
EXTENSION_NOINLINE double conv_4326_900913_x(const double x)
EXTENSION_NOINLINE double Log10__(const float x)
EXTENSION_NOINLINE bool is_point_size_in_merc_view(const double lon, const double lat, const double meters, const double min_lon, const double max_lon, const double min_lat, const double max_lat)
EXTENSION_NOINLINE int16_t Truncate__1(const int16_t x, const int32_t y)
EXTENSION_NOINLINE double Acos(const double x)