OmniSciDB  471d68cefb
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
ExtensionFunctionsGeo.hpp
Go to the documentation of this file.
2 
3 #ifndef __CUDACC__
4 #include "Shared/likely.h"
5 #else
6 #define LIKELY
7 #define UNLIKELY
8 #endif
9 
10 // #define DEBUG_STMT(x) x;
11 #define DEBUG_STMT(x)
12 
13 // Adjustable tolerance, determined by compression mode.
14 // The criteria is to still recognize a compressed+decompressed number.
15 // For example 1.0 longitude compressed with GEOINT32 and then decompressed
16 // comes back as 0.99999994086101651, which is still within GEOINT32
17 // tolerance val 0.0000001
18 DEVICE ALWAYS_INLINE double tol(int32_t ic) {
19  if (ic == COMPRESSION_GEOINT32) {
20  return TOLERANCE_GEOINT32;
21  }
22  return TOLERANCE_DEFAULT;
23 }
24 
25 // Combine tolerances for two-component calculations
26 DEVICE ALWAYS_INLINE double tol(const int32_t ic1, const int32_t ic2) {
27  return fmax(tol(ic1), tol(ic2));
28 }
29 
30 DEVICE ALWAYS_INLINE bool tol_zero(const double x,
31  const double tolerance = TOLERANCE_DEFAULT) {
32  return (-tolerance <= x) && (x <= tolerance);
33 }
34 
35 DEVICE ALWAYS_INLINE bool tol_eq(const double x,
36  const double y,
37  const double tolerance = TOLERANCE_DEFAULT) {
38  auto diff = x - y;
39  return (-tolerance <= diff) && (diff <= tolerance);
40 }
41 
42 DEVICE ALWAYS_INLINE bool tol_le(const double x,
43  const double y,
44  const double tolerance = TOLERANCE_DEFAULT) {
45  return x <= (y + tolerance);
46 }
47 
48 DEVICE ALWAYS_INLINE bool tol_ge(const double x,
49  const double y,
50  const double tolerance = TOLERANCE_DEFAULT) {
51  return (x + tolerance) >= y;
52 }
53 
54 DEVICE ALWAYS_INLINE double decompress_coord(const int8_t* data,
55  const int32_t index,
56  const int32_t ic,
57  const bool x) {
58  if (ic == COMPRESSION_GEOINT32) {
59  auto compressed_coords = reinterpret_cast<const int32_t*>(data);
60  auto compressed_coord = compressed_coords[index];
61  if (x) {
63  } else {
65  }
66  }
67  auto double_coords = reinterpret_cast<const double*>(data);
68  return double_coords[index];
69 }
70 
71 DEVICE ALWAYS_INLINE int32_t compression_unit_size(const int32_t ic) {
72  if (ic == COMPRESSION_GEOINT32) {
73  return 4;
74  }
75  return 8;
76 }
77 
78 DEVICE ALWAYS_INLINE double transform_coord(const double coord,
79  const int32_t isr,
80  const int32_t osr,
81  const bool x) {
82  if (isr == 4326) {
83  if (osr == 900913) {
84  // WGS 84 --> Web Mercator
85  if (x) {
86  return conv_4326_900913_x(coord);
87  } else {
88  return conv_4326_900913_y(coord);
89  }
90  }
91  }
92  return coord;
93 }
94 
95 DEVICE ALWAYS_INLINE double transform_coord_x(const double coord,
96  const int32_t isr,
97  const int32_t osr) {
98  if (isr == 4326) {
99  if (osr == 900913) {
100  // WGS 84 --> Web Mercator
101  return conv_4326_900913_x(coord);
102  }
103  }
104  return coord;
105 }
106 
107 DEVICE ALWAYS_INLINE double transform_coord_y(const double coord,
108  const int32_t isr,
109  const int32_t osr) {
110  if (isr == 4326) {
111  if (osr == 900913) {
112  // WGS 84 --> Web Mercator
113  return conv_4326_900913_y(coord);
114  }
115  }
116  return coord;
117 }
118 
119 // X coord accessor handling on-the-fly decommpression and transforms
120 DEVICE ALWAYS_INLINE double coord_x(const int8_t* data,
121  const int32_t index,
122  const int32_t ic,
123  const int32_t isr,
124  const int32_t osr) {
125  auto decompressed_coord_x = decompress_coord(data, index, ic, true);
126  auto decompressed_transformed_coord_x =
127  transform_coord(decompressed_coord_x, isr, osr, true);
128  return decompressed_transformed_coord_x;
129 }
130 
131 // Y coord accessor handling on-the-fly decommpression and transforms
132 DEVICE ALWAYS_INLINE double coord_y(const int8_t* data,
133  const int32_t index,
134  const int32_t ic,
135  const int32_t isr,
136  const int32_t osr) {
137  auto decompressed_coord_y = decompress_coord(data, index, ic, false);
138  auto decompressed_transformed_coord_y =
139  transform_coord(decompressed_coord_y, isr, osr, false);
140  return decompressed_transformed_coord_y;
141 }
142 
143 // coord accessor for compressed coords at location index
144 DEVICE ALWAYS_INLINE int32_t compressed_coord(const int8_t* data, const int32_t index) {
145  const auto compressed_coords = reinterpret_cast<const int32_t*>(data);
146  return compressed_coords[index];
147 }
148 
149 // Cartesian distance between points, squared
151  double p1y,
152  double p2x,
153  double p2y) {
154  auto x = p1x - p2x;
155  auto y = p1y - p2y;
156  auto x2 = x * x;
157  auto y2 = y * y;
158  auto d2 = x2 + y2;
160  return 0.0;
161  }
162  return d2;
163 }
164 
165 // Cartesian distance between points
167  double p1y,
168  double p2x,
169  double p2y) {
170  auto d2 = distance_point_point_squared(p1x, p1y, p2x, p2y);
171  return sqrt(d2);
172 }
173 
174 // Cartesian distance between a point and a line segment
175 DEVICE
176 double distance_point_line(double px,
177  double py,
178  double l1x,
179  double l1y,
180  double l2x,
181  double l2y) {
182  double length = distance_point_point(l1x, l1y, l2x, l2y);
183  if (tol_zero(length)) {
184  return distance_point_point(px, py, l1x, l1y);
185  }
186 
187  // Find projection of point P onto the line segment AB:
188  // Line containing that segment: A + k * (B - A)
189  // Projection of point P onto the line touches it at
190  // k = dot(P-A,B-A) / length^2
191  // AB segment is represented by k = [0,1]
192  // Clamping k to [0,1] will give the shortest distance from P to AB segment
193  double dotprod = (px - l1x) * (l2x - l1x) + (py - l1y) * (l2y - l1y);
194  double k = dotprod / (length * length);
195  k = fmax(0.0, fmin(1.0, k));
196  double projx = l1x + k * (l2x - l1x);
197  double projy = l1y + k * (l2y - l1y);
198  return distance_point_point(px, py, projx, projy);
199 }
200 
201 // Cartesian distance between a point and a line segment
202 DEVICE
204  double py,
205  double l1x,
206  double l1y,
207  double l2x,
208  double l2y) {
209  double length = distance_point_point_squared(l1x, l1y, l2x, l2y);
210  if (tol_zero(length, TOLERANCE_DEFAULT_SQUARED)) {
211  return distance_point_point_squared(px, py, l1x, l1y);
212  }
213 
214  // Find projection of point P onto the line segment AB:
215  // Line containing that segment: A + k * (B - A)
216  // Projection of point P onto the line touches it at
217  // k = dot(P-A,B-A) / length^2
218  // AB segment is represented by k = [0,1]
219  // Clamping k to [0,1] will give the shortest distance from P to AB segment
220  double dotprod = (px - l1x) * (l2x - l1x) + (py - l1y) * (l2y - l1y);
221  double k = dotprod / (length * length);
222  k = fmax(0.0, fmin(1.0, k));
223  double projx = l1x + k * (l2x - l1x);
224  double projy = l1y + k * (l2y - l1y);
225  return distance_point_point_squared(px, py, projx, projy);
226 }
227 
228 // Given three colinear points p, q, r, the function checks if
229 // point q lies on line segment 'pr'
231  double py,
232  double qx,
233  double qy,
234  double rx,
235  double ry) {
236  return (tol_le(qx, fmax(px, rx)) && tol_ge(qx, fmin(px, rx)) &&
237  tol_le(qy, fmax(py, ry)) && tol_ge(qy, fmin(py, ry)));
238 }
239 
240 DEVICE ALWAYS_INLINE int16_t
241 orientation(double px, double py, double qx, double qy, double rx, double ry) {
242  auto val = ((qy - py) * (rx - qx) - (qx - px) * (ry - qy));
243  if (tol_zero(val)) {
244  return 0; // Points p, q and r are colinear
245  }
246  if (val > 0.0) {
247  return 1; // Clockwise point orientation
248  }
249  return 2; // Counterclockwise point orientation
250 }
251 
252 // Cartesian intersection of two line segments l11-l12 and l21-l22
253 DEVICE
254 bool line_intersects_line(double l11x,
255  double l11y,
256  double l12x,
257  double l12y,
258  double l21x,
259  double l21y,
260  double l22x,
261  double l22y) {
262  auto o1 = orientation(l11x, l11y, l12x, l12y, l21x, l21y);
263  auto o2 = orientation(l11x, l11y, l12x, l12y, l22x, l22y);
264  auto o3 = orientation(l21x, l21y, l22x, l22y, l11x, l11y);
265  auto o4 = orientation(l21x, l21y, l22x, l22y, l12x, l12y);
266 
267  // General case
268  if (o1 != o2 && o3 != o4) {
269  return true;
270  }
271 
272  // Special Cases
273  // l11, l12 and l21 are colinear and l21 lies on segment l11-l12
274  if (o1 == 0 && on_segment(l11x, l11y, l21x, l21y, l12x, l12y)) {
275  return true;
276  }
277 
278  // l11, l12 and l21 are colinear and l22 lies on segment l11-l12
279  if (o2 == 0 && on_segment(l11x, l11y, l22x, l22y, l12x, l12y)) {
280  return true;
281  }
282 
283  // l21, l22 and l11 are colinear and l11 lies on segment l21-l22
284  if (o3 == 0 && on_segment(l21x, l21y, l11x, l11y, l22x, l22y)) {
285  return true;
286  }
287 
288  // l21, l22 and l12 are colinear and l12 lies on segment l21-l22
289  if (o4 == 0 && on_segment(l21x, l21y, l12x, l12y, l22x, l22y)) {
290  return true;
291  }
292 
293  return false;
294 }
295 
296 DEVICE
298  int32_t lnum_coords,
299  double l1x,
300  double l1y,
301  double l2x,
302  double l2y,
303  int32_t ic1,
304  int32_t isr1,
305  int32_t osr) {
306  double e1x = coord_x(l, 0, ic1, isr1, osr);
307  double e1y = coord_y(l, 1, ic1, isr1, osr);
308  for (int64_t i = 2; i < lnum_coords; i += 2) {
309  double e2x = coord_x(l, i, ic1, isr1, osr);
310  double e2y = coord_y(l, i + 1, ic1, isr1, osr);
311  if (line_intersects_line(e1x, e1y, e2x, e2y, l1x, l1y, l2x, l2y)) {
312  return true;
313  }
314  e1x = e2x;
315  e1y = e2y;
316  }
317  return false;
318 }
319 
320 DEVICE
321 bool ring_intersects_line(int8_t* ring,
322  int32_t ring_num_coords,
323  double l1x,
324  double l1y,
325  double l2x,
326  double l2y,
327  int32_t ic1,
328  int32_t isr1,
329  int32_t osr) {
330  double e1x = coord_x(ring, ring_num_coords - 2, ic1, isr1, osr);
331  double e1y = coord_y(ring, ring_num_coords - 1, ic1, isr1, osr);
332  double e2x = coord_x(ring, 0, ic1, isr1, osr);
333  double e2y = coord_y(ring, 1, ic1, isr1, osr);
334  if (line_intersects_line(e1x, e1y, e2x, e2y, l1x, l1y, l2x, l2y)) {
335  return true;
336  }
338  ring, ring_num_coords, l1x, l1y, l2x, l2y, ic1, isr1, osr);
339 }
340 
341 DEVICE
343  int32_t lnum_coords,
344  double l1x,
345  double l1y,
346  double l2x,
347  double l2y,
348  int32_t ic1,
349  int32_t isr1,
350  int32_t osr) {
351  double e1x = coord_x(l, 0, ic1, isr1, osr);
352  double e1y = coord_y(l, 1, ic1, isr1, osr);
353  for (int64_t i = 2; i < lnum_coords; i += 2) {
354  double e2x = coord_x(l, i, ic1, isr1, osr);
355  double e2y = coord_y(l, i + 1, ic1, isr1, osr);
356  if (line_intersects_line(e1x, e1y, e2x, e2y, l1x, l1y, l2x, l2y)) {
357  return true;
358  }
359  e1x = e2x;
360  e1y = e2y;
361  }
362  return false;
363 }
364 
365 // Cartesian distance between two line segments l11-l12 and l21-l22
366 DEVICE
367 double distance_line_line(double l11x,
368  double l11y,
369  double l12x,
370  double l12y,
371  double l21x,
372  double l21y,
373  double l22x,
374  double l22y) {
375  if (line_intersects_line(l11x, l11y, l12x, l12y, l21x, l21y, l22x, l22y)) {
376  return 0.0;
377  }
378  double dist12 = fmin(distance_point_line(l11x, l11y, l21x, l21y, l22x, l22y),
379  distance_point_line(l12x, l12y, l21x, l21y, l22x, l22y));
380  double dist21 = fmin(distance_point_line(l21x, l21y, l11x, l11y, l12x, l12y),
381  distance_point_line(l22x, l22y, l11x, l11y, l12x, l12y));
382  return fmin(dist12, dist21);
383 }
384 
385 // Cartesian squared distance between two line segments l11-l12 and l21-l22
386 DEVICE
387 double distance_line_line_squared(double l11x,
388  double l11y,
389  double l12x,
390  double l12y,
391  double l21x,
392  double l21y,
393  double l22x,
394  double l22y) {
395  if (line_intersects_line(l11x, l11y, l12x, l12y, l21x, l21y, l22x, l22y)) {
396  return 0.0;
397  }
398  double dist12 = fmin(distance_point_line_squared(l11x, l11y, l21x, l21y, l22x, l22y),
399  distance_point_line_squared(l12x, l12y, l21x, l21y, l22x, l22y));
400  double dist21 = fmin(distance_point_line_squared(l21x, l21y, l11x, l11y, l12x, l12y),
401  distance_point_line_squared(l22x, l22y, l11x, l11y, l12x, l12y));
402  return fmin(dist12, dist21);
403 }
404 
405 DEVICE
406 double distance_ring_linestring(int8_t* ring,
407  int32_t ring_num_coords,
408  int8_t* l,
409  int32_t lnum_coords,
410  int32_t ic1,
411  int32_t isr1,
412  int32_t ic2,
413  int32_t isr2,
414  int32_t osr,
415  double threshold) {
416  double min_distance = 0.0;
417 
418  double re1x = coord_x(ring, ring_num_coords - 2, ic1, isr1, osr);
419  double re1y = coord_y(ring, ring_num_coords - 1, ic1, isr1, osr);
420  for (auto i = 0; i < ring_num_coords; i += 2) {
421  double re2x = coord_x(ring, i, ic1, isr1, osr);
422  double re2y = coord_y(ring, i + 1, ic1, isr1, osr);
423 
424  double le1x = coord_x(l, 0, ic2, isr2, osr);
425  double le1y = coord_y(l, 1, ic2, isr2, osr);
426  for (auto j = 2; j < lnum_coords; j += 2) {
427  double le2x = coord_x(l, j, ic2, isr2, osr);
428  double le2y = coord_y(l, j + 1, ic2, isr2, osr);
429 
430  auto distance = distance_line_line(re1x, re1y, re2x, re2y, le1x, le1y, le2x, le2y);
431  if ((i == 0 && j == 2) || min_distance > distance) {
432  min_distance = distance;
433  if (tol_zero(min_distance)) {
434  return 0.0;
435  }
436  if (min_distance <= threshold) {
437  return min_distance;
438  }
439  }
440  le1x = le2x;
441  le1y = le2y;
442  }
443  re1x = re2x;
444  re1y = re2y;
445  }
446 
447  return min_distance;
448 }
449 
450 DEVICE
451 double distance_ring_ring(int8_t* ring1,
452  int32_t ring1_num_coords,
453  int8_t* ring2,
454  int32_t ring2_num_coords,
455  int32_t ic1,
456  int32_t isr1,
457  int32_t ic2,
458  int32_t isr2,
459  int32_t osr,
460  double threshold) {
461  double min_distance = 0.0;
462 
463  double e11x = coord_x(ring1, ring1_num_coords - 2, ic1, isr1, osr);
464  double e11y = coord_y(ring1, ring1_num_coords - 1, ic1, isr1, osr);
465  for (auto i = 0; i < ring1_num_coords; i += 2) {
466  double e12x = coord_x(ring1, i, ic1, isr1, osr);
467  double e12y = coord_y(ring1, i + 1, ic1, isr1, osr);
468 
469  double e21x = coord_x(ring2, ring2_num_coords - 2, ic2, isr2, osr);
470  double e21y = coord_y(ring2, ring2_num_coords - 1, ic2, isr2, osr);
471  for (auto j = 0; j < ring2_num_coords; j += 2) {
472  double e22x = coord_x(ring2, j, ic2, isr2, osr);
473  double e22y = coord_y(ring2, j + 1, ic2, isr2, osr);
474 
475  auto distance = distance_line_line(e11x, e11y, e12x, e12y, e21x, e21y, e22x, e22y);
476  if ((i == 0 && j == 0) || min_distance > distance) {
477  min_distance = distance;
478  if (tol_zero(min_distance)) {
479  return 0.0;
480  }
481  if (min_distance <= threshold) {
482  return min_distance;
483  }
484  }
485  e21x = e22x;
486  e21y = e22y;
487  }
488  e11x = e12x;
489  e11y = e12y;
490  }
491 
492  return min_distance;
493 }
494 
499 template <typename T>
501 is_left(const T lx0, const T ly0, const T lx1, const T ly1, const T px, const T py) {
502  return ((lx1 - lx0) * (py - ly0) - (px - lx0) * (ly1 - ly0));
503 }
504 
505 template <typename T>
507  const T tolerance = TOLERANCE_DEFAULT) {
508  if constexpr (!std::is_floating_point<T>::value) {
509  return x == 0; // assume integer 0s are always 0
510  }
511  return (-tolerance <= x) && (x <= tolerance);
512 }
513 
515 
529 template <typename T, EdgeBehavior TEdgeBehavior>
531  const int32_t poly_num_coords,
532  const T px,
533  const T py,
534  const int32_t ic1,
535  const int32_t isr1,
536  const int32_t osr) {
537  int32_t wn = 0;
538 
539  constexpr bool include_point_on_edge =
540  TEdgeBehavior == EdgeBehavior::kIncludePointOnEdge;
541 
542  auto get_x_coord = [=](const auto data, const auto index) -> T {
543  if constexpr (std::is_floating_point<T>::value) {
544  return coord_x(data, index, ic1, isr1, osr);
545  } else {
546  return compressed_coord(data, index);
547  }
548  return T{}; // https://stackoverflow.com/a/64561686/2700898
549  };
550 
551  auto get_y_coord = [=](const auto data, const auto index) -> T {
552  if constexpr (std::is_floating_point<T>::value) {
553  return coord_y(data, index, ic1, isr1, osr);
554  } else {
555  return compressed_coord(data, index);
556  }
557  return T{}; // https://stackoverflow.com/a/64561686/2700898
558  };
559 
560  DEBUG_STMT(printf("Poly num coords: %d\n", poly_num_coords));
561  const int64_t num_edges = poly_num_coords / 2;
562 
563  int64_t e0_index = 0;
564  T e0x = get_x_coord(poly, e0_index * 2);
565  T e0y = get_y_coord(poly, (e0_index * 2) + 1);
566 
567  int64_t e1_index;
568  T e1x, e1y;
569 
570  for (int64_t edge = 1; edge <= num_edges; edge++) {
571  DEBUG_STMT(printf("Edge %ld\n", edge));
572  // build the edge
573  e1_index = edge % num_edges;
574  e1x = get_x_coord(poly, e1_index * 2);
575  e1y = get_y_coord(poly, (e1_index * 2) + 1);
576 
577  DEBUG_STMT(printf("edge 0: %ld : %f, %f\n", e0_index, (double)e0x, (double)e0y));
578  DEBUG_STMT(printf("edge 1: %ld : %f, %f\n", e1_index, (double)e1x, (double)e1y));
579 
580  auto get_epsilon = [=]() -> T {
581  if constexpr (std::is_floating_point<T>::value) {
582  const T edge_vec_magnitude =
583  (e1x - e0x) * (e1x - e0x) + (e1y - e0y) * (e1y - e0y);
584  return 0.003 * edge_vec_magnitude;
585  } else {
586  return T(0);
587  }
588  return T{}; // https://stackoverflow.com/a/64561686/2700898
589  };
590 
591  const T epsilon = get_epsilon();
592  DEBUG_STMT(printf("using epsilon value %e\n", (double)epsilon));
593 
594  if constexpr (include_point_on_edge) {
595  const T xp = (px - e1x) * (e0y - e1y) - (py - e1y) * (e0x - e1x);
596  const T pt_vec_magnitude = (px - e0x) * (px - e0x) + (py - e0y) * (py - e0y);
597  const T edge_vec_magnitude = (e1x - e0x) * (e1x - e0x) + (e1y - e0y) * (e1y - e0y);
598  if (tol_zero_template(xp, epsilon) && (pt_vec_magnitude <= edge_vec_magnitude)) {
599  DEBUG_STMT(printf("point is on edge: %e && %e <= %e\n",
600  (double)xp,
601  (double)pt_vec_magnitude,
602  (double)edge_vec_magnitude));
603  return include_point_on_edge;
604  }
605  }
606 
607  if (e0y <= py) {
608  if (e1y > py) {
609  // upward crossing
610  DEBUG_STMT(printf("upward crossing\n"));
611  const auto is_left_val = is_left(e0x, e0y, e1x, e1y, px, py);
612  DEBUG_STMT(printf("Left val %f and tol zero left val %d\n",
613  (double)is_left_val,
614  tol_zero_template(is_left_val, epsilon)));
615  if (UNLIKELY(tol_zero_template(is_left_val, epsilon))) {
616  // p is on the edge
617  return include_point_on_edge;
618  } else if (is_left_val > T(0)) { // p left of edge
619  // valid up intersection
620  DEBUG_STMT(printf("++wn\n"));
621  ++wn;
622  }
623  }
624  } else if (e1y <= py) {
625  // downward crossing
626  DEBUG_STMT(printf("downward crossing\n"));
627  const auto is_left_val = is_left(e0x, e0y, e1x, e1y, px, py);
628  DEBUG_STMT(printf("Left val %f and tol zero left val %d\n",
629  (double)is_left_val,
630  tol_zero_template(is_left_val, epsilon)));
631  if (UNLIKELY(tol_zero_template(is_left_val, epsilon))) {
632  // p is on the edge
633  return include_point_on_edge;
634  } else if (is_left_val < T(0)) { // p right of edge
635  // valid down intersection
636  DEBUG_STMT(printf("--wn\n"));
637  --wn;
638  }
639  }
640 
641  e0_index = e1_index;
642  e0x = e1x;
643  e0y = e1y;
644  }
645  DEBUG_STMT(printf("wn: %d\n", wn));
646  // wn == 0 --> P is outside the polygon
647  return !(wn == 0);
648 }
649 
650 // Checks if a simple polygon (no holes) contains a point.
651 //
652 // Poly coords are extracted from raw data, based on compression (ic1) and input/output
653 // SRIDs (isr1/osr).
654 //
655 // Shoot a ray from point P to the right, register intersections with any of polygon's
656 // edges. Each intersection means entrance into or exit from the polygon. Odd number of
657 // intersections means the polygon does contain P. Account for special cases: touch+cross,
658 // touch+leave, touch+overlay+cross, touch+overlay+leave, on edge, etc.
659 //
660 // Secondary ray is shot from point P down for simple redundancy, to reduce main probe's
661 // chance of error. No intersections means P is outside, irrespective of main probe's
662 // result.
663 //
665  const int32_t poly_num_coords,
666  const double px,
667  const double py,
668  const int32_t ic1,
669  const int32_t isr1,
670  const int32_t osr) {
671  bool result = false;
672  int xray_touch = 0;
673  bool horizontal_edge = false;
674  bool yray_intersects = false;
675 
676  double e1x = coord_x(poly, poly_num_coords - 2, ic1, isr1, osr);
677  double e1y = coord_y(poly, poly_num_coords - 1, ic1, isr1, osr);
678  for (int64_t i = 0; i < poly_num_coords; i += 2) {
679  double e2x = coord_x(poly, i, ic1, isr1, osr);
680  double e2y = coord_y(poly, i + 1, ic1, isr1, osr);
681 
682  // Check if point sits on an edge.
683  if (tol_zero(distance_point_line(px, py, e1x, e1y, e2x, e2y))) {
684  return true;
685  }
686 
687  // Before flipping the switch, check if xray hit a horizontal edge
688  // - If an edge lays on the xray, one of the previous edges touched it
689  // so while moving horizontally we're in 'xray_touch' state
690  // - Last edge that touched xray at (e2x,e2y) didn't register intersection
691  // - Next edge that diverges from xray at (e1,e1y) will register intersection
692  // - Can have several horizontal edges, one after the other, keep moving though
693  // in 'xray_touch' state without flipping the switch
694  horizontal_edge = (xray_touch != 0) && tol_eq(py, e1y) && tol_eq(py, e2y);
695 
696  // Main probe: xray
697  // Overshoot the xray to detect an intersection if there is one.
698  double xray = fmax(e2x, e1x) + 1.0;
699  if (px <= xray && // Only check for intersection if the edge is on the right
700  !horizontal_edge && // Keep moving through horizontal edges
701  line_intersects_line(px, // xray shooting from point p to the right
702  py,
703  xray,
704  py,
705  e1x, // polygon edge
706  e1y,
707  e2x,
708  e2y)) {
709  // Register intersection
710  result = !result;
711 
712  // Adjust for special cases
713  if (xray_touch == 0) {
714  if (tol_zero(distance_point_line(e2x, e2y, px, py, xray + 1.0, py))) {
715  // Xray goes through the edge's second vertex, unregister intersection -
716  // that vertex will be crossed again when we look at the following edge(s)
717  result = !result;
718  // Enter the xray-touch state:
719  // (1) - xray was touched by the edge from above, (-1) from below
720  xray_touch = (e1y > py) ? 1 : -1;
721  }
722  } else {
723  // Previous edge touched the xray, intersection hasn't been registered,
724  // it has to be registered now if this edge continues across the xray.
725  if (xray_touch > 0) {
726  // Previous edge touched the xray from above
727  if (e2y <= py) {
728  // Current edge crosses under xray: intersection is already registered
729  } else {
730  // Current edge just touched the xray and pulled up: unregister intersection
731  result = !result;
732  }
733  } else {
734  // Previous edge touched the xray from below
735  if (e2y > py) {
736  // Current edge crosses over xray: intersection is already registered
737  } else {
738  // Current edge just touched the xray and pulled down: unregister intersection
739  result = !result;
740  }
741  }
742  // Exit the xray-touch state
743  xray_touch = 0;
744  }
745  }
746 
747  // Redundancy: vertical yray down
748  // Main probe xray may hit multiple complex fragments which increases a chance of
749  // error. Perform a simple secondary check for edge intersections to see if point is
750  // outside.
751  if (!yray_intersects) { // Continue checking on yray until intersection is found
752  double yray = fmin(e2y, e1y) - 1.0;
753  if (yray <= py) { // Only check for yray intersection if point P is above the edge
754  yray_intersects = line_intersects_line(px, // yray shooting from point P down
755  py,
756  px,
757  yray,
758  e1x, // polygon edge
759  e1y,
760  e2x,
761  e2y);
762  }
763  }
764 
765  // Advance to the next vertex
766  e1x = e2x;
767  e1y = e2y;
768  }
769  if (!yray_intersects) {
770  // yray has zero intersections - point is outside the polygon
771  return false;
772  }
773  // Otherwise rely on the main probe
774  return result;
775 }
776 
777 // Returns true if simple polygon (no holes) contains a linestring
778 DEVICE
779 bool polygon_contains_linestring(int8_t* poly,
780  int32_t poly_num_coords,
781  int8_t* l,
782  int64_t lnum_coords,
783  int32_t ic1,
784  int32_t isr1,
785  int32_t ic2,
786  int32_t isr2,
787  int32_t osr) {
788  // Check that the first point is in the polygon
789  double l1x = coord_x(l, 0, ic2, isr2, osr);
790  double l1y = coord_y(l, 1, ic2, isr2, osr);
791  if (!polygon_contains_point(poly, poly_num_coords, l1x, l1y, ic1, isr1, osr)) {
792  return false;
793  }
794 
795  // Go through line segments and check if there are no intersections with poly edges,
796  // i.e. linestring doesn't escape
797  for (int32_t i = 2; i < lnum_coords; i += 2) {
798  double l2x = coord_x(l, i, ic2, isr2, osr);
799  double l2y = coord_y(l, i + 1, ic2, isr2, osr);
800  if (ring_intersects_line(poly, poly_num_coords, l1x, l1y, l2x, l2y, ic1, isr1, osr)) {
801  return false;
802  }
803  l1x = l2x;
804  l1y = l2y;
805  }
806  return true;
807 }
808 
809 DEVICE ALWAYS_INLINE bool box_contains_point(const double* bounds,
810  const int64_t bounds_size,
811  const double px,
812  const double py) {
813  // TODO: mixed spatial references
814  return (tol_ge(px, bounds[0]) && tol_ge(py, bounds[1]) && tol_le(px, bounds[2]) &&
815  tol_le(py, bounds[3]));
816 }
817 
819  int64_t bounds_size,
820  double px,
821  double py) {
822  // TODO: mixed spatial references
823  return box_contains_point(bounds, bounds_size, px, py);
824 }
825 
826 DEVICE ALWAYS_INLINE bool box_contains_box(double* bounds1,
827  int64_t bounds1_size,
828  double* bounds2,
829  int64_t bounds2_size) {
830  // TODO: mixed spatial references
831  return (
833  bounds1, bounds1_size, bounds2[0], bounds2[1]) && // box1 <- box2: xmin, ymin
835  bounds1, bounds1_size, bounds2[2], bounds2[3])); // box1 <- box2: xmax, ymax
836 }
837 
839  int64_t bounds1_size,
840  double* bounds2,
841  int64_t bounds2_size) {
842  // TODO: mixed spatial references
843  return (
845  bounds1, bounds1_size, bounds2[0], bounds2[1]) || // box1 <- box2: xmin, ymin
847  bounds1, bounds1_size, bounds2[2], bounds2[3]) || // box1 <- box2: xmax, ymax
849  bounds1, bounds1_size, bounds2[0], bounds2[3]) || // box1 <- box2: xmin, ymax
851  bounds1, bounds1_size, bounds2[2], bounds2[1])); // box1 <- box2: xmax, ymin
852 }
853 
854 DEVICE ALWAYS_INLINE bool box_overlaps_box(double* bounds1,
855  int64_t bounds1_size,
856  double* bounds2,
857  int64_t bounds2_size) {
858  // TODO: tolerance
859  // TODO: mixed spatial references
860  if (bounds1[2] < bounds2[0] || // box1 is left of box2: box1.xmax < box2.xmin
861  bounds1[0] > bounds2[2] || // box1 is right of box2: box1.xmin > box2.xmax
862  bounds1[3] < bounds2[1] || // box1 is below box2: box1.ymax < box2.ymin
863  bounds1[1] > bounds2[3]) { // box1 is above box2: box1.ymin > box1.ymax
864  return false;
865  }
866  return true;
867 }
868 
870  int64_t p1size,
871  int32_t ic1,
872  int32_t isr1,
873  double* bounds2,
874  int64_t bounds2_size,
875  int32_t isr2,
876  int32_t osr,
877  double distance) {
878  // TODO: tolerance
879 
880  // Point has to be uncompressed and transformed to output SR.
881  // Bounding box has to be transformed to output SR.
882  auto px = coord_x(p1, 0, ic1, isr1, osr);
883  auto py = coord_y(p1, 1, ic1, isr1, osr);
884  if (
885  // point is left of box2: px < box2.xmin - distance
886  px < transform_coord_x(bounds2[0], isr2, osr) - distance ||
887  // point is right of box2: px > box2.xmax + distance
888  px > transform_coord_x(bounds2[2], isr2, osr) + distance ||
889  // point is below box2: py < box2.ymin - distance
890  py < transform_coord_y(bounds2[1], isr2, osr) - distance ||
891  // point is above box2: py > box2.ymax + distance
892  py > transform_coord_y(bounds2[3], isr2, osr) + distance) {
893  return false;
894  }
895  return true;
896 }
897 
898 DEVICE ALWAYS_INLINE bool box_dwithin_box(double* bounds1,
899  int64_t bounds1_size,
900  int32_t isr1,
901  double* bounds2,
902  int64_t bounds2_size,
903  int32_t isr2,
904  int32_t osr,
905  double distance) {
906  // TODO: tolerance
907 
908  // Bounds come in corresponding input spatial references.
909  // For example, here the first bounding box may come in 4326, second in 900913:
910  // ST_DWithin(ST_Transform(linestring4326, 900913), linestring900913, 10)
911  // Distance is given in osr spatial reference units, in this case 10 meters.
912  // Bounds should be converted to osr before calculations, if isr != osr.
913 
914  // TODO: revise all other functions that process bounds
915 
916  if (
917  // box1 is left of box2: box1.xmax + distance < box2.xmin
918  transform_coord_x(bounds1[2], isr1, osr) + distance <
919  transform_coord_x(bounds2[0], isr2, osr) ||
920  // box1 is right of box2: box1.xmin - distance > box2.xmax
921  transform_coord_x(bounds1[0], isr1, osr) - distance >
922  transform_coord_x(bounds2[2], isr2, osr) ||
923  // box1 is below box2: box1.ymax + distance < box2.ymin
924  transform_coord_y(bounds1[3], isr1, osr) + distance <
925  transform_coord_y(bounds2[1], isr2, osr) ||
926  // box1 is above box2: box1.ymin - distance > box1.ymax
927  transform_coord_y(bounds1[1], isr1, osr) - distance >
928  transform_coord_y(bounds2[3], isr2, osr)) {
929  return false;
930  }
931  return true;
932 }
933 
935 double ST_X_Point(int8_t* p, int64_t psize, int32_t ic, int32_t isr, int32_t osr) {
936  return coord_x(p, 0, ic, isr, osr);
937 }
938 
940 double ST_Y_Point(int8_t* p, int64_t psize, int32_t ic, int32_t isr, int32_t osr) {
941  return coord_y(p, 1, ic, isr, osr);
942 }
943 
945 double ST_XMin(int8_t* coords, int64_t size, int32_t ic, int32_t isr, int32_t osr) {
946  auto num_coords = size / compression_unit_size(ic);
947  double xmin = 0.0;
948  for (int32_t i = 0; i < num_coords; i += 2) {
949  double x = coord_x(coords, i, ic, isr, osr);
950  if (i == 0 || x < xmin) {
951  xmin = x;
952  }
953  }
954  return xmin;
955 }
956 
958 double ST_YMin(int8_t* coords, int64_t size, int32_t ic, int32_t isr, int32_t osr) {
959  auto num_coords = size / compression_unit_size(ic);
960  double ymin = 0.0;
961  for (int32_t i = 1; i < num_coords; i += 2) {
962  double y = coord_y(coords, i, ic, isr, osr);
963  if (i == 1 || y < ymin) {
964  ymin = y;
965  }
966  }
967  return ymin;
968 }
969 
971 double ST_XMax(int8_t* coords, int64_t size, int32_t ic, int32_t isr, int32_t osr) {
972  auto num_coords = size / compression_unit_size(ic);
973  double xmax = 0.0;
974  for (int32_t i = 0; i < num_coords; i += 2) {
975  double x = coord_x(coords, i, ic, isr, osr);
976  if (i == 0 || x > xmax) {
977  xmax = x;
978  }
979  }
980  return xmax;
981 }
982 
984 double ST_YMax(int8_t* coords, int64_t size, int32_t ic, int32_t isr, int32_t osr) {
985  auto num_coords = size / compression_unit_size(ic);
986  double ymax = 0.0;
987  for (int32_t i = 1; i < num_coords; i += 2) {
988  double y = coord_y(coords, i, ic, isr, osr);
989  if (i == 1 || y > ymax) {
990  ymax = y;
991  }
992  }
993  return ymax;
994 }
995 
997 double ST_XMin_Bounds(double* bounds, int64_t size, int32_t isr, int32_t osr) {
998  return transform_coord(bounds[0], isr, osr, true);
999 }
1000 
1002 double ST_YMin_Bounds(double* bounds, int64_t size, int32_t isr, int32_t osr) {
1003  return transform_coord(bounds[1], isr, osr, false);
1004 }
1005 
1007 double ST_XMax_Bounds(double* bounds, int64_t size, int32_t isr, int32_t osr) {
1008  return transform_coord(bounds[2], isr, osr, true);
1009 }
1010 
1012 double ST_YMax_Bounds(double* bounds, int64_t size, int32_t isr, int32_t osr) {
1013  return transform_coord(bounds[3], isr, osr, false);
1014 }
1015 
1016 //
1017 // ST_Length
1018 //
1019 
1021  int32_t lsize,
1022  int32_t ic,
1023  int32_t isr,
1024  int32_t osr,
1025  bool geodesic,
1026  bool check_closed) {
1027  auto l_num_coords = lsize / compression_unit_size(ic);
1028 
1029  double length = 0.0;
1030 
1031  double l0x = coord_x(l, 0, ic, isr, osr);
1032  double l0y = coord_y(l, 1, ic, isr, osr);
1033  double l2x = l0x;
1034  double l2y = l0y;
1035  for (int32_t i = 2; i < l_num_coords; i += 2) {
1036  double l1x = l2x;
1037  double l1y = l2y;
1038  l2x = coord_x(l, i, ic, isr, osr);
1039  l2y = coord_y(l, i + 1, ic, isr, osr);
1040  double ldist = geodesic ? distance_in_meters(l1x, l1y, l2x, l2y)
1041  : distance_point_point(l1x, l1y, l2x, l2y);
1042  length += ldist;
1043  }
1044  if (check_closed) {
1045  double ldist = geodesic ? distance_in_meters(l2x, l2y, l0x, l0y)
1046  : distance_point_point(l2x, l2y, l0x, l0y);
1047  length += ldist;
1048  }
1049  return length;
1050 }
1051 
1053 double ST_Length_LineString(int8_t* coords,
1054  int64_t coords_sz,
1055  int32_t ic,
1056  int32_t isr,
1057  int32_t osr) {
1058  return length_linestring(coords, coords_sz, ic, isr, osr, false, false);
1059 }
1060 
1062 double ST_Length_LineString_Geodesic(int8_t* coords,
1063  int64_t coords_sz,
1064  int32_t ic,
1065  int32_t isr,
1066  int32_t osr) {
1067  return length_linestring(coords, coords_sz, ic, isr, osr, true, false);
1068 }
1069 
1070 //
1071 // ST_Perimeter
1072 //
1073 
1075 double ST_Perimeter_Polygon(int8_t* poly,
1076  int32_t polysize,
1077  int8_t* poly_ring_sizes,
1078  int32_t poly_num_rings,
1079  int32_t ic,
1080  int32_t isr,
1081  int32_t osr) {
1082  if (poly_num_rings <= 0) {
1083  return 0.0;
1084  }
1085 
1086  auto exterior_ring_num_coords = poly_ring_sizes[0] * 2;
1087  auto exterior_ring_coords_size = exterior_ring_num_coords * compression_unit_size(ic);
1088 
1089  return length_linestring(poly, exterior_ring_coords_size, ic, isr, osr, false, true);
1090 }
1091 
1093 double ST_Perimeter_Polygon_Geodesic(int8_t* poly,
1094  int32_t polysize,
1095  int8_t* poly_ring_sizes_in,
1096  int32_t poly_num_rings,
1097  int32_t ic,
1098  int32_t isr,
1099  int32_t osr) {
1100  if (poly_num_rings <= 0) {
1101  return 0.0;
1102  }
1103 
1104  auto poly_ring_sizes = reinterpret_cast<int32_t*>(poly_ring_sizes_in);
1105  auto exterior_ring_num_coords = poly_ring_sizes[0] * 2;
1106  auto exterior_ring_coords_size = exterior_ring_num_coords * compression_unit_size(ic);
1107 
1108  return length_linestring(poly, exterior_ring_coords_size, ic, isr, osr, true, true);
1109 }
1110 
1111 DEVICE ALWAYS_INLINE double perimeter_multipolygon(int8_t* mpoly_coords,
1112  int32_t mpoly_coords_size,
1113  int8_t* mpoly_ring_sizes_in,
1114  int32_t mpoly_num_rings,
1115  int8_t* mpoly_poly_sizes,
1116  int32_t mpoly_num_polys,
1117  int32_t ic,
1118  int32_t isr,
1119  int32_t osr,
1120  bool geodesic) {
1121  if (mpoly_num_polys <= 0 || mpoly_num_rings <= 0) {
1122  return 0.0;
1123  }
1124 
1125  double perimeter = 0.0;
1126 
1127  auto mpoly_ring_sizes = reinterpret_cast<int32_t*>(mpoly_ring_sizes_in);
1128 
1129  // Set specific poly pointers as we move through the coords/ringsizes/polyrings arrays.
1130  auto next_poly_coords = mpoly_coords;
1131  auto next_poly_ring_sizes = mpoly_ring_sizes;
1132 
1133  for (auto poly = 0; poly < mpoly_num_polys; poly++) {
1134  auto poly_coords = next_poly_coords;
1135  auto poly_ring_sizes = next_poly_ring_sizes;
1136  auto poly_num_rings = reinterpret_cast<int32_t*>(mpoly_poly_sizes)[poly];
1137  // Count number of coords in all of poly's rings, advance ring size pointer.
1138  int32_t poly_num_coords = 0;
1139  for (auto ring = 0; ring < poly_num_rings; ring++) {
1140  poly_num_coords += 2 * *next_poly_ring_sizes++;
1141  }
1142  auto poly_coords_size = poly_num_coords * compression_unit_size(ic);
1143  next_poly_coords += poly_coords_size;
1144 
1145  auto exterior_ring_num_coords = poly_ring_sizes[0] * 2;
1146  auto exterior_ring_coords_size = exterior_ring_num_coords * compression_unit_size(ic);
1147 
1148  perimeter += length_linestring(
1149  poly_coords, exterior_ring_coords_size, ic, isr, osr, geodesic, true);
1150  }
1151 
1152  return perimeter;
1153 }
1154 
1156 double ST_Perimeter_MultiPolygon(int8_t* mpoly_coords,
1157  int32_t mpoly_coords_size,
1158  int8_t* mpoly_ring_sizes,
1159  int32_t mpoly_num_rings,
1160  int8_t* mpoly_poly_sizes,
1161  int32_t mpoly_num_polys,
1162  int32_t ic,
1163  int32_t isr,
1164  int32_t osr) {
1165  return perimeter_multipolygon(mpoly_coords,
1166  mpoly_coords_size,
1167  mpoly_ring_sizes,
1168  mpoly_num_rings,
1169  mpoly_poly_sizes,
1170  mpoly_num_polys,
1171  ic,
1172  isr,
1173  osr,
1174  false);
1175 }
1176 
1178 double ST_Perimeter_MultiPolygon_Geodesic(int8_t* mpoly_coords,
1179  int32_t mpoly_coords_size,
1180  int8_t* mpoly_ring_sizes,
1181  int32_t mpoly_num_rings,
1182  int8_t* mpoly_poly_sizes,
1183  int32_t mpoly_num_polys,
1184  int32_t ic,
1185  int32_t isr,
1186  int32_t osr) {
1187  return perimeter_multipolygon(mpoly_coords,
1188  mpoly_coords_size,
1189  mpoly_ring_sizes,
1190  mpoly_num_rings,
1191  mpoly_poly_sizes,
1192  mpoly_num_polys,
1193  ic,
1194  isr,
1195  osr,
1196  true);
1197 }
1198 
1199 //
1200 // ST_Area
1201 //
1202 
1204  double y1,
1205  double x2,
1206  double y2,
1207  double x3,
1208  double y3) {
1209  return (x1 * y2 - x2 * y1 + x3 * y1 - x1 * y3 + x2 * y3 - x3 * y2) / 2.0;
1210 }
1211 
1212 DEVICE ALWAYS_INLINE double area_ring(int8_t* ring,
1213  int64_t ringsize,
1214  int32_t ic,
1215  int32_t isr,
1216  int32_t osr) {
1217  auto ring_num_coords = ringsize / compression_unit_size(ic);
1218 
1219  if (ring_num_coords < 6) {
1220  return 0.0;
1221  }
1222 
1223  double area = 0.0;
1224 
1225  double x1 = coord_x(ring, 0, ic, isr, osr);
1226  double y1 = coord_y(ring, 1, ic, isr, osr);
1227  double x2 = coord_x(ring, 2, ic, isr, osr);
1228  double y2 = coord_y(ring, 3, ic, isr, osr);
1229  for (int32_t i = 4; i < ring_num_coords; i += 2) {
1230  double x3 = coord_x(ring, i, ic, isr, osr);
1231  double y3 = coord_y(ring, i + 1, ic, isr, osr);
1232  area += area_triangle(x1, y1, x2, y2, x3, y3);
1233  x2 = x3;
1234  y2 = y3;
1235  }
1236  return area;
1237 }
1238 
1239 DEVICE ALWAYS_INLINE double area_polygon(int8_t* poly_coords,
1240  int32_t poly_coords_size,
1241  int8_t* poly_ring_sizes_in,
1242  int32_t poly_num_rings,
1243  int32_t ic,
1244  int32_t isr,
1245  int32_t osr) {
1246  if (poly_num_rings <= 0) {
1247  return 0.0;
1248  }
1249 
1250  double area = 0.0;
1251  auto ring_coords = poly_coords;
1252  auto poly_ring_sizes = reinterpret_cast<int32_t*>(poly_ring_sizes_in);
1253 
1254  // Add up the areas of all rings.
1255  // External ring is CCW, open - positive area.
1256  // Internal rings (holes) are CW, open - negative areas.
1257  for (auto r = 0; r < poly_num_rings; r++) {
1258  auto ring_coords_size = poly_ring_sizes[r] * 2 * compression_unit_size(ic);
1259  area += area_ring(ring_coords, ring_coords_size, ic, isr, osr);
1260  // Advance to the next ring.
1261  ring_coords += ring_coords_size;
1262  }
1263  return area;
1264 }
1265 
1267 double ST_Area_Polygon(int8_t* poly_coords,
1268  int32_t poly_coords_size,
1269  int8_t* poly_ring_sizes,
1270  int32_t poly_num_rings,
1271  int32_t ic,
1272  int32_t isr,
1273  int32_t osr) {
1274  return area_polygon(
1275  poly_coords, poly_coords_size, poly_ring_sizes, poly_num_rings, ic, isr, osr);
1276 }
1277 
1279 double ST_Area_MultiPolygon(int8_t* mpoly_coords,
1280  int32_t mpoly_coords_size,
1281  int8_t* mpoly_ring_sizes,
1282  int32_t mpoly_num_rings,
1283  int8_t* mpoly_poly_sizes_in,
1284  int32_t mpoly_num_polys,
1285  int32_t ic,
1286  int32_t isr,
1287  int32_t osr) {
1288  if (mpoly_num_rings <= 0 || mpoly_num_polys <= 0) {
1289  return 0.0;
1290  }
1291 
1292  double area = 0.0;
1293 
1294  auto mpoly_poly_sizes = reinterpret_cast<int32_t*>(mpoly_poly_sizes_in);
1295 
1296  // Set specific poly pointers as we move through the coords/ringsizes/polyrings arrays.
1297  auto next_poly_coords = mpoly_coords;
1298  auto next_poly_ring_sizes = reinterpret_cast<int32_t*>(mpoly_ring_sizes);
1299 
1300  for (auto poly = 0; poly < mpoly_num_polys; poly++) {
1301  auto poly_coords = next_poly_coords;
1302  auto poly_ring_sizes = next_poly_ring_sizes;
1303  auto poly_num_rings = mpoly_poly_sizes[poly];
1304  // Count number of coords in all of poly's rings, advance ring size pointer.
1305  int32_t poly_num_coords = 0;
1306  for (auto ring = 0; ring < poly_num_rings; ring++) {
1307  poly_num_coords += 2 * *next_poly_ring_sizes++;
1308  }
1309  auto poly_coords_size = poly_num_coords * compression_unit_size(ic);
1310  next_poly_coords += poly_coords_size;
1311 
1312  area += area_polygon(poly_coords,
1313  poly_coords_size,
1314  reinterpret_cast<int8_t*>(poly_ring_sizes),
1315  poly_num_rings,
1316  ic,
1317  isr,
1318  osr);
1319  }
1320  return area;
1321 }
1322 
1323 //
1324 // ST_Centroid
1325 //
1326 
1327 // Point centroid
1328 // - Single point - drop, otherwise calculate average coordinates for all points
1329 // LineString centroid
1330 // - take midpoint of each line segment
1331 // - calculate average coordinates of all midpoints weighted by segment length
1332 // Polygon, MultiPolygon centroid, holes..
1333 // - weighted sum of centroids of triangles coming out of area decomposition
1334 //
1335 // On zero area fall back to linestring centroid
1336 // On zero length fall back to point centroid
1337 
1339 void ST_Centroid_Point(int8_t* p,
1340  int32_t psize,
1341  int32_t ic,
1342  int32_t isr,
1343  int32_t osr,
1344  double* point_centroid) {
1345  point_centroid[0] = coord_y(p, 0, ic, isr, osr);
1346  point_centroid[1] = coord_x(p, 1, ic, isr, osr);
1347 }
1348 
1350  double y1,
1351  double x2,
1352  double y2,
1353  double* length,
1354  double* linestring_centroid_sum) {
1355  double ldist = distance_point_point(x1, y1, x2, y2);
1356  *length += ldist;
1357  double segment_midpoint_x = (x1 + x2) / 2.0;
1358  double segment_midpoint_y = (y1 + y2) / 2.0;
1359  linestring_centroid_sum[0] += ldist * segment_midpoint_x;
1360  linestring_centroid_sum[1] += ldist * segment_midpoint_y;
1361  return true;
1362 }
1363 
1365  int64_t lsize,
1366  int32_t ic,
1367  int32_t isr,
1368  int32_t osr,
1369  bool closed,
1370  double* total_length,
1371  double* linestring_centroid_sum,
1372  int64_t* num_points,
1373  double* point_centroid_sum) {
1374  auto l_num_coords = lsize / compression_unit_size(ic);
1375  double length = 0.0;
1376  double l0x = coord_x(l, 0, ic, isr, osr);
1377  double l0y = coord_y(l, 1, ic, isr, osr);
1378  double l2x = l0x;
1379  double l2y = l0y;
1380  for (int32_t i = 2; i < l_num_coords; i += 2) {
1381  double l1x = l2x;
1382  double l1y = l2y;
1383  l2x = coord_x(l, i, ic, isr, osr);
1384  l2y = coord_y(l, i + 1, ic, isr, osr);
1385  centroid_add_segment(l1x, l1y, l2x, l2y, &length, linestring_centroid_sum);
1386  }
1387  if (l_num_coords > 4 && closed) {
1388  // Also add the closing segment between the last and the first points
1389  centroid_add_segment(l2x, l2y, l0x, l0y, &length, linestring_centroid_sum);
1390  }
1391  *total_length += length;
1392  if (length == 0.0 && l_num_coords > 0) {
1393  *num_points += 1;
1394  point_centroid_sum[0] += l0x;
1395  point_centroid_sum[1] += l0y;
1396  }
1397  return true;
1398 }
1399 
1401 void ST_Centroid_LineString(int8_t* coords,
1402  int32_t coords_sz,
1403  int32_t ic,
1404  int32_t isr,
1405  int32_t osr,
1406  double* linestring_centroid) {
1407  double length = 0.0;
1408  double linestring_centroid_sum[2] = {0.0, 0.0};
1409  int64_t num_points = 0;
1410  double point_centroid_sum[2] = {0.0, 0.0};
1411  centroid_add_linestring(coords,
1412  coords_sz,
1413  ic,
1414  isr,
1415  osr,
1416  false, // not closed
1417  &length,
1418  &linestring_centroid_sum[0],
1419  &num_points,
1420  &point_centroid_sum[0]);
1421  if (length > 0) {
1422  linestring_centroid[0] = linestring_centroid_sum[0] / length;
1423  linestring_centroid[1] = linestring_centroid_sum[1] / length;
1424  } else if (num_points > 0) {
1425  linestring_centroid[0] = point_centroid_sum[0] / num_points;
1426  linestring_centroid[1] = point_centroid_sum[1] / num_points;
1427  }
1428 }
1429 
1431  double y1,
1432  double x2,
1433  double y2,
1434  double x3,
1435  double y3,
1436  double sign,
1437  double* total_area2,
1438  double* cg3) {
1439  double cx = x1 + x2 + x3;
1440  double cy = y1 + y2 + y3;
1441  double area2 = x1 * y2 - x2 * y1 + x3 * y1 - x1 * y3 + x2 * y3 - x3 * y2;
1442  cg3[0] += sign * area2 * cx;
1443  cg3[1] += sign * area2 * cy;
1444  *total_area2 += sign * area2;
1445  return true;
1446 }
1447 
1449  int64_t ringsize,
1450  int32_t ic,
1451  int32_t isr,
1452  int32_t osr,
1453  double sign,
1454  double* total_area2,
1455  double* cg3,
1456  double* total_length,
1457  double* linestring_centroid_sum,
1458  int64_t* num_points,
1459  double* point_centroid_sum) {
1460  auto ring_num_coords = ringsize / compression_unit_size(ic);
1461 
1462  if (ring_num_coords < 6) {
1463  return false;
1464  }
1465 
1466  double x1 = coord_x(ring, 0, ic, isr, osr);
1467  double y1 = coord_y(ring, 1, ic, isr, osr);
1468  double x2 = coord_x(ring, 2, ic, isr, osr);
1469  double y2 = coord_y(ring, 3, ic, isr, osr);
1470  for (int32_t i = 4; i < ring_num_coords; i += 2) {
1471  double x3 = coord_x(ring, i, ic, isr, osr);
1472  double y3 = coord_y(ring, i + 1, ic, isr, osr);
1473  centroid_add_triangle(x1, y1, x2, y2, x3, y3, sign, total_area2, cg3);
1474  x2 = x3;
1475  y2 = y3;
1476  }
1477 
1479  ringsize,
1480  ic,
1481  isr,
1482  osr,
1483  true, // closed
1484  total_length,
1485  linestring_centroid_sum,
1486  num_points,
1487  point_centroid_sum);
1488  return true;
1489 }
1490 
1491 DEVICE ALWAYS_INLINE bool centroid_add_polygon(int8_t* poly_coords,
1492  int64_t poly_coords_size,
1493  int32_t* poly_ring_sizes,
1494  int64_t poly_num_rings,
1495  int32_t ic,
1496  int32_t isr,
1497  int32_t osr,
1498  double* total_area2,
1499  double* cg3,
1500  double* total_length,
1501  double* linestring_centroid_sum,
1502  int64_t* num_points,
1503  double* point_centroid_sum) {
1504  if (poly_num_rings <= 0) {
1505  return false;
1506  }
1507 
1508  auto ring_coords = poly_coords;
1509 
1510  for (auto r = 0; r < poly_num_rings; r++) {
1511  auto ring_coords_size = poly_ring_sizes[r] * 2 * compression_unit_size(ic);
1512  // Shape - positive area, holes - negative areas
1513  double sign = (r == 0) ? 1.0 : -1.0;
1514  centroid_add_ring(ring_coords,
1515  ring_coords_size,
1516  ic,
1517  isr,
1518  osr,
1519  sign,
1520  total_area2,
1521  cg3,
1522  total_length,
1523  linestring_centroid_sum,
1524  num_points,
1525  point_centroid_sum);
1526  // Advance to the next ring.
1527  ring_coords += ring_coords_size;
1528  }
1529  return true;
1530 }
1531 
1533 void ST_Centroid_Polygon(int8_t* poly_coords,
1534  int32_t poly_coords_size,
1535  int32_t* poly_ring_sizes,
1536  int32_t poly_num_rings,
1537  int32_t ic,
1538  int32_t isr,
1539  int32_t osr,
1540  double* poly_centroid) {
1541  if (poly_num_rings <= 0) {
1542  poly_centroid[0] = 0.0;
1543  poly_centroid[1] = 0.0;
1544  }
1545  double total_area2 = 0.0;
1546  double cg3[2] = {0.0, 0.0};
1547  double total_length = 0.0;
1548  double linestring_centroid_sum[2] = {0.0, 0.0};
1549  int64_t num_points = 0;
1550  double point_centroid_sum[2] = {0.0, 0.0};
1551  centroid_add_polygon(poly_coords,
1552  poly_coords_size,
1553  poly_ring_sizes,
1554  poly_num_rings,
1555  ic,
1556  isr,
1557  osr,
1558  &total_area2,
1559  &cg3[0],
1560  &total_length,
1561  &linestring_centroid_sum[0],
1562  &num_points,
1563  &point_centroid_sum[0]);
1564 
1565  if (total_area2 != 0.0) {
1566  poly_centroid[0] = cg3[0] / 3 / total_area2;
1567  poly_centroid[1] = cg3[1] / 3 / total_area2;
1568  } else if (total_length > 0.0) {
1569  // zero-area polygon, fall back to linear centroid
1570  poly_centroid[0] = linestring_centroid_sum[0] / total_length;
1571  poly_centroid[1] = linestring_centroid_sum[1] / total_length;
1572  } else if (num_points > 0) {
1573  // zero-area zero-length polygon, fall further back to point centroid
1574  poly_centroid[0] = point_centroid_sum[0] / num_points;
1575  poly_centroid[1] = point_centroid_sum[1] / num_points;
1576  } else {
1577  poly_centroid[0] = coord_x(poly_coords, 0, ic, isr, osr);
1578  poly_centroid[1] = coord_y(poly_coords, 1, ic, isr, osr);
1579  }
1580 }
1581 
1583 void ST_Centroid_MultiPolygon(int8_t* mpoly_coords,
1584  int32_t mpoly_coords_size,
1585  int32_t* mpoly_ring_sizes,
1586  int32_t mpoly_num_rings,
1587  int32_t* mpoly_poly_sizes,
1588  int32_t mpoly_num_polys,
1589  int32_t ic,
1590  int32_t isr,
1591  int32_t osr,
1592  double* mpoly_centroid) {
1593  if (mpoly_num_rings <= 0 || mpoly_num_polys <= 0) {
1594  mpoly_centroid[0] = 0.0;
1595  mpoly_centroid[1] = 0.0;
1596  }
1597 
1598  double total_area2 = 0.0;
1599  double cg3[2] = {0.0, 0.0};
1600  double total_length = 0.0;
1601  double linestring_centroid_sum[2] = {0.0, 0.0};
1602  int64_t num_points = 0;
1603  double point_centroid_sum[2] = {0.0, 0.0};
1604 
1605  // Set specific poly pointers as we move through the coords/ringsizes/polyrings arrays.
1606  auto next_poly_coords = mpoly_coords;
1607  auto next_poly_ring_sizes = mpoly_ring_sizes;
1608 
1609  for (auto poly = 0; poly < mpoly_num_polys; poly++) {
1610  auto poly_coords = next_poly_coords;
1611  auto poly_ring_sizes = next_poly_ring_sizes;
1612  auto poly_num_rings = mpoly_poly_sizes[poly];
1613  // Count number of coords in all of poly's rings, advance ring size pointer.
1614  int32_t poly_num_coords = 0;
1615  for (auto ring = 0; ring < poly_num_rings; ring++) {
1616  poly_num_coords += 2 * *next_poly_ring_sizes++;
1617  }
1618  auto poly_coords_size = poly_num_coords * compression_unit_size(ic);
1619  next_poly_coords += poly_coords_size;
1620 
1621  centroid_add_polygon(poly_coords,
1622  poly_coords_size,
1623  poly_ring_sizes,
1624  poly_num_rings,
1625  ic,
1626  isr,
1627  osr,
1628  &total_area2,
1629  &cg3[0],
1630  &total_length,
1631  &linestring_centroid_sum[0],
1632  &num_points,
1633  &point_centroid_sum[0]);
1634  }
1635 
1636  if (total_area2 != 0.0) {
1637  mpoly_centroid[0] = cg3[0] / 3 / total_area2;
1638  mpoly_centroid[1] = cg3[1] / 3 / total_area2;
1639  } else if (total_length > 0.0) {
1640  // zero-area multipolygon, fall back to linear centroid
1641  mpoly_centroid[0] = linestring_centroid_sum[0] / total_length;
1642  mpoly_centroid[1] = linestring_centroid_sum[1] / total_length;
1643  } else if (num_points > 0) {
1644  // zero-area zero-length multipolygon, fall further back to point centroid
1645  mpoly_centroid[0] = point_centroid_sum[0] / num_points;
1646  mpoly_centroid[1] = point_centroid_sum[1] / num_points;
1647  } else {
1648  mpoly_centroid[0] = coord_x(mpoly_coords, 0, ic, isr, osr);
1649  mpoly_centroid[1] = coord_y(mpoly_coords, 1, ic, isr, osr);
1650  }
1651 }
1652 
1653 //
1654 // ST_Distance
1655 //
1656 
1658 double ST_Distance_Point_Point(int8_t* p1,
1659  int64_t p1size,
1660  int8_t* p2,
1661  int64_t p2size,
1662  int32_t ic1,
1663  int32_t isr1,
1664  int32_t ic2,
1665  int32_t isr2,
1666  int32_t osr) {
1667  double p1x = coord_x(p1, 0, ic1, isr1, osr);
1668  double p1y = coord_y(p1, 1, ic1, isr1, osr);
1669  double p2x = coord_x(p2, 0, ic2, isr2, osr);
1670  double p2y = coord_y(p2, 1, ic2, isr2, osr);
1671  return distance_point_point(p1x, p1y, p2x, p2y);
1672 }
1673 
1676  int64_t p1size,
1677  int8_t* p2,
1678  int64_t p2size,
1679  int32_t ic1,
1680  int32_t isr1,
1681  int32_t ic2,
1682  int32_t isr2,
1683  int32_t osr) {
1684  double p1x = coord_x(p1, 0, ic1, isr1, osr);
1685  double p1y = coord_y(p1, 1, ic1, isr1, osr);
1686  double p2x = coord_x(p2, 0, ic2, isr2, osr);
1687  double p2y = coord_y(p2, 1, ic2, isr2, osr);
1688  return distance_point_point_squared(p1x, p1y, p2x, p2y);
1689 }
1690 
1693  int64_t p1size,
1694  int8_t* p2,
1695  int64_t p2size,
1696  int32_t ic1,
1697  int32_t isr1,
1698  int32_t ic2,
1699  int32_t isr2,
1700  int32_t osr) {
1701  double p1x = coord_x(p1, 0, ic1, 4326, 4326);
1702  double p1y = coord_y(p1, 1, ic1, 4326, 4326);
1703  double p2x = coord_x(p2, 0, ic2, 4326, 4326);
1704  double p2y = coord_y(p2, 1, ic2, 4326, 4326);
1705  return distance_in_meters(p1x, p1y, p2x, p2y);
1706 }
1707 
1710  int64_t psize,
1711  int8_t* l,
1712  int64_t lsize,
1713  int32_t ic1,
1714  int32_t isr1,
1715  int32_t ic2,
1716  int32_t isr2,
1717  int32_t osr) {
1718  // Currently only statically indexed LineString is supported
1719  double px = coord_x(p, 0, ic1, 4326, 4326);
1720  double py = coord_y(p, 1, ic1, 4326, 4326);
1721  const auto lpoints = lsize / (2 * compression_unit_size(ic2));
1722  double lx = coord_x(l, 2 * (lpoints - 1), ic2, 4326, 4326);
1723  double ly = coord_y(l, 2 * (lpoints - 1) + 1, ic2, 4326, 4326);
1724  return distance_in_meters(px, py, lx, ly);
1725 }
1726 
1729  int64_t lsize,
1730  int8_t* p,
1731  int64_t psize,
1732  int32_t ic1,
1733  int32_t isr1,
1734  int32_t ic2,
1735  int32_t isr2,
1736  int32_t osr) {
1737  // Currently only statically indexed LineString is supported
1739  p, psize, l, lsize, ic2, isr2, ic1, isr1, osr);
1740 }
1741 
1743  int64_t psize,
1744  int8_t* l,
1745  int64_t lsize,
1746  int32_t ic1,
1747  int32_t isr1,
1748  int32_t ic2,
1749  int32_t isr2,
1750  int32_t osr,
1751  bool check_closed,
1752  double threshold) {
1753  double px = coord_x(p, 0, ic1, isr1, osr);
1754  double py = coord_y(p, 1, ic1, isr1, osr);
1755 
1756  auto l_num_coords = lsize / compression_unit_size(ic2);
1757 
1758  double l1x = coord_x(l, 0, ic2, isr2, osr);
1759  double l1y = coord_y(l, 1, ic2, isr2, osr);
1760  double l2x = coord_x(l, 2, ic2, isr2, osr);
1761  double l2y = coord_y(l, 3, ic2, isr2, osr);
1762 
1763  double dist = distance_point_line(px, py, l1x, l1y, l2x, l2y);
1764  for (int32_t i = 4; i < l_num_coords; i += 2) {
1765  l1x = l2x; // advance one point
1766  l1y = l2y;
1767  l2x = coord_x(l, i, ic2, isr2, osr);
1768  l2y = coord_y(l, i + 1, ic2, isr2, osr);
1769  double ldist = distance_point_line(px, py, l1x, l1y, l2x, l2y);
1770  if (dist > ldist) {
1771  dist = ldist;
1772  }
1773  if (dist <= threshold) {
1774  return dist;
1775  }
1776  }
1777  if (l_num_coords > 4 && check_closed) {
1778  // Also check distance to the closing edge between the first and the last points
1779  l1x = coord_x(l, 0, ic2, isr2, osr);
1780  l1y = coord_y(l, 1, ic2, isr2, osr);
1781  double ldist = distance_point_line(px, py, l1x, l1y, l2x, l2y);
1782  if (dist > ldist) {
1783  dist = ldist;
1784  }
1785  }
1786  return dist;
1787 }
1788 
1791  int64_t psize,
1792  int8_t* l,
1793  int64_t lsize,
1794  int32_t ic1,
1795  int32_t isr1,
1796  int32_t ic2,
1797  int32_t isr2,
1798  int32_t osr,
1799  double threshold) {
1801  p, psize, l, lsize, ic1, isr1, ic2, isr2, osr, true, threshold);
1802 }
1803 
1806  int64_t psize,
1807  int8_t* l,
1808  int64_t lsize,
1809  int32_t ic1,
1810  int32_t isr1,
1811  int32_t ic2,
1812  int32_t isr2,
1813  int32_t osr,
1814  double threshold) {
1816  p, psize, l, lsize, ic1, isr1, ic2, isr2, osr, false, threshold);
1817 }
1818 
1820 double ST_Distance_Point_Polygon(int8_t* p,
1821  int64_t psize,
1822  int8_t* poly,
1823  int64_t polysize,
1824  int32_t* poly_ring_sizes,
1825  int64_t poly_num_rings,
1826  int32_t ic1,
1827  int32_t isr1,
1828  int32_t ic2,
1829  int32_t isr2,
1830  int32_t osr,
1831  double threshold) {
1832  auto exterior_ring_num_coords = polysize / compression_unit_size(ic2);
1833  if (poly_num_rings > 0) {
1834  exterior_ring_num_coords = poly_ring_sizes[0] * 2;
1835  }
1836  auto exterior_ring_coords_size = exterior_ring_num_coords * compression_unit_size(ic2);
1837 
1838  double px = coord_x(p, 0, ic1, isr1, osr);
1839  double py = coord_y(p, 1, ic1, isr1, osr);
1840  if (!polygon_contains_point(poly, exterior_ring_num_coords, px, py, ic2, isr2, osr)) {
1841  // Outside the exterior ring
1843  p, psize, poly, exterior_ring_coords_size, ic1, isr1, ic2, isr2, osr, threshold);
1844  }
1845  // Inside exterior ring
1846  // Advance to first interior ring
1847  poly += exterior_ring_coords_size;
1848  // Check if one of the polygon's holes contains that point
1849  for (auto r = 1; r < poly_num_rings; r++) {
1850  auto interior_ring_num_coords = poly_ring_sizes[r] * 2;
1851  auto interior_ring_coords_size =
1852  interior_ring_num_coords * compression_unit_size(ic2);
1853  if (polygon_contains_point(poly, interior_ring_num_coords, px, py, ic2, isr2, osr)) {
1854  // Inside an interior ring
1856  psize,
1857  poly,
1858  interior_ring_coords_size,
1859  ic1,
1860  isr1,
1861  ic2,
1862  isr2,
1863  osr,
1864  threshold);
1865  }
1866  poly += interior_ring_coords_size;
1867  }
1868  return 0.0;
1869 }
1870 
1873  int64_t psize,
1874  int8_t* mpoly_coords,
1875  int64_t mpoly_coords_size,
1876  int32_t* mpoly_ring_sizes,
1877  int64_t mpoly_num_rings,
1878  int32_t* mpoly_poly_sizes,
1879  int64_t mpoly_num_polys,
1880  int32_t ic1,
1881  int32_t isr1,
1882  int32_t ic2,
1883  int32_t isr2,
1884  int32_t osr,
1885  double threshold) {
1886  if (mpoly_num_polys <= 0) {
1887  return 0.0;
1888  }
1889  double min_distance = 0.0;
1890 
1891  // Set specific poly pointers as we move through the coords/ringsizes/polyrings arrays.
1892  auto next_poly_coords = mpoly_coords;
1893  auto next_poly_ring_sizes = mpoly_ring_sizes;
1894 
1895  for (auto poly = 0; poly < mpoly_num_polys; poly++) {
1896  auto poly_coords = next_poly_coords;
1897  auto poly_ring_sizes = next_poly_ring_sizes;
1898  auto poly_num_rings = mpoly_poly_sizes[poly];
1899  // Count number of coords in all of poly's rings, advance ring size pointer.
1900  int32_t poly_num_coords = 0;
1901  for (auto ring = 0; ring < poly_num_rings; ring++) {
1902  poly_num_coords += 2 * *next_poly_ring_sizes++;
1903  }
1904  auto poly_coords_size = poly_num_coords * compression_unit_size(ic2);
1905  next_poly_coords += poly_coords_size;
1906  double distance = ST_Distance_Point_Polygon(p,
1907  psize,
1908  poly_coords,
1909  poly_coords_size,
1910  poly_ring_sizes,
1911  poly_num_rings,
1912  ic1,
1913  isr1,
1914  ic2,
1915  isr2,
1916  osr,
1917  threshold);
1918  if (poly == 0 || min_distance > distance) {
1919  min_distance = distance;
1920  if (tol_zero(min_distance)) {
1921  min_distance = 0.0;
1922  break;
1923  }
1924  if (min_distance <= threshold) {
1925  break;
1926  }
1927  }
1928  }
1929 
1930  return min_distance;
1931 }
1932 
1935  int64_t lsize,
1936  int8_t* p,
1937  int64_t psize,
1938  int32_t ic1,
1939  int32_t isr1,
1940  int32_t ic2,
1941  int32_t isr2,
1942  int32_t osr,
1943  double threshold) {
1945  p, psize, l, lsize, ic2, isr2, ic1, isr1, osr, threshold);
1946 }
1947 
1950  int64_t l1size,
1951  int8_t* l2,
1952  int64_t l2size,
1953  int32_t ic1,
1954  int32_t isr1,
1955  int32_t ic2,
1956  int32_t isr2,
1957  int32_t osr,
1958  double threshold) {
1959  auto l1_num_coords = l1size / compression_unit_size(ic1);
1960  auto l2_num_coords = l2size / compression_unit_size(ic2);
1961 
1962  double threshold_squared = threshold * threshold;
1963  double dist_squared = 0.0;
1964  double l11x = coord_x(l1, 0, ic1, isr1, osr);
1965  double l11y = coord_y(l1, 1, ic1, isr1, osr);
1966  for (int32_t i1 = 2; i1 < l1_num_coords; i1 += 2) {
1967  double l12x = coord_x(l1, i1, ic1, isr1, osr);
1968  double l12y = coord_y(l1, i1 + 1, ic1, isr1, osr);
1969 
1970  double l21x = coord_x(l2, 0, ic2, isr2, osr);
1971  double l21y = coord_y(l2, 1, ic2, isr2, osr);
1972  for (int32_t i2 = 2; i2 < l2_num_coords; i2 += 2) {
1973  double l22x = coord_x(l2, i2, ic2, isr2, osr);
1974  double l22y = coord_y(l2, i2 + 1, ic2, isr2, osr);
1975 
1976  // double ldist_squared =
1977  // distance_line_line_squared(l11x, l11y, l12x, l12y, l21x, l21y, l22x, l22y);
1978  // TODO: fix distance_line_line_squared
1979  double ldist = distance_line_line(l11x, l11y, l12x, l12y, l21x, l21y, l22x, l22y);
1980  double ldist_squared = ldist * ldist;
1981 
1982  if (i1 == 2 && i2 == 2) {
1983  dist_squared = ldist_squared; // initialize dist with distance between the first
1984  // two segments
1985  } else if (dist_squared > ldist_squared) {
1986  dist_squared = ldist_squared;
1987  }
1988  if (tol_zero(dist_squared, TOLERANCE_DEFAULT_SQUARED)) {
1989  return 0.0; // Segments touch, short-circuit
1990  }
1991  if (dist_squared <= threshold_squared) {
1992  // If threashold_squared is defined and the calculated dist_squared dips under it,
1993  // short-circuit and return it
1994  return sqrt(dist_squared);
1995  }
1996 
1997  l21x = l22x; // advance to the next point on l2
1998  l21y = l22y;
1999  }
2000 
2001  l11x = l12x; // advance to the next point on l1
2002  l11y = l12y;
2003  }
2004  return sqrt(dist_squared);
2005 }
2006 
2009  int64_t lsize,
2010  int8_t* poly_coords,
2011  int64_t poly_coords_size,
2012  int32_t* poly_ring_sizes,
2013  int64_t poly_num_rings,
2014  int32_t ic1,
2015  int32_t isr1,
2016  int32_t ic2,
2017  int32_t isr2,
2018  int32_t osr,
2019  double threshold) {
2020  auto lnum_coords = lsize / compression_unit_size(ic1);
2021  auto psize = 2 * compression_unit_size(ic1);
2022  auto min_distance =
2023  ST_Distance_Point_Polygon(l, // pointer to start of linestring for first point
2024  psize,
2025  poly_coords,
2026  poly_coords_size,
2027  poly_ring_sizes,
2028  poly_num_rings,
2029  ic1,
2030  isr1,
2031  ic2,
2032  isr2,
2033  osr,
2034  threshold);
2035  if (tol_zero(min_distance)) {
2036  // Linestring's first point is inside the poly
2037  return 0.0;
2038  }
2039  if (min_distance <= threshold) {
2040  return min_distance;
2041  }
2042 
2043  // Otherwise, linestring's first point is outside the external ring or inside
2044  // an internal ring. Measure minimum distance between linestring segments and
2045  // poly rings. Crossing a ring zeroes the distance and causes an early return.
2046  auto poly_ring_coords = poly_coords;
2047  for (auto r = 0; r < poly_num_rings; r++) {
2048  int64_t poly_ring_num_coords = poly_ring_sizes[r] * 2;
2049 
2050  auto distance = distance_ring_linestring(poly_ring_coords,
2051  poly_ring_num_coords,
2052  l,
2053  lnum_coords,
2054  ic2,
2055  isr2,
2056  ic1,
2057  isr1,
2058  osr,
2059  threshold);
2060  if (min_distance > distance) {
2061  min_distance = distance;
2062  if (tol_zero(min_distance)) {
2063  return 0.0;
2064  }
2065  if (min_distance <= threshold) {
2066  return min_distance;
2067  }
2068  }
2069 
2070  poly_ring_coords += poly_ring_num_coords * compression_unit_size(ic2);
2071  }
2072 
2073  return min_distance;
2074 }
2075 
2078  int64_t lsize,
2079  int8_t* mpoly_coords,
2080  int64_t mpoly_coords_size,
2081  int32_t* mpoly_ring_sizes,
2082  int64_t mpoly_num_rings,
2083  int32_t* mpoly_poly_sizes,
2084  int64_t mpoly_num_polys,
2085  int32_t ic1,
2086  int32_t isr1,
2087  int32_t ic2,
2088  int32_t isr2,
2089  int32_t osr,
2090  double threshold) {
2091  // TODO: revisit implementation, cover all cases
2092  double min_distance = 0.0;
2093 
2094  // Set specific poly pointers as we move through the coords/ringsizes/polyrings arrays.
2095  auto next_poly_coords = mpoly_coords;
2096  auto next_poly_ring_sizes = mpoly_ring_sizes;
2097 
2098  for (auto poly = 0; poly < mpoly_num_polys; poly++) {
2099  auto poly_coords = next_poly_coords;
2100  auto poly_ring_sizes = next_poly_ring_sizes;
2101  auto poly_num_rings = mpoly_poly_sizes[poly];
2102  // Count number of coords in all of poly's rings, advance ring size pointer.
2103  int32_t poly_num_coords = 0;
2104  for (auto ring = 0; ring < poly_num_rings; ring++) {
2105  poly_num_coords += 2 * *next_poly_ring_sizes++;
2106  }
2107  auto poly_coords_size = poly_num_coords * compression_unit_size(ic2);
2108  next_poly_coords += poly_coords_size;
2109  double distance = ST_Distance_LineString_Polygon(l,
2110  lsize,
2111  poly_coords,
2112  poly_coords_size,
2113  poly_ring_sizes,
2114  poly_num_rings,
2115  ic1,
2116  isr1,
2117  ic2,
2118  isr2,
2119  osr,
2120  threshold);
2121  if (poly == 0 || min_distance > distance) {
2122  min_distance = distance;
2123  if (tol_zero(min_distance)) {
2124  min_distance = 0.0;
2125  break;
2126  }
2127  if (min_distance <= threshold) {
2128  return min_distance;
2129  }
2130  }
2131  }
2132 
2133  return min_distance;
2134 }
2135 
2137 double ST_Distance_Polygon_Point(int8_t* poly_coords,
2138  int64_t poly_coords_size,
2139  int32_t* poly_ring_sizes,
2140  int64_t poly_num_rings,
2141  int8_t* p,
2142  int64_t psize,
2143  int32_t ic1,
2144  int32_t isr1,
2145  int32_t ic2,
2146  int32_t isr2,
2147  int32_t osr,
2148  double threshold) {
2149  return ST_Distance_Point_Polygon(p,
2150  psize,
2151  poly_coords,
2152  poly_coords_size,
2153  poly_ring_sizes,
2154  poly_num_rings,
2155  ic2,
2156  isr2,
2157  ic1,
2158  isr1,
2159  osr,
2160  threshold);
2161 }
2162 
2164 double ST_Distance_Polygon_LineString(int8_t* poly_coords,
2165  int64_t poly_coords_size,
2166  int32_t* poly_ring_sizes,
2167  int64_t poly_num_rings,
2168  int8_t* l,
2169  int64_t lsize,
2170  int32_t ic1,
2171  int32_t isr1,
2172  int32_t ic2,
2173  int32_t isr2,
2174  int32_t osr,
2175  double threshold) {
2177  lsize,
2178  poly_coords,
2179  poly_coords_size,
2180  poly_ring_sizes,
2181  poly_num_rings,
2182  ic2,
2183  isr2,
2184  ic1,
2185  isr2,
2186  osr,
2187  threshold);
2188 }
2189 
2191 double ST_Distance_Polygon_Polygon(int8_t* poly1_coords,
2192  int64_t poly1_coords_size,
2193  int32_t* poly1_ring_sizes,
2194  int64_t poly1_num_rings,
2195  int8_t* poly2_coords,
2196  int64_t poly2_coords_size,
2197  int32_t* poly2_ring_sizes,
2198  int64_t poly2_num_rings,
2199  int32_t ic1,
2200  int32_t isr1,
2201  int32_t ic2,
2202  int32_t isr2,
2203  int32_t osr,
2204  double threshold) {
2205  // Check if poly1 contains the first point of poly2's shape, i.e. the external ring
2206  auto poly2_first_point_coords = poly2_coords;
2207  auto poly2_first_point_coords_size = compression_unit_size(ic2) * 2;
2208  auto min_distance = ST_Distance_Polygon_Point(poly1_coords,
2209  poly1_coords_size,
2210  poly1_ring_sizes,
2211  poly1_num_rings,
2212  poly2_first_point_coords,
2213  poly2_first_point_coords_size,
2214  ic1,
2215  isr1,
2216  ic2,
2217  isr2,
2218  osr,
2219  threshold);
2220  if (tol_zero(min_distance)) {
2221  // Polygons overlap
2222  return 0.0;
2223  }
2224  if (min_distance <= threshold) {
2225  return min_distance;
2226  }
2227 
2228  // Poly2's first point is either outside poly1's external ring or inside one of the
2229  // internal rings. Measure the smallest distance between a poly1 ring (external or
2230  // internal) and a poly2 ring (external or internal). If poly2 is completely outside
2231  // poly1, then the min distance would be between poly1's and poly2's external rings. If
2232  // poly2 is completely inside one of poly1 internal rings then the min distance would be
2233  // between that poly1 internal ring and poly2's external ring. If poly1 is completely
2234  // inside one of poly2 internal rings, min distance is between that internal ring and
2235  // poly1's external ring. In each case other rings don't get in the way. Any ring
2236  // intersection means zero distance - short-circuit and return.
2237 
2238  auto poly1_ring_coords = poly1_coords;
2239  for (auto r1 = 0; r1 < poly1_num_rings; r1++) {
2240  int64_t poly1_ring_num_coords = poly1_ring_sizes[r1] * 2;
2241 
2242  auto poly2_ring_coords = poly2_coords;
2243  for (auto r2 = 0; r2 < poly2_num_rings; r2++) {
2244  int64_t poly2_ring_num_coords = poly2_ring_sizes[r2] * 2;
2245 
2246  auto distance = distance_ring_ring(poly1_ring_coords,
2247  poly1_ring_num_coords,
2248  poly2_ring_coords,
2249  poly2_ring_num_coords,
2250  ic1,
2251  isr1,
2252  ic2,
2253  isr2,
2254  osr,
2255  threshold);
2256  if (min_distance > distance) {
2257  min_distance = distance;
2258  if (tol_zero(min_distance)) {
2259  return 0.0;
2260  }
2261  if (min_distance <= threshold) {
2262  return min_distance;
2263  }
2264  }
2265 
2266  poly2_ring_coords += poly2_ring_num_coords * compression_unit_size(ic2);
2267  }
2268 
2269  poly1_ring_coords += poly1_ring_num_coords * compression_unit_size(ic1);
2270  }
2271 
2272  return min_distance;
2273 }
2274 
2276 double ST_Distance_Polygon_MultiPolygon(int8_t* poly1_coords,
2277  int64_t poly1_coords_size,
2278  int32_t* poly1_ring_sizes,
2279  int64_t poly1_num_rings,
2280  int8_t* mpoly_coords,
2281  int64_t mpoly_coords_size,
2282  int32_t* mpoly_ring_sizes,
2283  int64_t mpoly_num_rings,
2284  int32_t* mpoly_poly_sizes,
2285  int64_t mpoly_num_polys,
2286  int32_t ic1,
2287  int32_t isr1,
2288  int32_t ic2,
2289  int32_t isr2,
2290  int32_t osr,
2291  double threshold) {
2292  double min_distance = 0.0;
2293 
2294  // Set specific poly pointers as we move through the coords/ringsizes/polyrings arrays.
2295  auto next_poly_coords = mpoly_coords;
2296  auto next_poly_ring_sizes = mpoly_ring_sizes;
2297 
2298  for (auto poly = 0; poly < mpoly_num_polys; poly++) {
2299  auto poly_coords = next_poly_coords;
2300  auto poly_ring_sizes = next_poly_ring_sizes;
2301  auto poly_num_rings = mpoly_poly_sizes[poly];
2302  // Count number of coords in all of poly's rings, advance ring size pointer.
2303  int32_t poly_num_coords = 0;
2304  for (auto ring = 0; ring < poly_num_rings; ring++) {
2305  poly_num_coords += 2 * *next_poly_ring_sizes++;
2306  }
2307  auto poly_coords_size = poly_num_coords * compression_unit_size(ic2);
2308  next_poly_coords += poly_coords_size;
2309  double distance = ST_Distance_Polygon_Polygon(poly1_coords,
2310  poly1_coords_size,
2311  poly1_ring_sizes,
2312  poly1_num_rings,
2313  poly_coords,
2314  poly_coords_size,
2315  poly_ring_sizes,
2316  poly_num_rings,
2317  ic1,
2318  isr1,
2319  ic2,
2320  isr2,
2321  osr,
2322  threshold);
2323  if (poly == 0 || min_distance > distance) {
2324  min_distance = distance;
2325  if (tol_zero(min_distance)) {
2326  min_distance = 0.0;
2327  break;
2328  }
2329  if (min_distance <= threshold) {
2330  break;
2331  }
2332  }
2333  }
2334 
2335  return min_distance;
2336 }
2337 
2339 double ST_Distance_MultiPolygon_Point(int8_t* mpoly_coords,
2340  int64_t mpoly_coords_size,
2341  int32_t* mpoly_ring_sizes,
2342  int64_t mpoly_num_rings,
2343  int32_t* mpoly_poly_sizes,
2344  int64_t mpoly_num_polys,
2345  int8_t* p,
2346  int64_t psize,
2347  int32_t ic1,
2348  int32_t isr1,
2349  int32_t ic2,
2350  int32_t isr2,
2351  int32_t osr,
2352  double threshold) {
2354  psize,
2355  mpoly_coords,
2356  mpoly_coords_size,
2357  mpoly_ring_sizes,
2358  mpoly_num_rings,
2359  mpoly_poly_sizes,
2360  mpoly_num_polys,
2361  ic2,
2362  isr2,
2363  ic1,
2364  isr1,
2365  osr,
2366  threshold);
2367 }
2368 
2370 double ST_Distance_MultiPolygon_LineString(int8_t* mpoly_coords,
2371  int64_t mpoly_coords_size,
2372  int32_t* mpoly_ring_sizes,
2373  int64_t mpoly_num_rings,
2374  int32_t* mpoly_poly_sizes,
2375  int64_t mpoly_num_polys,
2376  int8_t* l,
2377  int64_t lsize,
2378  int32_t ic1,
2379  int32_t isr1,
2380  int32_t ic2,
2381  int32_t isr2,
2382  int32_t osr,
2383  double threshold) {
2385  lsize,
2386  mpoly_coords,
2387  mpoly_coords_size,
2388  mpoly_ring_sizes,
2389  mpoly_num_rings,
2390  mpoly_poly_sizes,
2391  mpoly_num_polys,
2392  ic2,
2393  isr2,
2394  ic1,
2395  isr1,
2396  osr,
2397  threshold);
2398 }
2399 
2401 double ST_Distance_MultiPolygon_Polygon(int8_t* mpoly_coords,
2402  int64_t mpoly_coords_size,
2403  int32_t* mpoly_ring_sizes,
2404  int64_t mpoly_num_rings,
2405  int32_t* mpoly_poly_sizes,
2406  int64_t mpoly_num_polys,
2407  int8_t* poly1_coords,
2408  int64_t poly1_coords_size,
2409  int32_t* poly1_ring_sizes,
2410  int64_t poly1_num_rings,
2411  int32_t ic1,
2412  int32_t isr1,
2413  int32_t ic2,
2414  int32_t isr2,
2415  int32_t osr,
2416  double threshold) {
2417  return ST_Distance_Polygon_MultiPolygon(poly1_coords,
2418  poly1_coords_size,
2419  poly1_ring_sizes,
2420  poly1_num_rings,
2421  mpoly_coords,
2422  mpoly_coords_size,
2423  mpoly_ring_sizes,
2424  mpoly_num_rings,
2425  mpoly_poly_sizes,
2426  mpoly_num_polys,
2427  ic2,
2428  isr2,
2429  ic1,
2430  isr1,
2431  osr,
2432  threshold);
2433 }
2434 
2436 double ST_Distance_MultiPolygon_MultiPolygon(int8_t* mpoly1_coords,
2437  int64_t mpoly1_coords_size,
2438  int32_t* mpoly1_ring_sizes,
2439  int64_t mpoly1_num_rings,
2440  int32_t* mpoly1_poly_sizes,
2441  int64_t mpoly1_num_polys,
2442  int8_t* mpoly2_coords,
2443  int64_t mpoly2_coords_size,
2444  int32_t* mpoly2_ring_sizes,
2445  int64_t mpoly2_num_rings,
2446  int32_t* mpoly2_poly_sizes,
2447  int64_t mpoly2_num_polys,
2448  int32_t ic1,
2449  int32_t isr1,
2450  int32_t ic2,
2451  int32_t isr2,
2452  int32_t osr,
2453  double threshold) {
2454  double min_distance = 0.0;
2455 
2456  // Set specific poly pointers as we move through mpoly1's coords/ringsizes/polyrings
2457  // arrays.
2458  auto next_poly_coords = mpoly1_coords;
2459  auto next_poly_ring_sizes = mpoly1_ring_sizes;
2460 
2461  for (auto poly = 0; poly < mpoly1_num_polys; poly++) {
2462  auto poly_coords = next_poly_coords;
2463  auto poly_ring_sizes = next_poly_ring_sizes;
2464  auto poly_num_rings = mpoly1_poly_sizes[poly];
2465  // Count number of coords in all of poly's rings, advance ring size pointer.
2466  int32_t poly_num_coords = 0;
2467  for (auto ring = 0; ring < poly_num_rings; ring++) {
2468  poly_num_coords += 2 * *next_poly_ring_sizes++;
2469  }
2470  auto poly_coords_size = poly_num_coords * compression_unit_size(ic1);
2471  next_poly_coords += poly_coords_size;
2472  double distance = ST_Distance_Polygon_MultiPolygon(poly_coords,
2473  poly_coords_size,
2474  poly_ring_sizes,
2475  poly_num_rings,
2476  mpoly2_coords,
2477  mpoly2_coords_size,
2478  mpoly2_ring_sizes,
2479  mpoly2_num_rings,
2480  mpoly2_poly_sizes,
2481  mpoly2_num_polys,
2482  ic1,
2483  isr1,
2484  ic2,
2485  isr2,
2486  osr,
2487  threshold);
2488  if (poly == 0 || min_distance > distance) {
2489  min_distance = distance;
2490  if (tol_zero(min_distance)) {
2491  min_distance = 0.0;
2492  break;
2493  }
2494  if (min_distance <= threshold) {
2495  break;
2496  }
2497  }
2498  }
2499 
2500  return min_distance;
2501 }
2502 
2503 //
2504 // ST_DWithin
2505 //
2506 
2508 bool ST_DWithin_Point_Point(int8_t* p1,
2509  int64_t p1size,
2510  int8_t* p2,
2511  int64_t p2size,
2512  int32_t ic1,
2513  int32_t isr1,
2514  int32_t ic2,
2515  int32_t isr2,
2516  int32_t osr,
2517  double distance_within) {
2519  p1, p1size, p2, p2size, ic1, isr1, ic2, isr2, osr) <=
2520  distance_within * distance_within;
2521 }
2522 
2525  int64_t p1size,
2526  int8_t* p2,
2527  int64_t p2size,
2528  int32_t ic1,
2529  int32_t isr1,
2530  int32_t ic2,
2531  int32_t isr2,
2532  int32_t osr,
2533  double distance_within) {
2534  auto dist_meters =
2535  ST_Distance_Point_Point_Geodesic(p1, p1size, p2, p2size, ic1, isr1, ic2, isr2, osr);
2536  return (dist_meters <= distance_within);
2537 }
2538 
2541  int64_t p1size,
2542  int8_t* l2,
2543  int64_t l2size,
2544  double* l2bounds,
2545  int64_t l2bounds_size,
2546  int32_t ic1,
2547  int32_t isr1,
2548  int32_t ic2,
2549  int32_t isr2,
2550  int32_t osr,
2551  double distance_within) {
2552  if (l2bounds) {
2553  // Bounding boxes need to be transformed to output SR before proximity check
2554  if (!point_dwithin_box(
2555  p1, p1size, ic1, isr1, l2bounds, l2bounds_size, isr2, osr, distance_within)) {
2556  return false;
2557  }
2558  }
2559 
2560  // May need to adjust the threshold by TOLERANCE_DEFAULT
2561  const double threshold = distance_within;
2563  p1, p1size, l2, l2size, ic1, isr1, ic2, isr2, osr, threshold) <=
2564  distance_within;
2565 }
2566 
2569  int64_t psize,
2570  int8_t* poly_coords,
2571  int64_t poly_coords_size,
2572  int32_t* poly_ring_sizes,
2573  int64_t poly_num_rings,
2574  double* poly_bounds,
2575  int64_t poly_bounds_size,
2576  int32_t ic1,
2577  int32_t isr1,
2578  int32_t ic2,
2579  int32_t isr2,
2580  int32_t osr,
2581  double distance_within) {
2582  if (poly_bounds) {
2583  if (!point_dwithin_box(p,
2584  psize,
2585  ic1,
2586  isr1,
2587  poly_bounds,
2588  poly_bounds_size,
2589  isr2,
2590  osr,
2591  distance_within)) {
2592  return false;
2593  }
2594  }
2595 
2596  // May need to adjust the threshold by TOLERANCE_DEFAULT
2597  const double threshold = distance_within;
2598  return ST_Distance_Point_Polygon(p,
2599  psize,
2600  poly_coords,
2601  poly_coords_size,
2602  poly_ring_sizes,
2603  poly_num_rings,
2604  ic1,
2605  isr1,
2606  ic2,
2607  isr2,
2608  osr,
2609  threshold) <= distance_within;
2610 }
2611 
2614  int64_t psize,
2615  int8_t* mpoly_coords,
2616  int64_t mpoly_coords_size,
2617  int32_t* mpoly_ring_sizes,
2618  int64_t mpoly_num_rings,
2619  int32_t* mpoly_poly_sizes,
2620  int64_t mpoly_num_polys,
2621  double* mpoly_bounds,
2622  int64_t mpoly_bounds_size,
2623  int32_t ic1,
2624  int32_t isr1,
2625  int32_t ic2,
2626  int32_t isr2,
2627  int32_t osr,
2628  double distance_within) {
2629  if (mpoly_bounds) {
2630  if (!point_dwithin_box(p,
2631  psize,
2632  ic1,
2633  isr1,
2634  mpoly_bounds,
2635  mpoly_bounds_size,
2636  isr2,
2637  osr,
2638  distance_within)) {
2639  return false;
2640  }
2641  }
2642 
2643  // May need to adjust the threshold by TOLERANCE_DEFAULT
2644  const double threshold = distance_within;
2646  psize,
2647  mpoly_coords,
2648  mpoly_coords_size,
2649  mpoly_ring_sizes,
2650  mpoly_num_rings,
2651  mpoly_poly_sizes,
2652  mpoly_num_polys,
2653  ic1,
2654  isr1,
2655  ic2,
2656  isr2,
2657  osr,
2658  threshold) <= distance_within;
2659 }
2660 
2663  int64_t l1size,
2664  double* l1bounds,
2665  int64_t l1bounds_size,
2666  int8_t* l2,
2667  int64_t l2size,
2668  double* l2bounds,
2669  int64_t l2bounds_size,
2670  int32_t ic1,
2671  int32_t isr1,
2672  int32_t ic2,
2673  int32_t isr2,
2674  int32_t osr,
2675  double distance_within) {
2676  if (l1bounds && l2bounds) {
2677  // Bounding boxes need to be transformed to output SR before proximity check
2678  if (!box_dwithin_box(l1bounds,
2679  l1bounds_size,
2680  isr1,
2681  l2bounds,
2682  l2bounds_size,
2683  isr2,
2684  osr,
2685  distance_within)) {
2686  return false;
2687  }
2688  }
2689 
2690  // May need to adjust the threshold by TOLERANCE_DEFAULT
2691  const double threshold = distance_within;
2693  l1, l1size, l2, l2size, ic1, isr1, ic2, isr2, osr, threshold) <=
2694  distance_within;
2695 }
2696 
2699  int64_t l1size,
2700  double* l1bounds,
2701  int64_t l1bounds_size,
2702  int8_t* poly_coords,
2703  int64_t poly_coords_size,
2704  int32_t* poly_ring_sizes,
2705  int64_t poly_num_rings,
2706  double* poly_bounds,
2707  int64_t poly_bounds_size,
2708  int32_t ic1,
2709  int32_t isr1,
2710  int32_t ic2,
2711  int32_t isr2,
2712  int32_t osr,
2713  double distance_within) {
2714  if (l1bounds && poly_bounds) {
2715  // Bounding boxes need to be transformed to output SR before proximity check
2716  if (!box_dwithin_box(l1bounds,
2717  l1bounds_size,
2718  isr1,
2719  poly_bounds,
2720  poly_bounds_size,
2721  isr2,
2722  osr,
2723  distance_within)) {
2724  return false;
2725  }
2726  }
2727 
2728  // May need to adjust the threshold by TOLERANCE_DEFAULT
2729  const double threshold = distance_within;
2731  l1size,
2732  poly_coords,
2733  poly_coords_size,
2734  poly_ring_sizes,
2735  poly_num_rings,
2736  ic1,
2737  isr1,
2738  ic2,
2739  isr2,
2740  osr,
2741  threshold) <= distance_within;
2742 }
2743 
2746  int64_t l1size,
2747  double* l1bounds,
2748  int64_t l1bounds_size,
2749  int8_t* mpoly_coords,
2750  int64_t mpoly_coords_size,
2751  int32_t* mpoly_ring_sizes,
2752  int64_t mpoly_num_rings,
2753  int32_t* mpoly_poly_sizes,
2754  int64_t mpoly_num_polys,
2755  double* mpoly_bounds,
2756  int64_t mpoly_bounds_size,
2757  int32_t ic1,
2758  int32_t isr1,
2759  int32_t ic2,
2760  int32_t isr2,
2761  int32_t osr,
2762  double distance_within) {
2763  if (l1bounds && mpoly_bounds) {
2764  // Bounding boxes need to be transformed to output SR before proximity check
2765  if (!box_dwithin_box(l1bounds,
2766  l1bounds_size,
2767  isr1,
2768  mpoly_bounds,
2769  mpoly_bounds_size,
2770  isr2,
2771  osr,
2772  distance_within)) {
2773  return false;
2774  }
2775  }
2776 
2777  // May need to adjust the threshold by TOLERANCE_DEFAULT
2778  const double threshold = distance_within;
2780  l1size,
2781  mpoly_coords,
2782  mpoly_coords_size,
2783  mpoly_ring_sizes,
2784  mpoly_num_rings,
2785  mpoly_poly_sizes,
2786  mpoly_num_polys,
2787  ic1,
2788  isr1,
2789  ic2,
2790  isr2,
2791  osr,
2792  threshold) <= distance_within;
2793 }
2794 
2796 bool ST_DWithin_Polygon_Polygon(int8_t* poly1_coords,
2797  int64_t poly1_coords_size,
2798  int32_t* poly1_ring_sizes,
2799  int64_t poly1_num_rings,
2800  double* poly1_bounds,
2801  int64_t poly1_bounds_size,
2802  int8_t* poly2_coords,
2803  int64_t poly2_coords_size,
2804  int32_t* poly2_ring_sizes,
2805  int64_t poly2_num_rings,
2806  double* poly2_bounds,
2807  int64_t poly2_bounds_size,
2808  int32_t ic1,
2809  int32_t isr1,
2810  int32_t ic2,
2811  int32_t isr2,
2812  int32_t osr,
2813  double distance_within) {
2814  if (poly1_bounds && poly2_bounds) {
2815  // Bounding boxes need to be transformed to output SR before proximity check
2816  if (!box_dwithin_box(poly1_bounds,
2817  poly1_bounds_size,
2818  isr1,
2819  poly2_bounds,
2820  poly2_bounds_size,
2821  isr2,
2822  osr,
2823  distance_within)) {
2824  return false;
2825  }
2826  }
2827 
2828  // May need to adjust the threshold by TOLERANCE_DEFAULT
2829  const double threshold = distance_within;
2830  return ST_Distance_Polygon_Polygon(poly1_coords,
2831  poly1_coords_size,
2832  poly1_ring_sizes,
2833  poly1_num_rings,
2834  poly2_coords,
2835  poly2_coords_size,
2836  poly2_ring_sizes,
2837  poly2_num_rings,
2838  ic1,
2839  isr1,
2840  ic2,
2841  isr2,
2842  osr,
2843  threshold) <= distance_within;
2844 }
2845 
2847 bool ST_DWithin_Polygon_MultiPolygon(int8_t* poly_coords,
2848  int64_t poly_coords_size,
2849  int32_t* poly_ring_sizes,
2850  int64_t poly_num_rings,
2851  double* poly_bounds,
2852  int64_t poly_bounds_size,
2853  int8_t* mpoly_coords,
2854  int64_t mpoly_coords_size,
2855  int32_t* mpoly_ring_sizes,
2856  int64_t mpoly_num_rings,
2857  int32_t* mpoly_poly_sizes,
2858  int64_t mpoly_num_polys,
2859  double* mpoly_bounds,
2860  int64_t mpoly_bounds_size,
2861  int32_t ic1,
2862  int32_t isr1,
2863  int32_t ic2,
2864  int32_t isr2,
2865  int32_t osr,
2866  double distance_within) {
2867  if (poly_bounds && mpoly_bounds) {
2868  // Bounding boxes need to be transformed to output SR before proximity check
2869  if (!box_dwithin_box(poly_bounds,
2870  poly_bounds_size,
2871  isr1,
2872  mpoly_bounds,
2873  mpoly_bounds_size,
2874  isr2,
2875  osr,
2876  distance_within)) {
2877  return false;
2878  }
2879  }
2880 
2881  // May need to adjust the threshold by TOLERANCE_DEFAULT
2882  const double threshold = distance_within;
2883  return ST_Distance_Polygon_MultiPolygon(poly_coords,
2884  poly_coords_size,
2885  poly_ring_sizes,
2886  poly_num_rings,
2887  mpoly_coords,
2888  mpoly_coords_size,
2889  mpoly_ring_sizes,
2890  mpoly_num_rings,
2891  mpoly_poly_sizes,
2892  mpoly_num_polys,
2893  ic1,
2894  isr1,
2895  ic2,
2896  isr2,
2897  osr,
2898  threshold) <= distance_within;
2899 }
2900 
2902 bool ST_DWithin_MultiPolygon_MultiPolygon(int8_t* mpoly1_coords,
2903  int64_t mpoly1_coords_size,
2904  int32_t* mpoly1_ring_sizes,
2905  int64_t mpoly1_num_rings,
2906  int32_t* mpoly1_poly_sizes,
2907  int64_t mpoly1_num_polys,
2908  double* mpoly1_bounds,
2909  int64_t mpoly1_bounds_size,
2910  int8_t* mpoly2_coords,
2911  int64_t mpoly2_coords_size,
2912  int32_t* mpoly2_ring_sizes,
2913  int64_t mpoly2_num_rings,
2914  int32_t* mpoly2_poly_sizes,
2915  int64_t mpoly2_num_polys,
2916  double* mpoly2_bounds,
2917  int64_t mpoly2_bounds_size,
2918  int32_t ic1,
2919  int32_t isr1,
2920  int32_t ic2,
2921  int32_t isr2,
2922  int32_t osr,
2923  double distance_within) {
2924  if (mpoly1_bounds && mpoly2_bounds) {
2925  // Bounding boxes need to be transformed to output SR before proximity check
2926  if (!box_dwithin_box(mpoly1_bounds,
2927  mpoly1_bounds_size,
2928  isr1,
2929  mpoly2_bounds,
2930  mpoly2_bounds_size,
2931  isr2,
2932  osr,
2933  distance_within)) {
2934  return false;
2935  }
2936  }
2937 
2938  // May need to adjust the threshold by TOLERANCE_DEFAULT
2939  const double threshold = distance_within;
2940  return ST_Distance_MultiPolygon_MultiPolygon(mpoly1_coords,
2941  mpoly1_coords_size,
2942  mpoly1_ring_sizes,
2943  mpoly1_num_rings,
2944  mpoly1_poly_sizes,
2945  mpoly1_num_polys,
2946  mpoly2_coords,
2947  mpoly2_coords_size,
2948  mpoly2_ring_sizes,
2949  mpoly2_num_rings,
2950  mpoly2_poly_sizes,
2951  mpoly2_num_polys,
2952  ic1,
2953  isr1,
2954  ic2,
2955  isr2,
2956  osr,
2957  threshold) <= distance_within;
2958 }
2959 
2960 //
2961 // ST_MaxDistance
2962 //
2963 
2964 // Max cartesian distance between a point and a line segment
2965 DEVICE
2966 double max_distance_point_line(double px,
2967  double py,
2968  double l1x,
2969  double l1y,
2970  double l2x,
2971  double l2y) {
2972  // TODO: switch to squared distances
2973  double length1 = distance_point_point(px, py, l1x, l1y);
2974  double length2 = distance_point_point(px, py, l2x, l2y);
2975  if (length1 > length2) {
2976  return length1;
2977  }
2978  return length2;
2979 }
2980 
2982  int64_t psize,
2983  int8_t* l,
2984  int64_t lsize,
2985  int32_t ic1,
2986  int32_t isr1,
2987  int32_t ic2,
2988  int32_t isr2,
2989  int32_t osr,
2990  bool check_closed) {
2991  // TODO: switch to squared distances
2992  double px = coord_x(p, 0, ic1, isr1, osr);
2993  double py = coord_y(p, 1, ic1, isr1, osr);
2994 
2995  auto l_num_coords = lsize / compression_unit_size(ic2);
2996 
2997  double l1x = coord_x(l, 0, ic2, isr2, osr);
2998  double l1y = coord_y(l, 1, ic2, isr2, osr);
2999  double l2x = coord_x(l, 2, ic2, isr2, osr);
3000  double l2y = coord_y(l, 3, ic2, isr2, osr);
3001 
3002  double max_dist = max_distance_point_line(px, py, l1x, l1y, l2x, l2y);
3003  for (int32_t i = 4; i < l_num_coords; i += 2) {
3004  l1x = l2x; // advance one point
3005  l1y = l2y;
3006  l2x = coord_x(l, i, ic2, isr2, osr);
3007  l2y = coord_y(l, i + 1, ic2, isr2, osr);
3008  double ldist = max_distance_point_line(px, py, l1x, l1y, l2x, l2y);
3009  if (max_dist < ldist) {
3010  max_dist = ldist;
3011  }
3012  }
3013  if (l_num_coords > 4 && check_closed) {
3014  // Also check distance to the closing edge between the first and the last points
3015  l1x = coord_x(l, 0, ic2, isr2, osr);
3016  l1y = coord_y(l, 1, ic2, isr2, osr);
3017  double ldist = max_distance_point_line(px, py, l1x, l1y, l2x, l2y);
3018  if (max_dist < ldist) {
3019  max_dist = ldist;
3020  }
3021  }
3022  return max_dist;
3023 }
3024 
3027  int64_t psize,
3028  int8_t* l,
3029  int64_t lsize,
3030  int32_t ic1,
3031  int32_t isr1,
3032  int32_t ic2,
3033  int32_t isr2,
3034  int32_t osr) {
3036  p, psize, l, lsize, ic1, isr1, ic2, isr2, osr, false);
3037 }
3038 
3041  int64_t lsize,
3042  int8_t* p,
3043  int64_t psize,
3044  int32_t ic1,
3045  int32_t isr1,
3046  int32_t ic2,
3047  int32_t isr2,
3048  int32_t osr) {
3050  p, psize, l, lsize, ic2, isr2, ic1, isr1, osr, false);
3051 }
3052 
3053 // TODO: add ST_MaxDistance_LineString_LineString (with short-circuit threshold)
3054 
3055 //
3056 // ST_Contains
3057 //
3058 
3060 bool ST_Contains_Point_Point(int8_t* p1,
3061  int64_t p1size,
3062  int8_t* p2,
3063  int64_t p2size,
3064  int32_t ic1,
3065  int32_t isr1,
3066  int32_t ic2,
3067  int32_t isr2,
3068  int32_t osr) {
3069  double p1x = coord_x(p1, 0, ic1, isr1, osr);
3070  double p1y = coord_y(p1, 1, ic1, isr1, osr);
3071  double p2x = coord_x(p2, 0, ic2, isr2, osr);
3072  double p2y = coord_y(p2, 1, ic2, isr2, osr);
3073  double tolerance = tol(ic1, ic2);
3074  return tol_eq(p1x, p2x, tolerance) && tol_eq(p1y, p2y, tolerance);
3075 }
3076 
3079  int64_t psize,
3080  int8_t* l,
3081  int64_t lsize,
3082  double* lbounds,
3083  int64_t lbounds_size,
3084  int32_t ic1,
3085  int32_t isr1,
3086  int32_t ic2,
3087  int32_t isr2,
3088  int32_t osr) {
3089  double px = coord_x(p, 0, ic1, isr1, osr);
3090  double py = coord_y(p, 1, ic1, isr1, osr);
3091 
3092  if (lbounds) {
3093  if (tol_eq(px, lbounds[0]) && tol_eq(py, lbounds[1]) && tol_eq(px, lbounds[2]) &&
3094  tol_eq(py, lbounds[3])) {
3095  return true;
3096  }
3097  }
3098 
3099  auto l_num_coords = lsize / compression_unit_size(ic2);
3100  for (int i = 0; i < l_num_coords; i += 2) {
3101  double lx = coord_x(l, i, ic2, isr2, osr);
3102  double ly = coord_y(l, i + 1, ic2, isr2, osr);
3103  if (tol_eq(px, lx) && tol_eq(py, ly)) {
3104  continue;
3105  }
3106  return false;
3107  }
3108  return true;
3109 }
3110 
3113  int64_t psize,
3114  int8_t* poly_coords,
3115  int64_t poly_coords_size,
3116  int32_t* poly_ring_sizes,
3117  int64_t poly_num_rings,
3118  double* poly_bounds,
3119  int64_t poly_bounds_size,
3120  int32_t ic1,
3121  int32_t isr1,
3122  int32_t ic2,
3123  int32_t isr2,
3124  int32_t osr) {
3125  auto exterior_ring_num_coords = poly_coords_size / compression_unit_size(ic2);
3126  if (poly_num_rings > 0) {
3127  exterior_ring_num_coords = poly_ring_sizes[0] * 2;
3128  }
3129  auto exterior_ring_coords_size = exterior_ring_num_coords * compression_unit_size(ic2);
3130 
3132  psize,
3133  poly_coords,
3134  exterior_ring_coords_size,
3135  poly_bounds,
3136  poly_bounds_size,
3137  ic1,
3138  isr1,
3139  ic2,
3140  isr2,
3141  osr);
3142 }
3143 
3146  int64_t lsize,
3147  double* lbounds,
3148  int64_t lbounds_size,
3149  int8_t* p,
3150  int64_t psize,
3151  int32_t ic1,
3152  int32_t isr1,
3153  int32_t ic2,
3154  int32_t isr2,
3155  int32_t osr) {
3156  return tol_zero(
3157  ST_Distance_Point_LineString(p, psize, l, lsize, ic2, isr2, ic1, isr1, osr, 0.0));
3158 }
3159 
3162  int64_t lsize,
3163  double* lbounds,
3164  int64_t lbounds_size,
3165  int8_t* poly_coords,
3166  int64_t poly_coords_size,
3167  int32_t* poly_ring_sizes,
3168  int64_t poly_num_rings,
3169  double* poly_bounds,
3170  int64_t poly_bounds_size,
3171  int32_t ic1,
3172  int32_t isr1,
3173  int32_t ic2,
3174  int32_t isr2,
3175  int32_t osr) {
3176  // TODO
3177  return false;
3178 }
3179 
3180 template <typename T, EdgeBehavior TEdgeBehavior>
3181 DEVICE ALWAYS_INLINE bool Contains_Polygon_Point_Impl(const int8_t* poly_coords,
3182  const int64_t poly_coords_size,
3183  const int32_t* poly_ring_sizes,
3184  const int64_t poly_num_rings,
3185  const double* poly_bounds,
3186  const int64_t poly_bounds_size,
3187  const int8_t* p,
3188  const int64_t psize,
3189  const int32_t ic1,
3190  const int32_t isr1,
3191  const int32_t ic2,
3192  const int32_t isr2,
3193  const int32_t osr) {
3194  if (poly_bounds) {
3195  // TODO: codegen
3196  if (!box_contains_point(poly_bounds,
3197  poly_bounds_size,
3198  coord_x(p, 0, ic2, isr2, osr),
3199  coord_y(p, 1, ic2, isr2, osr))) {
3200  DEBUG_STMT(printf("Bounding box does not contain point, exiting.\n"));
3201  return false;
3202  }
3203  }
3204 
3205  auto get_x_coord = [=]() -> T {
3206  if constexpr (std::is_floating_point<T>::value) {
3207  return coord_x(p, 0, ic2, isr2, osr);
3208  } else {
3209  return compressed_coord(p, 0);
3210  }
3211  return T{}; // https://stackoverflow.com/a/64561686/2700898
3212  };
3213 
3214  auto get_y_coord = [=]() -> T {
3215  if constexpr (std::is_floating_point<T>::value) {
3216  return coord_y(p, 1, ic2, isr2, osr);
3217  } else {
3218  return compressed_coord(p, 1);
3219  }
3220  return T{}; // https://stackoverflow.com/a/64561686/2700898
3221  };
3222 
3223  const T px = get_x_coord();
3224  const T py = get_y_coord();
3225  DEBUG_STMT(printf("pt: %f, %f\n", (double)px, (double)py));
3226 
3227  const auto poly_num_coords = poly_coords_size / compression_unit_size(ic1);
3228  auto exterior_ring_num_coords =
3229  poly_num_rings > 0 ? poly_ring_sizes[0] * 2 : poly_num_coords;
3230 
3231  auto poly = poly_coords;
3232  if (point_in_polygon_winding_number<T, TEdgeBehavior>(
3233  poly, exterior_ring_num_coords, px, py, ic1, isr1, osr)) {
3234  // Inside exterior ring
3235  poly += exterior_ring_num_coords * compression_unit_size(ic1);
3236  // Check that none of the polygon's holes contain that point
3237  for (auto r = 1; r < poly_num_rings; r++) {
3238  const int64_t interior_ring_num_coords = poly_ring_sizes[r] * 2;
3239  if (point_in_polygon_winding_number<T, TEdgeBehavior>(
3240  poly, interior_ring_num_coords, px, py, ic1, isr1, osr)) {
3241  return false;
3242  }
3243  poly += interior_ring_num_coords * compression_unit_size(ic1);
3244  }
3245  return true;
3246  }
3247  return false;
3248 }
3249 
3250 EXTENSION_NOINLINE bool ST_Contains_Polygon_Point(const int8_t* poly_coords,
3251  const int64_t poly_coords_size,
3252  const int32_t* poly_ring_sizes,
3253  const int64_t poly_num_rings,
3254  const double* poly_bounds,
3255  const int64_t poly_bounds_size,
3256  const int8_t* p,
3257  const int64_t psize,
3258  const int32_t ic1,
3259  const int32_t isr1,
3260  const int32_t ic2,
3261  const int32_t isr2,
3262  const int32_t osr) {
3263  return Contains_Polygon_Point_Impl<double, EdgeBehavior::kExcludePointOnEdge>(
3264  poly_coords,
3265  poly_coords_size,
3266  poly_ring_sizes,
3267  poly_num_rings,
3268  poly_bounds,
3269  poly_bounds_size,
3270  p,
3271  psize,
3272  ic1,
3273  isr1,
3274  ic2,
3275  isr2,
3276  osr);
3277 }
3278 
3279 EXTENSION_NOINLINE bool ST_cContains_Polygon_Point(const int8_t* poly_coords,
3280  const int64_t poly_coords_size,
3281  const int32_t* poly_ring_sizes,
3282  const int64_t poly_num_rings,
3283  const double* poly_bounds,
3284  const int64_t poly_bounds_size,
3285  const int8_t* p,
3286  const int64_t psize,
3287  const int32_t ic1,
3288  const int32_t isr1,
3289  const int32_t ic2,
3290  const int32_t isr2,
3291  const int32_t osr) {
3292  return Contains_Polygon_Point_Impl<int64_t, EdgeBehavior::kExcludePointOnEdge>(
3293  poly_coords,
3294  poly_coords_size,
3295  poly_ring_sizes,
3296  poly_num_rings,
3297  poly_bounds,
3298  poly_bounds_size,
3299  p,
3300  psize,
3301  ic1,
3302  isr1,
3303  ic2,
3304  isr2,
3305  osr);
3306 }
3307 
3309 bool ST_Contains_Polygon_LineString(int8_t* poly_coords,
3310  int64_t poly_coords_size,
3311  int32_t* poly_ring_sizes,
3312  int64_t poly_num_rings,
3313  double* poly_bounds,
3314  int64_t poly_bounds_size,
3315  int8_t* l,
3316  int64_t lsize,
3317  double* lbounds,
3318  int64_t lbounds_size,
3319  int32_t ic1,
3320  int32_t isr1,
3321  int32_t ic2,
3322  int32_t isr2,
3323  int32_t osr) {
3324  if (poly_num_rings > 1) {
3325  return false; // TODO: support polygons with interior rings
3326  }
3327 
3328  // Bail out if poly bounding box doesn't contain linestring bounding box
3329  if (poly_bounds && lbounds) {
3330  if (!box_contains_box(poly_bounds, poly_bounds_size, lbounds, lbounds_size)) {
3331  return false;
3332  }
3333  }
3334 
3335  const auto poly_num_coords = poly_coords_size / compression_unit_size(ic1);
3336  const auto lnum_coords = lsize / compression_unit_size(ic2);
3337 
3339  poly_coords, poly_num_coords, l, lnum_coords, ic1, isr1, ic2, isr2, osr);
3340 }
3341 
3343 bool ST_Contains_Polygon_Polygon(int8_t* poly1_coords,
3344  int64_t poly1_coords_size,
3345  int32_t* poly1_ring_sizes,
3346  int64_t poly1_num_rings,
3347  double* poly1_bounds,
3348  int64_t poly1_bounds_size,
3349  int8_t* poly2_coords,
3350  int64_t poly2_coords_size,
3351  int32_t* poly2_ring_sizes,
3352  int64_t poly2_num_rings,
3353  double* poly2_bounds,
3354  int64_t poly2_bounds_size,
3355  int32_t ic1,
3356  int32_t isr1,
3357  int32_t ic2,
3358  int32_t isr2,
3359  int32_t osr) {
3360  // TODO: needs to be extended, cover more cases
3361  // Right now only checking if simple poly1 (no holes) contains poly2's exterior shape
3362  if (poly1_num_rings > 1) {
3363  return false; // TODO: support polygons with interior rings
3364  }
3365 
3366  if (poly1_bounds && poly2_bounds) {
3367  if (!box_contains_box(
3368  poly1_bounds, poly1_bounds_size, poly2_bounds, poly2_bounds_size)) {
3369  return false;
3370  }
3371  }
3372 
3373  int64_t poly2_exterior_ring_coords_size = poly2_coords_size;
3374  if (poly2_num_rings > 0) {
3375  poly2_exterior_ring_coords_size =
3376  2 * poly2_ring_sizes[0] * compression_unit_size(ic2);
3377  }
3378  return ST_Contains_Polygon_LineString(poly1_coords,
3379  poly1_coords_size,
3380  poly1_ring_sizes,
3381  poly1_num_rings,
3382  poly1_bounds,
3383  poly1_bounds_size,
3384  poly2_coords,
3385  poly2_exterior_ring_coords_size,
3386  poly2_bounds,
3387  poly2_bounds_size,
3388  ic1,
3389  isr1,
3390  ic2,
3391  isr2,
3392  osr);
3393 }
3394 
3395 template <typename T, EdgeBehavior TEdgeBehavior>
3397  const int8_t* mpoly_coords,
3398  const int64_t mpoly_coords_size,
3399  const int32_t* mpoly_ring_sizes,
3400  const int64_t mpoly_num_rings,
3401  const int32_t* mpoly_poly_sizes,
3402  const int64_t mpoly_num_polys,
3403  const double* mpoly_bounds,
3404  const int64_t mpoly_bounds_size,
3405  const int8_t* p,
3406  const int64_t psize,
3407  const int32_t ic1,
3408  const int32_t isr1,
3409  const int32_t ic2,
3410  const int32_t isr2,
3411  const int32_t osr) {
3412  if (mpoly_num_polys <= 0) {
3413  return false;
3414  }
3415 
3416  // TODO: mpoly_bounds could contain individual bounding boxes too:
3417  // first two points - box for the entire multipolygon, then a pair for each polygon
3418  if (mpoly_bounds) {
3419  if (!box_contains_point(mpoly_bounds,
3420  mpoly_bounds_size,
3421  coord_x(p, 0, ic2, isr2, osr),
3422  coord_y(p, 1, ic2, isr2, osr))) {
3423  return false;
3424  }
3425  }
3426 
3427  // Set specific poly pointers as we move through the coords/ringsizes/polyrings
3428  // arrays.
3429  auto next_poly_coords = mpoly_coords;
3430  auto next_poly_ring_sizes = mpoly_ring_sizes;
3431 
3432  for (auto poly = 0; poly < mpoly_num_polys; poly++) {
3433  const auto poly_coords = next_poly_coords;
3434  const auto poly_ring_sizes = next_poly_ring_sizes;
3435  const auto poly_num_rings = mpoly_poly_sizes[poly];
3436  // Count number of coords in all of poly's rings, advance ring size pointer.
3437  int32_t poly_num_coords = 0;
3438  for (auto ring = 0; ring < poly_num_rings; ring++) {
3439  poly_num_coords += 2 * *next_poly_ring_sizes++;
3440  }
3441  const auto poly_coords_size = poly_num_coords * compression_unit_size(ic1);
3442  next_poly_coords += poly_coords_size;
3443  // TODO: pass individual bounding boxes for each polygon
3444  if (Contains_Polygon_Point_Impl<T, TEdgeBehavior>(poly_coords,
3445  poly_coords_size,
3446  poly_ring_sizes,
3447  poly_num_rings,
3448  nullptr,
3449  0,
3450  p,
3451  psize,
3452  ic1,
3453  isr1,
3454  ic2,
3455  isr2,
3456  osr)) {
3457  return true;
3458  }
3459  }
3460 
3461  return false;
3462 }
3463 
3465 bool ST_Contains_MultiPolygon_Point(int8_t* mpoly_coords,
3466  int64_t mpoly_coords_size,
3467  int32_t* mpoly_ring_sizes,
3468  int64_t mpoly_num_rings,
3469  int32_t* mpoly_poly_sizes,
3470  int64_t mpoly_num_polys,
3471  double* mpoly_bounds,
3472  int64_t mpoly_bounds_size,
3473  int8_t* p,
3474  int64_t psize,
3475  int32_t ic1,
3476  int32_t isr1,
3477  int32_t ic2,
3478  int32_t isr2,
3479  int32_t osr) {
3480  return Contains_MultiPolygon_Point_Impl<double, EdgeBehavior::kExcludePointOnEdge>(
3481  mpoly_coords,
3482  mpoly_coords_size,
3483  mpoly_ring_sizes,
3484  mpoly_num_rings,
3485  mpoly_poly_sizes,
3486  mpoly_num_polys,
3487  mpoly_bounds,
3488  mpoly_bounds_size,
3489  p,
3490  psize,
3491  ic1,
3492  isr1,
3493  ic2,
3494  isr2,
3495  osr);
3496 }
3497 
3499  int64_t mpoly_coords_size,
3500  int32_t* mpoly_ring_sizes,
3501  int64_t mpoly_num_rings,
3502  int32_t* mpoly_poly_sizes,
3503  int64_t mpoly_num_polys,
3504  double* mpoly_bounds,
3505  int64_t mpoly_bounds_size,
3506  int8_t* p,
3507  int64_t psize,
3508  int32_t ic1,
3509  int32_t isr1,
3510  int32_t ic2,
3511  int32_t isr2,
3512  int32_t osr) {
3513  return Contains_MultiPolygon_Point_Impl<int64_t, EdgeBehavior::kExcludePointOnEdge>(
3514  mpoly_coords,
3515  mpoly_coords_size,
3516  mpoly_ring_sizes,
3517  mpoly_num_rings,
3518  mpoly_poly_sizes,
3519  mpoly_num_polys,
3520  mpoly_bounds,
3521  mpoly_bounds_size,
3522  p,
3523  psize,
3524  ic1,
3525  isr1,
3526  ic2,
3527  isr2,
3528  osr);
3529 }
3530 
3532 bool ST_Contains_MultiPolygon_LineString(int8_t* mpoly_coords,
3533  int64_t mpoly_coords_size,
3534  int32_t* mpoly_ring_sizes,
3535  int64_t mpoly_num_rings,
3536  int32_t* mpoly_poly_sizes,
3537  int64_t mpoly_num_polys,
3538  double* mpoly_bounds,
3539  int64_t mpoly_bounds_size,
3540  int8_t* l,
3541  int64_t lsize,
3542  double* lbounds,
3543  int64_t lbounds_size,
3544  int32_t ic1,
3545  int32_t isr1,
3546  int32_t ic2,
3547  int32_t isr2,
3548  int32_t osr) {
3549  if (mpoly_num_polys <= 0) {
3550  return false;
3551  }
3552 
3553  if (mpoly_bounds && lbounds) {
3554  if (!box_contains_box(mpoly_bounds, mpoly_bounds_size, lbounds, lbounds_size)) {
3555  return false;
3556  }
3557  }
3558 
3559  // Set specific poly pointers as we move through the coords/ringsizes/polyrings arrays.
3560  auto next_poly_coords = mpoly_coords;
3561  auto next_poly_ring_sizes = mpoly_ring_sizes;
3562 
3563  for (auto poly = 0; poly < mpoly_num_polys; poly++) {
3564  auto poly_coords = next_poly_coords;
3565  auto poly_ring_sizes = next_poly_ring_sizes;
3566  auto poly_num_rings = mpoly_poly_sizes[poly];
3567  // Count number of coords in all of poly's rings, advance ring size pointer.
3568  int32_t poly_num_coords = 0;
3569  for (auto ring = 0; ring < poly_num_rings; ring++) {
3570  poly_num_coords += 2 * *next_poly_ring_sizes++;
3571  }
3572  auto poly_coords_size = poly_num_coords * compression_unit_size(ic1);
3573  next_poly_coords += poly_coords_size;
3574 
3575  if (ST_Contains_Polygon_LineString(poly_coords,
3576  poly_coords_size,
3577  poly_ring_sizes,
3578  poly_num_rings,
3579  nullptr,
3580  0,
3581  l,
3582  lsize,
3583  nullptr,
3584  0,
3585  ic1,
3586  isr1,
3587  ic2,
3588  isr2,
3589  osr)) {
3590  return true;
3591  }
3592  }
3593 
3594  return false;
3595 }
3596 
3597 //
3598 // ST_Intersects
3599 //
3600 
3603  int64_t p1size,
3604  int8_t* p2,
3605  int64_t p2size,
3606  int32_t ic1,
3607  int32_t isr1,
3608  int32_t ic2,
3609  int32_t isr2,
3610  int32_t osr) {
3611  return tol_zero(
3612  ST_Distance_Point_Point(p1, p1size, p2, p2size, ic1, isr1, ic2, isr2, osr));
3613 }
3614 
3617  int64_t psize,
3618  int8_t* l,
3619  int64_t lsize,
3620  double* lbounds,
3621  int64_t lbounds_size,
3622  int32_t ic1,
3623  int32_t isr1,
3624  int32_t ic2,
3625  int32_t isr2,
3626  int32_t osr) {
3627  double px = coord_x(p, 0, ic1, isr1, osr);
3628  double py = coord_y(p, 1, ic1, isr1, osr);
3629 
3630  if (lbounds) {
3631  if (!box_contains_point(lbounds, lbounds_size, px, py)) {
3632  return false;
3633  }
3634  }
3635  return tol_zero(
3636  ST_Distance_Point_LineString(p, psize, l, lsize, ic1, isr1, ic2, isr2, osr, 0.0));
3637 }
3638 
3641  int64_t psize,
3642  int8_t* poly,
3643  int64_t polysize,
3644  int32_t* poly_ring_sizes,
3645  int64_t poly_num_rings,
3646  double* poly_bounds,
3647  int64_t poly_bounds_size,
3648  int32_t ic1,
3649  int32_t isr1,
3650  int32_t ic2,
3651  int32_t isr2,
3652  int32_t osr) {
3653  return Contains_Polygon_Point_Impl<double, EdgeBehavior::kIncludePointOnEdge>(
3654  poly,
3655  polysize,
3656  poly_ring_sizes,
3657  poly_num_rings,
3658  poly_bounds,
3659  poly_bounds_size,
3660  p,
3661  psize,
3662  ic2,
3663  isr2,
3664  ic1,
3665  isr1,
3666  osr);
3667 }
3668 
3671  int64_t psize,
3672  int8_t* mpoly_coords,
3673  int64_t mpoly_coords_size,
3674  int32_t* mpoly_ring_sizes,
3675  int64_t mpoly_num_rings,
3676  int32_t* mpoly_poly_sizes,
3677  int64_t mpoly_num_polys,
3678  double* mpoly_bounds,
3679  int64_t mpoly_bounds_size,
3680  int32_t ic1,
3681  int32_t isr1,
3682  int32_t ic2,
3683  int32_t isr2,
3684  int32_t osr) {
3685  return ST_Contains_MultiPolygon_Point(mpoly_coords,
3686  mpoly_coords_size,
3687  mpoly_ring_sizes,
3688  mpoly_num_rings,
3689  mpoly_poly_sizes,
3690  mpoly_num_polys,
3691  mpoly_bounds,
3692  mpoly_bounds_size,
3693  p,
3694  psize,
3695  ic2,
3696  isr2,
3697  ic1,
3698  isr1,
3699  osr);
3700 }
3701 
3704  int64_t lsize,
3705  double* lbounds,
3706  int64_t lbounds_size,
3707  int8_t* p,
3708  int64_t psize,
3709  int32_t ic1,
3710  int32_t isr1,
3711  int32_t ic2,
3712  int32_t isr2,
3713  int32_t osr) {
3715  p, psize, l, lsize, lbounds, lbounds_size, ic2, isr2, ic1, isr1, osr);
3716 }
3717 
3720  int64_t l1size,
3721  double* l1bounds,
3722  int64_t l1bounds_size,
3723  int8_t* l2,
3724  int64_t l2size,
3725  double* l2bounds,
3726  int64_t l2bounds_size,
3727  int32_t ic1,
3728  int32_t isr1,
3729  int32_t ic2,
3730  int32_t isr2,
3731  int32_t osr) {
3732  if (l1bounds && l2bounds) {
3733  if (!box_overlaps_box(l1bounds, l1bounds_size, l2bounds, l2bounds_size)) {
3734  return false;
3735  }
3736  }
3737 
3739  l1, l1size, l2, l2size, ic1, isr1, ic2, isr2, osr, 0.0));
3740 }
3741 
3744  int64_t lsize,
3745  double* lbounds,
3746  int64_t lbounds_size,
3747  int8_t* poly,
3748  int64_t polysize,
3749  int32_t* poly_ring_sizes,
3750  int64_t poly_num_rings,
3751  double* poly_bounds,
3752  int64_t poly_bounds_size,
3753  int32_t ic1,
3754  int32_t isr1,
3755  int32_t ic2,
3756  int32_t isr2,
3757  int32_t osr) {
3758  if (lbounds && poly_bounds) {
3759  if (!box_overlaps_box(lbounds, lbounds_size, poly_bounds, poly_bounds_size)) {
3760  return false;
3761  }
3762  }
3763 
3764  // Check for spatial intersection.
3765  // One way to do that would be to start with linestring's first point, if it's inside
3766  // the polygon - it means intersection. Otherwise follow the linestring, segment by
3767  // segment, checking for intersections with polygon rings, bail as soon as we cross into
3768  // the polygon.
3769 
3770  // Or, alternatively, just measure the distance:
3772  lsize,
3773  poly,
3774  polysize,
3775  poly_ring_sizes,
3776  poly_num_rings,
3777  ic1,
3778  isr1,
3779  ic2,
3780  isr2,
3781  osr,
3782  0.0));
3783 }
3784 
3787  int64_t lsize,
3788  double* lbounds,
3789  int64_t lbounds_size,
3790  int8_t* mpoly_coords,
3791  int64_t mpoly_coords_size,
3792  int32_t* mpoly_ring_sizes,
3793  int64_t mpoly_num_rings,
3794  int32_t* mpoly_poly_sizes,
3795  int64_t mpoly_num_polys,
3796  double* mpoly_bounds,
3797  int64_t mpoly_bounds_size,
3798  int32_t ic1,
3799  int32_t isr1,
3800  int32_t ic2,
3801  int32_t isr2,
3802  int32_t osr) {
3803  if (lbounds && mpoly_bounds) {
3804  if (!box_overlaps_box(lbounds, lbounds_size, mpoly_bounds, mpoly_bounds_size)) {
3805  return false;
3806  }
3807  }
3808 
3809  // Check for spatial intersection.
3810  // One way to do that would be to start with linestring's first point, if it's inside
3811  // any of the polygons - it means intersection. Otherwise follow the linestring, segment
3812  // by segment, checking for intersections with polygon shapes/holes, bail as soon as we
3813  // cross into a polygon.
3814 
3815  // Or, alternatively, just measure the distance:
3817  lsize,
3818  mpoly_coords,
3819  mpoly_coords_size,
3820  mpoly_ring_sizes,
3821  mpoly_num_rings,
3822  mpoly_poly_sizes,
3823  mpoly_num_polys,
3824  ic1,
3825  isr1,
3826  ic2,
3827  isr2,
3828  osr,
3829  0.0));
3830 }
3831 
3834  int64_t polysize,
3835  int32_t* poly_ring_sizes,
3836  int64_t poly_num_rings,
3837  double* poly_bounds,
3838  int64_t poly_bounds_size,
3839  int8_t* p,
3840  int64_t psize,
3841  int32_t ic1,
3842  int32_t isr1,
3843  int32_t ic2,
3844  int32_t isr2,
3845  int32_t osr) {
3846  return Contains_Polygon_Point_Impl<double, EdgeBehavior::kIncludePointOnEdge>(
3847  poly,
3848  polysize,
3849  poly_ring_sizes,
3850  poly_num_rings,
3851  poly_bounds,
3852  poly_bounds_size,
3853  p,
3854  psize,
3855  ic1,
3856  isr1,
3857  ic2,
3858  isr2,
3859  osr);
3860 }
3861 
3864  int64_t polysize,
3865  int32_t* poly_ring_sizes,
3866  int64_t poly_num_rings,
3867  double* poly_bounds,
3868  int64_t poly_bounds_size,
3869  int8_t* l,
3870  int64_t lsize,
3871  double* lbounds,
3872  int64_t lbounds_size,
3873  int32_t ic1,
3874  int32_t isr1,
3875  int32_t ic2,
3876  int32_t isr2,
3877  int32_t osr) {
3879  lsize,
3880  lbounds,
3881  lbounds_size,
3882  poly,
3883  polysize,
3884  poly_ring_sizes,
3885  poly_num_rings,
3886  poly_bounds,
3887  poly_bounds_size,
3888  ic2,
3889  isr2,
3890  ic1,
3891  isr1,
3892  osr);
3893 }
3894 
3896 bool ST_Intersects_Polygon_Polygon(int8_t* poly1_coords,
3897  int64_t poly1_coords_size,
3898  int32_t* poly1_ring_sizes,
3899  int64_t poly1_num_rings,
3900  double* poly1_bounds,
3901  int64_t poly1_bounds_size,
3902  int8_t* poly2_coords,
3903  int64_t poly2_coords_size,
3904  int32_t* poly2_ring_sizes,
3905  int64_t poly2_num_rings,
3906  double* poly2_bounds,
3907  int64_t poly2_bounds_size,
3908  int32_t ic1,
3909  int32_t isr1,
3910  int32_t ic2,
3911  int32_t isr2,
3912  int32_t osr) {
3913  if (poly1_bounds && poly2_bounds) {
3914  if (!box_overlaps_box(
3915  poly1_bounds, poly1_bounds_size, poly2_bounds, poly2_bounds_size)) {
3916  return false;
3917  }
3918  }
3919 
3920  return tol_zero(ST_Distance_Polygon_Polygon(poly1_coords,
3921  poly1_coords_size,
3922  poly1_ring_sizes,
3923  poly1_num_rings,
3924  poly2_coords,
3925  poly2_coords_size,
3926  poly2_ring_sizes,
3927  poly2_num_rings,
3928  ic1,
3929  isr1,
3930  ic2,
3931  isr2,
3932  osr,
3933  0.0));
3934 }
3935 
3937 bool ST_Intersects_Polygon_MultiPolygon(int8_t* poly_coords,
3938  int64_t poly_coords_size,
3939  int32_t* poly_ring_sizes,
3940  int64_t poly_num_rings,
3941  double* poly_bounds,
3942  int64_t poly_bounds_size,
3943  int8_t* mpoly_coords,
3944  int64_t mpoly_coords_size,
3945  int32_t* mpoly_ring_sizes,
3946  int64_t mpoly_num_rings,
3947  int32_t* mpoly_poly_sizes,
3948  int64_t mpoly_num_polys,
3949  double* mpoly_bounds,
3950  int64_t mpoly_bounds_size,
3951  int32_t ic1,
3952  int32_t isr1,
3953  int32_t ic2,
3954  int32_t isr2,
3955  int32_t osr) {
3956  if (poly_bounds && mpoly_bounds) {
3957  if (!box_overlaps_box(
3958  poly_bounds, poly_bounds_size, mpoly_bounds, mpoly_bounds_size)) {
3959  return false;
3960  }
3961  }
3962 
3963  return tol_zero(ST_Distance_Polygon_MultiPolygon(poly_coords,
3964  poly_coords_size,
3965  poly_ring_sizes,
3966  poly_num_rings,
3967  mpoly_coords,
3968  mpoly_coords_size,
3969  mpoly_ring_sizes,
3970  mpoly_num_rings,
3971  mpoly_poly_sizes,
3972  mpoly_num_polys,
3973  ic1,
3974  isr1,
3975  ic2,
3976  isr2,
3977  osr,
3978  0.0));
3979 }
3980 
3982 bool ST_Intersects_MultiPolygon_Point(int8_t* mpoly_coords,
3983  int64_t mpoly_coords_size,
3984  int32_t* mpoly_ring_sizes,
3985  int64_t mpoly_num_rings,
3986  int32_t* mpoly_poly_sizes,
3987  int64_t mpoly_num_polys,
3988  double* mpoly_bounds,
3989  int64_t mpoly_bounds_size,
3990  int8_t* p,
3991  int64_t psize,
3992  int32_t ic1,
3993  int32_t isr1,
3994  int32_t ic2,
3995  int32_t isr2,
3996  int32_t osr) {
3997  return Contains_MultiPolygon_Point_Impl<double, EdgeBehavior::kIncludePointOnEdge>(
3998  mpoly_coords,
3999  mpoly_coords_size,
4000  mpoly_ring_sizes,
4001  mpoly_num_rings,
4002  mpoly_poly_sizes,
4003  mpoly_num_polys,
4004  mpoly_bounds,
4005  mpoly_bounds_size,
4006  p,
4007  psize,
4008  ic1,
4009  isr1,
4010  ic2,
4011  isr2,
4012  osr);
4013 }
4014 
4016 bool ST_Intersects_MultiPolygon_LineString(int8_t* mpoly_coords,
4017  int64_t mpoly_coords_size,
4018  int32_t* mpoly_ring_sizes,
4019  int64_t mpoly_num_rings,
4020  int32_t* mpoly_poly_sizes,
4021  int64_t mpoly_num_polys,
4022  double* mpoly_bounds,
4023  int64_t mpoly_bounds_size,
4024  int8_t* l,
4025  int64_t lsize,
4026  double* lbounds,
4027  int64_t lbounds_size,
4028  int32_t ic1,
4029  int32_t isr1,
4030  int32_t ic2,
4031  int32_t isr2,
4032  int32_t osr) {
4034  lsize,
4035  lbounds,
4036  lbounds_size,
4037  mpoly_coords,
4038  mpoly_coords_size,
4039  mpoly_ring_sizes,
4040  mpoly_num_rings,
4041  mpoly_poly_sizes,
4042  mpoly_num_polys,
4043  mpoly_bounds,
4044  mpoly_bounds_size,
4045  ic2,
4046  isr2,
4047  ic1,
4048  isr1,
4049  osr);
4050 }
4051 
4053 bool ST_Intersects_MultiPolygon_Polygon(int8_t* mpoly_coords,
4054  int64_t mpoly_coords_size,
4055  int32_t* mpoly_ring_sizes,
4056  int64_t mpoly_num_rings,
4057  int32_t* mpoly_poly_sizes,
4058  int64_t mpoly_num_polys,
4059  double* mpoly_bounds,
4060  int64_t mpoly_bounds_size,
4061  int8_t* poly_coords,
4062  int64_t poly_coords_size,
4063  int32_t* poly_ring_sizes,
4064  int64_t poly_num_rings,
4065  double* poly_bounds,
4066  int64_t poly_bounds_size,
4067  int32_t ic1,
4068  int32_t isr1,
4069  int32_t ic2,
4070  int32_t isr2,
4071  int32_t osr) {
4072  return ST_Intersects_Polygon_MultiPolygon(poly_coords,
4073  poly_coords_size,
4074  poly_ring_sizes,
4075  poly_num_rings,
4076  poly_bounds,
4077  poly_bounds_size,
4078  mpoly_coords,
4079  mpoly_coords_size,
4080  mpoly_ring_sizes,
4081  mpoly_num_rings,
4082  mpoly_poly_sizes,
4083  mpoly_num_polys,
4084  mpoly_bounds,
4085  mpoly_bounds_size,
4086  ic2,
4087  isr2,
4088  ic1,
4089  isr1,
4090  osr);
4091 }
4092 
4094 bool ST_Intersects_MultiPolygon_MultiPolygon(int8_t* mpoly1_coords,
4095  int64_t mpoly1_coords_size,
4096  int32_t* mpoly1_ring_sizes,
4097  int64_t mpoly1_num_rings,
4098  int32_t* mpoly1_poly_sizes,
4099  int64_t mpoly1_num_polys,
4100  double* mpoly1_bounds,
4101  int64_t mpoly1_bounds_size,
4102  int8_t* mpoly2_coords,
4103  int64_t mpoly2_coords_size,
4104  int32_t* mpoly2_ring_sizes,
4105  int64_t mpoly2_num_rings,
4106  int32_t* mpoly2_poly_sizes,
4107  int64_t mpoly2_num_polys,
4108  double* mpoly2_bounds,
4109  int64_t mpoly2_bounds_size,
4110  int32_t ic1,
4111  int32_t isr1,
4112  int32_t ic2,
4113  int32_t isr2,
4114  int32_t osr) {
4115  if (mpoly1_bounds && mpoly2_bounds) {
4116  if (!box_overlaps_box(
4117  mpoly1_bounds, mpoly1_bounds_size, mpoly2_bounds, mpoly2_bounds_size)) {
4118  return false;
4119  }
4120  }
4121 
4122  return tol_zero(ST_Distance_MultiPolygon_MultiPolygon(mpoly1_coords,
4123  mpoly1_coords_size,
4124  mpoly1_ring_sizes,
4125  mpoly1_num_rings,
4126  mpoly1_poly_sizes,
4127  mpoly1_num_polys,
4128  mpoly2_coords,
4129  mpoly2_coords_size,
4130  mpoly2_ring_sizes,
4131  mpoly2_num_rings,
4132  mpoly2_poly_sizes,
4133  mpoly2_num_polys,
4134  ic1,
4135  isr1,
4136  ic2,
4137  isr2,
4138  osr,
4139  0.0));
4140 }
4141 
4142 //
4143 // Accessors for poly bounds and render group for in-situ poly render queries
4144 //
4145 // The MapD_* varieties are deprecated and renamed to "OmniSci_Geo*"
4146 // There may be some clients out there who are playing with the MapD_* so leaving
4147 // them for backwards compatibility.
4148 //
4149 
4151 int64_t OmniSci_Geo_PolyBoundsPtr(double* bounds, int64_t size) {
4152  return reinterpret_cast<int64_t>(bounds);
4153 }
4154 
4156 int32_t OmniSci_Geo_PolyRenderGroup(int32_t render_group) {
4157  return render_group;
4158 }
4159 
4161 int64_t MapD_GeoPolyBoundsPtr(double* bounds, int64_t size) {
4162  return OmniSci_Geo_PolyBoundsPtr(bounds, size);
4163 }
4164 
4166 int32_t MapD_GeoPolyRenderGroup(int32_t render_group) {
4167  return OmniSci_Geo_PolyRenderGroup(render_group);
4168 }
4169 
4171 double convert_meters_to_pixel_width(const double meters,
4172  int8_t* p,
4173  const int64_t psize,
4174  const int32_t ic,
4175  const int32_t isr,
4176  const int32_t osr,
4177  const double min_lon,
4178  const double max_lon,
4179  const int32_t img_width,
4180  const double min_width) {
4181  const double const1 = 0.017453292519943295769236907684886;
4182  const double const2 = 6372797.560856;
4183  const auto lon = decompress_coord(p, 0, ic, true);
4184  const auto lat = decompress_coord(p, 1, ic, false);
4185  double t1 = sinf(meters / (2.0 * const2));
4186  double t2 = cosf(const1 * lat);
4187  const double newlon = lon - (2.0 * asinf(t1 / t2)) / const1;
4188  t1 = transform_coord(lon, isr, osr, true);
4189  t2 = transform_coord(newlon, isr, osr, true);
4190  const double min_domain_x = transform_coord(min_lon, isr, osr, true);
4191  const double max_domain_x = transform_coord(max_lon, isr, osr, true);
4192  const double domain_diff = max_domain_x - min_domain_x;
4193  t1 = ((t1 - min_domain_x) / domain_diff) * static_cast<double>(img_width);
4194  t2 = ((t2 - min_domain_x) / domain_diff) * static_cast<double>(img_width);
4195 
4196  // TODO(croot): need to account for edge cases, such as getting close to the poles.
4197  const double sz = fabs(t1 - t2);
4198  return (sz < min_width ? min_width : sz);
4199 }
4200 
4202 double convert_meters_to_pixel_height(const double meters,
4203  int8_t* p,
4204  const int64_t psize,
4205  const int32_t ic,
4206  const int32_t isr,
4207  const int32_t osr,
4208  const double min_lat,
4209  const double max_lat,
4210  const int32_t img_height,
4211  const double min_height) {
4212  const double const1 = 0.017453292519943295769236907684886;
4213  const double const2 = 6372797.560856;
4214  const auto lat = decompress_coord(p, 1, ic, false);
4215  const double latdiff = meters / (const1 * const2);
4216  const double newlat =
4217  (lat < 0) ? lat + latdiff : lat - latdiff; // assumes a lat range of [-90, 90]
4218  double t1 = transform_coord(lat, isr, osr, false);
4219  double t2 = transform_coord(newlat, isr, osr, false);
4220  const double min_domain_y = transform_coord(min_lat, isr, osr, false);
4221  const double max_domain_y = transform_coord(max_lat, isr, osr, false);
4222  const double domain_diff = max_domain_y - min_domain_y;
4223  t1 = ((t1 - min_domain_y) / domain_diff) * static_cast<double>(img_height);
4224  t2 = ((t2 - min_domain_y) / domain_diff) * static_cast<double>(img_height);
4225 
4226  // TODO(croot): need to account for edge cases, such as getting close to the poles.
4227  const double sz = fabs(t1 - t2);
4228  return (sz < min_height ? min_height : sz);
4229 }
4230 
4232  const int64_t psize,
4233  const int32_t ic,
4234  const double min_lon,
4235  const double max_lon,
4236  const double min_lat,
4237  const double max_lat) {
4238  const auto lon = decompress_coord(p, 0, ic, true);
4239  const auto lat = decompress_coord(p, 1, ic, false);
4240  return !(lon < min_lon || lon > max_lon || lat < min_lat || lat > max_lat);
4241 }
4242 
4244  const int64_t psize,
4245  const int32_t ic,
4246  const double meters,
4247  const double min_lon,
4248  const double max_lon,
4249  const double min_lat,
4250  const double max_lat) {
4251  const double const1 = 0.017453292519943295769236907684886;
4252  const double const2 = 6372797.560856;
4253  const auto lon = decompress_coord(p, 0, ic, true);
4254  const auto lat = decompress_coord(p, 1, ic, false);
4255  const double latdiff = meters / (const1 * const2);
4256  const double t1 = sinf(meters / (2.0 * const2));
4257  const double t2 = cosf(const1 * lat);
4258  const double londiff = (2.0 * asinf(t1 / t2)) / const1;
4259  return !(lon + londiff < min_lon || lon - londiff > max_lon ||
4260  lat + latdiff < min_lat || lat - latdiff > max_lat);
4261 }
EXTENSION_NOINLINE double ST_Perimeter_MultiPolygon(int8_t *mpoly_coords, int32_t mpoly_coords_size, int8_t *mpoly_ring_sizes, int32_t mpoly_num_rings, int8_t *mpoly_poly_sizes, int32_t mpoly_num_polys, int32_t ic, int32_t isr, int32_t osr)
EXTENSION_NOINLINE double convert_meters_to_pixel_width(const double meters, int8_t *p, const int64_t psize, const int32_t ic, const int32_t isr, const int32_t osr, const double min_lon, const double max_lon, const int32_t img_width, const double min_width)
EXTENSION_INLINE bool ST_Intersects_Point_Polygon(int8_t *p, int64_t psize, int8_t *poly, int64_t polysize, int32_t *poly_ring_sizes, int64_t poly_num_rings, double *poly_bounds, int64_t poly_bounds_size, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr)
DEVICE ALWAYS_INLINE double length_linestring(int8_t *l, int32_t lsize, int32_t ic, int32_t isr, int32_t osr, bool geodesic, bool check_closed)
EXTENSION_NOINLINE double ST_Distance_LineString_LineString(int8_t *l1, int64_t l1size, int8_t *l2, int64_t l2size, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr, double threshold)
EXTENSION_NOINLINE double ST_Distance_Point_Point(int8_t *p1, int64_t p1size, int8_t *p2, int64_t p2size, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr)
EXTENSION_INLINE double ST_Distance_LineString_Point_Geodesic(int8_t *l, int64_t lsize, int8_t *p, int64_t psize, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr)
EXTENSION_NOINLINE double ST_X_Point(int8_t *p, int64_t psize, int32_t ic, int32_t isr, int32_t osr)
DEVICE ALWAYS_INLINE bool on_segment(double px, double py, double qx, double qy, double rx, double ry)
DEVICE ALWAYS_INLINE double area_polygon(int8_t *poly_coords, int32_t poly_coords_size, int8_t *poly_ring_sizes_in, int32_t poly_num_rings, int32_t ic, int32_t isr, int32_t osr)
EXTENSION_INLINE bool ST_DWithin_MultiPolygon_MultiPolygon(int8_t *mpoly1_coords, int64_t mpoly1_coords_size, int32_t *mpoly1_ring_sizes, int64_t mpoly1_num_rings, int32_t *mpoly1_poly_sizes, int64_t mpoly1_num_polys, double *mpoly1_bounds, int64_t mpoly1_bounds_size, int8_t *mpoly2_coords, int64_t mpoly2_coords_size, int32_t *mpoly2_ring_sizes, int64_t mpoly2_num_rings, int32_t *mpoly2_poly_sizes, int64_t mpoly2_num_polys, double *mpoly2_bounds, int64_t mpoly2_bounds_size, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr, double distance_within)
DEVICE ALWAYS_INLINE double transform_coord(const double coord, const int32_t isr, const int32_t osr, const bool x)
EXTENSION_INLINE bool ST_Intersects_MultiPolygon_Polygon(int8_t *mpoly_coords, int64_t mpoly_coords_size, int32_t *mpoly_ring_sizes, int64_t mpoly_num_rings, int32_t *mpoly_poly_sizes, int64_t mpoly_num_polys, double *mpoly_bounds, int64_t mpoly_bounds_size, int8_t *poly_coords, int64_t poly_coords_size, int32_t *poly_ring_sizes, int64_t poly_num_rings, double *poly_bounds, int64_t poly_bounds_size, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr)
DEVICE ALWAYS_INLINE bool box_dwithin_box(double *bounds1, int64_t bounds1_size, int32_t isr1, double *bounds2, int64_t bounds2_size, int32_t isr2, int32_t osr, double distance)
EXTENSION_NOINLINE bool ST_Contains_Point_Point(int8_t *p1, int64_t p1size, int8_t *p2, int64_t p2size, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr)
EXTENSION_NOINLINE double ST_Distance_Point_LineString_Geodesic(int8_t *p, int64_t psize, int8_t *l, int64_t lsize, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr)
EXTENSION_INLINE int32_t OmniSci_Geo_PolyRenderGroup(int32_t render_group)
DEVICE ALWAYS_INLINE bool centroid_add_triangle(double x1, double y1, double x2, double y2, double x3, double y3, double sign, double *total_area2, double *cg3)
EXTENSION_NOINLINE double ST_XMax(int8_t *coords, int64_t size, int32_t ic, int32_t isr, int32_t osr)
DEVICE ALWAYS_INLINE double max_distance_point_linestring(int8_t *p, int64_t psize, int8_t *l, int64_t lsize, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr, bool check_closed)
EXTENSION_NOINLINE bool ST_Intersects_LineString_MultiPolygon(int8_t *l, int64_t lsize, double *lbounds, int64_t lbounds_size, int8_t *mpoly_coords, int64_t mpoly_coords_size, int32_t *mpoly_ring_sizes, int64_t mpoly_num_rings, int32_t *mpoly_poly_sizes, int64_t mpoly_num_polys, double *mpoly_bounds, int64_t mpoly_bounds_size, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr)
EXTENSION_NOINLINE double ST_Distance_Point_LineString(int8_t *p, int64_t psize, int8_t *l, int64_t lsize, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr, double threshold)
EXTENSION_INLINE bool ST_DWithin_Point_Polygon(int8_t *p, int64_t psize, int8_t *poly_coords, int64_t poly_coords_size, int32_t *poly_ring_sizes, int64_t poly_num_rings, double *poly_bounds, int64_t poly_bounds_size, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr, double distance_within)
EXTENSION_NOINLINE bool ST_Intersects_LineString_Polygon(int8_t *l, int64_t lsize, double *lbounds, int64_t lbounds_size, int8_t *poly, int64_t polysize, int32_t *poly_ring_sizes, int64_t poly_num_rings, double *poly_bounds, int64_t poly_bounds_size, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr)
DEVICE ALWAYS_INLINE double perimeter_multipolygon(int8_t *mpoly_coords, int32_t mpoly_coords_size, int8_t *mpoly_ring_sizes_in, int32_t mpoly_num_rings, int8_t *mpoly_poly_sizes, int32_t mpoly_num_polys, int32_t ic, int32_t isr, int32_t osr, bool geodesic)
EXTENSION_NOINLINE double ST_Distance_Point_MultiPolygon(int8_t *p, int64_t psize, int8_t *mpoly_coords, int64_t mpoly_coords_size, int32_t *mpoly_ring_sizes, int64_t mpoly_num_rings, int32_t *mpoly_poly_sizes, int64_t mpoly_num_polys, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr, double threshold)
EXTENSION_NOINLINE bool ST_Intersects_Polygon_Polygon(int8_t *poly1_coords, int64_t poly1_coords_size, int32_t *poly1_ring_sizes, int64_t poly1_num_rings, double *poly1_bounds, int64_t poly1_bounds_size, int8_t *poly2_coords, int64_t poly2_coords_size, int32_t *poly2_ring_sizes, int64_t poly2_num_rings, double *poly2_bounds, int64_t poly2_bounds_size, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr)
EXTENSION_NOINLINE bool is_point_size_in_view(int8_t *p, const int64_t psize, const int32_t ic, const double meters, const double min_lon, const double max_lon, const double min_lat, const double max_lat)
DEVICE double distance_point_line(double px, double py, double l1x, double l1y, double l2x, double l2y)
EXTENSION_INLINE double ST_YMin_Bounds(double *bounds, int64_t size, int32_t isr, int32_t osr)
EXTENSION_NOINLINE double ST_YMax(int8_t *coords, int64_t size, int32_t ic, int32_t isr, int32_t osr)
DEVICE ALWAYS_INLINE int32_t compressed_coord(const int8_t *data, const int32_t index)
EXTENSION_INLINE double ST_XMin_Bounds(double *bounds, int64_t size, int32_t isr, int32_t osr)
EXTENSION_NOINLINE bool ST_cContains_Polygon_Point(const int8_t *poly_coords, const int64_t poly_coords_size, const int32_t *poly_ring_sizes, const int64_t poly_num_rings, const double *poly_bounds, const int64_t poly_bounds_size, const int8_t *p, const int64_t psize, const int32_t ic1, const int32_t isr1, const int32_t ic2, const int32_t isr2, const int32_t osr)
DEVICE ALWAYS_INLINE T is_left(const T lx0, const T ly0, const T lx1, const T ly1, const T px, const T py)
DEVICE ALWAYS_INLINE bool Contains_MultiPolygon_Point_Impl(const int8_t *mpoly_coords, const int64_t mpoly_coords_size, const int32_t *mpoly_ring_sizes, const int64_t mpoly_num_rings, const int32_t *mpoly_poly_sizes, const int64_t mpoly_num_polys, const double *mpoly_bounds, const int64_t mpoly_bounds_size, const int8_t *p, const int64_t psize, const int32_t ic1, const int32_t isr1, const int32_t ic2, const int32_t isr2, const int32_t osr)
DEVICE double decompress_lattitude_coord_geoint32(const int32_t compressed)
DEVICE ALWAYS_INLINE double area_triangle(double x1, double y1, double x2, double y2, double x3, double y3)
EXTENSION_INLINE bool ST_DWithin_Polygon_Polygon(int8_t *poly1_coords, int64_t poly1_coords_size, int32_t *poly1_ring_sizes, int64_t poly1_num_rings, double *poly1_bounds, int64_t poly1_bounds_size, int8_t *poly2_coords, int64_t poly2_coords_size, int32_t *poly2_ring_sizes, int64_t poly2_num_rings, double *poly2_bounds, int64_t poly2_bounds_size, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr, double distance_within)
EXTENSION_NOINLINE double ST_XMin(int8_t *coords, int64_t size, int32_t ic, int32_t isr, int32_t osr)
EXTENSION_INLINE bool ST_Intersects_Polygon_Point(int8_t *poly, int64_t polysize, int32_t *poly_ring_sizes, int64_t poly_num_rings, double *poly_bounds, int64_t poly_bounds_size, int8_t *p, int64_t psize, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr)
EXTENSION_NOINLINE double ST_Distance_Point_Polygon(int8_t *p, int64_t psize, int8_t *poly, int64_t polysize, int32_t *poly_ring_sizes, int64_t poly_num_rings, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr, double threshold)
EXTENSION_NOINLINE double ST_Distance_Point_Point_Squared(int8_t *p1, int64_t p1size, int8_t *p2, int64_t p2size, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr)
EXTENSION_INLINE bool ST_DWithin_Point_Point(int8_t *p1, int64_t p1size, int8_t *p2, int64_t p2size, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr, double distance_within)
EXTENSION_NOINLINE double ST_Y_Point(int8_t *p, int64_t psize, int32_t ic, int32_t isr, int32_t osr)
EXTENSION_NOINLINE double convert_meters_to_pixel_height(const double meters, int8_t *p, const int64_t psize, const int32_t ic, const int32_t isr, const int32_t osr, const double min_lat, const double max_lat, const int32_t img_height, const double min_height)
#define TOLERANCE_GEOINT32
DEVICE double max_distance_point_line(double px, double py, double l1x, double l1y, double l2x, double l2y)
EXTENSION_INLINE int32_t MapD_GeoPolyRenderGroup(int32_t render_group)
EXTENSION_NOINLINE bool Point_Overlaps_Box(double *bounds, int64_t bounds_size, double px, double py)
EXTENSION_INLINE int64_t OmniSci_Geo_PolyBoundsPtr(double *bounds, int64_t size)
#define DEVICE
#define EXTENSION_NOINLINE
Definition: OmniSciTypes.h:28
EXTENSION_NOINLINE double ST_Distance_Polygon_MultiPolygon(int8_t *poly1_coords, int64_t poly1_coords_size, int32_t *poly1_ring_sizes, int64_t poly1_num_rings, int8_t *mpoly_coords, int64_t mpoly_coords_size, int32_t *mpoly_ring_sizes, int64_t mpoly_num_rings, int32_t *mpoly_poly_sizes, int64_t mpoly_num_polys, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr, double threshold)
EXTENSION_NOINLINE void ST_Centroid_Polygon(int8_t *poly_coords, int32_t poly_coords_size, int32_t *poly_ring_sizes, int32_t poly_num_rings, int32_t ic, int32_t isr, int32_t osr, double *poly_centroid)
EXTENSION_NOINLINE double ST_Distance_Polygon_Polygon(int8_t *poly1_coords, int64_t poly1_coords_size, int32_t *poly1_ring_sizes, int64_t poly1_num_rings, int8_t *poly2_coords, int64_t poly2_coords_size, int32_t *poly2_ring_sizes, int64_t poly2_num_rings, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr, double threshold)
EXTENSION_NOINLINE bool ST_Intersects_Point_LineString(int8_t *p, int64_t psize, int8_t *l, int64_t lsize, double *lbounds, int64_t lbounds_size, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr)
EXTENSION_INLINE double ST_Distance_MultiPolygon_Point(int8_t *mpoly_coords, int64_t mpoly_coords_size, int32_t *mpoly_ring_sizes, int64_t mpoly_num_rings, int32_t *mpoly_poly_sizes, int64_t mpoly_num_polys, int8_t *p, int64_t psize, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr, double threshold)
EXTENSION_NOINLINE double ST_MaxDistance_Point_LineString(int8_t *p, int64_t psize, int8_t *l, int64_t lsize, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr)
EXTENSION_NOINLINE bool ST_Contains_Polygon_Point(const int8_t *poly_coords, const int64_t poly_coords_size, const int32_t *poly_ring_sizes, const int64_t poly_num_rings, const double *poly_bounds, const int64_t poly_bounds_size, const int8_t *p, const int64_t psize, const int32_t ic1, const int32_t isr1, const int32_t ic2, const int32_t isr2, const int32_t osr)
DEVICE ALWAYS_INLINE int16_t orientation(double px, double py, double qx, double qy, double rx, double ry)
DEVICE double distance_line_line_squared(double l11x, double l11y, double l12x, double l12y, double l21x, double l21y, double l22x, double l22y)
EXTENSION_INLINE double ST_Distance_MultiPolygon_LineString(int8_t *mpoly_coords, int64_t mpoly_coords_size, int32_t *mpoly_ring_sizes, int64_t mpoly_num_rings, int32_t *mpoly_poly_sizes, int64_t mpoly_num_polys, int8_t *l, int64_t lsize, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr, double threshold)
EXTENSION_NOINLINE double ST_YMin(int8_t *coords, int64_t size, int32_t ic, int32_t isr, int32_t osr)
DEVICE ALWAYS_INLINE double area_ring(int8_t *ring, int64_t ringsize, int32_t ic, int32_t isr, int32_t osr)
EXTENSION_INLINE bool ST_DWithin_Polygon_MultiPolygon(int8_t *poly_coords, int64_t poly_coords_size, int32_t *poly_ring_sizes, int64_t poly_num_rings, double *poly_bounds, int64_t poly_bounds_size, int8_t *mpoly_coords, int64_t mpoly_coords_size, int32_t *mpoly_ring_sizes, int64_t mpoly_num_rings, int32_t *mpoly_poly_sizes, int64_t mpoly_num_polys, double *mpoly_bounds, int64_t mpoly_bounds_size, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr, double distance_within)
EXTENSION_INLINE double ST_Distance_LineString_Point(int8_t *l, int64_t lsize, int8_t *p, int64_t psize, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr, double threshold)
EXTENSION_INLINE bool ST_Contains_LineString_Point(int8_t *l, int64_t lsize, double *lbounds, int64_t lbounds_size, int8_t *p, int64_t psize, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr)
EXTENSION_NOINLINE bool ST_cContains_MultiPolygon_Point(int8_t *mpoly_coords, int64_t mpoly_coords_size, int32_t *mpoly_ring_sizes, int64_t mpoly_num_rings, int32_t *mpoly_poly_sizes, int64_t mpoly_num_polys, double *mpoly_bounds, int64_t mpoly_bounds_size, int8_t *p, int64_t psize, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr)
DEVICE ALWAYS_INLINE bool box_overlaps_box(double *bounds1, int64_t bounds1_size, double *bounds2, int64_t bounds2_size)
#define DEBUG_STMT(x)
EXTENSION_INLINE int64_t MapD_GeoPolyBoundsPtr(double *bounds, int64_t size)
EXTENSION_NOINLINE void ST_Centroid_Point(int8_t *p, int32_t psize, int32_t ic, int32_t isr, int32_t osr, double *point_centroid)
EXTENSION_NOINLINE void ST_Centroid_MultiPolygon(int8_t *mpoly_coords, int32_t mpoly_coords_size, int32_t *mpoly_ring_sizes, int32_t mpoly_num_rings, int32_t *mpoly_poly_sizes, int32_t mpoly_num_polys, int32_t ic, int32_t isr, int32_t osr, double *mpoly_centroid)
EXTENSION_INLINE bool ST_Intersects_Point_MultiPolygon(int8_t *p, int64_t psize, int8_t *mpoly_coords, int64_t mpoly_coords_size, int32_t *mpoly_ring_sizes, int64_t mpoly_num_rings, int32_t *mpoly_poly_sizes, int64_t mpoly_num_polys, double *mpoly_bounds, int64_t mpoly_bounds_size, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr)
EXTENSION_INLINE double ST_Distance_MultiPolygon_Polygon(int8_t *mpoly_coords, int64_t mpoly_coords_size, int32_t *mpoly_ring_sizes, int64_t mpoly_num_rings, int32_t *mpoly_poly_sizes, int64_t mpoly_num_polys, int8_t *poly1_coords, int64_t poly1_coords_size, int32_t *poly1_ring_sizes, int64_t poly1_num_rings, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr, double threshold)
DEVICE bool linestring_intersects_line(int8_t *l, int32_t lnum_coords, double l1x, double l1y, double l2x, double l2y, int32_t ic1, int32_t isr1, int32_t osr)
#define TOLERANCE_DEFAULT
EXTENSION_NOINLINE double ST_Distance_LineString_Polygon(int8_t *l, int64_t lsize, int8_t *poly_coords, int64_t poly_coords_size, int32_t *poly_ring_sizes, int64_t poly_num_rings, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr, double threshold)
DEVICE bool line_intersects_line(double l11x, double l11y, double l12x, double l12y, double l21x, double l21y, double l22x, double l22y)
EXTENSION_NOINLINE bool ST_Contains_Polygon_LineString(int8_t *poly_coords, int64_t poly_coords_size, int32_t *poly_ring_sizes, int64_t poly_num_rings, double *poly_bounds, int64_t poly_bounds_size, int8_t *l, int64_t lsize, double *lbounds, int64_t lbounds_size, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr)
EXTENSION_NOINLINE bool ST_Intersects_Polygon_MultiPolygon(int8_t *poly_coords, int64_t poly_coords_size, int32_t *poly_ring_sizes, int64_t poly_num_rings, double *poly_bounds, int64_t poly_bounds_size, int8_t *mpoly_coords, int64_t mpoly_coords_size, int32_t *mpoly_ring_sizes, int64_t mpoly_num_rings, int32_t *mpoly_poly_sizes, int64_t mpoly_num_polys, double *mpoly_bounds, int64_t mpoly_bounds_size, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr)
EXTENSION_NOINLINE double ST_Distance_LineString_MultiPolygon(int8_t *l, int64_t lsize, int8_t *mpoly_coords, int64_t mpoly_coords_size, int32_t *mpoly_ring_sizes, int64_t mpoly_num_rings, int32_t *mpoly_poly_sizes, int64_t mpoly_num_polys, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr, double threshold)
EXTENSION_NOINLINE double ST_MaxDistance_LineString_Point(int8_t *l, int64_t lsize, int8_t *p, int64_t psize, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr)
DEVICE ALWAYS_INLINE bool tol_zero_template(const T x, const T tolerance=TOLERANCE_DEFAULT)
EXTENSION_NOINLINE bool is_point_in_view(int8_t *p, const int64_t psize, const int32_t ic, const double min_lon, const double max_lon, const double min_lat, const double max_lat)
EXTENSION_INLINE bool ST_Intersects_Polygon_LineString(int8_t *poly, int64_t polysize, int32_t *poly_ring_sizes, int64_t poly_num_rings, double *poly_bounds, int64_t poly_bounds_size, int8_t *l, int64_t lsize, double *lbounds, int64_t lbounds_size, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr)
#define UNLIKELY(x)
Definition: likely.h:25
EXTENSION_INLINE bool ST_DWithin_Point_LineString(int8_t *p1, int64_t p1size, int8_t *l2, int64_t l2size, double *l2bounds, int64_t l2bounds_size, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr, double distance_within)
EXTENSION_NOINLINE bool ST_Contains_MultiPolygon_LineString(int8_t *mpoly_coords, int64_t mpoly_coords_size, int32_t *mpoly_ring_sizes, int64_t mpoly_num_rings, int32_t *mpoly_poly_sizes, int64_t mpoly_num_polys, double *mpoly_bounds, int64_t mpoly_bounds_size, int8_t *l, int64_t lsize, double *lbounds, int64_t lbounds_size, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr)
DEVICE ALWAYS_INLINE double decompress_coord(const int8_t *data, const int32_t index, const int32_t ic, const bool x)
DEVICE ALWAYS_INLINE bool centroid_add_polygon(int8_t *poly_coords, int64_t poly_coords_size, int32_t *poly_ring_sizes, int64_t poly_num_rings, int32_t ic, int32_t isr, int32_t osr, double *total_area2, double *cg3, double *total_length, double *linestring_centroid_sum, int64_t *num_points, double *point_centroid_sum)
EXTENSION_NOINLINE bool ST_Contains_Point_Polygon(int8_t *p, int64_t psize, int8_t *poly_coords, int64_t poly_coords_size, int32_t *poly_ring_sizes, int64_t poly_num_rings, double *poly_bounds, int64_t poly_bounds_size, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr)
EXTENSION_NOINLINE double ST_Distance_Point_Point_Geodesic(int8_t *p1, int64_t p1size, int8_t *p2, int64_t p2size, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr)
DEVICE ALWAYS_INLINE bool box_contains_point(const double *bounds, const int64_t bounds_size, const double px, const double py)
DEVICE double distance_ring_ring(int8_t *ring1, int32_t ring1_num_coords, int8_t *ring2, int32_t ring2_num_coords, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr, double threshold)
EXTENSION_NOINLINE double ST_Length_LineString(int8_t *coords, int64_t coords_sz, int32_t ic, int32_t isr, int32_t osr)
DEVICE double decompress_longitude_coord_geoint32(const int32_t compressed)
EXTENSION_INLINE bool ST_Intersects_Point_Point(int8_t *p1, int64_t p1size, int8_t *p2, int64_t p2size, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr)
EXTENSION_NOINLINE double ST_Perimeter_MultiPolygon_Geodesic(int8_t *mpoly_coords, int32_t mpoly_coords_size, int8_t *mpoly_ring_sizes, int32_t mpoly_num_rings, int8_t *mpoly_poly_sizes, int32_t mpoly_num_polys, int32_t ic, int32_t isr, int32_t osr)
DEVICE ALWAYS_INLINE bool tol_le(const double x, const double y, const double tolerance=TOLERANCE_DEFAULT)
DEVICE ALWAYS_INLINE bool centroid_add_segment(double x1, double y1, double x2, double y2, double *length, double *linestring_centroid_sum)
EXTENSION_INLINE bool ST_DWithin_LineString_LineString(int8_t *l1, int64_t l1size, double *l1bounds, int64_t l1bounds_size, int8_t *l2, int64_t l2size, double *l2bounds, int64_t l2bounds_size, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr, double distance_within)
DEVICE ALWAYS_INLINE bool tol_zero(const double x, const double tolerance=TOLERANCE_DEFAULT)
DEVICE ALWAYS_INLINE bool centroid_add_linestring(int8_t *l, int64_t lsize, int32_t ic, int32_t isr, int32_t osr, bool closed, double *total_length, double *linestring_centroid_sum, int64_t *num_points, double *point_centroid_sum)
EXTENSION_NOINLINE bool ST_Intersects_LineString_Linestring(int8_t *l1, int64_t l1size, double *l1bounds, int64_t l1bounds_size, int8_t *l2, int64_t l2size, double *l2bounds, int64_t l2bounds_size, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr)
DEVICE ALWAYS_INLINE double distance_point_point_squared(double p1x, double p1y, double p2x, double p2y)
DEVICE ALWAYS_INLINE double distance_point_linestring(int8_t *p, int64_t psize, int8_t *l, int64_t lsize, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr, bool check_closed, double threshold)
EXTENSION_NOINLINE double ST_Area_MultiPolygon(int8_t *mpoly_coords, int32_t mpoly_coords_size, int8_t *mpoly_ring_sizes, int32_t mpoly_num_rings, int8_t *mpoly_poly_sizes_in, int32_t mpoly_num_polys, int32_t ic, int32_t isr, int32_t osr)
DEVICE bool linestring_intersects_linestring(int8_t *l, int32_t lnum_coords, double l1x, double l1y, double l2x, double l2y, int32_t ic1, int32_t isr1, int32_t osr)
EXTENSION_INLINE bool ST_Intersects_MultiPolygon_MultiPolygon(int8_t *mpoly1_coords, int64_t mpoly1_coords_size, int32_t *mpoly1_ring_sizes, int64_t mpoly1_num_rings, int32_t *mpoly1_poly_sizes, int64_t mpoly1_num_polys, double *mpoly1_bounds, int64_t mpoly1_bounds_size, int8_t *mpoly2_coords, int64_t mpoly2_coords_size, int32_t *mpoly2_ring_sizes, int64_t mpoly2_num_rings, int32_t *mpoly2_poly_sizes, int64_t mpoly2_num_polys, double *mpoly2_bounds, int64_t mpoly2_bounds_size, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr)
EXTENSION_NOINLINE double ST_Perimeter_Polygon(int8_t *poly, int32_t polysize, int8_t *poly_ring_sizes, int32_t poly_num_rings, int32_t ic, int32_t isr, int32_t osr)
EXTENSION_INLINE bool ST_Intersects_MultiPolygon_LineString(int8_t *mpoly_coords, int64_t mpoly_coords_size, int32_t *mpoly_ring_sizes, int64_t mpoly_num_rings, int32_t *mpoly_poly_sizes, int64_t mpoly_num_polys, double *mpoly_bounds, int64_t mpoly_bounds_size, int8_t *l, int64_t lsize, double *lbounds, int64_t lbounds_size, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr)
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 ST_Length_LineString_Geodesic(int8_t *coords, int64_t coords_sz, int32_t ic, int32_t isr, int32_t osr)
DEVICE ALWAYS_INLINE bool tol_ge(const double x, const double y, const double tolerance=TOLERANCE_DEFAULT)
DEVICE ALWAYS_INLINE bool polygon_contains_point(const int8_t *poly, const int32_t poly_num_coords, const double px, const double py, const int32_t ic1, const int32_t isr1, const int32_t osr)
DEVICE ALWAYS_INLINE int32_t compression_unit_size(const int32_t ic)
DEVICE ALWAYS_INLINE double transform_coord_y(const double coord, const int32_t isr, const int32_t osr)
#define TOLERANCE_DEFAULT_SQUARED
EXTENSION_INLINE bool ST_Intersects_LineString_Point(int8_t *l, int64_t lsize, double *lbounds, int64_t lbounds_size, int8_t *p, int64_t psize, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr)
EXTENSION_NOINLINE double ST_Distance_MultiPolygon_MultiPolygon(int8_t *mpoly1_coords, int64_t mpoly1_coords_size, int32_t *mpoly1_ring_sizes, int64_t mpoly1_num_rings, int32_t *mpoly1_poly_sizes, int64_t mpoly1_num_polys, int8_t *mpoly2_coords, int64_t mpoly2_coords_size, int32_t *mpoly2_ring_sizes, int64_t mpoly2_num_rings, int32_t *mpoly2_poly_sizes, int64_t mpoly2_num_polys, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr, double threshold)
EXTENSION_NOINLINE double ST_Distance_Point_ClosedLineString(int8_t *p, int64_t psize, int8_t *l, int64_t lsize, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr, double threshold)
#define EXTENSION_INLINE
Definition: OmniSciTypes.h:27
EXTENSION_NOINLINE double ST_Area_Polygon(int8_t *poly_coords, int32_t poly_coords_size, int8_t *poly_ring_sizes, int32_t poly_num_rings, int32_t ic, int32_t isr, int32_t osr)
EXTENSION_NOINLINE bool ST_Contains_Polygon_Polygon(int8_t *poly1_coords, int64_t poly1_coords_size, int32_t *poly1_ring_sizes, int64_t poly1_num_rings, double *poly1_bounds, int64_t poly1_bounds_size, int8_t *poly2_coords, int64_t poly2_coords_size, int32_t *poly2_ring_sizes, int64_t poly2_num_rings, double *poly2_bounds, int64_t poly2_bounds_size, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr)
DEVICE ALWAYS_INLINE bool point_in_polygon_winding_number(const int8_t *poly, const int32_t poly_num_coords, const T px, const T py, const int32_t ic1, const int32_t isr1, const int32_t osr)
EXTENSION_INLINE double ST_Distance_Polygon_Point(int8_t *poly_coords, int64_t poly_coords_size, int32_t *poly_ring_sizes, int64_t poly_num_rings, int8_t *p, int64_t psize, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr, double threshold)
DEVICE double distance_point_line_squared(double px, double py, double l1x, double l1y, double l2x, double l2y)
#define COMPRESSION_GEOINT32
EXTENSION_NOINLINE double ST_Perimeter_Polygon_Geodesic(int8_t *poly, int32_t polysize, int8_t *poly_ring_sizes_in, int32_t poly_num_rings, int32_t ic, int32_t isr, int32_t osr)
DEVICE ALWAYS_INLINE double coord_y(const int8_t *data, const int32_t index, const int32_t ic, const int32_t isr, const int32_t osr)
DEVICE ALWAYS_INLINE double tol(int32_t ic)
EXTENSION_INLINE bool ST_DWithin_LineString_Polygon(int8_t *l1, int64_t l1size, double *l1bounds, int64_t l1bounds_size, int8_t *poly_coords, int64_t poly_coords_size, int32_t *poly_ring_sizes, int64_t poly_num_rings, double *poly_bounds, int64_t poly_bounds_size, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr, double distance_within)
DEVICE double distance_line_line(double l11x, double l11y, double l12x, double l12y, double l21x, double l21y, double l22x, double l22y)
DEVICE bool polygon_contains_linestring(int8_t *poly, int32_t poly_num_coords, int8_t *l, int64_t lnum_coords, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr)
EXTENSION_NOINLINE bool ST_Contains_LineString_Polygon(int8_t *l, int64_t lsize, double *lbounds, int64_t lbounds_size, int8_t *poly_coords, int64_t poly_coords_size, int32_t *poly_ring_sizes, int64_t poly_num_rings, double *poly_bounds, int64_t poly_bounds_size, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr)
DEVICE ALWAYS_INLINE bool Contains_Polygon_Point_Impl(const int8_t *poly_coords, const int64_t poly_coords_size, const int32_t *poly_ring_sizes, const int64_t poly_num_rings, const double *poly_bounds, const int64_t poly_bounds_size, const int8_t *p, const int64_t psize, const int32_t ic1, const int32_t isr1, const int32_t ic2, const int32_t isr2, const int32_t osr)
EXTENSION_NOINLINE bool ST_Contains_MultiPolygon_Point(int8_t *mpoly_coords, int64_t mpoly_coords_size, int32_t *mpoly_ring_sizes, int64_t mpoly_num_rings, int32_t *mpoly_poly_sizes, int64_t mpoly_num_polys, double *mpoly_bounds, int64_t mpoly_bounds_size, int8_t *p, int64_t psize, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr)
EXTENSION_INLINE bool ST_Intersects_MultiPolygon_Point(int8_t *mpoly_coords, int64_t mpoly_coords_size, int32_t *mpoly_ring_sizes, int64_t mpoly_num_rings, int32_t *mpoly_poly_sizes, int64_t mpoly_num_polys, double *mpoly_bounds, int64_t mpoly_bounds_size, int8_t *p, int64_t psize, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr)
DEVICE bool ring_intersects_line(int8_t *ring, int32_t ring_num_coords, double l1x, double l1y, double l2x, double l2y, int32_t ic1, int32_t isr1, int32_t osr)
DEVICE double distance_ring_linestring(int8_t *ring, int32_t ring_num_coords, int8_t *l, int32_t lnum_coords, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr, double threshold)
#define ALWAYS_INLINE
EXTENSION_NOINLINE double conv_4326_900913_y(const double y)
EXTENSION_INLINE bool ST_DWithin_Point_Point_Geodesic(int8_t *p1, int64_t p1size, int8_t *p2, int64_t p2size, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr, double distance_within)
EXTENSION_NOINLINE double conv_4326_900913_x(const double x)
DEVICE ALWAYS_INLINE bool tol_eq(const double x, const double y, const double tolerance=TOLERANCE_DEFAULT)
DEVICE ALWAYS_INLINE bool box_contains_box_vertex(double *bounds1, int64_t bounds1_size, double *bounds2, int64_t bounds2_size)
EXTENSION_INLINE double ST_Distance_Polygon_LineString(int8_t *poly_coords, int64_t poly_coords_size, int32_t *poly_ring_sizes, int64_t poly_num_rings, int8_t *l, int64_t lsize, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr, double threshold)
DEVICE ALWAYS_INLINE bool box_contains_box(double *bounds1, int64_t bounds1_size, double *bounds2, int64_t bounds2_size)
DEVICE ALWAYS_INLINE bool point_dwithin_box(int8_t *p1, int64_t p1size, int32_t ic1, int32_t isr1, double *bounds2, int64_t bounds2_size, int32_t isr2, int32_t osr, double distance)
DEVICE ALWAYS_INLINE double coord_x(const int8_t *data, const int32_t index, const int32_t ic, const int32_t isr, const int32_t osr)
EXTENSION_INLINE double ST_YMax_Bounds(double *bounds, int64_t size, int32_t isr, int32_t osr)
EXTENSION_NOINLINE void ST_Centroid_LineString(int8_t *coords, int32_t coords_sz, int32_t ic, int32_t isr, int32_t osr, double *linestring_centroid)
EXTENSION_INLINE double ST_XMax_Bounds(double *bounds, int64_t size, int32_t isr, int32_t osr)
DEVICE ALWAYS_INLINE double transform_coord_x(const double coord, const int32_t isr, const int32_t osr)
EXTENSION_INLINE bool ST_DWithin_LineString_MultiPolygon(int8_t *l1, int64_t l1size, double *l1bounds, int64_t l1bounds_size, int8_t *mpoly_coords, int64_t mpoly_coords_size, int32_t *mpoly_ring_sizes, int64_t mpoly_num_rings, int32_t *mpoly_poly_sizes, int64_t mpoly_num_polys, double *mpoly_bounds, int64_t mpoly_bounds_size, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr, double distance_within)
EXTENSION_INLINE bool ST_DWithin_Point_MultiPolygon(int8_t *p, int64_t psize, int8_t *mpoly_coords, int64_t mpoly_coords_size, int32_t *mpoly_ring_sizes, int64_t mpoly_num_rings, int32_t *mpoly_poly_sizes, int64_t mpoly_num_polys, double *mpoly_bounds, int64_t mpoly_bounds_size, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr, double distance_within)
EXTENSION_NOINLINE bool ST_Contains_Point_LineString(int8_t *p, int64_t psize, int8_t *l, int64_t lsize, double *lbounds, int64_t lbounds_size, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr)
DEVICE ALWAYS_INLINE double distance_point_point(double p1x, double p1y, double p2x, double p2y)
DEVICE ALWAYS_INLINE bool centroid_add_ring(int8_t *ring, int64_t ringsize, int32_t ic, int32_t isr, int32_t osr, double sign, double *total_area2, double *cg3, double *total_length, double *linestring_centroid_sum, int64_t *num_points, double *point_centroid_sum)