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