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