OmniSciDB  a667adc9c8
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros 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_X_LineString(int8_t* l,
946  int64_t lsize,
947  int32_t lindex,
948  int32_t ic,
949  int32_t isr,
950  int32_t osr) {
951  auto l_num_points = lsize / (2 * compression_unit_size(ic));
952  if (lindex < 0 || lindex > l_num_points) {
953  lindex = l_num_points; // Endpoint
954  }
955  return coord_x(l, 2 * (lindex - 1), ic, isr, osr);
956 }
957 
959 double ST_Y_LineString(int8_t* l,
960  int64_t lsize,
961  int32_t lindex,
962  int32_t ic,
963  int32_t isr,
964  int32_t osr) {
965  auto l_num_points = lsize / (2 * compression_unit_size(ic));
966  if (lindex < 0 || lindex > l_num_points) {
967  lindex = l_num_points; // Endpoint
968  }
969  return coord_y(l, 2 * (lindex - 1) + 1, ic, isr, osr);
970 }
971 
973 double ST_XMin(int8_t* coords, int64_t size, int32_t ic, int32_t isr, int32_t osr) {
974  auto num_coords = size / compression_unit_size(ic);
975  double xmin = 0.0;
976  for (int32_t i = 0; i < num_coords; i += 2) {
977  double x = coord_x(coords, i, ic, isr, osr);
978  if (i == 0 || x < xmin) {
979  xmin = x;
980  }
981  }
982  return xmin;
983 }
984 
986 double ST_YMin(int8_t* coords, int64_t size, int32_t ic, int32_t isr, int32_t osr) {
987  auto num_coords = size / compression_unit_size(ic);
988  double ymin = 0.0;
989  for (int32_t i = 1; i < num_coords; i += 2) {
990  double y = coord_y(coords, i, ic, isr, osr);
991  if (i == 1 || y < ymin) {
992  ymin = y;
993  }
994  }
995  return ymin;
996 }
997 
999 double ST_XMax(int8_t* coords, int64_t size, int32_t ic, int32_t isr, int32_t osr) {
1000  auto num_coords = size / compression_unit_size(ic);
1001  double xmax = 0.0;
1002  for (int32_t i = 0; i < num_coords; i += 2) {
1003  double x = coord_x(coords, i, ic, isr, osr);
1004  if (i == 0 || x > xmax) {
1005  xmax = x;
1006  }
1007  }
1008  return xmax;
1009 }
1010 
1012 double ST_YMax(int8_t* coords, int64_t size, int32_t ic, int32_t isr, int32_t osr) {
1013  auto num_coords = size / compression_unit_size(ic);
1014  double ymax = 0.0;
1015  for (int32_t i = 1; i < num_coords; i += 2) {
1016  double y = coord_y(coords, i, ic, isr, osr);
1017  if (i == 1 || y > ymax) {
1018  ymax = y;
1019  }
1020  }
1021  return ymax;
1022 }
1023 
1025 double ST_XMin_Bounds(double* bounds, int64_t size, int32_t isr, int32_t osr) {
1026  return transform_coord(bounds[0], isr, osr, true);
1027 }
1028 
1030 double ST_YMin_Bounds(double* bounds, int64_t size, int32_t isr, int32_t osr) {
1031  return transform_coord(bounds[1], isr, osr, false);
1032 }
1033 
1035 double ST_XMax_Bounds(double* bounds, int64_t size, int32_t isr, int32_t osr) {
1036  return transform_coord(bounds[2], isr, osr, true);
1037 }
1038 
1040 double ST_YMax_Bounds(double* bounds, int64_t size, int32_t isr, int32_t osr) {
1041  return transform_coord(bounds[3], isr, osr, false);
1042 }
1043 
1044 //
1045 // ST_Length
1046 //
1047 
1049  int64_t lsize,
1050  int32_t ic,
1051  int32_t isr,
1052  int32_t osr,
1053  bool geodesic,
1054  bool check_closed) {
1055  auto l_num_coords = lsize / compression_unit_size(ic);
1056 
1057  double length = 0.0;
1058 
1059  double l0x = coord_x(l, 0, ic, isr, osr);
1060  double l0y = coord_y(l, 1, ic, isr, osr);
1061  double l2x = l0x;
1062  double l2y = l0y;
1063  for (int32_t i = 2; i < l_num_coords; i += 2) {
1064  double l1x = l2x;
1065  double l1y = l2y;
1066  l2x = coord_x(l, i, ic, isr, osr);
1067  l2y = coord_y(l, i + 1, ic, isr, osr);
1068  double ldist = geodesic ? distance_in_meters(l1x, l1y, l2x, l2y)
1069  : distance_point_point(l1x, l1y, l2x, l2y);
1070  length += ldist;
1071  }
1072  if (check_closed) {
1073  double ldist = geodesic ? distance_in_meters(l2x, l2y, l0x, l0y)
1074  : distance_point_point(l2x, l2y, l0x, l0y);
1075  length += ldist;
1076  }
1077  return length;
1078 }
1079 
1081 double ST_Length_LineString(int8_t* coords,
1082  int64_t coords_sz,
1083  int32_t ic,
1084  int32_t isr,
1085  int32_t osr) {
1086  return length_linestring(coords, coords_sz, ic, isr, osr, false, false);
1087 }
1088 
1090 double ST_Length_LineString_Geodesic(int8_t* coords,
1091  int64_t coords_sz,
1092  int32_t ic,
1093  int32_t isr,
1094  int32_t osr) {
1095  return length_linestring(coords, coords_sz, ic, isr, osr, true, false);
1096 }
1097 
1098 //
1099 // ST_Perimeter
1100 //
1101 
1103 double ST_Perimeter_Polygon(int8_t* poly,
1104  int64_t polysize,
1105  int32_t* poly_ring_sizes,
1106  int64_t poly_num_rings,
1107  int32_t ic,
1108  int32_t isr,
1109  int32_t osr) {
1110  if (poly_num_rings <= 0) {
1111  return 0.0;
1112  }
1113 
1114  auto exterior_ring_num_coords = poly_ring_sizes[0] * 2;
1115  auto exterior_ring_coords_size = exterior_ring_num_coords * compression_unit_size(ic);
1116 
1117  return length_linestring(poly, exterior_ring_coords_size, ic, isr, osr, false, true);
1118 }
1119 
1121 double ST_Perimeter_Polygon_Geodesic(int8_t* poly,
1122  int64_t polysize,
1123  int32_t* poly_ring_sizes,
1124  int64_t poly_num_rings,
1125  int32_t ic,
1126  int32_t isr,
1127  int32_t osr) {
1128  if (poly_num_rings <= 0) {
1129  return 0.0;
1130  }
1131 
1132  auto exterior_ring_num_coords = poly_ring_sizes[0] * 2;
1133  auto exterior_ring_coords_size = exterior_ring_num_coords * compression_unit_size(ic);
1134 
1135  return length_linestring(poly, exterior_ring_coords_size, ic, isr, osr, true, true);
1136 }
1137 
1138 DEVICE ALWAYS_INLINE double perimeter_multipolygon(int8_t* mpoly_coords,
1139  int64_t mpoly_coords_size,
1140  int32_t* mpoly_ring_sizes,
1141  int64_t mpoly_num_rings,
1142  int32_t* mpoly_poly_sizes,
1143  int64_t mpoly_num_polys,
1144  int32_t ic,
1145  int32_t isr,
1146  int32_t osr,
1147  bool geodesic) {
1148  if (mpoly_num_polys <= 0 || mpoly_num_rings <= 0) {
1149  return 0.0;
1150  }
1151 
1152  double perimeter = 0.0;
1153 
1154  // Set specific poly pointers as we move through the coords/ringsizes/polyrings arrays.
1155  auto next_poly_coords = mpoly_coords;
1156  auto next_poly_ring_sizes = mpoly_ring_sizes;
1157 
1158  for (auto poly = 0; poly < mpoly_num_polys; poly++) {
1159  auto poly_coords = next_poly_coords;
1160  auto poly_ring_sizes = next_poly_ring_sizes;
1161  auto poly_num_rings = mpoly_poly_sizes[poly];
1162  // Count number of coords in all of poly's rings, advance ring size pointer.
1163  int32_t poly_num_coords = 0;
1164  for (auto ring = 0; ring < poly_num_rings; ring++) {
1165  poly_num_coords += 2 * *next_poly_ring_sizes++;
1166  }
1167  auto poly_coords_size = poly_num_coords * compression_unit_size(ic);
1168  next_poly_coords += poly_coords_size;
1169 
1170  auto exterior_ring_num_coords = poly_ring_sizes[0] * 2;
1171  auto exterior_ring_coords_size = exterior_ring_num_coords * compression_unit_size(ic);
1172 
1173  perimeter += length_linestring(
1174  poly_coords, exterior_ring_coords_size, ic, isr, osr, geodesic, true);
1175  }
1176 
1177  return perimeter;
1178 }
1179 
1181 double ST_Perimeter_MultiPolygon(int8_t* mpoly_coords,
1182  int64_t mpoly_coords_size,
1183  int32_t* mpoly_ring_sizes,
1184  int64_t mpoly_num_rings,
1185  int32_t* mpoly_poly_sizes,
1186  int64_t mpoly_num_polys,
1187  int32_t ic,
1188  int32_t isr,
1189  int32_t osr) {
1190  return perimeter_multipolygon(mpoly_coords,
1191  mpoly_coords_size,
1192  mpoly_ring_sizes,
1193  mpoly_num_rings,
1194  mpoly_poly_sizes,
1195  mpoly_num_polys,
1196  ic,
1197  isr,
1198  osr,
1199  false);
1200 }
1201 
1203 double ST_Perimeter_MultiPolygon_Geodesic(int8_t* mpoly_coords,
1204  int64_t mpoly_coords_size,
1205  int32_t* mpoly_ring_sizes,
1206  int64_t mpoly_num_rings,
1207  int32_t* mpoly_poly_sizes,
1208  int64_t mpoly_num_polys,
1209  int32_t ic,
1210  int32_t isr,
1211  int32_t osr) {
1212  return perimeter_multipolygon(mpoly_coords,
1213  mpoly_coords_size,
1214  mpoly_ring_sizes,
1215  mpoly_num_rings,
1216  mpoly_poly_sizes,
1217  mpoly_num_polys,
1218  ic,
1219  isr,
1220  osr,
1221  true);
1222 }
1223 
1224 //
1225 // ST_Area
1226 //
1227 
1229  double y1,
1230  double x2,
1231  double y2,
1232  double x3,
1233  double y3) {
1234  return (x1 * y2 - x2 * y1 + x3 * y1 - x1 * y3 + x2 * y3 - x3 * y2) / 2.0;
1235 }
1236 
1237 DEVICE ALWAYS_INLINE double area_ring(int8_t* ring,
1238  int64_t ringsize,
1239  int32_t ic,
1240  int32_t isr,
1241  int32_t osr) {
1242  auto ring_num_coords = ringsize / compression_unit_size(ic);
1243 
1244  if (ring_num_coords < 6) {
1245  return 0.0;
1246  }
1247 
1248  double area = 0.0;
1249 
1250  double x1 = coord_x(ring, 0, ic, isr, osr);
1251  double y1 = coord_y(ring, 1, ic, isr, osr);
1252  double x2 = coord_x(ring, 2, ic, isr, osr);
1253  double y2 = coord_y(ring, 3, ic, isr, osr);
1254  for (int32_t i = 4; i < ring_num_coords; i += 2) {
1255  double x3 = coord_x(ring, i, ic, isr, osr);
1256  double y3 = coord_y(ring, i + 1, ic, isr, osr);
1257  area += area_triangle(x1, y1, x2, y2, x3, y3);
1258  x2 = x3;
1259  y2 = y3;
1260  }
1261  return area;
1262 }
1263 
1264 DEVICE ALWAYS_INLINE double area_polygon(int8_t* poly_coords,
1265  int64_t poly_coords_size,
1266  int32_t* poly_ring_sizes,
1267  int64_t poly_num_rings,
1268  int32_t ic,
1269  int32_t isr,
1270  int32_t osr) {
1271  if (poly_num_rings <= 0) {
1272  return 0.0;
1273  }
1274 
1275  double area = 0.0;
1276  auto ring_coords = poly_coords;
1277 
1278  // Add up the areas of all rings.
1279  // External ring is CCW, open - positive area.
1280  // Internal rings (holes) are CW, open - negative areas.
1281  for (auto r = 0; r < poly_num_rings; r++) {
1282  auto ring_coords_size = poly_ring_sizes[r] * 2 * compression_unit_size(ic);
1283  area += area_ring(ring_coords, ring_coords_size, ic, isr, osr);
1284  // Advance to the next ring.
1285  ring_coords += ring_coords_size;
1286  }
1287  return area;
1288 }
1289 
1291 double ST_Area_Polygon(int8_t* poly_coords,
1292  int64_t poly_coords_size,
1293  int32_t* poly_ring_sizes,
1294  int64_t poly_num_rings,
1295  int32_t ic,
1296  int32_t isr,
1297  int32_t osr) {
1298  return area_polygon(
1299  poly_coords, poly_coords_size, poly_ring_sizes, poly_num_rings, ic, isr, osr);
1300 }
1301 
1303 double ST_Area_Polygon_Geodesic(int8_t* poly_coords,
1304  int64_t poly_coords_size,
1305  int32_t* poly_ring_sizes,
1306  int64_t poly_num_rings,
1307  int32_t ic,
1308  int32_t isr,
1309  int32_t osr) {
1310  return ST_Area_Polygon(
1311  poly_coords, poly_coords_size, poly_ring_sizes, poly_num_rings, ic, isr, osr);
1312 }
1313 
1315 double ST_Area_MultiPolygon(int8_t* mpoly_coords,
1316  int64_t mpoly_coords_size,
1317  int32_t* mpoly_ring_sizes,
1318  int64_t mpoly_num_rings,
1319  int32_t* mpoly_poly_sizes,
1320  int64_t mpoly_num_polys,
1321  int32_t ic,
1322  int32_t isr,
1323  int32_t osr) {
1324  if (mpoly_num_rings <= 0 || mpoly_num_polys <= 0) {
1325  return 0.0;
1326  }
1327 
1328  double area = 0.0;
1329 
1330  // Set specific poly pointers as we move through the coords/ringsizes/polyrings arrays.
1331  auto next_poly_coords = mpoly_coords;
1332  auto next_poly_ring_sizes = mpoly_ring_sizes;
1333 
1334  for (auto poly = 0; poly < mpoly_num_polys; poly++) {
1335  auto poly_coords = next_poly_coords;
1336  auto poly_ring_sizes = next_poly_ring_sizes;
1337  auto poly_num_rings = mpoly_poly_sizes[poly];
1338  // Count number of coords in all of poly's rings, advance ring size pointer.
1339  int32_t poly_num_coords = 0;
1340  for (auto ring = 0; ring < poly_num_rings; ring++) {
1341  poly_num_coords += 2 * *next_poly_ring_sizes++;
1342  }
1343  auto poly_coords_size = poly_num_coords * compression_unit_size(ic);
1344  next_poly_coords += poly_coords_size;
1345 
1346  area += area_polygon(
1347  poly_coords, poly_coords_size, poly_ring_sizes, poly_num_rings, ic, isr, osr);
1348  }
1349  return area;
1350 }
1351 
1353 double ST_Area_MultiPolygon_Geodesic(int8_t* mpoly_coords,
1354  int64_t mpoly_coords_size,
1355  int32_t* mpoly_ring_sizes,
1356  int64_t mpoly_num_rings,
1357  int32_t* mpoly_poly_sizes,
1358  int64_t mpoly_num_polys,
1359  int32_t ic,
1360  int32_t isr,
1361  int32_t osr) {
1362  return ST_Area_MultiPolygon(mpoly_coords,
1363  mpoly_coords_size,
1364  mpoly_ring_sizes,
1365  mpoly_num_rings,
1366  mpoly_poly_sizes,
1367  mpoly_num_polys,
1368  ic,
1369  isr,
1370  osr);
1371 }
1372 
1373 //
1374 // ST_Centroid
1375 //
1376 
1377 // Point centroid
1378 // - calculate average coordinates for all points
1379 // LineString centroid
1380 // - take midpoint of each line segment
1381 // - calculate average coordinates of all midpoints weighted by segment length
1382 // Polygon, MultiPolygon centroid, holes..
1383 // - weighted sum of centroids of triangles coming out of area decomposition
1384 //
1385 // On zero area fall back to linestring centroid
1386 // On zero length fall back to point centroid
1387 
1389 double ST_Centroid_Point(int8_t* p,
1390  int64_t psize,
1391  int32_t ic,
1392  int32_t isr,
1393  int32_t osr,
1394  bool ycoord) {
1395  if (ycoord) {
1396  return coord_y(p, 1, ic, isr, osr);
1397  }
1398  return coord_x(p, 0, ic, isr, osr);
1399 }
1400 
1402  double y1,
1403  double x2,
1404  double y2,
1405  double* length,
1406  double* linestring_centroid_sum) {
1407  double ldist = distance_point_point(x1, y1, x2, y2);
1408  *length += ldist;
1409  double segment_midpoint_x = (x1 + x2) / 2.0;
1410  double segment_midpoint_y = (y1 + y2) / 2.0;
1411  linestring_centroid_sum[0] += ldist * segment_midpoint_x;
1412  linestring_centroid_sum[1] += ldist * segment_midpoint_y;
1413  return true;
1414 }
1415 
1417  int64_t lsize,
1418  int32_t ic,
1419  int32_t isr,
1420  int32_t osr,
1421  bool closed,
1422  double* total_length,
1423  double* linestring_centroid_sum,
1424  int64_t* num_points,
1425  double* point_centroid_sum) {
1426  auto l_num_coords = lsize / compression_unit_size(ic);
1427  double length = 0.0;
1428  double l0x = coord_x(l, 0, ic, isr, osr);
1429  double l0y = coord_y(l, 1, ic, isr, osr);
1430  double l2x = l0x;
1431  double l2y = l0y;
1432  for (int32_t i = 2; i < l_num_coords; i += 2) {
1433  double l1x = l2x;
1434  double l1y = l2y;
1435  l2x = coord_x(l, i, ic, isr, osr);
1436  l2y = coord_y(l, i + 1, ic, isr, osr);
1437  centroid_add_segment(l1x, l1y, l2x, l2y, &length, linestring_centroid_sum);
1438  }
1439  if (l_num_coords > 4 && closed) {
1440  // Also add the closing segment between the last and the first points
1441  centroid_add_segment(l2x, l2y, l0x, l0y, &length, linestring_centroid_sum);
1442  }
1443  *total_length += length;
1444  if (length == 0.0 && l_num_coords > 0) {
1445  *num_points += 1;
1446  point_centroid_sum[0] += l0x;
1447  point_centroid_sum[1] += l0y;
1448  }
1449  return true;
1450 }
1451 
1453 double ST_Centroid_LineString(int8_t* coords,
1454  int64_t coords_sz,
1455  int32_t ic,
1456  int32_t isr,
1457  int32_t osr,
1458  bool ycoord) {
1459  double length = 0.0;
1460  double linestring_centroid_sum[2] = {0.0, 0.0};
1461  int64_t num_points = 0;
1462  double point_centroid_sum[2] = {0.0, 0.0};
1463  centroid_add_linestring(coords,
1464  coords_sz,
1465  ic,
1466  isr,
1467  osr,
1468  false, // not closed
1469  &length,
1470  &linestring_centroid_sum[0],
1471  &num_points,
1472  &point_centroid_sum[0]);
1473  double linestring_centroid[2] = {0.0, 0.0};
1474  if (length > 0) {
1475  linestring_centroid[0] = linestring_centroid_sum[0] / length;
1476  linestring_centroid[1] = linestring_centroid_sum[1] / length;
1477  } else if (num_points > 0) {
1478  linestring_centroid[0] = point_centroid_sum[0] / num_points;
1479  linestring_centroid[1] = point_centroid_sum[1] / num_points;
1480  }
1481  if (ycoord) {
1482  return linestring_centroid[1];
1483  }
1484  return linestring_centroid[0];
1485 }
1486 
1488  double y1,
1489  double x2,
1490  double y2,
1491  double x3,
1492  double y3,
1493  double sign,
1494  double* total_area2,
1495  double* cg3) {
1496  double cx = x1 + x2 + x3;
1497  double cy = y1 + y2 + y3;
1498  double area2 = x1 * y2 - x2 * y1 + x3 * y1 - x1 * y3 + x2 * y3 - x3 * y2;
1499  cg3[0] += sign * area2 * cx;
1500  cg3[1] += sign * area2 * cy;
1501  *total_area2 += sign * area2;
1502  return true;
1503 }
1504 
1506  int64_t ringsize,
1507  int32_t ic,
1508  int32_t isr,
1509  int32_t osr,
1510  double sign,
1511  double* total_area2,
1512  double* cg3,
1513  double* total_length,
1514  double* linestring_centroid_sum,
1515  int64_t* num_points,
1516  double* point_centroid_sum) {
1517  auto ring_num_coords = ringsize / compression_unit_size(ic);
1518 
1519  if (ring_num_coords < 6) {
1520  return false;
1521  }
1522 
1523  double x1 = coord_x(ring, 0, ic, isr, osr);
1524  double y1 = coord_y(ring, 1, ic, isr, osr);
1525  double x2 = coord_x(ring, 2, ic, isr, osr);
1526  double y2 = coord_y(ring, 3, ic, isr, osr);
1527  for (int32_t i = 4; i < ring_num_coords; i += 2) {
1528  double x3 = coord_x(ring, i, ic, isr, osr);
1529  double y3 = coord_y(ring, i + 1, ic, isr, osr);
1530  centroid_add_triangle(x1, y1, x2, y2, x3, y3, sign, total_area2, cg3);
1531  x2 = x3;
1532  y2 = y3;
1533  }
1534 
1536  ringsize,
1537  ic,
1538  isr,
1539  osr,
1540  true, // closed
1541  total_length,
1542  linestring_centroid_sum,
1543  num_points,
1544  point_centroid_sum);
1545  return true;
1546 }
1547 
1548 DEVICE ALWAYS_INLINE bool centroid_add_polygon(int8_t* poly_coords,
1549  int64_t poly_coords_size,
1550  int32_t* poly_ring_sizes,
1551  int64_t poly_num_rings,
1552  int32_t ic,
1553  int32_t isr,
1554  int32_t osr,
1555  double* total_area2,
1556  double* cg3,
1557  double* total_length,
1558  double* linestring_centroid_sum,
1559  int64_t* num_points,
1560  double* point_centroid_sum) {
1561  if (poly_num_rings <= 0) {
1562  return false;
1563  }
1564 
1565  auto ring_coords = poly_coords;
1566 
1567  for (auto r = 0; r < poly_num_rings; r++) {
1568  auto ring_coords_size = poly_ring_sizes[r] * 2 * compression_unit_size(ic);
1569  // Shape - positive area, holes - negative areas
1570  double sign = (r == 0) ? 1.0 : -1.0;
1571  centroid_add_ring(ring_coords,
1572  ring_coords_size,
1573  ic,
1574  isr,
1575  osr,
1576  sign,
1577  total_area2,
1578  cg3,
1579  total_length,
1580  linestring_centroid_sum,
1581  num_points,
1582  point_centroid_sum);
1583  // Advance to the next ring.
1584  ring_coords += ring_coords_size;
1585  }
1586  return true;
1587 }
1588 
1590 double ST_Centroid_Polygon(int8_t* poly_coords,
1591  int64_t poly_coords_size,
1592  int32_t* poly_ring_sizes,
1593  int64_t poly_num_rings,
1594  int32_t ic,
1595  int32_t isr,
1596  int32_t osr,
1597  bool ycoord) {
1598  if (poly_num_rings <= 0) {
1599  return 0.0;
1600  }
1601  double total_area2 = 0.0;
1602  double cg3[2] = {0.0, 0.0};
1603  double total_length = 0.0;
1604  double linestring_centroid_sum[2] = {0.0, 0.0};
1605  int64_t num_points = 0;
1606  double point_centroid_sum[2] = {0.0, 0.0};
1607  centroid_add_polygon(poly_coords,
1608  poly_coords_size,
1609  poly_ring_sizes,
1610  poly_num_rings,
1611  ic,
1612  isr,
1613  osr,
1614  &total_area2,
1615  &cg3[0],
1616  &total_length,
1617  &linestring_centroid_sum[0],
1618  &num_points,
1619  &point_centroid_sum[0]);
1620 
1621  double x1 = coord_x(poly_coords, 0, ic, isr, osr);
1622  double y1 = coord_y(poly_coords, 1, ic, isr, osr);
1623  double poly_centroid[2] = {x1, y1};
1624  if (total_area2 != 0.0) {
1625  poly_centroid[0] = cg3[0] / 3 / total_area2;
1626  poly_centroid[1] = cg3[1] / 3 / total_area2;
1627  } else if (total_length > 0.0) {
1628  // zero-area polygon, fall back to linear centroid
1629  poly_centroid[0] = linestring_centroid_sum[0] / total_length;
1630  poly_centroid[1] = linestring_centroid_sum[1] / total_length;
1631  } else if (num_points > 0) {
1632  // zero-area zero-length polygon, fall further back to point centroid
1633  poly_centroid[0] = point_centroid_sum[0] / num_points;
1634  poly_centroid[1] = point_centroid_sum[1] / num_points;
1635  }
1636 
1637  if (ycoord) {
1638  return poly_centroid[1];
1639  }
1640  return poly_centroid[0];
1641 }
1642 
1644 double ST_Centroid_MultiPolygon(int8_t* mpoly_coords,
1645  int64_t mpoly_coords_size,
1646  int32_t* mpoly_ring_sizes,
1647  int64_t mpoly_num_rings,
1648  int32_t* mpoly_poly_sizes,
1649  int64_t mpoly_num_polys,
1650  int32_t ic,
1651  int32_t isr,
1652  int32_t osr,
1653  bool ycoord) {
1654  if (mpoly_num_rings <= 0 || mpoly_num_polys <= 0) {
1655  return 0.0;
1656  }
1657 
1658  double total_area2 = 0.0;
1659  double cg3[2] = {0.0, 0.0};
1660  double total_length = 0.0;
1661  double linestring_centroid_sum[2] = {0.0, 0.0};
1662  int64_t num_points = 0;
1663  double point_centroid_sum[2] = {0.0, 0.0};
1664 
1665  // Set specific poly pointers as we move through the coords/ringsizes/polyrings arrays.
1666  auto next_poly_coords = mpoly_coords;
1667  auto next_poly_ring_sizes = mpoly_ring_sizes;
1668 
1669  for (auto poly = 0; poly < mpoly_num_polys; poly++) {
1670  auto poly_coords = next_poly_coords;
1671  auto poly_ring_sizes = next_poly_ring_sizes;
1672  auto poly_num_rings = mpoly_poly_sizes[poly];
1673  // Count number of coords in all of poly's rings, advance ring size pointer.
1674  int32_t poly_num_coords = 0;
1675  for (auto ring = 0; ring < poly_num_rings; ring++) {
1676  poly_num_coords += 2 * *next_poly_ring_sizes++;
1677  }
1678  auto poly_coords_size = poly_num_coords * compression_unit_size(ic);
1679  next_poly_coords += poly_coords_size;
1680 
1681  centroid_add_polygon(poly_coords,
1682  poly_coords_size,
1683  poly_ring_sizes,
1684  poly_num_rings,
1685  ic,
1686  isr,
1687  osr,
1688  &total_area2,
1689  &cg3[0],
1690  &total_length,
1691  &linestring_centroid_sum[0],
1692  &num_points,
1693  &point_centroid_sum[0]);
1694  }
1695 
1696  double x1 = coord_x(mpoly_coords, 0, ic, isr, osr);
1697  double y1 = coord_y(mpoly_coords, 1, ic, isr, osr);
1698  double mpoly_centroid[2] = {x1, y1};
1699  if (total_area2 != 0.0) {
1700  mpoly_centroid[0] = cg3[0] / 3 / total_area2;
1701  mpoly_centroid[1] = cg3[1] / 3 / total_area2;
1702  } else if (total_length > 0.0) {
1703  // zero-area multipolygon, fall back to linear centroid
1704  mpoly_centroid[0] = linestring_centroid_sum[0] / total_length;
1705  mpoly_centroid[1] = linestring_centroid_sum[1] / total_length;
1706  } else if (num_points > 0) {
1707  // zero-area zero-length multipolygon, fall further back to point centroid
1708  mpoly_centroid[0] = point_centroid_sum[0] / num_points;
1709  mpoly_centroid[1] = point_centroid_sum[1] / num_points;
1710  }
1711 
1712  if (ycoord) {
1713  return mpoly_centroid[1];
1714  }
1715  return mpoly_centroid[0];
1716 }
1717 
1719 int32_t ST_NPoints(int8_t* coords, int64_t coords_sz, int32_t ic) {
1720  auto num_pts = coords_sz / compression_unit_size(ic);
1721  return static_cast<int32_t>(num_pts / 2);
1722 }
1723 
1725 int32_t ST_NRings(int32_t* poly_ring_sizes, int64_t poly_num_rings) {
1726  return static_cast<int32_t>(poly_num_rings);
1727 }
1728 
1729 //
1730 // ST_Distance
1731 //
1732 
1734 double ST_Distance_Point_Point(int8_t* p1,
1735  int64_t p1size,
1736  int8_t* p2,
1737  int64_t p2size,
1738  int32_t ic1,
1739  int32_t isr1,
1740  int32_t ic2,
1741  int32_t isr2,
1742  int32_t osr) {
1743  double p1x = coord_x(p1, 0, ic1, isr1, osr);
1744  double p1y = coord_y(p1, 1, ic1, isr1, osr);
1745  double p2x = coord_x(p2, 0, ic2, isr2, osr);
1746  double p2y = coord_y(p2, 1, ic2, isr2, osr);
1747  return distance_point_point(p1x, p1y, p2x, p2y);
1748 }
1749 
1752  int64_t p1size,
1753  int8_t* p2,
1754  int64_t p2size,
1755  int32_t ic1,
1756  int32_t isr1,
1757  int32_t ic2,
1758  int32_t isr2,
1759  int32_t osr) {
1760  double p1x = coord_x(p1, 0, ic1, isr1, osr);
1761  double p1y = coord_y(p1, 1, ic1, isr1, osr);
1762  double p2x = coord_x(p2, 0, ic2, isr2, osr);
1763  double p2y = coord_y(p2, 1, ic2, isr2, osr);
1764  return distance_point_point_squared(p1x, p1y, p2x, p2y);
1765 }
1766 
1769  int64_t p1size,
1770  int8_t* p2,
1771  int64_t p2size,
1772  int32_t ic1,
1773  int32_t isr1,
1774  int32_t ic2,
1775  int32_t isr2,
1776  int32_t osr) {
1777  double p1x = coord_x(p1, 0, ic1, 4326, 4326);
1778  double p1y = coord_y(p1, 1, ic1, 4326, 4326);
1779  double p2x = coord_x(p2, 0, ic2, 4326, 4326);
1780  double p2y = coord_y(p2, 1, ic2, 4326, 4326);
1781  return distance_in_meters(p1x, p1y, p2x, p2y);
1782 }
1783 
1786  int64_t psize,
1787  int8_t* l,
1788  int64_t lsize,
1789  int32_t lindex,
1790  int32_t ic1,
1791  int32_t isr1,
1792  int32_t ic2,
1793  int32_t isr2,
1794  int32_t osr) {
1795  // Currently only statically indexed LineString is supported
1796  double px = coord_x(p, 0, ic1, 4326, 4326);
1797  double py = coord_y(p, 1, ic1, 4326, 4326);
1798  auto lpoints = lsize / (2 * compression_unit_size(ic2));
1799  if (lindex < 0 || lindex > lpoints) {
1800  lindex = lpoints; // Endpoint
1801  }
1802  double lx = coord_x(l, 2 * (lindex - 1), ic2, 4326, 4326);
1803  double ly = coord_y(l, 2 * (lindex - 1) + 1, ic2, 4326, 4326);
1804  return distance_in_meters(px, py, lx, ly);
1805 }
1806 
1809  int64_t lsize,
1810  int32_t lindex,
1811  int8_t* p,
1812  int64_t psize,
1813  int32_t ic1,
1814  int32_t isr1,
1815  int32_t ic2,
1816  int32_t isr2,
1817  int32_t osr) {
1818  // Currently only statically indexed LineString is supported
1820  p, psize, l, lsize, lindex, ic2, isr2, ic1, isr1, osr);
1821 }
1822 
1825  int64_t l1size,
1826  int32_t l1index,
1827  int8_t* l2,
1828  int64_t l2size,
1829  int32_t l2index,
1830  int32_t ic1,
1831  int32_t isr1,
1832  int32_t ic2,
1833  int32_t isr2,
1834  int32_t osr) {
1835  // Currently only statically indexed LineStrings are supported
1836  auto l1points = l1size / (2 * compression_unit_size(ic1));
1837  if (l1index < 0 || l1index > l1points) {
1838  l1index = l1points; // Endpoint
1839  }
1840  double l1x = coord_x(l1, 2 * (l1index - 1), ic1, 4326, 4326);
1841  double l1y = coord_y(l1, 2 * (l1index - 1) + 1, ic1, 4326, 4326);
1842  auto l2points = l2size / (2 * compression_unit_size(ic2));
1843  if (l2index < 0 || l2index > l2points) {
1844  l2index = l2points; // Endpoint
1845  }
1846  double l2x = coord_x(l2, 2 * (l2index - 1), ic2, 4326, 4326);
1847  double l2y = coord_y(l2, 2 * (l2index - 1) + 1, ic2, 4326, 4326);
1848  return distance_in_meters(l1x, l1y, l2x, l2y);
1849 }
1850 
1852  int64_t psize,
1853  int8_t* l,
1854  int64_t lsize,
1855  int32_t lindex,
1856  int32_t ic1,
1857  int32_t isr1,
1858  int32_t ic2,
1859  int32_t isr2,
1860  int32_t osr,
1861  bool check_closed,
1862  double threshold) {
1863  double px = coord_x(p, 0, ic1, isr1, osr);
1864  double py = coord_y(p, 1, ic1, isr1, osr);
1865 
1866  auto l_num_coords = lsize / compression_unit_size(ic2);
1867  auto l_num_points = l_num_coords / 2;
1868  if (lindex != 0) { // Statically indexed linestring
1869  if (lindex < 0 || lindex > l_num_points) {
1870  lindex = l_num_points; // Endpoint
1871  }
1872  double lx = coord_x(l, 2 * (lindex - 1), ic2, isr2, osr);
1873  double ly = coord_y(l, 2 * (lindex - 1) + 1, ic2, isr2, osr);
1874  return distance_point_point(px, py, lx, ly);
1875  }
1876 
1877  double l1x = coord_x(l, 0, ic2, isr2, osr);
1878  double l1y = coord_y(l, 1, ic2, isr2, osr);
1879  double l2x = coord_x(l, 2, ic2, isr2, osr);
1880  double l2y = coord_y(l, 3, ic2, isr2, osr);
1881 
1882  double dist = distance_point_line(px, py, l1x, l1y, l2x, l2y);
1883  for (int32_t i = 4; i < l_num_coords; i += 2) {
1884  l1x = l2x; // advance one point
1885  l1y = l2y;
1886  l2x = coord_x(l, i, ic2, isr2, osr);
1887  l2y = coord_y(l, i + 1, ic2, isr2, osr);
1888  double ldist = distance_point_line(px, py, l1x, l1y, l2x, l2y);
1889  if (dist > ldist) {
1890  dist = ldist;
1891  }
1892  if (dist <= threshold) {
1893  return dist;
1894  }
1895  }
1896  if (l_num_coords > 4 && check_closed) {
1897  // Also check distance to the closing edge between the first and the last points
1898  l1x = coord_x(l, 0, ic2, isr2, osr);
1899  l1y = coord_y(l, 1, ic2, isr2, osr);
1900  double ldist = distance_point_line(px, py, l1x, l1y, l2x, l2y);
1901  if (dist > ldist) {
1902  dist = ldist;
1903  }
1904  }
1905  return dist;
1906 }
1907 
1910  int64_t psize,
1911  int8_t* l,
1912  int64_t lsize,
1913  int32_t lindex,
1914  int32_t ic1,
1915  int32_t isr1,
1916  int32_t ic2,
1917  int32_t isr2,
1918  int32_t osr,
1919  double threshold) {
1921  p, psize, l, lsize, lindex, ic1, isr1, ic2, isr2, osr, true, threshold);
1922 }
1923 
1926  int64_t psize,
1927  int8_t* l,
1928  int64_t lsize,
1929  int32_t lindex,
1930  int32_t ic1,
1931  int32_t isr1,
1932  int32_t ic2,
1933  int32_t isr2,
1934  int32_t osr,
1935  double threshold) {
1936  if (lindex != 0) { // Statically indexed linestring
1937  auto l_num_coords = lsize / compression_unit_size(ic2);
1938  auto l_num_points = l_num_coords / 2;
1939  if (lindex < 0 || lindex > l_num_points) {
1940  lindex = l_num_points; // Endpoint
1941  }
1942  double px = coord_x(p, 0, ic1, isr1, osr);
1943  double py = coord_y(p, 1, ic1, isr1, osr);
1944  double lx = coord_x(l, 2 * (lindex - 1), ic2, isr2, osr);
1945  double ly = coord_y(l, 2 * (lindex - 1) + 1, ic2, isr2, osr);
1946  return distance_point_point(px, py, lx, ly);
1947  }
1948 
1950  p, psize, l, lsize, lindex, ic1, isr1, ic2, isr2, osr, false, threshold);
1951 }
1952 
1954 double ST_Distance_Point_Polygon(int8_t* p,
1955  int64_t psize,
1956  int8_t* poly,
1957  int64_t polysize,
1958  int32_t* poly_ring_sizes,
1959  int64_t poly_num_rings,
1960  int32_t ic1,
1961  int32_t isr1,
1962  int32_t ic2,
1963  int32_t isr2,
1964  int32_t osr,
1965  double threshold) {
1966  auto exterior_ring_num_coords = polysize / compression_unit_size(ic2);
1967  if (poly_num_rings > 0) {
1968  exterior_ring_num_coords = poly_ring_sizes[0] * 2;
1969  }
1970  auto exterior_ring_coords_size = exterior_ring_num_coords * compression_unit_size(ic2);
1971 
1972  double px = coord_x(p, 0, ic1, isr1, osr);
1973  double py = coord_y(p, 1, ic1, isr1, osr);
1974  if (!polygon_contains_point(poly, exterior_ring_num_coords, px, py, ic2, isr2, osr)) {
1975  // Outside the exterior ring
1977  psize,
1978  poly,
1979  exterior_ring_coords_size,
1980  0,
1981  ic1,
1982  isr1,
1983  ic2,
1984  isr2,
1985  osr,
1986  threshold);
1987  }
1988  // Inside exterior ring
1989  // Advance to first interior ring
1990  poly += exterior_ring_coords_size;
1991  // Check if one of the polygon's holes contains that point
1992  for (auto r = 1; r < poly_num_rings; r++) {
1993  auto interior_ring_num_coords = poly_ring_sizes[r] * 2;
1994  auto interior_ring_coords_size =
1995  interior_ring_num_coords * compression_unit_size(ic2);
1996  if (polygon_contains_point(poly, interior_ring_num_coords, px, py, ic2, isr2, osr)) {
1997  // Inside an interior ring
1999  psize,
2000  poly,
2001  interior_ring_coords_size,
2002  0,
2003  ic1,
2004  isr1,
2005  ic2,
2006  isr2,
2007  osr,
2008  threshold);
2009  }
2010  poly += interior_ring_coords_size;
2011  }
2012  return 0.0;
2013 }
2014 
2017  int64_t psize,
2018  int8_t* mpoly_coords,
2019  int64_t mpoly_coords_size,
2020  int32_t* mpoly_ring_sizes,
2021  int64_t mpoly_num_rings,
2022  int32_t* mpoly_poly_sizes,
2023  int64_t mpoly_num_polys,
2024  int32_t ic1,
2025  int32_t isr1,
2026  int32_t ic2,
2027  int32_t isr2,
2028  int32_t osr,
2029  double threshold) {
2030  if (mpoly_num_polys <= 0) {
2031  return 0.0;
2032  }
2033  double min_distance = 0.0;
2034 
2035  // Set specific poly pointers as we move through the coords/ringsizes/polyrings arrays.
2036  auto next_poly_coords = mpoly_coords;
2037  auto next_poly_ring_sizes = mpoly_ring_sizes;
2038 
2039  for (auto poly = 0; poly < mpoly_num_polys; poly++) {
2040  auto poly_coords = next_poly_coords;
2041  auto poly_ring_sizes = next_poly_ring_sizes;
2042  auto poly_num_rings = mpoly_poly_sizes[poly];
2043  // Count number of coords in all of poly's rings, advance ring size pointer.
2044  int32_t poly_num_coords = 0;
2045  for (auto ring = 0; ring < poly_num_rings; ring++) {
2046  poly_num_coords += 2 * *next_poly_ring_sizes++;
2047  }
2048  auto poly_coords_size = poly_num_coords * compression_unit_size(ic2);
2049  next_poly_coords += poly_coords_size;
2050  double distance = ST_Distance_Point_Polygon(p,
2051  psize,
2052  poly_coords,
2053  poly_coords_size,
2054  poly_ring_sizes,
2055  poly_num_rings,
2056  ic1,
2057  isr1,
2058  ic2,
2059  isr2,
2060  osr,
2061  threshold);
2062  if (poly == 0 || min_distance > distance) {
2063  min_distance = distance;
2064  if (tol_zero(min_distance)) {
2065  min_distance = 0.0;
2066  break;
2067  }
2068  if (min_distance <= threshold) {
2069  break;
2070  }
2071  }
2072  }
2073 
2074  return min_distance;
2075 }
2076 
2079  int64_t lsize,
2080  int32_t lindex,
2081  int8_t* p,
2082  int64_t psize,
2083  int32_t ic1,
2084  int32_t isr1,
2085  int32_t ic2,
2086  int32_t isr2,
2087  int32_t osr,
2088  double threshold) {
2090  p, psize, l, lsize, lindex, ic2, isr2, ic1, isr1, osr, threshold);
2091 }
2092 
2095  int64_t l1size,
2096  int32_t l1index,
2097  int8_t* l2,
2098  int64_t l2size,
2099  int32_t l2index,
2100  int32_t ic1,
2101  int32_t isr1,
2102  int32_t ic2,
2103  int32_t isr2,
2104  int32_t osr,
2105  double threshold) {
2106  auto l1_num_coords = l1size / compression_unit_size(ic1);
2107  if (l1index != 0) {
2108  // l1 is a statically indexed linestring
2109  auto l1_num_points = l1_num_coords / 2;
2110  if (l1index < 0 || l1index > l1_num_points) {
2111  l1index = l1_num_points;
2112  }
2113  int8_t* p = l1 + 2 * (l1index - 1) * compression_unit_size(ic1);
2114  int64_t psize = 2 * compression_unit_size(ic1);
2116  p, psize, l2, l2size, l2index, ic1, isr1, ic2, isr2, osr, threshold);
2117  }
2118 
2119  auto l2_num_coords = l2size / compression_unit_size(ic2);
2120  if (l2index != 0) {
2121  // l2 is a statically indexed linestring
2122  auto l2_num_points = l2_num_coords / 2;
2123  if (l2index < 0 || l2index > l2_num_points) {
2124  l2index = l2_num_points;
2125  }
2126  int8_t* p = l2 + 2 * (l2index - 1) * compression_unit_size(ic2);
2127  int64_t psize = 2 * compression_unit_size(ic2);
2129  p, psize, l1, l1size, l1index, ic2, isr2, ic1, isr1, osr, threshold);
2130  }
2131 
2132  double threshold_squared = threshold * threshold;
2133  double dist_squared = 0.0;
2134  double l11x = coord_x(l1, 0, ic1, isr1, osr);
2135  double l11y = coord_y(l1, 1, ic1, isr1, osr);
2136  for (int32_t i1 = 2; i1 < l1_num_coords; i1 += 2) {
2137  double l12x = coord_x(l1, i1, ic1, isr1, osr);
2138  double l12y = coord_y(l1, i1 + 1, ic1, isr1, osr);
2139 
2140  double l21x = coord_x(l2, 0, ic2, isr2, osr);
2141  double l21y = coord_y(l2, 1, ic2, isr2, osr);
2142  for (int32_t i2 = 2; i2 < l2_num_coords; i2 += 2) {
2143  double l22x = coord_x(l2, i2, ic2, isr2, osr);
2144  double l22y = coord_y(l2, i2 + 1, ic2, isr2, osr);
2145 
2146  // double ldist_squared =
2147  // distance_line_line_squared(l11x, l11y, l12x, l12y, l21x, l21y, l22x, l22y);
2148  // TODO: fix distance_line_line_squared
2149  double ldist = distance_line_line(l11x, l11y, l12x, l12y, l21x, l21y, l22x, l22y);
2150  double ldist_squared = ldist * ldist;
2151 
2152  if (i1 == 2 && i2 == 2) {
2153  dist_squared = ldist_squared; // initialize dist with distance between the first
2154  // two segments
2155  } else if (dist_squared > ldist_squared) {
2156  dist_squared = ldist_squared;
2157  }
2158  if (tol_zero(dist_squared, TOLERANCE_DEFAULT_SQUARED)) {
2159  return 0.0; // Segments touch, short-circuit
2160  }
2161  if (dist_squared <= threshold_squared) {
2162  // If threashold_squared is defined and the calculated dist_squared dips under it,
2163  // short-circuit and return it
2164  return sqrt(dist_squared);
2165  }
2166 
2167  l21x = l22x; // advance to the next point on l2
2168  l21y = l22y;
2169  }
2170 
2171  l11x = l12x; // advance to the next point on l1
2172  l11y = l12y;
2173  }
2174  return sqrt(dist_squared);
2175 }
2176 
2179  int64_t lsize,
2180  int32_t lindex,
2181  int8_t* poly_coords,
2182  int64_t poly_coords_size,
2183  int32_t* poly_ring_sizes,
2184  int64_t poly_num_rings,
2185  int32_t ic1,
2186  int32_t isr1,
2187  int32_t ic2,
2188  int32_t isr2,
2189  int32_t osr,
2190  double threshold) {
2191  auto lnum_coords = lsize / compression_unit_size(ic1);
2192  auto lnum_points = lnum_coords / 2;
2193  if (lindex < 0 || lindex > lnum_points) {
2194  lindex = lnum_points;
2195  }
2196  auto p = l + lindex * compression_unit_size(ic1);
2197  auto psize = 2 * compression_unit_size(ic1);
2198  auto min_distance = ST_Distance_Point_Polygon(p,
2199  psize,
2200  poly_coords,
2201  poly_coords_size,
2202  poly_ring_sizes,
2203  poly_num_rings,
2204  ic1,
2205  isr1,
2206  ic2,
2207  isr2,
2208  osr,
2209  threshold);
2210  if (lindex != 0) {
2211  // Statically indexed linestring: return distance from the indexed point to poly
2212  return min_distance;
2213  }
2214  if (tol_zero(min_distance)) {
2215  // Linestring's first point is inside the poly
2216  return 0.0;
2217  }
2218  if (min_distance <= threshold) {
2219  return min_distance;
2220  }
2221 
2222  // Otherwise, linestring's first point is outside the external ring or inside
2223  // an internal ring. Measure minimum distance between linestring segments and
2224  // poly rings. Crossing a ring zeroes the distance and causes an early return.
2225  auto poly_ring_coords = poly_coords;
2226  for (auto r = 0; r < poly_num_rings; r++) {
2227  int64_t poly_ring_num_coords = poly_ring_sizes[r] * 2;
2228 
2229  auto distance = distance_ring_linestring(poly_ring_coords,
2230  poly_ring_num_coords,
2231  l,
2232  lnum_coords,
2233  ic2,
2234  isr2,
2235  ic1,
2236  isr1,
2237  osr,
2238  threshold);
2239  if (min_distance > distance) {
2240  min_distance = distance;
2241  if (tol_zero(min_distance)) {
2242  return 0.0;
2243  }
2244  if (min_distance <= threshold) {
2245  return min_distance;
2246  }
2247  }
2248 
2249  poly_ring_coords += poly_ring_num_coords * compression_unit_size(ic2);
2250  }
2251 
2252  return min_distance;
2253 }
2254 
2257  int64_t lsize,
2258  int32_t lindex,
2259  int8_t* mpoly_coords,
2260  int64_t mpoly_coords_size,
2261  int32_t* mpoly_ring_sizes,
2262  int64_t mpoly_num_rings,
2263  int32_t* mpoly_poly_sizes,
2264  int64_t mpoly_num_polys,
2265  int32_t ic1,
2266  int32_t isr1,
2267  int32_t ic2,
2268  int32_t isr2,
2269  int32_t osr,
2270  double threshold) {
2271  // TODO: revisit implementation, cover all cases
2272 
2273  auto lnum_coords = lsize / compression_unit_size(ic1);
2274  auto lnum_points = lnum_coords / 2;
2275  if (lindex != 0) {
2276  // Statically indexed linestring
2277  if (lindex < 0 || lindex > lnum_points) {
2278  lindex = lnum_points;
2279  }
2280  auto p = l + lindex * compression_unit_size(ic1);
2281  auto psize = 2 * compression_unit_size(ic1);
2283  psize,
2284  mpoly_coords,
2285  mpoly_coords_size,
2286  mpoly_ring_sizes,
2287  mpoly_num_rings,
2288  mpoly_poly_sizes,
2289  mpoly_num_polys,
2290  ic1,
2291  isr1,
2292  ic2,
2293  isr2,
2294  osr,
2295  threshold);
2296  }
2297 
2298  double min_distance = 0.0;
2299 
2300  // Set specific poly pointers as we move through the coords/ringsizes/polyrings arrays.
2301  auto next_poly_coords = mpoly_coords;
2302  auto next_poly_ring_sizes = mpoly_ring_sizes;
2303 
2304  for (auto poly = 0; poly < mpoly_num_polys; poly++) {
2305  auto poly_coords = next_poly_coords;
2306  auto poly_ring_sizes = next_poly_ring_sizes;
2307  auto poly_num_rings = mpoly_poly_sizes[poly];
2308  // Count number of coords in all of poly's rings, advance ring size pointer.
2309  int32_t poly_num_coords = 0;
2310  for (auto ring = 0; ring < poly_num_rings; ring++) {
2311  poly_num_coords += 2 * *next_poly_ring_sizes++;
2312  }
2313  auto poly_coords_size = poly_num_coords * compression_unit_size(ic2);
2314  next_poly_coords += poly_coords_size;
2315  double distance = ST_Distance_LineString_Polygon(l,
2316  lsize,
2317  lindex,
2318  poly_coords,
2319  poly_coords_size,
2320  poly_ring_sizes,
2321  poly_num_rings,
2322  ic1,
2323  isr1,
2324  ic2,
2325  isr2,
2326  osr,
2327  threshold);
2328  if (poly == 0 || min_distance > distance) {
2329  min_distance = distance;
2330  if (tol_zero(min_distance)) {
2331  min_distance = 0.0;
2332  break;
2333  }
2334  if (min_distance <= threshold) {
2335  return min_distance;
2336  }
2337  }
2338  }
2339 
2340  return min_distance;
2341 }
2342 
2344 double ST_Distance_Polygon_Point(int8_t* poly_coords,
2345  int64_t poly_coords_size,
2346  int32_t* poly_ring_sizes,
2347  int64_t poly_num_rings,
2348  int8_t* p,
2349  int64_t psize,
2350  int32_t ic1,
2351  int32_t isr1,
2352  int32_t ic2,
2353  int32_t isr2,
2354  int32_t osr,
2355  double threshold) {
2356  return ST_Distance_Point_Polygon(p,
2357  psize,
2358  poly_coords,
2359  poly_coords_size,
2360  poly_ring_sizes,
2361  poly_num_rings,
2362  ic2,
2363  isr2,
2364  ic1,
2365  isr1,
2366  osr,
2367  threshold);
2368 }
2369 
2371 double ST_Distance_Polygon_LineString(int8_t* poly_coords,
2372  int64_t poly_coords_size,
2373  int32_t* poly_ring_sizes,
2374  int64_t poly_num_rings,
2375  int8_t* l,
2376  int64_t lsize,
2377  int32_t li,
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  li,
2387  poly_coords,
2388  poly_coords_size,
2389  poly_ring_sizes,
2390  poly_num_rings,
2391  ic2,
2392  isr2,
2393  ic1,
2394  isr2,
2395  osr,
2396  threshold);
2397 }
2398 
2400 double ST_Distance_Polygon_Polygon(int8_t* poly1_coords,
2401  int64_t poly1_coords_size,
2402  int32_t* poly1_ring_sizes,
2403  int64_t poly1_num_rings,
2404  int8_t* poly2_coords,
2405  int64_t poly2_coords_size,
2406  int32_t* poly2_ring_sizes,
2407  int64_t poly2_num_rings,
2408  int32_t ic1,
2409  int32_t isr1,
2410  int32_t ic2,
2411  int32_t isr2,
2412  int32_t osr,
2413  double threshold) {
2414  // Check if poly1 contains the first point of poly2's shape, i.e. the external ring
2415  auto poly2_first_point_coords = poly2_coords;
2416  auto poly2_first_point_coords_size = compression_unit_size(ic2) * 2;
2417  auto min_distance = ST_Distance_Polygon_Point(poly1_coords,
2418  poly1_coords_size,
2419  poly1_ring_sizes,
2420  poly1_num_rings,
2421  poly2_first_point_coords,
2422  poly2_first_point_coords_size,
2423  ic1,
2424  isr1,
2425  ic2,
2426  isr2,
2427  osr,
2428  threshold);
2429  if (tol_zero(min_distance)) {
2430  // Polygons overlap
2431  return 0.0;
2432  }
2433  if (min_distance <= threshold) {
2434  return min_distance;
2435  }
2436 
2437  // Poly2's first point is either outside poly1's external ring or inside one of the
2438  // internal rings. Measure the smallest distance between a poly1 ring (external or
2439  // internal) and a poly2 ring (external or internal). If poly2 is completely outside
2440  // poly1, then the min distance would be between poly1's and poly2's external rings. If
2441  // poly2 is completely inside one of poly1 internal rings then the min distance would be
2442  // between that poly1 internal ring and poly2's external ring. If poly1 is completely
2443  // inside one of poly2 internal rings, min distance is between that internal ring and
2444  // poly1's external ring. In each case other rings don't get in the way. Any ring
2445  // intersection means zero distance - short-circuit and return.
2446 
2447  auto poly1_ring_coords = poly1_coords;
2448  for (auto r1 = 0; r1 < poly1_num_rings; r1++) {
2449  int64_t poly1_ring_num_coords = poly1_ring_sizes[r1] * 2;
2450 
2451  auto poly2_ring_coords = poly2_coords;
2452  for (auto r2 = 0; r2 < poly2_num_rings; r2++) {
2453  int64_t poly2_ring_num_coords = poly2_ring_sizes[r2] * 2;
2454 
2455  auto distance = distance_ring_ring(poly1_ring_coords,
2456  poly1_ring_num_coords,
2457  poly2_ring_coords,
2458  poly2_ring_num_coords,
2459  ic1,
2460  isr1,
2461  ic2,
2462  isr2,
2463  osr,
2464  threshold);
2465  if (min_distance > distance) {
2466  min_distance = distance;
2467  if (tol_zero(min_distance)) {
2468  return 0.0;
2469  }
2470  if (min_distance <= threshold) {
2471  return min_distance;
2472  }
2473  }
2474 
2475  poly2_ring_coords += poly2_ring_num_coords * compression_unit_size(ic2);
2476  }
2477 
2478  poly1_ring_coords += poly1_ring_num_coords * compression_unit_size(ic1);
2479  }
2480 
2481  return min_distance;
2482 }
2483 
2485 double ST_Distance_Polygon_MultiPolygon(int8_t* poly1_coords,
2486  int64_t poly1_coords_size,
2487  int32_t* poly1_ring_sizes,
2488  int64_t poly1_num_rings,
2489  int8_t* mpoly_coords,
2490  int64_t mpoly_coords_size,
2491  int32_t* mpoly_ring_sizes,
2492  int64_t mpoly_num_rings,
2493  int32_t* mpoly_poly_sizes,
2494  int64_t mpoly_num_polys,
2495  int32_t ic1,
2496  int32_t isr1,
2497  int32_t ic2,
2498  int32_t isr2,
2499  int32_t osr,
2500  double threshold) {
2501  double min_distance = 0.0;
2502 
2503  // Set specific poly pointers as we move through the coords/ringsizes/polyrings arrays.
2504  auto next_poly_coords = mpoly_coords;
2505  auto next_poly_ring_sizes = mpoly_ring_sizes;
2506 
2507  for (auto poly = 0; poly < mpoly_num_polys; poly++) {
2508  auto poly_coords = next_poly_coords;
2509  auto poly_ring_sizes = next_poly_ring_sizes;
2510  auto poly_num_rings = mpoly_poly_sizes[poly];
2511  // Count number of coords in all of poly's rings, advance ring size pointer.
2512  int32_t poly_num_coords = 0;
2513  for (auto ring = 0; ring < poly_num_rings; ring++) {
2514  poly_num_coords += 2 * *next_poly_ring_sizes++;
2515  }
2516  auto poly_coords_size = poly_num_coords * compression_unit_size(ic2);
2517  next_poly_coords += poly_coords_size;
2518  double distance = ST_Distance_Polygon_Polygon(poly1_coords,
2519  poly1_coords_size,
2520  poly1_ring_sizes,
2521  poly1_num_rings,
2522  poly_coords,
2523  poly_coords_size,
2524  poly_ring_sizes,
2525  poly_num_rings,
2526  ic1,
2527  isr1,
2528  ic2,
2529  isr2,
2530  osr,
2531  threshold);
2532  if (poly == 0 || min_distance > distance) {
2533  min_distance = distance;
2534  if (tol_zero(min_distance)) {
2535  min_distance = 0.0;
2536  break;
2537  }
2538  if (min_distance <= threshold) {
2539  break;
2540  }
2541  }
2542  }
2543 
2544  return min_distance;
2545 }
2546 
2548 double ST_Distance_MultiPolygon_Point(int8_t* mpoly_coords,
2549  int64_t mpoly_coords_size,
2550  int32_t* mpoly_ring_sizes,
2551  int64_t mpoly_num_rings,
2552  int32_t* mpoly_poly_sizes,
2553  int64_t mpoly_num_polys,
2554  int8_t* p,
2555  int64_t psize,
2556  int32_t ic1,
2557  int32_t isr1,
2558  int32_t ic2,
2559  int32_t isr2,
2560  int32_t osr,
2561  double threshold) {
2563  psize,
2564  mpoly_coords,
2565  mpoly_coords_size,
2566  mpoly_ring_sizes,
2567  mpoly_num_rings,
2568  mpoly_poly_sizes,
2569  mpoly_num_polys,
2570  ic2,
2571  isr2,
2572  ic1,
2573  isr1,
2574  osr,
2575  threshold);
2576 }
2577 
2579 double ST_Distance_MultiPolygon_LineString(int8_t* mpoly_coords,
2580  int64_t mpoly_coords_size,
2581  int32_t* mpoly_ring_sizes,
2582  int64_t mpoly_num_rings,
2583  int32_t* mpoly_poly_sizes,
2584  int64_t mpoly_num_polys,
2585  int8_t* l,
2586  int64_t lsize,
2587  int32_t lindex,
2588  int32_t ic1,
2589  int32_t isr1,
2590  int32_t ic2,
2591  int32_t isr2,
2592  int32_t osr,
2593  double threshold) {
2595  lsize,
2596  lindex,
2597  mpoly_coords,
2598  mpoly_coords_size,
2599  mpoly_ring_sizes,
2600  mpoly_num_rings,
2601  mpoly_poly_sizes,
2602  mpoly_num_polys,
2603  ic2,
2604  isr2,
2605  ic1,
2606  isr1,
2607  osr,
2608  threshold);
2609 }
2610 
2612 double ST_Distance_MultiPolygon_Polygon(int8_t* mpoly_coords,
2613  int64_t mpoly_coords_size,
2614  int32_t* mpoly_ring_sizes,
2615  int64_t mpoly_num_rings,
2616  int32_t* mpoly_poly_sizes,
2617  int64_t mpoly_num_polys,
2618  int8_t* poly1_coords,
2619  int64_t poly1_coords_size,
2620  int32_t* poly1_ring_sizes,
2621  int64_t poly1_num_rings,
2622  int32_t ic1,
2623  int32_t isr1,
2624  int32_t ic2,
2625  int32_t isr2,
2626  int32_t osr,
2627  double threshold) {
2628  return ST_Distance_Polygon_MultiPolygon(poly1_coords,
2629  poly1_coords_size,
2630  poly1_ring_sizes,
2631  poly1_num_rings,
2632  mpoly_coords,
2633  mpoly_coords_size,
2634  mpoly_ring_sizes,
2635  mpoly_num_rings,
2636  mpoly_poly_sizes,
2637  mpoly_num_polys,
2638  ic2,
2639  isr2,
2640  ic1,
2641  isr1,
2642  osr,
2643  threshold);
2644 }
2645 
2647 double ST_Distance_MultiPolygon_MultiPolygon(int8_t* mpoly1_coords,
2648  int64_t mpoly1_coords_size,
2649  int32_t* mpoly1_ring_sizes,
2650  int64_t mpoly1_num_rings,
2651  int32_t* mpoly1_poly_sizes,
2652  int64_t mpoly1_num_polys,
2653  int8_t* mpoly2_coords,
2654  int64_t mpoly2_coords_size,
2655  int32_t* mpoly2_ring_sizes,
2656  int64_t mpoly2_num_rings,
2657  int32_t* mpoly2_poly_sizes,
2658  int64_t mpoly2_num_polys,
2659  int32_t ic1,
2660  int32_t isr1,
2661  int32_t ic2,
2662  int32_t isr2,
2663  int32_t osr,
2664  double threshold) {
2665  double min_distance = 0.0;
2666 
2667  // Set specific poly pointers as we move through mpoly1's coords/ringsizes/polyrings
2668  // arrays.
2669  auto next_poly_coords = mpoly1_coords;
2670  auto next_poly_ring_sizes = mpoly1_ring_sizes;
2671 
2672  for (auto poly = 0; poly < mpoly1_num_polys; poly++) {
2673  auto poly_coords = next_poly_coords;
2674  auto poly_ring_sizes = next_poly_ring_sizes;
2675  auto poly_num_rings = mpoly1_poly_sizes[poly];
2676  // Count number of coords in all of poly's rings, advance ring size pointer.
2677  int32_t poly_num_coords = 0;
2678  for (auto ring = 0; ring < poly_num_rings; ring++) {
2679  poly_num_coords += 2 * *next_poly_ring_sizes++;
2680  }
2681  auto poly_coords_size = poly_num_coords * compression_unit_size(ic1);
2682  next_poly_coords += poly_coords_size;
2683  double distance = ST_Distance_Polygon_MultiPolygon(poly_coords,
2684  poly_coords_size,
2685  poly_ring_sizes,
2686  poly_num_rings,
2687  mpoly2_coords,
2688  mpoly2_coords_size,
2689  mpoly2_ring_sizes,
2690  mpoly2_num_rings,
2691  mpoly2_poly_sizes,
2692  mpoly2_num_polys,
2693  ic1,
2694  isr1,
2695  ic2,
2696  isr2,
2697  osr,
2698  threshold);
2699  if (poly == 0 || min_distance > distance) {
2700  min_distance = distance;
2701  if (tol_zero(min_distance)) {
2702  min_distance = 0.0;
2703  break;
2704  }
2705  if (min_distance <= threshold) {
2706  break;
2707  }
2708  }
2709  }
2710 
2711  return min_distance;
2712 }
2713 
2714 //
2715 // ST_DWithin
2716 //
2717 
2719 bool ST_DWithin_Point_Point(int8_t* p1,
2720  int64_t p1size,
2721  int8_t* p2,
2722  int64_t p2size,
2723  int32_t ic1,
2724  int32_t isr1,
2725  int32_t ic2,
2726  int32_t isr2,
2727  int32_t osr,
2728  double distance_within) {
2730  p1, p1size, p2, p2size, ic1, isr1, ic2, isr2, osr) <=
2731  distance_within * distance_within;
2732 }
2733 
2736  int64_t p1size,
2737  int8_t* l2,
2738  int64_t l2size,
2739  double* l2bounds,
2740  int64_t l2bounds_size,
2741  int32_t l2index,
2742  int32_t ic1,
2743  int32_t isr1,
2744  int32_t ic2,
2745  int32_t isr2,
2746  int32_t osr,
2747  double distance_within) {
2748  if (l2bounds) {
2749  // Bounding boxes need to be transformed to output SR before proximity check
2750  if (!point_dwithin_box(
2751  p1, p1size, ic1, isr1, l2bounds, l2bounds_size, isr2, osr, distance_within)) {
2752  return false;
2753  }
2754  }
2755 
2756  // May need to adjust the threshold by TOLERANCE_DEFAULT
2757  const double threshold = distance_within;
2759  p1, p1size, l2, l2size, l2index, ic1, isr1, ic2, isr2, osr, threshold) <=
2760  distance_within;
2761 }
2762 
2765  int64_t psize,
2766  int8_t* poly_coords,
2767  int64_t poly_coords_size,
2768  int32_t* poly_ring_sizes,
2769  int64_t poly_num_rings,
2770  double* poly_bounds,
2771  int64_t poly_bounds_size,
2772  int32_t ic1,
2773  int32_t isr1,
2774  int32_t ic2,
2775  int32_t isr2,
2776  int32_t osr,
2777  double distance_within) {
2778  if (poly_bounds) {
2779  if (!point_dwithin_box(p,
2780  psize,
2781  ic1,
2782  isr1,
2783  poly_bounds,
2784  poly_bounds_size,
2785  isr2,
2786  osr,
2787  distance_within)) {
2788  return false;
2789  }
2790  }
2791 
2792  // May need to adjust the threshold by TOLERANCE_DEFAULT
2793  const double threshold = distance_within;
2794  return ST_Distance_Point_Polygon(p,
2795  psize,
2796  poly_coords,
2797  poly_coords_size,
2798  poly_ring_sizes,
2799  poly_num_rings,
2800  ic1,
2801  isr1,
2802  ic2,
2803  isr2,
2804  osr,
2805  threshold) <= distance_within;
2806 }
2807 
2810  int64_t psize,
2811  int8_t* mpoly_coords,
2812  int64_t mpoly_coords_size,
2813  int32_t* mpoly_ring_sizes,
2814  int64_t mpoly_num_rings,
2815  int32_t* mpoly_poly_sizes,
2816  int64_t mpoly_num_polys,
2817  double* mpoly_bounds,
2818  int64_t mpoly_bounds_size,
2819  int32_t ic1,
2820  int32_t isr1,
2821  int32_t ic2,
2822  int32_t isr2,
2823  int32_t osr,
2824  double distance_within) {
2825  if (mpoly_bounds) {
2826  if (!point_dwithin_box(p,
2827  psize,
2828  ic1,
2829  isr1,
2830  mpoly_bounds,
2831  mpoly_bounds_size,
2832  isr2,
2833  osr,
2834  distance_within)) {
2835  return false;
2836  }
2837  }
2838 
2839  // May need to adjust the threshold by TOLERANCE_DEFAULT
2840  const double threshold = distance_within;
2842  psize,
2843  mpoly_coords,
2844  mpoly_coords_size,
2845  mpoly_ring_sizes,
2846  mpoly_num_rings,
2847  mpoly_poly_sizes,
2848  mpoly_num_polys,
2849  ic1,
2850  isr1,
2851  ic2,
2852  isr2,
2853  osr,
2854  threshold) <= distance_within;
2855 }
2856 
2859  int64_t l1size,
2860  double* l1bounds,
2861  int64_t l1bounds_size,
2862  int32_t l1index,
2863  int8_t* l2,
2864  int64_t l2size,
2865  double* l2bounds,
2866  int64_t l2bounds_size,
2867  int32_t l2index,
2868  int32_t ic1,
2869  int32_t isr1,
2870  int32_t ic2,
2871  int32_t isr2,
2872  int32_t osr,
2873  double distance_within) {
2874  if (l1bounds && l2bounds) {
2875  // Bounding boxes need to be transformed to output SR before proximity check
2876  if (!box_dwithin_box(l1bounds,
2877  l1bounds_size,
2878  isr1,
2879  l2bounds,
2880  l2bounds_size,
2881  isr2,
2882  osr,
2883  distance_within)) {
2884  return false;
2885  }
2886  }
2887 
2888  // May need to adjust the threshold by TOLERANCE_DEFAULT
2889  const double threshold = distance_within;
2891  l1size,
2892  l1index,
2893  l2,
2894  l2size,
2895  l2index,
2896  ic1,
2897  isr1,
2898  ic2,
2899  isr2,
2900  osr,
2901  threshold) <= distance_within;
2902 }
2903 
2906  int64_t l1size,
2907  double* l1bounds,
2908  int64_t l1bounds_size,
2909  int32_t l1index,
2910  int8_t* poly_coords,
2911  int64_t poly_coords_size,
2912  int32_t* poly_ring_sizes,
2913  int64_t poly_num_rings,
2914  double* poly_bounds,
2915  int64_t poly_bounds_size,
2916  int32_t ic1,
2917  int32_t isr1,
2918  int32_t ic2,
2919  int32_t isr2,
2920  int32_t osr,
2921  double distance_within) {
2922  if (l1bounds && poly_bounds) {
2923  // Bounding boxes need to be transformed to output SR before proximity check
2924  if (!box_dwithin_box(l1bounds,
2925  l1bounds_size,
2926  isr1,
2927  poly_bounds,
2928  poly_bounds_size,
2929  isr2,
2930  osr,
2931  distance_within)) {
2932  return false;
2933  }
2934  }
2935 
2936  // May need to adjust the threshold by TOLERANCE_DEFAULT
2937  const double threshold = distance_within;
2939  l1size,
2940  l1index,
2941  poly_coords,
2942  poly_coords_size,
2943  poly_ring_sizes,
2944  poly_num_rings,
2945  ic1,
2946  isr1,
2947  ic2,
2948  isr2,
2949  osr,
2950  threshold) <= distance_within;
2951 }
2952 
2955  int64_t l1size,
2956  double* l1bounds,
2957  int64_t l1bounds_size,
2958  int32_t l1index,
2959  int8_t* mpoly_coords,
2960  int64_t mpoly_coords_size,
2961  int32_t* mpoly_ring_sizes,
2962  int64_t mpoly_num_rings,
2963  int32_t* mpoly_poly_sizes,
2964  int64_t mpoly_num_polys,
2965  double* mpoly_bounds,
2966  int64_t mpoly_bounds_size,
2967  int32_t ic1,
2968  int32_t isr1,
2969  int32_t ic2,
2970  int32_t isr2,
2971  int32_t osr,
2972  double distance_within) {
2973  if (l1bounds && mpoly_bounds) {
2974  // Bounding boxes need to be transformed to output SR before proximity check
2975  if (!box_dwithin_box(l1bounds,
2976  l1bounds_size,
2977  isr1,
2978  mpoly_bounds,
2979  mpoly_bounds_size,
2980  isr2,
2981  osr,
2982  distance_within)) {
2983  return false;
2984  }
2985  }
2986 
2987  // May need to adjust the threshold by TOLERANCE_DEFAULT
2988  const double threshold = distance_within;
2990  l1size,
2991  l1index,
2992  mpoly_coords,
2993  mpoly_coords_size,
2994  mpoly_ring_sizes,
2995  mpoly_num_rings,
2996  mpoly_poly_sizes,
2997  mpoly_num_polys,
2998  ic1,
2999  isr1,
3000  ic2,
3001  isr2,
3002  osr,
3003  threshold) <= distance_within;
3004 }
3005 
3007 bool ST_DWithin_Polygon_Polygon(int8_t* poly1_coords,
3008  int64_t poly1_coords_size,
3009  int32_t* poly1_ring_sizes,
3010  int64_t poly1_num_rings,
3011  double* poly1_bounds,
3012  int64_t poly1_bounds_size,
3013  int8_t* poly2_coords,
3014  int64_t poly2_coords_size,
3015  int32_t* poly2_ring_sizes,
3016  int64_t poly2_num_rings,
3017  double* poly2_bounds,
3018  int64_t poly2_bounds_size,
3019  int32_t ic1,
3020  int32_t isr1,
3021  int32_t ic2,
3022  int32_t isr2,
3023  int32_t osr,
3024  double distance_within) {
3025  if (poly1_bounds && poly2_bounds) {
3026  // Bounding boxes need to be transformed to output SR before proximity check
3027  if (!box_dwithin_box(poly1_bounds,
3028  poly1_bounds_size,
3029  isr1,
3030  poly2_bounds,
3031  poly2_bounds_size,
3032  isr2,
3033  osr,
3034  distance_within)) {
3035  return false;
3036  }
3037  }
3038 
3039  // May need to adjust the threshold by TOLERANCE_DEFAULT
3040  const double threshold = distance_within;
3041  return ST_Distance_Polygon_Polygon(poly1_coords,
3042  poly1_coords_size,
3043  poly1_ring_sizes,
3044  poly1_num_rings,
3045  poly2_coords,
3046  poly2_coords_size,
3047  poly2_ring_sizes,
3048  poly2_num_rings,
3049  ic1,
3050  isr1,
3051  ic2,
3052  isr2,
3053  osr,
3054  threshold) <= distance_within;
3055 }
3056 
3058 bool ST_DWithin_Polygon_MultiPolygon(int8_t* poly_coords,
3059  int64_t poly_coords_size,
3060  int32_t* poly_ring_sizes,
3061  int64_t poly_num_rings,
3062  double* poly_bounds,
3063  int64_t poly_bounds_size,
3064  int8_t* mpoly_coords,
3065  int64_t mpoly_coords_size,
3066  int32_t* mpoly_ring_sizes,
3067  int64_t mpoly_num_rings,
3068  int32_t* mpoly_poly_sizes,
3069  int64_t mpoly_num_polys,
3070  double* mpoly_bounds,
3071  int64_t mpoly_bounds_size,
3072  int32_t ic1,
3073  int32_t isr1,
3074  int32_t ic2,
3075  int32_t isr2,
3076  int32_t osr,
3077  double distance_within) {
3078  if (poly_bounds && mpoly_bounds) {
3079  // Bounding boxes need to be transformed to output SR before proximity check
3080  if (!box_dwithin_box(poly_bounds,
3081  poly_bounds_size,
3082  isr1,
3083  mpoly_bounds,
3084  mpoly_bounds_size,
3085  isr2,
3086  osr,
3087  distance_within)) {
3088  return false;
3089  }
3090  }
3091 
3092  // May need to adjust the threshold by TOLERANCE_DEFAULT
3093  const double threshold = distance_within;
3094  return ST_Distance_Polygon_MultiPolygon(poly_coords,
3095  poly_coords_size,
3096  poly_ring_sizes,
3097  poly_num_rings,
3098  mpoly_coords,
3099  mpoly_coords_size,
3100  mpoly_ring_sizes,
3101  mpoly_num_rings,
3102  mpoly_poly_sizes,
3103  mpoly_num_polys,
3104  ic1,
3105  isr1,
3106  ic2,
3107  isr2,
3108  osr,
3109  threshold) <= distance_within;
3110 }
3111 
3113 bool ST_DWithin_MultiPolygon_MultiPolygon(int8_t* mpoly1_coords,
3114  int64_t mpoly1_coords_size,
3115  int32_t* mpoly1_ring_sizes,
3116  int64_t mpoly1_num_rings,
3117  int32_t* mpoly1_poly_sizes,
3118  int64_t mpoly1_num_polys,
3119  double* mpoly1_bounds,
3120  int64_t mpoly1_bounds_size,
3121  int8_t* mpoly2_coords,
3122  int64_t mpoly2_coords_size,
3123  int32_t* mpoly2_ring_sizes,
3124  int64_t mpoly2_num_rings,
3125  int32_t* mpoly2_poly_sizes,
3126  int64_t mpoly2_num_polys,
3127  double* mpoly2_bounds,
3128  int64_t mpoly2_bounds_size,
3129  int32_t ic1,
3130  int32_t isr1,
3131  int32_t ic2,
3132  int32_t isr2,
3133  int32_t osr,
3134  double distance_within) {
3135  if (mpoly1_bounds && mpoly2_bounds) {
3136  // Bounding boxes need to be transformed to output SR before proximity check
3137  if (!box_dwithin_box(mpoly1_bounds,
3138  mpoly1_bounds_size,
3139  isr1,
3140  mpoly2_bounds,
3141  mpoly2_bounds_size,
3142  isr2,
3143  osr,
3144  distance_within)) {
3145  return false;
3146  }
3147  }
3148 
3149  // May need to adjust the threshold by TOLERANCE_DEFAULT
3150  const double threshold = distance_within;
3151  return ST_Distance_MultiPolygon_MultiPolygon(mpoly1_coords,
3152  mpoly1_coords_size,
3153  mpoly1_ring_sizes,
3154  mpoly1_num_rings,
3155  mpoly1_poly_sizes,
3156  mpoly1_num_polys,
3157  mpoly2_coords,
3158  mpoly2_coords_size,
3159  mpoly2_ring_sizes,
3160  mpoly2_num_rings,
3161  mpoly2_poly_sizes,
3162  mpoly2_num_polys,
3163  ic1,
3164  isr1,
3165  ic2,
3166  isr2,
3167  osr,
3168  threshold) <= distance_within;
3169 }
3170 
3171 //
3172 // ST_MaxDistance
3173 //
3174 
3175 // Max cartesian distance between a point and a line segment
3176 DEVICE
3177 double max_distance_point_line(double px,
3178  double py,
3179  double l1x,
3180  double l1y,
3181  double l2x,
3182  double l2y) {
3183  // TODO: switch to squared distances
3184  double length1 = distance_point_point(px, py, l1x, l1y);
3185  double length2 = distance_point_point(px, py, l2x, l2y);
3186  if (length1 > length2) {
3187  return length1;
3188  }
3189  return length2;
3190 }
3191 
3193  int64_t psize,
3194  int8_t* l,
3195  int64_t lsize,
3196  int32_t lindex,
3197  int32_t ic1,
3198  int32_t isr1,
3199  int32_t ic2,
3200  int32_t isr2,
3201  int32_t osr,
3202  bool check_closed) {
3203  // TODO: switch to squared distances
3204  double px = coord_x(p, 0, ic1, isr1, osr);
3205  double py = coord_y(p, 1, ic1, isr1, osr);
3206 
3207  auto l_num_coords = lsize / compression_unit_size(ic2);
3208  auto l_num_points = l_num_coords / 2;
3209  if (lindex != 0) { // Statically indexed linestring
3210  if (lindex < 0 || lindex > l_num_points) {
3211  lindex = l_num_points; // Endpoint
3212  }
3213  double lx = coord_x(l, 2 * (lindex - 1), ic2, isr2, osr);
3214  double ly = coord_y(l, 2 * (lindex - 1) + 1, ic2, isr2, osr);
3215  return distance_point_point(px, py, lx, ly);
3216  }
3217 
3218  double l1x = coord_x(l, 0, ic2, isr2, osr);
3219  double l1y = coord_y(l, 1, ic2, isr2, osr);
3220  double l2x = coord_x(l, 2, ic2, isr2, osr);
3221  double l2y = coord_y(l, 3, ic2, isr2, osr);
3222 
3223  double max_dist = max_distance_point_line(px, py, l1x, l1y, l2x, l2y);
3224  for (int32_t i = 4; i < l_num_coords; i += 2) {
3225  l1x = l2x; // advance one point
3226  l1y = l2y;
3227  l2x = coord_x(l, i, ic2, isr2, osr);
3228  l2y = coord_y(l, i + 1, ic2, isr2, osr);
3229  double ldist = max_distance_point_line(px, py, l1x, l1y, l2x, l2y);
3230  if (max_dist < ldist) {
3231  max_dist = ldist;
3232  }
3233  }
3234  if (l_num_coords > 4 && check_closed) {
3235  // Also check distance to the closing edge between the first and the last points
3236  l1x = coord_x(l, 0, ic2, isr2, osr);
3237  l1y = coord_y(l, 1, ic2, isr2, osr);
3238  double ldist = max_distance_point_line(px, py, l1x, l1y, l2x, l2y);
3239  if (max_dist < ldist) {
3240  max_dist = ldist;
3241  }
3242  }
3243  return max_dist;
3244 }
3245 
3248  int64_t psize,
3249  int8_t* l,
3250  int64_t lsize,
3251  int32_t lindex,
3252  int32_t ic1,
3253  int32_t isr1,
3254  int32_t ic2,
3255  int32_t isr2,
3256  int32_t osr) {
3258  p, psize, l, lsize, lindex, ic1, isr1, ic2, isr2, osr, false);
3259 }
3260 
3263  int64_t lsize,
3264  int32_t lindex,
3265  int8_t* p,
3266  int64_t psize,
3267  int32_t ic1,
3268  int32_t isr1,
3269  int32_t ic2,
3270  int32_t isr2,
3271  int32_t osr) {
3273  p, psize, l, lsize, lindex, ic2, isr2, ic1, isr1, osr, false);
3274 }
3275 
3276 // TODO: add ST_MaxDistance_LineString_LineString (with short-circuit threshold)
3277 
3278 //
3279 // ST_Contains
3280 //
3281 
3283 bool ST_Contains_Point_Point(int8_t* p1,
3284  int64_t p1size,
3285  int8_t* p2,
3286  int64_t p2size,
3287  int32_t ic1,
3288  int32_t isr1,
3289  int32_t ic2,
3290  int32_t isr2,
3291  int32_t osr) {
3292  double p1x = coord_x(p1, 0, ic1, isr1, osr);
3293  double p1y = coord_y(p1, 1, ic1, isr1, osr);
3294  double p2x = coord_x(p2, 0, ic2, isr2, osr);
3295  double p2y = coord_y(p2, 1, ic2, isr2, osr);
3296  double tolerance = tol(ic1, ic2);
3297  return tol_eq(p1x, p2x, tolerance) && tol_eq(p1y, p2y, tolerance);
3298 }
3299 
3302  int64_t psize,
3303  int8_t* l,
3304  int64_t lsize,
3305  double* lbounds,
3306  int64_t lbounds_size,
3307  int32_t li,
3308  int32_t ic1,
3309  int32_t isr1,
3310  int32_t ic2,
3311  int32_t isr2,
3312  int32_t osr) {
3313  double px = coord_x(p, 0, ic1, isr1, osr);
3314  double py = coord_y(p, 1, ic1, isr1, osr);
3315 
3316  if (lbounds) {
3317  if (tol_eq(px, lbounds[0]) && tol_eq(py, lbounds[1]) && tol_eq(px, lbounds[2]) &&
3318  tol_eq(py, lbounds[3])) {
3319  return true;
3320  }
3321  }
3322 
3323  auto l_num_coords = lsize / compression_unit_size(ic2);
3324  for (int i = 0; i < l_num_coords; i += 2) {
3325  double lx = coord_x(l, i, ic2, isr2, osr);
3326  double ly = coord_y(l, i + 1, ic2, isr2, osr);
3327  if (tol_eq(px, lx) && tol_eq(py, ly)) {
3328  continue;
3329  }
3330  return false;
3331  }
3332  return true;
3333 }
3334 
3337  int64_t psize,
3338  int8_t* poly_coords,
3339  int64_t poly_coords_size,
3340  int32_t* poly_ring_sizes,
3341  int64_t poly_num_rings,
3342  double* poly_bounds,
3343  int64_t poly_bounds_size,
3344  int32_t ic1,
3345  int32_t isr1,
3346  int32_t ic2,
3347  int32_t isr2,
3348  int32_t osr) {
3349  auto exterior_ring_num_coords = poly_coords_size / compression_unit_size(ic2);
3350  if (poly_num_rings > 0) {
3351  exterior_ring_num_coords = poly_ring_sizes[0] * 2;
3352  }
3353  auto exterior_ring_coords_size = exterior_ring_num_coords * compression_unit_size(ic2);
3354 
3356  psize,
3357  poly_coords,
3358  exterior_ring_coords_size,
3359  poly_bounds,
3360  poly_bounds_size,
3361  0,
3362  ic1,
3363  isr1,
3364  ic2,
3365  isr2,
3366  osr);
3367 }
3368 
3371  int64_t lsize,
3372  double* lbounds,
3373  int64_t lbounds_size,
3374  int32_t li,
3375  int8_t* p,
3376  int64_t psize,
3377  int32_t ic1,
3378  int32_t isr1,
3379  int32_t ic2,
3380  int32_t isr2,
3381  int32_t osr) {
3383  p, psize, l, lsize, li, ic2, isr2, ic1, isr1, osr, 0.0));
3384 }
3385 
3388  int64_t l1size,
3389  double* l1bounds,
3390  int64_t l1bounds_size,
3391  int32_t l1i,
3392  int8_t* l2,
3393  int64_t l2size,
3394  double* l2bounds,
3395  int64_t l2bounds_size,
3396  int32_t l2i,
3397  int32_t ic1,
3398  int32_t isr1,
3399  int32_t ic2,
3400  int32_t isr2,
3401  int32_t osr) {
3402  if (l1i != 0 || l2i != 0) {
3403  // At least one linestring is indexed, can rely on distance
3405  l1, l1size, l1i, l2, l2size, l2i, ic1, isr1, ic2, isr2, osr, 0.0));
3406  }
3407 
3408  // TODO: sublinestring
3409  // For each line segment in l2 check if there is a segment in l1
3410  // that it's colinear with and both l2 vertices are on l1 segment.
3411  // Bail if any line segment deviates from the path.
3412  return false;
3413 }
3414 
3417  int64_t lsize,
3418  double* lbounds,
3419  int64_t lbounds_size,
3420  int32_t li,
3421  int8_t* poly_coords,
3422  int64_t poly_coords_size,
3423  int32_t* poly_ring_sizes,
3424  int64_t poly_num_rings,
3425  double* poly_bounds,
3426  int64_t poly_bounds_size,
3427  int32_t ic1,
3428  int32_t isr1,
3429  int32_t ic2,
3430  int32_t isr2,
3431  int32_t osr) {
3432  // TODO
3433  return false;
3434 }
3435 
3436 template <typename T, EdgeBehavior TEdgeBehavior>
3437 DEVICE ALWAYS_INLINE bool Contains_Polygon_Point_Impl(const int8_t* poly_coords,
3438  const int64_t poly_coords_size,
3439  const int32_t* poly_ring_sizes,
3440  const int64_t poly_num_rings,
3441  const double* poly_bounds,
3442  const int64_t poly_bounds_size,
3443  const int8_t* p,
3444  const int64_t psize,
3445  const int32_t ic1,
3446  const int32_t isr1,
3447  const int32_t ic2,
3448  const int32_t isr2,
3449  const int32_t osr) {
3450  if (poly_bounds) {
3451  // TODO: codegen
3452  if (!box_contains_point(poly_bounds,
3453  poly_bounds_size,
3454  coord_x(p, 0, ic2, isr2, osr),
3455  coord_y(p, 1, ic2, isr2, osr))) {
3456  DEBUG_STMT(printf("Bounding box does not contain point, exiting.\n"));
3457  return false;
3458  }
3459  }
3460 
3461  auto get_x_coord = [=]() -> T {
3462  if constexpr (std::is_floating_point<T>::value) {
3463  return coord_x(p, 0, ic2, isr2, osr);
3464  } else {
3465  return compressed_coord(p, 0);
3466  }
3467  return T{}; // https://stackoverflow.com/a/64561686/2700898
3468  };
3469 
3470  auto get_y_coord = [=]() -> T {
3471  if constexpr (std::is_floating_point<T>::value) {
3472  return coord_y(p, 1, ic2, isr2, osr);
3473  } else {
3474  return compressed_coord(p, 1);
3475  }
3476  return T{}; // https://stackoverflow.com/a/64561686/2700898
3477  };
3478 
3479  const T px = get_x_coord();
3480  const T py = get_y_coord();
3481  DEBUG_STMT(printf("pt: %f, %f\n", (double)px, (double)py));
3482 
3483  const auto poly_num_coords = poly_coords_size / compression_unit_size(ic1);
3484  auto exterior_ring_num_coords =
3485  poly_num_rings > 0 ? poly_ring_sizes[0] * 2 : poly_num_coords;
3486 
3487  auto poly = poly_coords;
3488  if (point_in_polygon_winding_number<T, TEdgeBehavior>(
3489  poly, exterior_ring_num_coords, px, py, ic1, isr1, osr)) {
3490  // Inside exterior ring
3491  poly += exterior_ring_num_coords * compression_unit_size(ic1);
3492  // Check that none of the polygon's holes contain that point
3493  for (auto r = 1; r < poly_num_rings; r++) {
3494  const int64_t interior_ring_num_coords = poly_ring_sizes[r] * 2;
3495  if (point_in_polygon_winding_number<T, TEdgeBehavior>(
3496  poly, interior_ring_num_coords, px, py, ic1, isr1, osr)) {
3497  return false;
3498  }
3499  poly += interior_ring_num_coords * compression_unit_size(ic1);
3500  }
3501  return true;
3502  }
3503  return false;
3504 }
3505 
3506 EXTENSION_NOINLINE bool ST_Contains_Polygon_Point(const int8_t* poly_coords,
3507  const int64_t poly_coords_size,
3508  const int32_t* poly_ring_sizes,
3509  const int64_t poly_num_rings,
3510  const double* poly_bounds,
3511  const int64_t poly_bounds_size,
3512  const int8_t* p,
3513  const int64_t psize,
3514  const int32_t ic1,
3515  const int32_t isr1,
3516  const int32_t ic2,
3517  const int32_t isr2,
3518  const int32_t osr) {
3519  return Contains_Polygon_Point_Impl<double, EdgeBehavior::kExcludePointOnEdge>(
3520  poly_coords,
3521  poly_coords_size,
3522  poly_ring_sizes,
3523  poly_num_rings,
3524  poly_bounds,
3525  poly_bounds_size,
3526  p,
3527  psize,
3528  ic1,
3529  isr1,
3530  ic2,
3531  isr2,
3532  osr);
3533 }
3534 
3535 EXTENSION_NOINLINE bool ST_cContains_Polygon_Point(const int8_t* poly_coords,
3536  const int64_t poly_coords_size,
3537  const int32_t* poly_ring_sizes,
3538  const int64_t poly_num_rings,
3539  const double* poly_bounds,
3540  const int64_t poly_bounds_size,
3541  const int8_t* p,
3542  const int64_t psize,
3543  const int32_t ic1,
3544  const int32_t isr1,
3545  const int32_t ic2,
3546  const int32_t isr2,
3547  const int32_t osr) {
3548  return Contains_Polygon_Point_Impl<int64_t, EdgeBehavior::kExcludePointOnEdge>(
3549  poly_coords,
3550  poly_coords_size,
3551  poly_ring_sizes,
3552  poly_num_rings,
3553  poly_bounds,
3554  poly_bounds_size,
3555  p,
3556  psize,
3557  ic1,
3558  isr1,
3559  ic2,
3560  isr2,
3561  osr);
3562 }
3563 
3565 bool ST_Contains_Polygon_LineString(int8_t* poly_coords,
3566  int64_t poly_coords_size,
3567  int32_t* poly_ring_sizes,
3568  int64_t poly_num_rings,
3569  double* poly_bounds,
3570  int64_t poly_bounds_size,
3571  int8_t* l,
3572  int64_t lsize,
3573  double* lbounds,
3574  int64_t lbounds_size,
3575  int32_t li,
3576  int32_t ic1,
3577  int32_t isr1,
3578  int32_t ic2,
3579  int32_t isr2,
3580  int32_t osr) {
3581  if (poly_num_rings > 1) {
3582  return false; // TODO: support polygons with interior rings
3583  }
3584 
3585  auto poly_num_coords = poly_coords_size / compression_unit_size(ic1);
3586  auto lnum_coords = lsize / compression_unit_size(ic2);
3587  auto lnum_points = lnum_coords / 2;
3588  if (li != 0) {
3589  // Statically indexed linestring
3590  if (li < 0 || li > lnum_points) {
3591  li = lnum_points;
3592  }
3593  double lx = coord_x(l, 2 * (li - 1), ic2, isr2, osr);
3594  double ly = coord_y(l, 2 * (li - 1) + 1, ic2, isr2, osr);
3595 
3596  if (poly_bounds) {
3597  if (!box_contains_point(poly_bounds, poly_bounds_size, lx, ly)) {
3598  return false;
3599  }
3600  }
3601  // TODO: should be exclude
3602  return point_in_polygon_winding_number<double, EdgeBehavior::kIncludePointOnEdge>(
3603  poly_coords, poly_num_coords, lx, ly, ic1, isr1, osr);
3604  }
3605 
3606  // Bail out if poly bounding box doesn't contain linestring bounding box
3607  if (poly_bounds && lbounds) {
3608  if (!box_contains_box(poly_bounds, poly_bounds_size, lbounds, lbounds_size)) {
3609  return false;
3610  }
3611  }
3612 
3614  poly_coords, poly_num_coords, l, lnum_coords, ic1, isr1, ic2, isr2, osr);
3615 }
3616 
3618 bool ST_Contains_Polygon_Polygon(int8_t* poly1_coords,
3619  int64_t poly1_coords_size,
3620  int32_t* poly1_ring_sizes,
3621  int64_t poly1_num_rings,
3622  double* poly1_bounds,
3623  int64_t poly1_bounds_size,
3624  int8_t* poly2_coords,
3625  int64_t poly2_coords_size,
3626  int32_t* poly2_ring_sizes,
3627  int64_t poly2_num_rings,
3628  double* poly2_bounds,
3629  int64_t poly2_bounds_size,
3630  int32_t ic1,
3631  int32_t isr1,
3632  int32_t ic2,
3633  int32_t isr2,
3634  int32_t osr) {
3635  // TODO: needs to be extended, cover more cases
3636  // Right now only checking if simple poly1 (no holes) contains poly2's exterior shape
3637  if (poly1_num_rings > 1) {
3638  return false; // TODO: support polygons with interior rings
3639  }
3640 
3641  if (poly1_bounds && poly2_bounds) {
3642  if (!box_contains_box(
3643  poly1_bounds, poly1_bounds_size, poly2_bounds, poly2_bounds_size)) {
3644  return false;
3645  }
3646  }
3647 
3648  int64_t poly2_exterior_ring_coords_size = poly2_coords_size;
3649  if (poly2_num_rings > 0) {
3650  poly2_exterior_ring_coords_size =
3651  2 * poly2_ring_sizes[0] * compression_unit_size(ic2);
3652  }
3653  return ST_Contains_Polygon_LineString(poly1_coords,
3654  poly1_coords_size,
3655  poly1_ring_sizes,
3656  poly1_num_rings,
3657  poly1_bounds,
3658  poly1_bounds_size,
3659  poly2_coords,
3660  poly2_exterior_ring_coords_size,
3661  poly2_bounds,
3662  poly2_bounds_size,
3663  0,
3664  ic1,
3665  isr1,
3666  ic2,
3667  isr2,
3668  osr);
3669 }
3670 
3671 template <typename T, EdgeBehavior TEdgeBehavior>
3673  const int8_t* mpoly_coords,
3674  const int64_t mpoly_coords_size,
3675  const int32_t* mpoly_ring_sizes,
3676  const int64_t mpoly_num_rings,
3677  const int32_t* mpoly_poly_sizes,
3678  const int64_t mpoly_num_polys,
3679  const double* mpoly_bounds,
3680  const int64_t mpoly_bounds_size,
3681  const int8_t* p,
3682  const int64_t psize,
3683  const int32_t ic1,
3684  const int32_t isr1,
3685  const int32_t ic2,
3686  const int32_t isr2,
3687  const int32_t osr) {
3688  if (mpoly_num_polys <= 0) {
3689  return false;
3690  }
3691 
3692  // TODO: mpoly_bounds could contain individual bounding boxes too:
3693  // first two points - box for the entire multipolygon, then a pair for each polygon
3694  if (mpoly_bounds) {
3695  if (!box_contains_point(mpoly_bounds,
3696  mpoly_bounds_size,
3697  coord_x(p, 0, ic2, isr2, osr),
3698  coord_y(p, 1, ic2, isr2, osr))) {
3699  return false;
3700  }
3701  }
3702 
3703  // Set specific poly pointers as we move through the coords/ringsizes/polyrings
3704  // arrays.
3705  auto next_poly_coords = mpoly_coords;
3706  auto next_poly_ring_sizes = mpoly_ring_sizes;
3707 
3708  for (auto poly = 0; poly < mpoly_num_polys; poly++) {
3709  const auto poly_coords = next_poly_coords;
3710  const auto poly_ring_sizes = next_poly_ring_sizes;
3711  const auto poly_num_rings = mpoly_poly_sizes[poly];
3712  // Count number of coords in all of poly's rings, advance ring size pointer.
3713  int32_t poly_num_coords = 0;
3714  for (auto ring = 0; ring < poly_num_rings; ring++) {
3715  poly_num_coords += 2 * *next_poly_ring_sizes++;
3716  }
3717  const auto poly_coords_size = poly_num_coords * compression_unit_size(ic1);
3718  next_poly_coords += poly_coords_size;
3719  // TODO: pass individual bounding boxes for each polygon
3720  if (Contains_Polygon_Point_Impl<T, TEdgeBehavior>(poly_coords,
3721  poly_coords_size,
3722  poly_ring_sizes,
3723  poly_num_rings,
3724  nullptr,
3725  0,
3726  p,
3727  psize,
3728  ic1,
3729  isr1,
3730  ic2,
3731  isr2,
3732  osr)) {
3733  return true;
3734  }
3735  }
3736 
3737  return false;
3738 }
3739 
3741 bool ST_Contains_MultiPolygon_Point(int8_t* mpoly_coords,
3742  int64_t mpoly_coords_size,
3743  int32_t* mpoly_ring_sizes,
3744  int64_t mpoly_num_rings,
3745  int32_t* mpoly_poly_sizes,
3746  int64_t mpoly_num_polys,
3747  double* mpoly_bounds,
3748  int64_t mpoly_bounds_size,
3749  int8_t* p,
3750  int64_t psize,
3751  int32_t ic1,
3752  int32_t isr1,
3753  int32_t ic2,
3754  int32_t isr2,
3755  int32_t osr) {
3756  return Contains_MultiPolygon_Point_Impl<double, EdgeBehavior::kExcludePointOnEdge>(
3757  mpoly_coords,
3758  mpoly_coords_size,
3759  mpoly_ring_sizes,
3760  mpoly_num_rings,
3761  mpoly_poly_sizes,
3762  mpoly_num_polys,
3763  mpoly_bounds,
3764  mpoly_bounds_size,
3765  p,
3766  psize,
3767  ic1,
3768  isr1,
3769  ic2,
3770  isr2,
3771  osr);
3772 }
3773 
3775  int64_t mpoly_coords_size,
3776  int32_t* mpoly_ring_sizes,
3777  int64_t mpoly_num_rings,
3778  int32_t* mpoly_poly_sizes,
3779  int64_t mpoly_num_polys,
3780  double* mpoly_bounds,
3781  int64_t mpoly_bounds_size,
3782  int8_t* p,
3783  int64_t psize,
3784  int32_t ic1,
3785  int32_t isr1,
3786  int32_t ic2,
3787  int32_t isr2,
3788  int32_t osr) {
3789  return Contains_MultiPolygon_Point_Impl<int64_t, EdgeBehavior::kExcludePointOnEdge>(
3790  mpoly_coords,
3791  mpoly_coords_size,
3792  mpoly_ring_sizes,
3793  mpoly_num_rings,
3794  mpoly_poly_sizes,
3795  mpoly_num_polys,
3796  mpoly_bounds,
3797  mpoly_bounds_size,
3798  p,
3799  psize,
3800  ic1,
3801  isr1,
3802  ic2,
3803  isr2,
3804  osr);
3805 }
3806 
3808 bool ST_Contains_MultiPolygon_LineString(int8_t* mpoly_coords,
3809  int64_t mpoly_coords_size,
3810  int32_t* mpoly_ring_sizes,
3811  int64_t mpoly_num_rings,
3812  int32_t* mpoly_poly_sizes,
3813  int64_t mpoly_num_polys,
3814  double* mpoly_bounds,
3815  int64_t mpoly_bounds_size,
3816  int8_t* l,
3817  int64_t lsize,
3818  double* lbounds,
3819  int64_t lbounds_size,
3820  int32_t li,
3821  int32_t ic1,
3822  int32_t isr1,
3823  int32_t ic2,
3824  int32_t isr2,
3825  int32_t osr) {
3826  if (mpoly_num_polys <= 0) {
3827  return false;
3828  }
3829 
3830  auto lnum_coords = lsize / compression_unit_size(ic2);
3831  auto lnum_points = lnum_coords / 2;
3832  if (li != 0) {
3833  // Statically indexed linestring
3834  if (li < 0 || li > lnum_points) {
3835  li = lnum_points;
3836  }
3837  double lx = coord_x(l, 2 * (li - 1), ic2, isr2, osr);
3838  double ly = coord_y(l, 2 * (li - 1) + 1, ic2, isr2, osr);
3839 
3840  if (mpoly_bounds) {
3841  if (!box_contains_point(mpoly_bounds, mpoly_bounds_size, lx, ly)) {
3842  return false;
3843  }
3844  }
3845  auto p = l + li * compression_unit_size(ic2);
3846  auto psize = 2 * compression_unit_size(ic2);
3847  return ST_Contains_MultiPolygon_Point(mpoly_coords,
3848  mpoly_coords_size,
3849  mpoly_ring_sizes,
3850  mpoly_num_rings,
3851  mpoly_poly_sizes,
3852  mpoly_num_polys,
3853  mpoly_bounds,
3854  mpoly_bounds_size,
3855  p,
3856  psize,
3857  ic1,
3858  isr1,
3859  ic2,
3860  isr2,
3861  osr);
3862  }
3863 
3864  if (mpoly_bounds && lbounds) {
3865  if (!box_contains_box(mpoly_bounds, mpoly_bounds_size, lbounds, lbounds_size)) {
3866  return false;
3867  }
3868  }
3869 
3870  // Set specific poly pointers as we move through the coords/ringsizes/polyrings arrays.
3871  auto next_poly_coords = mpoly_coords;
3872  auto next_poly_ring_sizes = mpoly_ring_sizes;
3873 
3874  for (auto poly = 0; poly < mpoly_num_polys; poly++) {
3875  auto poly_coords = next_poly_coords;
3876  auto poly_ring_sizes = next_poly_ring_sizes;
3877  auto poly_num_rings = mpoly_poly_sizes[poly];
3878  // Count number of coords in all of poly's rings, advance ring size pointer.
3879  int32_t poly_num_coords = 0;
3880  for (auto ring = 0; ring < poly_num_rings; ring++) {
3881  poly_num_coords += 2 * *next_poly_ring_sizes++;
3882  }
3883  auto poly_coords_size = poly_num_coords * compression_unit_size(ic1);
3884  next_poly_coords += poly_coords_size;
3885 
3886  if (ST_Contains_Polygon_LineString(poly_coords,
3887  poly_coords_size,
3888  poly_ring_sizes,
3889  poly_num_rings,
3890  nullptr,
3891  0,
3892  l,
3893  lsize,
3894  nullptr,
3895  0,
3896  li,
3897  ic1,
3898  isr1,
3899  ic2,
3900  isr2,
3901  osr)) {
3902  return true;
3903  }
3904  }
3905 
3906  return false;
3907 }
3908 
3909 //
3910 // ST_Intersects
3911 //
3912 
3915  int64_t p1size,
3916  int8_t* p2,
3917  int64_t p2size,
3918  int32_t ic1,
3919  int32_t isr1,
3920  int32_t ic2,
3921  int32_t isr2,
3922  int32_t osr) {
3923  return tol_zero(
3924  ST_Distance_Point_Point(p1, p1size, p2, p2size, ic1, isr1, ic2, isr2, osr));
3925 }
3926 
3929  int64_t psize,
3930  int8_t* l,
3931  int64_t lsize,
3932  double* lbounds,
3933  int64_t lbounds_size,
3934  int32_t li,
3935  int32_t ic1,
3936  int32_t isr1,
3937  int32_t ic2,
3938  int32_t isr2,
3939  int32_t osr) {
3940  double px = coord_x(p, 0, ic1, isr1, osr);
3941  double py = coord_y(p, 1, ic1, isr1, osr);
3942 
3943  auto lnum_coords = lsize / compression_unit_size(ic2);
3944  auto lnum_points = lnum_coords / 2;
3945  if (li != 0) {
3946  // Statically indexed linestring
3947  if (li < 0 || li > lnum_points) {
3948  li = lnum_points;
3949  }
3950  auto p2 = l + li * compression_unit_size(ic2);
3951  auto p2size = 2 * compression_unit_size(ic2);
3952  return tol_zero(
3953  ST_Distance_Point_Point(p2, p2size, p, psize, ic2, isr2, ic1, isr1, osr));
3954  }
3955 
3956  if (lbounds) {
3957  if (!box_contains_point(lbounds, lbounds_size, px, py)) {
3958  return false;
3959  }
3960  }
3962  p, psize, l, lsize, li, ic1, isr1, ic2, isr2, osr, 0.0));
3963 }
3964 
3967  int64_t psize,
3968  int8_t* poly,
3969  int64_t polysize,
3970  int32_t* poly_ring_sizes,
3971  int64_t poly_num_rings,
3972  double* poly_bounds,
3973  int64_t poly_bounds_size,
3974  int32_t ic1,
3975  int32_t isr1,
3976  int32_t ic2,
3977  int32_t isr2,
3978  int32_t osr) {
3979  return Contains_Polygon_Point_Impl<double, EdgeBehavior::kIncludePointOnEdge>(
3980  poly,
3981  polysize,
3982  poly_ring_sizes,
3983  poly_num_rings,
3984  poly_bounds,
3985  poly_bounds_size,
3986  p,
3987  psize,
3988  ic2,
3989  isr2,
3990  ic1,
3991  isr1,
3992  osr);
3993 }
3994 
3997  int64_t psize,
3998  int8_t* mpoly_coords,
3999  int64_t mpoly_coords_size,
4000  int32_t* mpoly_ring_sizes,
4001  int64_t mpoly_num_rings,
4002  int32_t* mpoly_poly_sizes,
4003  int64_t mpoly_num_polys,
4004  double* mpoly_bounds,
4005  int64_t mpoly_bounds_size,
4006  int32_t ic1,
4007  int32_t isr1,
4008  int32_t ic2,
4009  int32_t isr2,
4010  int32_t osr) {
4011  return ST_Contains_MultiPolygon_Point(mpoly_coords,
4012  mpoly_coords_size,
4013  mpoly_ring_sizes,
4014  mpoly_num_rings,
4015  mpoly_poly_sizes,
4016  mpoly_num_polys,
4017  mpoly_bounds,
4018  mpoly_bounds_size,
4019  p,
4020  psize,
4021  ic2,
4022  isr2,
4023  ic1,
4024  isr1,
4025  osr);
4026 }
4027 
4030  int64_t lsize,
4031  double* lbounds,
4032  int64_t lbounds_size,
4033  int32_t li,
4034  int8_t* p,
4035  int64_t psize,
4036  int32_t ic1,
4037  int32_t isr1,
4038  int32_t ic2,
4039  int32_t isr2,
4040  int32_t osr) {
4042  p, psize, l, lsize, lbounds, lbounds_size, li, ic2, isr2, ic1, isr1, osr);
4043 }
4044 
4047  int64_t l1size,
4048  double* l1bounds,
4049  int64_t l1bounds_size,
4050  int32_t l1i,
4051  int8_t* l2,
4052  int64_t l2size,
4053  double* l2bounds,
4054  int64_t l2bounds_size,
4055  int32_t l2i,
4056  int32_t ic1,
4057  int32_t isr1,
4058  int32_t ic2,
4059  int32_t isr2,
4060  int32_t osr) {
4061  auto l2num_coords = l2size / compression_unit_size(ic2);
4062  auto l2num_points = l2num_coords / 2;
4063  if (l2i != 0) {
4064  // Statically indexed linestring
4065  if (l2i < 0 || l2i > l2num_points) {
4066  l2i = l2num_points;
4067  }
4068  auto p2 = l2 + l2i * compression_unit_size(ic2);
4069  auto p2size = 2 * compression_unit_size(ic2);
4071  l1, l1size, l1bounds, l1bounds_size, l1i, p2, p2size, ic1, isr1, ic2, isr2, osr);
4072  }
4073  auto l1num_coords = l1size / compression_unit_size(ic1);
4074  auto l1num_points = l1num_coords / 2;
4075  if (l1i != 0) {
4076  // Statically indexed linestring
4077  if (l1i < 0 || l1i > l1num_points) {
4078  l1i = l1num_points;
4079  }
4080  auto p1 = l1 + l1i * compression_unit_size(ic1);
4081  auto p1size = 2 * compression_unit_size(ic1);
4083  l2, l2size, l2bounds, l2bounds_size, l2i, p1, p1size, ic2, isr2, ic1, isr1, osr);
4084  }
4085 
4086  if (l1bounds && l2bounds) {
4087  if (!box_overlaps_box(l1bounds, l1bounds_size, l2bounds, l2bounds_size)) {
4088  return false;
4089  }
4090  }
4091 
4093  l1, l1size, l1i, l2, l2size, l2i, ic1, isr1, ic2, isr2, osr, 0.0));
4094 }
4095 
4098  int64_t lsize,
4099  double* lbounds,
4100  int64_t lbounds_size,
4101  int32_t li,
4102  int8_t* poly,
4103  int64_t polysize,
4104  int32_t* poly_ring_sizes,
4105  int64_t poly_num_rings,
4106  double* poly_bounds,
4107  int64_t poly_bounds_size,
4108  int32_t ic1,
4109  int32_t isr1,
4110  int32_t ic2,
4111  int32_t isr2,
4112  int32_t osr) {
4113  auto lnum_coords = lsize / compression_unit_size(ic1);
4114  auto lnum_points = lnum_coords / 2;
4115  if (li != 0) {
4116  // Statically indexed linestring
4117  if (li < 0 || li > lnum_points) {
4118  li = lnum_points;
4119  }
4120  auto p = l + li * compression_unit_size(ic1);
4121  auto psize = 2 * compression_unit_size(ic1);
4122  return Contains_Polygon_Point_Impl<double, EdgeBehavior::kIncludePointOnEdge>(
4123  poly,
4124  polysize,
4125  poly_ring_sizes,
4126  poly_num_rings,
4127  poly_bounds,
4128  poly_bounds_size,
4129  p,
4130  psize,
4131  ic2,
4132  isr2,
4133  ic1,
4134  isr1,
4135  osr);
4136  }
4137 
4138  if (lbounds && poly_bounds) {
4139  if (!box_overlaps_box(lbounds, lbounds_size, poly_bounds, poly_bounds_size)) {
4140  return false;
4141  }
4142  }
4143 
4144  // Check for spatial intersection.
4145  // One way to do that would be to start with linestring's first point, if it's inside
4146  // the polygon - it means intersection. Otherwise follow the linestring, segment by
4147  // segment, checking for intersections with polygon rings, bail as soon as we cross into
4148  // the polygon.
4149 
4150  // Or, alternatively, just measure the distance:
4152  lsize,
4153  li,
4154  poly,
4155  polysize,
4156  poly_ring_sizes,
4157  poly_num_rings,
4158  ic1,
4159  isr1,
4160  ic2,
4161  isr2,
4162  osr,
4163  0.0));
4164 }
4165 
4168  int64_t lsize,
4169  double* lbounds,
4170  int64_t lbounds_size,
4171  int32_t li,
4172  int8_t* mpoly_coords,
4173  int64_t mpoly_coords_size,
4174  int32_t* mpoly_ring_sizes,
4175  int64_t mpoly_num_rings,
4176  int32_t* mpoly_poly_sizes,
4177  int64_t mpoly_num_polys,
4178  double* mpoly_bounds,
4179  int64_t mpoly_bounds_size,
4180  int32_t ic1,
4181  int32_t isr1,
4182  int32_t ic2,
4183  int32_t isr2,
4184  int32_t osr) {
4185  auto lnum_coords = lsize / compression_unit_size(ic1);
4186  auto lnum_points = lnum_coords / 2;
4187  if (li != 0) {
4188  // Statically indexed linestring
4189  if (li < 0 || li > lnum_points) {
4190  li = lnum_points;
4191  }
4192  auto p = l + li * compression_unit_size(ic1);
4193  auto psize = 2 * compression_unit_size(ic1);
4194  return ST_Contains_MultiPolygon_Point(mpoly_coords,
4195  mpoly_coords_size,
4196  mpoly_ring_sizes,
4197  mpoly_num_rings,
4198  mpoly_poly_sizes,
4199  mpoly_num_polys,
4200  mpoly_bounds,
4201  mpoly_bounds_size,
4202  p,
4203  psize,
4204  ic2,
4205  isr2,
4206  ic1,
4207  isr1,
4208  osr);
4209  }
4210 
4211  if (lbounds && mpoly_bounds) {
4212  if (!box_overlaps_box(lbounds, lbounds_size, mpoly_bounds, mpoly_bounds_size)) {
4213  return false;
4214  }
4215  }
4216 
4217  // Check for spatial intersection.
4218  // One way to do that would be to start with linestring's first point, if it's inside
4219  // any of the polygons - it means intersection. Otherwise follow the linestring, segment
4220  // by segment, checking for intersections with polygon shapes/holes, bail as soon as we
4221  // cross into a polygon.
4222 
4223  // Or, alternatively, just measure the distance:
4225  lsize,
4226  li,
4227  mpoly_coords,
4228  mpoly_coords_size,
4229  mpoly_ring_sizes,
4230  mpoly_num_rings,
4231  mpoly_poly_sizes,
4232  mpoly_num_polys,
4233  ic1,
4234  isr1,
4235  ic2,
4236  isr2,
4237  osr,
4238  0.0));
4239 }
4240 
4243  int64_t polysize,
4244  int32_t* poly_ring_sizes,
4245  int64_t poly_num_rings,
4246  double* poly_bounds,
4247  int64_t poly_bounds_size,
4248  int8_t* p,
4249  int64_t psize,
4250  int32_t ic1,
4251  int32_t isr1,
4252  int32_t ic2,
4253  int32_t isr2,
4254  int32_t osr) {
4255  return Contains_Polygon_Point_Impl<double, EdgeBehavior::kIncludePointOnEdge>(
4256  poly,
4257  polysize,
4258  poly_ring_sizes,
4259  poly_num_rings,
4260  poly_bounds,
4261  poly_bounds_size,
4262  p,
4263  psize,
4264  ic1,
4265  isr1,
4266  ic2,
4267  isr2,
4268  osr);
4269 }
4270 
4273  int64_t polysize,
4274  int32_t* poly_ring_sizes,
4275  int64_t poly_num_rings,
4276  double* poly_bounds,
4277  int64_t poly_bounds_size,
4278  int8_t* l,
4279  int64_t lsize,
4280  double* lbounds,
4281  int64_t lbounds_size,
4282  int32_t li,
4283  int32_t ic1,
4284  int32_t isr1,
4285  int32_t ic2,
4286  int32_t isr2,
4287  int32_t osr) {
4289  lsize,
4290  lbounds,
4291  lbounds_size,
4292  li,
4293  poly,
4294  polysize,
4295  poly_ring_sizes,
4296  poly_num_rings,
4297  poly_bounds,
4298  poly_bounds_size,
4299  ic2,
4300  isr2,
4301  ic1,
4302  isr1,
4303  osr);
4304 }
4305 
4307 bool ST_Intersects_Polygon_Polygon(int8_t* poly1_coords,
4308  int64_t poly1_coords_size,
4309  int32_t* poly1_ring_sizes,
4310  int64_t poly1_num_rings,
4311  double* poly1_bounds,
4312  int64_t poly1_bounds_size,
4313  int8_t* poly2_coords,
4314  int64_t poly2_coords_size,
4315  int32_t* poly2_ring_sizes,
4316  int64_t poly2_num_rings,
4317  double* poly2_bounds,
4318  int64_t poly2_bounds_size,
4319  int32_t ic1,
4320  int32_t isr1,
4321  int32_t ic2,
4322  int32_t isr2,
4323  int32_t osr) {
4324  if (poly1_bounds && poly2_bounds) {
4325  if (!box_overlaps_box(
4326  poly1_bounds, poly1_bounds_size, poly2_bounds, poly2_bounds_size)) {
4327  return false;
4328  }
4329  }
4330 
4331  return tol_zero(ST_Distance_Polygon_Polygon(poly1_coords,
4332  poly1_coords_size,
4333  poly1_ring_sizes,
4334  poly1_num_rings,
4335  poly2_coords,
4336  poly2_coords_size,
4337  poly2_ring_sizes,
4338  poly2_num_rings,
4339  ic1,
4340  isr1,
4341  ic2,
4342  isr2,
4343  osr,
4344  0.0));
4345 }
4346 
4348 bool ST_Intersects_Polygon_MultiPolygon(int8_t* poly_coords,
4349  int64_t poly_coords_size,
4350  int32_t* poly_ring_sizes,
4351  int64_t poly_num_rings,
4352  double* poly_bounds,
4353  int64_t poly_bounds_size,
4354  int8_t* mpoly_coords,
4355  int64_t mpoly_coords_size,
4356  int32_t* mpoly_ring_sizes,
4357  int64_t mpoly_num_rings,
4358  int32_t* mpoly_poly_sizes,
4359  int64_t mpoly_num_polys,
4360  double* mpoly_bounds,
4361  int64_t mpoly_bounds_size,
4362  int32_t ic1,
4363  int32_t isr1,
4364  int32_t ic2,
4365  int32_t isr2,
4366  int32_t osr) {
4367  if (poly_bounds && mpoly_bounds) {
4368  if (!box_overlaps_box(
4369  poly_bounds, poly_bounds_size, mpoly_bounds, mpoly_bounds_size)) {
4370  return false;
4371  }
4372  }
4373 
4374  return tol_zero(ST_Distance_Polygon_MultiPolygon(poly_coords,
4375  poly_coords_size,
4376  poly_ring_sizes,
4377  poly_num_rings,
4378  mpoly_coords,
4379  mpoly_coords_size,
4380  mpoly_ring_sizes,
4381  mpoly_num_rings,
4382  mpoly_poly_sizes,
4383  mpoly_num_polys,
4384  ic1,
4385  isr1,
4386  ic2,
4387  isr2,
4388  osr,
4389  0.0));
4390 }
4391 
4393 bool ST_Intersects_MultiPolygon_Point(int8_t* mpoly_coords,
4394  int64_t mpoly_coords_size,
4395  int32_t* mpoly_ring_sizes,
4396  int64_t mpoly_num_rings,
4397  int32_t* mpoly_poly_sizes,
4398  int64_t mpoly_num_polys,
4399  double* mpoly_bounds,
4400  int64_t mpoly_bounds_size,
4401  int8_t* p,
4402  int64_t psize,
4403  int32_t ic1,
4404  int32_t isr1,
4405  int32_t ic2,
4406  int32_t isr2,
4407  int32_t osr) {
4408  return Contains_MultiPolygon_Point_Impl<double, EdgeBehavior::kIncludePointOnEdge>(
4409  mpoly_coords,
4410  mpoly_coords_size,
4411  mpoly_ring_sizes,
4412  mpoly_num_rings,
4413  mpoly_poly_sizes,
4414  mpoly_num_polys,
4415  mpoly_bounds,
4416  mpoly_bounds_size,
4417  p,
4418  psize,
4419  ic1,
4420  isr1,
4421  ic2,
4422  isr2,
4423  osr);
4424 }
4425 
4427 bool ST_Intersects_MultiPolygon_LineString(int8_t* mpoly_coords,
4428  int64_t mpoly_coords_size,
4429  int32_t* mpoly_ring_sizes,
4430  int64_t mpoly_num_rings,
4431  int32_t* mpoly_poly_sizes,
4432  int64_t mpoly_num_polys,
4433  double* mpoly_bounds,
4434  int64_t mpoly_bounds_size,
4435  int8_t* l,
4436  int64_t lsize,
4437  double* lbounds,
4438  int64_t lbounds_size,
4439  int32_t li,
4440  int32_t ic1,
4441  int32_t isr1,
4442  int32_t ic2,
4443  int32_t isr2,
4444  int32_t osr) {
4446  lsize,
4447  lbounds,
4448  lbounds_size,
4449  li,
4450  mpoly_coords,
4451  mpoly_coords_size,
4452  mpoly_ring_sizes,
4453  mpoly_num_rings,
4454  mpoly_poly_sizes,
4455  mpoly_num_polys,
4456  mpoly_bounds,
4457  mpoly_bounds_size,
4458  ic2,
4459  isr2,
4460  ic1,
4461  isr1,
4462  osr);
4463 }
4464 
4466 bool ST_Intersects_MultiPolygon_Polygon(int8_t* mpoly_coords,
4467  int64_t mpoly_coords_size,
4468  int32_t* mpoly_ring_sizes,
4469  int64_t mpoly_num_rings,
4470  int32_t* mpoly_poly_sizes,
4471  int64_t mpoly_num_polys,
4472  double* mpoly_bounds,
4473  int64_t mpoly_bounds_size,
4474  int8_t* poly_coords,
4475  int64_t poly_coords_size,
4476  int32_t* poly_ring_sizes,
4477  int64_t poly_num_rings,
4478  double* poly_bounds,
4479  int64_t poly_bounds_size,
4480  int32_t ic1,
4481  int32_t isr1,
4482  int32_t ic2,
4483  int32_t isr2,
4484  int32_t osr) {
4485  return ST_Intersects_Polygon_MultiPolygon(poly_coords,
4486  poly_coords_size,
4487  poly_ring_sizes,
4488  poly_num_rings,
4489  poly_bounds,
4490  poly_bounds_size,
4491  mpoly_coords,
4492  mpoly_coords_size,
4493  mpoly_ring_sizes,
4494  mpoly_num_rings,
4495  mpoly_poly_sizes,
4496  mpoly_num_polys,
4497  mpoly_bounds,
4498  mpoly_bounds_size,
4499  ic2,
4500  isr2,
4501  ic1,
4502  isr1,
4503  osr);
4504 }
4505 
4507 bool ST_Intersects_MultiPolygon_MultiPolygon(int8_t* mpoly1_coords,
4508  int64_t mpoly1_coords_size,
4509  int32_t* mpoly1_ring_sizes,
4510  int64_t mpoly1_num_rings,
4511  int32_t* mpoly1_poly_sizes,
4512  int64_t mpoly1_num_polys,
4513  double* mpoly1_bounds,
4514  int64_t mpoly1_bounds_size,
4515  int8_t* mpoly2_coords,
4516  int64_t mpoly2_coords_size,
4517  int32_t* mpoly2_ring_sizes,
4518  int64_t mpoly2_num_rings,
4519  int32_t* mpoly2_poly_sizes,
4520  int64_t mpoly2_num_polys,
4521  double* mpoly2_bounds,
4522  int64_t mpoly2_bounds_size,
4523  int32_t ic1,
4524  int32_t isr1,
4525  int32_t ic2,
4526  int32_t isr2,
4527  int32_t osr) {
4528  if (mpoly1_bounds && mpoly2_bounds) {
4529  if (!box_overlaps_box(
4530  mpoly1_bounds, mpoly1_bounds_size, mpoly2_bounds, mpoly2_bounds_size)) {
4531  return false;
4532  }
4533  }
4534 
4535  return tol_zero(ST_Distance_MultiPolygon_MultiPolygon(mpoly1_coords,
4536  mpoly1_coords_size,
4537  mpoly1_ring_sizes,
4538  mpoly1_num_rings,
4539  mpoly1_poly_sizes,
4540  mpoly1_num_polys,
4541  mpoly2_coords,
4542  mpoly2_coords_size,
4543  mpoly2_ring_sizes,
4544  mpoly2_num_rings,
4545  mpoly2_poly_sizes,
4546  mpoly2_num_polys,
4547  ic1,
4548  isr1,
4549  ic2,
4550  isr2,
4551  osr,
4552  0.0));
4553 }
4554 
4555 //
4556 // Accessors for poly bounds and render group for in-situ poly render queries
4557 //
4558 // The MapD_* varieties are deprecated and renamed to "OmniSci_Geo*"
4559 // There may be some clients out there who are playing with the MapD_* so leaving
4560 // them for backwards compatibility.
4561 //
4562 
4564 int64_t OmniSci_Geo_PolyBoundsPtr(double* bounds, int64_t size) {
4565  return reinterpret_cast<int64_t>(bounds);
4566 }
4567 
4569 int32_t OmniSci_Geo_PolyRenderGroup(int32_t render_group) {
4570  return render_group;
4571 }
4572 
4574 int64_t MapD_GeoPolyBoundsPtr(double* bounds, int64_t size) {
4575  return OmniSci_Geo_PolyBoundsPtr(bounds, size);
4576 }
4577 
4579 int32_t MapD_GeoPolyRenderGroup(int32_t render_group) {
4580  return OmniSci_Geo_PolyRenderGroup(render_group);
4581 }
4582 
4584 double convert_meters_to_pixel_width(const double meters,
4585  int8_t* p,
4586  const int64_t psize,
4587  const int32_t ic,
4588  const int32_t isr,
4589  const int32_t osr,
4590  const double min_lon,
4591  const double max_lon,
4592  const int32_t img_width,
4593  const double min_width) {
4594  const double const1 = 0.017453292519943295769236907684886;
4595  const double const2 = 6372797.560856;
4596  const auto lon = decompress_coord(p, 0, ic, true);
4597  const auto lat = decompress_coord(p, 1, ic, false);
4598  double t1 = sinf(meters / (2.0 * const2));
4599  double t2 = cosf(const1 * lat);
4600  const double newlon = lon - (2.0 * asinf(t1 / t2)) / const1;
4601  t1 = transform_coord(lon, isr, osr, true);
4602  t2 = transform_coord(newlon, isr, osr, true);
4603  const double min_domain_x = transform_coord(min_lon, isr, osr, true);
4604  const double max_domain_x = transform_coord(max_lon, isr, osr, true);
4605  const double domain_diff = max_domain_x - min_domain_x;
4606  t1 = ((t1 - min_domain_x) / domain_diff) * static_cast<double>(img_width);
4607  t2 = ((t2 - min_domain_x) / domain_diff) * static_cast<double>(img_width);
4608 
4609  // TODO(croot): need to account for edge cases, such as getting close to the poles.
4610  const double sz = fabs(t1 - t2);
4611  return (sz < min_width ? min_width : sz);
4612 }
4613 
4615 double convert_meters_to_pixel_height(const double meters,
4616  int8_t* p,
4617  const int64_t psize,
4618  const int32_t ic,
4619  const int32_t isr,
4620  const int32_t osr,
4621  const double min_lat,
4622  const double max_lat,
4623  const int32_t img_height,
4624  const double min_height) {
4625  const double const1 = 0.017453292519943295769236907684886;
4626  const double const2 = 6372797.560856;
4627  const auto lat = decompress_coord(p, 1, ic, false);
4628  const double latdiff = meters / (const1 * const2);
4629  const double newlat =
4630  (lat < 0) ? lat + latdiff : lat - latdiff; // assumes a lat range of [-90, 90]
4631  double t1 = transform_coord(lat, isr, osr, false);
4632  double t2 = transform_coord(newlat, isr, osr, false);
4633  const double min_domain_y = transform_coord(min_lat, isr, osr, false);
4634  const double max_domain_y = transform_coord(max_lat, isr, osr, false);
4635  const double domain_diff = max_domain_y - min_domain_y;
4636  t1 = ((t1 - min_domain_y) / domain_diff) * static_cast<double>(img_height);
4637  t2 = ((t2 - min_domain_y) / domain_diff) * static_cast<double>(img_height);
4638 
4639  // TODO(croot): need to account for edge cases, such as getting close to the poles.
4640  const double sz = fabs(t1 - t2);
4641  return (sz < min_height ? min_height : sz);
4642 }
4643 
4645  const int64_t psize,
4646  const int32_t ic,
4647  const double min_lon,
4648  const double max_lon,
4649  const double min_lat,
4650  const double max_lat) {
4651  const auto lon = decompress_coord(p, 0, ic, true);
4652  const auto lat = decompress_coord(p, 1, ic, false);
4653  return !(lon < min_lon || lon > max_lon || lat < min_lat || lat > max_lat);
4654 }
4655 
4657  const int64_t psize,
4658  const int32_t ic,
4659  const double meters,
4660  const double min_lon,
4661  const double max_lon,
4662  const double min_lat,
4663  const double max_lat) {
4664  const double const1 = 0.017453292519943295769236907684886;
4665  const double const2 = 6372797.560856;
4666  const auto lon = decompress_coord(p, 0, ic, true);
4667  const auto lat = decompress_coord(p, 1, ic, false);
4668  const double latdiff = meters / (const1 * const2);
4669  const double t1 = sinf(meters / (2.0 * const2));
4670  const double t2 = cosf(const1 * lat);
4671  const double londiff = (2.0 * asinf(t1 / t2)) / const1;
4672  return !(lon + londiff < min_lon || lon - londiff > max_lon ||
4673  lat + latdiff < min_lat || lat - latdiff > max_lat);
4674 }
EXTENSION_NOINLINE double ST_Distance_LineString_LineString_Geodesic(int8_t *l1, int64_t l1size, int32_t l1index, int8_t *l2, int64_t l2size, int32_t l2index, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, 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, int64_t lsize, int32_t ic, int32_t isr, int32_t osr, bool geodesic, bool check_closed)
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_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)
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_NOINLINE double ST_Centroid_LineString(int8_t *coords, int64_t coords_sz, int32_t ic, int32_t isr, int32_t osr, bool ycoord)
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 double area_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)
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_Perimeter_Polygon(int8_t *poly, int64_t polysize, int32_t *poly_ring_sizes, int64_t poly_num_rings, int32_t ic, int32_t isr, int32_t osr)
EXTENSION_NOINLINE bool ST_Contains_LineString_Polygon(int8_t *l, int64_t lsize, double *lbounds, int64_t lbounds_size, int32_t li, 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_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)
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_MultiPolygon(int8_t *l, int64_t lsize, double *lbounds, int64_t lbounds_size, int32_t li, 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_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_NOINLINE double ST_X_LineString(int8_t *l, int64_t lsize, int32_t lindex, int32_t ic, int32_t isr, int32_t osr)
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 ALWAYS_INLINE double perimeter_multipolygon(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 ic, int32_t isr, int32_t osr, bool geodesic)
tuple r
Definition: test_fsi.py:16
DEVICE ALWAYS_INLINE double max_distance_point_linestring(int8_t *p, int64_t psize, int8_t *l, int64_t lsize, int32_t lindex, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr, bool check_closed)
DEVICE double decompress_lattitude_coord_geoint32(const int32_t compressed)
EXTENSION_NOINLINE double ST_Area_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)
EXTENSION_NOINLINE double ST_Distance_LineString_MultiPolygon(int8_t *l, int64_t lsize, int32_t lindex, 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_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 li, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr)
DEVICE ALWAYS_INLINE double area_triangle(double x1, double y1, double x2, double y2, double x3, double y3)
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 li, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr, double threshold)
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 double ST_Distance_LineString_Point(int8_t *l, int64_t lsize, int32_t lindex, 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 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 li, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, 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 ST_Perimeter_MultiPolygon_Geodesic(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 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
EXTENSION_INLINE bool ST_DWithin_LineString_MultiPolygon(int8_t *l1, int64_t l1size, double *l1bounds, int64_t l1bounds_size, int32_t l1index, 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 double ST_Y_LineString(int8_t *l, int64_t lsize, int32_t lindex, int32_t ic, int32_t isr, int32_t osr)
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)
EXTENSION_NOINLINE double ST_Distance_Point_LineString(int8_t *p, int64_t psize, int8_t *l, int64_t lsize, int32_t lindex, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr, double threshold)
#define DEVICE
EXTENSION_NOINLINE bool ST_Intersects_LineString_Linestring(int8_t *l1, int64_t l1size, double *l1bounds, int64_t l1bounds_size, int32_t l1i, int8_t *l2, int64_t l2size, double *l2bounds, int64_t l2bounds_size, int32_t l2i, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr)
#define EXTENSION_NOINLINE
Definition: OmniSciTypes.h:27
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 li, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr)
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 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_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 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_NOINLINE double ST_Centroid_MultiPolygon(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 ic, int32_t isr, int32_t osr, bool ycoord)
EXTENSION_NOINLINE double ST_YMin(int8_t *coords, int64_t size, int32_t ic, int32_t isr, int32_t osr)
EXTENSION_INLINE bool ST_Intersects_LineString_Point(int8_t *l, int64_t lsize, double *lbounds, int64_t lbounds_size, int32_t li, int8_t *p, int64_t psize, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, 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_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)
EXTENSION_INLINE double ST_Area_Polygon_Geodesic(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)
EXTENSION_NOINLINE bool ST_Intersects_LineString_Polygon(int8_t *l, int64_t lsize, double *lbounds, int64_t lbounds_size, int32_t li, 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 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_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)
EXTENSION_NOINLINE double ST_Centroid_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, bool ycoord)
#define TOLERANCE_DEFAULT
EXTENSION_NOINLINE double ST_Distance_LineString_LineString(int8_t *l1, int64_t l1size, int32_t l1index, int8_t *l2, int64_t l2size, int32_t l2index, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr, double threshold)
EXTENSION_INLINE double ST_Distance_LineString_Point_Geodesic(int8_t *l, int64_t lsize, int32_t lindex, int8_t *p, int64_t psize, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr)
DEVICE bool line_intersects_line(double l11x, double l11y, double l12x, double l12y, double l21x, double l21y, double l22x, double l22y)
EXTENSION_NOINLINE double ST_Area_MultiPolygon(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 ic, int32_t isr, 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)
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)
#define UNLIKELY(x)
Definition: likely.h:25
EXTENSION_INLINE bool ST_Contains_LineString_Point(int8_t *l, int64_t lsize, double *lbounds, int64_t lbounds_size, int32_t li, 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_Contains_LineString_LineString(int8_t *l1, int64_t l1size, double *l1bounds, int64_t l1bounds_size, int32_t l1i, int8_t *l2, int64_t l2size, double *l2bounds, int64_t l2bounds_size, int32_t l2i, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr)
EXTENSION_INLINE bool ST_DWithin_LineString_Polygon(int8_t *l1, int64_t l1size, double *l1bounds, int64_t l1bounds_size, int32_t l1index, 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 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)
EXTENSION_NOINLINE double ST_MaxDistance_Point_LineString(int8_t *p, int64_t psize, int8_t *l, int64_t lsize, int32_t lindex, 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_NOINLINE double ST_Perimeter_MultiPolygon(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 ic, int32_t isr, int32_t osr)
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)
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 double ST_Area_MultiPolygon_Geodesic(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 ic, int32_t isr, int32_t osr)
DEVICE ALWAYS_INLINE bool tol_zero(const double x, const double tolerance=TOLERANCE_DEFAULT)
EXTENSION_INLINE int32_t ST_NRings(int32_t *poly_ring_sizes, int64_t poly_num_rings)
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 double ST_Perimeter_Polygon_Geodesic(int8_t *poly, int64_t polysize, int32_t *poly_ring_sizes, int64_t poly_num_rings, int32_t ic, int32_t isr, int32_t osr)
DEVICE ALWAYS_INLINE double distance_point_point_squared(double p1x, double p1y, double p2x, double p2y)
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 lindex, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr, double threshold)
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 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)
EXTENSION_NOINLINE double ST_Distance_Point_LineString_Geodesic(int8_t *p, int64_t psize, int8_t *l, int64_t lsize, int32_t lindex, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr)
DEVICE ALWAYS_INLINE double distance_point_linestring(int8_t *p, int64_t psize, int8_t *l, int64_t lsize, int32_t lindex, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr, bool check_closed, double threshold)
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_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)
#define EXTENSION_INLINE
Definition: OmniSciTypes.h:26
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)
EXTENSION_INLINE bool ST_DWithin_LineString_LineString(int8_t *l1, int64_t l1size, double *l1bounds, int64_t l1bounds_size, int32_t l1index, int8_t *l2, int64_t l2size, double *l2bounds, int64_t l2bounds_size, int32_t l2index, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr, double distance_within)
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 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 l2index, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr, double distance_within)
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_MaxDistance_LineString_Point(int8_t *l, int64_t lsize, int32_t lindex, int8_t *p, int64_t psize, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, 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)
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 double ST_Centroid_Point(int8_t *p, int64_t psize, int32_t ic, int32_t isr, int32_t osr, bool ycoord)
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 double ST_Distance_LineString_Polygon(int8_t *l, int64_t lsize, int32_t lindex, 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)
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_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 li, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr)
EXTENSION_INLINE int32_t ST_NPoints(int8_t *coords, int64_t coords_sz, int32_t ic)
EXTENSION_NOINLINE double conv_4326_900913_y(const double y)
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)
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)
EXTENSION_NOINLINE double ST_Distance_Point_ClosedLineString(int8_t *p, int64_t psize, int8_t *l, int64_t lsize, int32_t lindex, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr, double threshold)
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_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_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 li, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr)
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 li, 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)