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