OmniSciDB  a987f07e93
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) {
30  }
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);
87  if constexpr (C == Coords::X) {
89  } else {
91  }
92  }
93  auto double_coords = reinterpret_cast<const double*>(data);
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
1179 double ST_Length_MultiLineString(int8_t* coords,
1180  int64_t coords_sz,
1181  int8_t* linestring_sizes_in,
1182  int64_t linestring_sizes_sz,
1183  int32_t ic,
1184  int32_t isr,
1185  int32_t osr) {
1186  if (linestring_sizes_sz <= 0) {
1187  return 0.0;
1188  }
1189  double mls_length = 0.0;
1190
1191  auto next_linestring_coords = coords;
1192
1193  for (auto l = 0; l < linestring_sizes_sz; l++) {
1194  auto linestring_coords = next_linestring_coords;
1195  auto linestring_sizes = reinterpret_cast<int32_t*>(linestring_sizes_in);
1196  auto linestring_num_points = linestring_sizes[l];
1197  auto linestring_num_coords = 2 * linestring_num_points;
1198  auto linestring_coords_size = linestring_num_coords * compression_unit_size(ic);
1199  next_linestring_coords += linestring_coords_size;
1200  double ls_length = length_linestring(
1201  linestring_coords, linestring_coords_size, ic, isr, osr, false, false);
1202  mls_length += ls_length;
1203  }
1204
1205  return mls_length;
1206 }
1207
1208 //
1209 // ST_Perimeter
1210 //
1211
1213 double ST_Perimeter_Polygon(int8_t* poly,
1214  int32_t polysize,
1215  int8_t* poly_ring_sizes,
1216  int32_t poly_num_rings,
1217  int32_t ic,
1218  int32_t isr,
1219  int32_t osr) {
1220  if (poly_num_rings <= 0) {
1221  return 0.0;
1222  }
1223
1224  auto exterior_ring_num_coords = poly_ring_sizes[0] * 2;
1225  auto exterior_ring_coords_size = exterior_ring_num_coords * compression_unit_size(ic);
1226
1227  return length_linestring(poly, exterior_ring_coords_size, ic, isr, osr, false, true);
1228 }
1229
1231 double ST_Perimeter_Polygon_Geodesic(int8_t* poly,
1232  int32_t polysize,
1233  int8_t* poly_ring_sizes_in,
1234  int32_t poly_num_rings,
1235  int32_t ic,
1236  int32_t isr,
1237  int32_t osr) {
1238  if (poly_num_rings <= 0) {
1239  return 0.0;
1240  }
1241
1242  auto poly_ring_sizes = reinterpret_cast<int32_t*>(poly_ring_sizes_in);
1243  auto exterior_ring_num_coords = poly_ring_sizes[0] * 2;
1244  auto exterior_ring_coords_size = exterior_ring_num_coords * compression_unit_size(ic);
1245
1246  return length_linestring(poly, exterior_ring_coords_size, ic, isr, osr, true, true);
1247 }
1248
1249 DEVICE ALWAYS_INLINE double perimeter_multipolygon(int8_t* mpoly_coords,
1250  int32_t mpoly_coords_size,
1251  int8_t* mpoly_ring_sizes_in,
1252  int32_t mpoly_num_rings,
1253  int8_t* mpoly_poly_sizes,
1254  int32_t mpoly_num_polys,
1255  int32_t ic,
1256  int32_t isr,
1257  int32_t osr,
1258  bool geodesic) {
1259  if (mpoly_num_polys <= 0 || mpoly_num_rings <= 0) {
1260  return 0.0;
1261  }
1262
1263  double perimeter = 0.0;
1264
1265  auto mpoly_ring_sizes = reinterpret_cast<int32_t*>(mpoly_ring_sizes_in);
1266
1267  // Set specific poly pointers as we move through the coords/ringsizes/polyrings arrays.
1268  auto next_poly_coords = mpoly_coords;
1269  auto next_poly_ring_sizes = mpoly_ring_sizes;
1270
1271  for (auto poly = 0; poly < mpoly_num_polys; poly++) {
1272  auto poly_coords = next_poly_coords;
1273  auto poly_ring_sizes = next_poly_ring_sizes;
1274  auto poly_num_rings = reinterpret_cast<int32_t*>(mpoly_poly_sizes)[poly];
1275  // Count number of coords in all of poly's rings, advance ring size pointer.
1276  int32_t poly_num_coords = 0;
1277  for (auto ring = 0; ring < poly_num_rings; ring++) {
1278  poly_num_coords += 2 * *next_poly_ring_sizes++;
1279  }
1280  auto poly_coords_size = poly_num_coords * compression_unit_size(ic);
1281  next_poly_coords += poly_coords_size;
1282
1283  auto exterior_ring_num_coords = poly_ring_sizes[0] * 2;
1284  auto exterior_ring_coords_size = exterior_ring_num_coords * compression_unit_size(ic);
1285
1286  perimeter += length_linestring(
1287  poly_coords, exterior_ring_coords_size, ic, isr, osr, geodesic, true);
1288  }
1289
1290  return perimeter;
1291 }
1292
1294 double ST_Perimeter_MultiPolygon(int8_t* mpoly_coords,
1295  int32_t mpoly_coords_size,
1296  int8_t* mpoly_ring_sizes,
1297  int32_t mpoly_num_rings,
1298  int8_t* mpoly_poly_sizes,
1299  int32_t mpoly_num_polys,
1300  int32_t ic,
1301  int32_t isr,
1302  int32_t osr) {
1303  return perimeter_multipolygon(mpoly_coords,
1304  mpoly_coords_size,
1305  mpoly_ring_sizes,
1306  mpoly_num_rings,
1307  mpoly_poly_sizes,
1308  mpoly_num_polys,
1309  ic,
1310  isr,
1311  osr,
1312  false);
1313 }
1314
1316 double ST_Perimeter_MultiPolygon_Geodesic(int8_t* mpoly_coords,
1317  int32_t mpoly_coords_size,
1318  int8_t* mpoly_ring_sizes,
1319  int32_t mpoly_num_rings,
1320  int8_t* mpoly_poly_sizes,
1321  int32_t mpoly_num_polys,
1322  int32_t ic,
1323  int32_t isr,
1324  int32_t osr) {
1325  return perimeter_multipolygon(mpoly_coords,
1326  mpoly_coords_size,
1327  mpoly_ring_sizes,
1328  mpoly_num_rings,
1329  mpoly_poly_sizes,
1330  mpoly_num_polys,
1331  ic,
1332  isr,
1333  osr,
1334  true);
1335 }
1336
1337 //
1338 // ST_Area
1339 //
1340
1342  double y1,
1343  double x2,
1344  double y2,
1345  double x3,
1346  double y3) {
1347  return (x1 * y2 - x2 * y1 + x3 * y1 - x1 * y3 + x2 * y3 - x3 * y2) / 2.0;
1348 }
1349
1350 DEVICE ALWAYS_INLINE double area_ring(int8_t* ring,
1351  int64_t ringsize,
1352  int32_t ic,
1353  int32_t isr,
1354  int32_t osr) {
1355  auto ring_num_coords = ringsize / compression_unit_size(ic);
1356
1357  if (ring_num_coords < 6) {
1358  return 0.0;
1359  }
1360
1361  double area = 0.0;
1362
1363  Point2D p1 = get_point(ring, 0, ic, isr, osr);
1364  Point2D p2 = get_point(ring, 2, ic, isr, osr);
1365  for (int32_t i = 4; i < ring_num_coords; i += 2) {
1366  Point2D p3 = get_point(ring, i, ic, isr, osr);
1367  area += area_triangle(p1.x, p1.y, p2.x, p2.y, p3.x, p3.y);
1368  p2 = p3;
1369  }
1370  return area;
1371 }
1372
1373 DEVICE ALWAYS_INLINE double area_polygon(int8_t* poly_coords,
1374  int32_t poly_coords_size,
1375  int8_t* poly_ring_sizes_in,
1376  int32_t poly_num_rings,
1377  int32_t ic,
1378  int32_t isr,
1379  int32_t osr) {
1380  if (poly_num_rings <= 0) {
1381  return 0.0;
1382  }
1383
1384  double area = 0.0;
1385  auto ring_coords = poly_coords;
1386  auto poly_ring_sizes = reinterpret_cast<int32_t*>(poly_ring_sizes_in);
1387
1388  // Add up the areas of all rings.
1389  // External ring is CCW, open - positive area.
1390  // Internal rings (holes) are CW, open - negative areas.
1391  for (auto r = 0; r < poly_num_rings; r++) {
1392  auto ring_coords_size = poly_ring_sizes[r] * 2 * compression_unit_size(ic);
1393  area += area_ring(ring_coords, ring_coords_size, ic, isr, osr);
1394  // Advance to the next ring.
1395  ring_coords += ring_coords_size;
1396  }
1397  return area;
1398 }
1399
1401 double ST_Area_Polygon(int8_t* poly_coords,
1402  int32_t poly_coords_size,
1403  int8_t* poly_ring_sizes,
1404  int32_t poly_num_rings,
1405  int32_t ic,
1406  int32_t isr,
1407  int32_t osr) {
1408  return area_polygon(
1409  poly_coords, poly_coords_size, poly_ring_sizes, poly_num_rings, ic, isr, osr);
1410 }
1411
1413 double ST_Area_MultiPolygon(int8_t* mpoly_coords,
1414  int32_t mpoly_coords_size,
1415  int8_t* mpoly_ring_sizes,
1416  int32_t mpoly_num_rings,
1417  int8_t* mpoly_poly_sizes_in,
1418  int32_t mpoly_num_polys,
1419  int32_t ic,
1420  int32_t isr,
1421  int32_t osr) {
1422  if (mpoly_num_rings <= 0 || mpoly_num_polys <= 0) {
1423  return 0.0;
1424  }
1425
1426  double area = 0.0;
1427
1428  auto mpoly_poly_sizes = reinterpret_cast<int32_t*>(mpoly_poly_sizes_in);
1429
1430  // Set specific poly pointers as we move through the coords/ringsizes/polyrings arrays.
1431  auto next_poly_coords = mpoly_coords;
1432  auto next_poly_ring_sizes = reinterpret_cast<int32_t*>(mpoly_ring_sizes);
1433
1434  for (auto poly = 0; poly < mpoly_num_polys; poly++) {
1435  auto poly_coords = next_poly_coords;
1436  auto poly_ring_sizes = next_poly_ring_sizes;
1437  auto poly_num_rings = mpoly_poly_sizes[poly];
1438  // Count number of coords in all of poly's rings, advance ring size pointer.
1439  int32_t poly_num_coords = 0;
1440  for (auto ring = 0; ring < poly_num_rings; ring++) {
1441  poly_num_coords += 2 * *next_poly_ring_sizes++;
1442  }
1443  auto poly_coords_size = poly_num_coords * compression_unit_size(ic);
1444  next_poly_coords += poly_coords_size;
1445
1446  area += area_polygon(poly_coords,
1447  poly_coords_size,
1448  reinterpret_cast<int8_t*>(poly_ring_sizes),
1449  poly_num_rings,
1450  ic,
1451  isr,
1452  osr);
1453  }
1454  return area;
1455 }
1456
1457 //
1458 // ST_Centroid
1459 //
1460
1461 // Point centroid
1462 // - Single point - drop, otherwise calculate average coordinates for all points
1463 // LineString centroid
1464 // - take midpoint of each line segment
1465 // - calculate average coordinates of all midpoints weighted by segment length
1466 // Polygon, MultiPolygon centroid, holes..
1467 // - weighted sum of centroids of triangles coming out of area decomposition
1468 //
1469 // On zero area fall back to linestring centroid
1470 // On zero length fall back to point centroid
1471
1473 void ST_Centroid_Point(int8_t* p,
1474  int32_t psize,
1475  int32_t ic,
1476  int32_t isr,
1477  int32_t osr,
1478  double* point_centroid) {
1479  Point2D const centroid = get_point(p, 0, ic, isr, osr);
1480  point_centroid[0] = centroid.x;
1481  point_centroid[1] = centroid.y;
1482 }
1483
1485 void ST_Centroid_MultiPoint(int8_t* mp,
1486  int32_t mpsize,
1487  int32_t ic,
1488  int32_t isr,
1489  int32_t osr,
1490  double* multipoint_centroid) {
1491  auto mp_num_coords = mpsize / compression_unit_size(ic);
1492  double x = 0.0;
1493  double y = 0.0;
1494  for (int32_t i = 0; i < mp_num_coords; i += 2) {
1495  Point2D mpp = get_point(mp, i, ic, isr, osr);
1496  x += mpp.x; // In case of overflows on large geometries,
1497  y += mpp.y; // move division into the loop
1498  }
1499  auto mp_num_points = mp_num_coords >> 1;
1500  multipoint_centroid[0] = x / mp_num_points;
1501  multipoint_centroid[1] = y / mp_num_points;
1502 }
1503
1505  double y1,
1506  double x2,
1507  double y2,
1508  double* length,
1509  double* linestring_centroid_sum) {
1510  double ldist = distance_point_point(x1, y1, x2, y2);
1511  *length += ldist;
1512  double segment_midpoint_x = (x1 + x2) / 2.0;
1513  double segment_midpoint_y = (y1 + y2) / 2.0;
1514  linestring_centroid_sum[0] += ldist * segment_midpoint_x;
1515  linestring_centroid_sum[1] += ldist * segment_midpoint_y;
1516  return true;
1517 }
1518
1520  int64_t lsize,
1521  int32_t ic,
1522  int32_t isr,
1523  int32_t osr,
1524  bool closed,
1525  double* total_length,
1526  double* linestring_centroid_sum,
1527  int64_t* num_points,
1528  double* point_centroid_sum) {
1529  auto l_num_coords = lsize / compression_unit_size(ic);
1530  double length = 0.0;
1531  Point2D const l0 = get_point(l, 0, ic, isr, osr);
1532  Point2D l2 = l0;
1533  for (int32_t i = 2; i < l_num_coords; i += 2) {
1534  Point2D const l1 = l2;
1535  l2 = get_point(l, i, ic, isr, osr);
1536  centroid_add_segment(l1.x, l1.y, l2.x, l2.y, &length, linestring_centroid_sum);
1537  }
1538  if (l_num_coords > 4 && closed) {
1539  // Also add the closing segment between the last and the first points
1540  centroid_add_segment(l2.x, l2.y, l0.x, l0.y, &length, linestring_centroid_sum);
1541  }
1542  *total_length += length;
1543  if (length == 0.0 && l_num_coords > 0) {
1544  *num_points += 1;
1545  point_centroid_sum[0] += l0.x;
1546  point_centroid_sum[1] += l0.y;
1547  }
1548  return true;
1549 }
1550
1552 void ST_Centroid_LineString(int8_t* coords,
1553  int32_t coords_sz,
1554  int32_t ic,
1555  int32_t isr,
1556  int32_t osr,
1557  double* linestring_centroid) {
1558  double length = 0.0;
1559  double linestring_centroid_sum[2] = {0.0, 0.0};
1560  int64_t num_points = 0;
1561  double point_centroid_sum[2] = {0.0, 0.0};
1563  coords_sz,
1564  ic,
1565  isr,
1566  osr,
1567  false, // not closed
1568  &length,
1569  &linestring_centroid_sum[0],
1570  &num_points,
1571  &point_centroid_sum[0]);
1572  if (length > 0) {
1573  linestring_centroid[0] = linestring_centroid_sum[0] / length;
1574  linestring_centroid[1] = linestring_centroid_sum[1] / length;
1575  } else if (num_points > 0) {
1576  linestring_centroid[0] = point_centroid_sum[0] / num_points;
1577  linestring_centroid[1] = point_centroid_sum[1] / num_points;
1578  }
1579 }
1580
1582  double y1,
1583  double x2,
1584  double y2,
1585  double x3,
1586  double y3,
1587  double sign,
1588  double* total_area2,
1589  double* cg3) {
1590  double cx = x1 + x2 + x3;
1591  double cy = y1 + y2 + y3;
1592  double area2 = x1 * y2 - x2 * y1 + x3 * y1 - x1 * y3 + x2 * y3 - x3 * y2;
1593  cg3[0] += sign * area2 * cx;
1594  cg3[1] += sign * area2 * cy;
1595  *total_area2 += sign * area2;
1596  return true;
1597 }
1598
1600  int64_t ringsize,
1601  int32_t ic,
1602  int32_t isr,
1603  int32_t osr,
1604  double sign,
1605  double* total_area2,
1606  double* cg3,
1607  double* total_length,
1608  double* linestring_centroid_sum,
1609  int64_t* num_points,
1610  double* point_centroid_sum) {
1611  auto ring_num_coords = ringsize / compression_unit_size(ic);
1612
1613  if (ring_num_coords < 6) {
1614  return false;
1615  }
1616
1617  Point2D p1 = get_point(ring, 0, ic, isr, osr);
1618  Point2D p2 = get_point(ring, 2, ic, isr, osr);
1619  for (int32_t i = 4; i < ring_num_coords; i += 2) {
1620  Point2D p3 = get_point(ring, i, ic, isr, osr);
1621  centroid_add_triangle(p1.x, p1.y, p2.x, p2.y, p3.x, p3.y, sign, total_area2, cg3);
1622  p2 = p3;
1623  }
1624
1626  ringsize,
1627  ic,
1628  isr,
1629  osr,
1630  true, // closed
1631  total_length,
1632  linestring_centroid_sum,
1633  num_points,
1634  point_centroid_sum);
1635  return true;
1636 }
1637
1638 DEVICE ALWAYS_INLINE bool centroid_add_polygon(int8_t* poly_coords,
1639  int64_t poly_coords_size,
1640  int32_t* poly_ring_sizes,
1641  int64_t poly_num_rings,
1642  int32_t ic,
1643  int32_t isr,
1644  int32_t osr,
1645  double* total_area2,
1646  double* cg3,
1647  double* total_length,
1648  double* linestring_centroid_sum,
1649  int64_t* num_points,
1650  double* point_centroid_sum) {
1651  if (poly_num_rings <= 0) {
1652  return false;
1653  }
1654
1655  auto ring_coords = poly_coords;
1656
1657  for (auto r = 0; r < poly_num_rings; r++) {
1658  auto ring_coords_size = poly_ring_sizes[r] * 2 * compression_unit_size(ic);
1659  // Shape - positive area, holes - negative areas
1660  double sign = (r == 0) ? 1.0 : -1.0;
1662  ring_coords_size,
1663  ic,
1664  isr,
1665  osr,
1666  sign,
1667  total_area2,
1668  cg3,
1669  total_length,
1670  linestring_centroid_sum,
1671  num_points,
1672  point_centroid_sum);
1673  // Advance to the next ring.
1674  ring_coords += ring_coords_size;
1675  }
1676  return true;
1677 }
1678
1680 void ST_Centroid_Polygon(int8_t* poly_coords,
1681  int32_t poly_coords_size,
1682  int32_t* poly_ring_sizes,
1683  int32_t poly_num_rings,
1684  int32_t ic,
1685  int32_t isr,
1686  int32_t osr,
1687  double* poly_centroid) {
1688  if (poly_num_rings <= 0) {
1689  poly_centroid[0] = 0.0;
1690  poly_centroid[1] = 0.0;
1691  }
1692  double total_area2 = 0.0;
1693  double cg3[2] = {0.0, 0.0};
1694  double total_length = 0.0;
1695  double linestring_centroid_sum[2] = {0.0, 0.0};
1696  int64_t num_points = 0;
1697  double point_centroid_sum[2] = {0.0, 0.0};
1699  poly_coords_size,
1700  poly_ring_sizes,
1701  poly_num_rings,
1702  ic,
1703  isr,
1704  osr,
1705  &total_area2,
1706  &cg3[0],
1707  &total_length,
1708  &linestring_centroid_sum[0],
1709  &num_points,
1710  &point_centroid_sum[0]);
1711
1712  if (total_area2 != 0.0) {
1713  poly_centroid[0] = cg3[0] / 3 / total_area2;
1714  poly_centroid[1] = cg3[1] / 3 / total_area2;
1715  } else if (total_length > 0.0) {
1716  // zero-area polygon, fall back to linear centroid
1717  poly_centroid[0] = linestring_centroid_sum[0] / total_length;
1718  poly_centroid[1] = linestring_centroid_sum[1] / total_length;
1719  } else if (num_points > 0) {
1720  // zero-area zero-length polygon, fall further back to point centroid
1721  poly_centroid[0] = point_centroid_sum[0] / num_points;
1722  poly_centroid[1] = point_centroid_sum[1] / num_points;
1723  } else {
1724  Point2D centroid = get_point(poly_coords, 0, ic, isr, osr);
1725  poly_centroid[0] = centroid.x;
1726  poly_centroid[1] = centroid.y;
1727  }
1728 }
1729
1731 void ST_Centroid_MultiPolygon(int8_t* mpoly_coords,
1732  int32_t mpoly_coords_size,
1733  int32_t* mpoly_ring_sizes,
1734  int32_t mpoly_num_rings,
1735  int32_t* mpoly_poly_sizes,
1736  int32_t mpoly_num_polys,
1737  int32_t ic,
1738  int32_t isr,
1739  int32_t osr,
1740  double* mpoly_centroid) {
1741  if (mpoly_num_rings <= 0 || mpoly_num_polys <= 0) {
1742  mpoly_centroid[0] = 0.0;
1743  mpoly_centroid[1] = 0.0;
1744  }
1745
1746  double total_area2 = 0.0;
1747  double cg3[2] = {0.0, 0.0};
1748  double total_length = 0.0;
1749  double linestring_centroid_sum[2] = {0.0, 0.0};
1750  int64_t num_points = 0;
1751  double point_centroid_sum[2] = {0.0, 0.0};
1752
1753  // Set specific poly pointers as we move through the coords/ringsizes/polyrings arrays.
1754  auto next_poly_coords = mpoly_coords;
1755  auto next_poly_ring_sizes = mpoly_ring_sizes;
1756
1757  for (auto poly = 0; poly < mpoly_num_polys; poly++) {
1758  auto poly_coords = next_poly_coords;
1759  auto poly_ring_sizes = next_poly_ring_sizes;
1760  auto poly_num_rings = mpoly_poly_sizes[poly];
1761  // Count number of coords in all of poly's rings, advance ring size pointer.
1762  int32_t poly_num_coords = 0;
1763  for (auto ring = 0; ring < poly_num_rings; ring++) {
1764  poly_num_coords += 2 * *next_poly_ring_sizes++;
1765  }
1766  auto poly_coords_size = poly_num_coords * compression_unit_size(ic);
1767  next_poly_coords += poly_coords_size;
1768
1770  poly_coords_size,
1771  poly_ring_sizes,
1772  poly_num_rings,
1773  ic,
1774  isr,
1775  osr,
1776  &total_area2,
1777  &cg3[0],
1778  &total_length,
1779  &linestring_centroid_sum[0],
1780  &num_points,
1781  &point_centroid_sum[0]);
1782  }
1783
1784  if (total_area2 != 0.0) {
1785  mpoly_centroid[0] = cg3[0] / 3 / total_area2;
1786  mpoly_centroid[1] = cg3[1] / 3 / total_area2;
1787  } else if (total_length > 0.0) {
1788  // zero-area multipolygon, fall back to linear centroid
1789  mpoly_centroid[0] = linestring_centroid_sum[0] / total_length;
1790  mpoly_centroid[1] = linestring_centroid_sum[1] / total_length;
1791  } else if (num_points > 0) {
1792  // zero-area zero-length multipolygon, fall further back to point centroid
1793  mpoly_centroid[0] = point_centroid_sum[0] / num_points;
1794  mpoly_centroid[1] = point_centroid_sum[1] / num_points;
1795  } else {
1796  Point2D centroid = get_point(mpoly_coords, 0, ic, isr, osr);
1797  mpoly_centroid[0] = centroid.x;
1798  mpoly_centroid[1] = centroid.y;
1799  }
1800 }
1801
1802 //
1803 // ST_Distance
1804 //
1805
1807 double ST_Distance_Point_Point(int8_t* p1,
1808  int64_t p1size,
1809  int8_t* p2,
1810  int64_t p2size,
1811  int32_t ic1,
1812  int32_t isr1,
1813  int32_t ic2,
1814  int32_t isr2,
1815  int32_t osr) {
1816  Point2D const pt1 = get_point(p1, 0, ic1, isr1, osr);
1817  Point2D const pt2 = get_point(p2, 0, ic2, isr2, osr);
1818  return distance_point_point(pt1.x, pt1.y, pt2.x, pt2.y);
1819 }
1820
1823  int64_t p1size,
1824  int8_t* p2,
1825  int64_t p2size,
1826  int32_t ic1,
1827  int32_t isr1,
1828  int32_t ic2,
1829  int32_t isr2,
1830  int32_t osr) {
1831  Point2D const pt1 = get_point(p1, 0, ic1, isr1, osr);
1832  Point2D const pt2 = get_point(p2, 0, ic2, isr2, osr);
1833  return distance_point_point_squared(pt1.x, pt1.y, pt2.x, pt2.y);
1834 }
1835
1838  int64_t p1size,
1839  int8_t* p2,
1840  int64_t p2size,
1841  int32_t ic1,
1842  int32_t isr1,
1843  int32_t ic2,
1844  int32_t isr2,
1845  int32_t osr) {
1846  Point2D const pt1 = get_point(p1, 0, ic1, 4326, 4326);
1847  Point2D const pt2 = get_point(p2, 0, ic2, 4326, 4326);
1848  return distance_in_meters(pt1.x, pt1.y, pt2.x, pt2.y);
1849 }
1850
1853  int64_t psize,
1854  int8_t* l,
1855  int64_t lsize,
1856  int32_t ic1,
1857  int32_t isr1,
1858  int32_t ic2,
1859  int32_t isr2,
1860  int32_t osr) {
1861  // Currently only statically indexed LineString is supported
1862  Point2D const pt = get_point(p, 0, ic1, 4326, 4326);
1863  const auto lpoints = lsize / (2 * compression_unit_size(ic2));
1864  Point2D const pl = get_point(l, 2 * (lpoints - 1), ic2, 4326, 4326);
1865  return distance_in_meters(pt.x, pt.y, pl.x, pl.y);
1866 }
1867
1870  int64_t lsize,
1871  int8_t* p,
1872  int64_t psize,
1873  int32_t ic1,
1874  int32_t isr1,
1875  int32_t ic2,
1876  int32_t isr2,
1877  int32_t osr) {
1878  // Currently only statically indexed LineString is supported
1880  p, psize, l, lsize, ic2, isr2, ic1, isr1, osr);
1881 }
1882
1885  int64_t psize,
1886  int8_t* mp,
1887  int64_t mpsize,
1888  int32_t ic1,
1889  int32_t isr1,
1890  int32_t ic2,
1891  int32_t isr2,
1892  int32_t osr,
1893  double threshold) {
1894  Point2D pt = get_point(p, 0, ic1, isr1, osr);
1895  auto mp_num_coords = mpsize / compression_unit_size(ic2);
1896  Point2D mpp = get_point(mp, 0, ic2, isr2, osr);
1897  double dist = distance_point_point(pt.x, pt.y, mpp.x, mpp.y);
1898  for (int32_t i = 2; i < mp_num_coords; i += 2) {
1899  mpp = get_point(mp, i, ic2, isr2, osr);
1900  double ldist = distance_point_point(pt.x, pt.y, mpp.x, mpp.y);
1901  if (dist > ldist) {
1902  dist = ldist;
1903  }
1904  if (dist <= threshold) {
1905  return dist;
1906  }
1907  }
1908  return dist;
1909 }
1910
1913  int64_t psize,
1914  int8_t* mp,
1915  int64_t mpsize,
1916  int32_t ic1,
1917  int32_t isr1,
1918  int32_t ic2,
1919  int32_t isr2,
1920  int32_t osr,
1921  double threshold) {
1922  Point2D pt = get_point(p, 0, ic1, isr1, osr);
1923  auto mp_num_coords = mpsize / compression_unit_size(ic2);
1924  Point2D mpp = get_point(mp, 0, ic2, isr2, osr);
1925  double dist = distance_point_point_squared(pt.x, pt.y, mpp.x, mpp.y);
1926  for (int32_t i = 2; i < mp_num_coords; i += 2) {
1927  mpp = get_point(mp, i, ic2, isr2, osr);
1928  double ldist = distance_point_point_squared(pt.x, pt.y, mpp.x, mpp.y);
1929  if (dist > ldist) {
1930  dist = ldist;
1931  }
1932  if (dist <= threshold) {
1933  return dist;
1934  }
1935  }
1936  return dist;
1937 }
1938
1940  int64_t psize,
1941  int8_t* l,
1942  int64_t lsize,
1943  int32_t ic1,
1944  int32_t isr1,
1945  int32_t ic2,
1946  int32_t isr2,
1947  int32_t osr,
1948  bool check_closed,
1949  double threshold) {
1950  Point2D pt = get_point(p, 0, ic1, isr1, osr);
1951
1952  auto l_num_coords = lsize / compression_unit_size(ic2);
1953
1954  Point2D l1 = get_point(l, 0, ic2, isr2, osr);
1955  Point2D l2 = get_point(l, 2, ic2, isr2, osr);
1956
1957  double dist = distance_point_line(pt.x, pt.y, l1.x, l1.y, l2.x, l2.y);
1958  for (int32_t i = 4; i < l_num_coords; i += 2) {
1959  l1 = l2; // advance one point
1960  l2 = get_point(l, i, ic2, isr2, osr);
1961  double ldist = distance_point_line(pt.x, pt.y, l1.x, l1.y, l2.x, l2.y);
1962  if (dist > ldist) {
1963  dist = ldist;
1964  }
1965  if (dist <= threshold) {
1966  return dist;
1967  }
1968  }
1969  if (l_num_coords > 4 && check_closed) {
1970  // Also check distance to the closing edge between the first and the last points
1971  l1 = get_point(l, 0, ic2, isr2, osr);
1972  double ldist = distance_point_line(pt.x, pt.y, l1.x, l1.y, l2.x, l2.y);
1973  if (dist > ldist) {
1974  dist = ldist;
1975  }
1976  }
1977  return dist;
1978 }
1979
1982  int64_t psize,
1983  int8_t* l,
1984  int64_t lsize,
1985  int32_t ic1,
1986  int32_t isr1,
1987  int32_t ic2,
1988  int32_t isr2,
1989  int32_t osr,
1990  double threshold) {
1992  p, psize, l, lsize, ic1, isr1, ic2, isr2, osr, true, threshold);
1993 }
1994
1997  int64_t psize,
1998  int8_t* l,
1999  int64_t lsize,
2000  int32_t ic1,
2001  int32_t isr1,
2002  int32_t ic2,
2003  int32_t isr2,
2004  int32_t osr,
2005  double threshold) {
2007  p, psize, l, lsize, ic1, isr1, ic2, isr2, osr, false, threshold);
2008 }
2009
2011  int64_t psize,
2012  int8_t* mls,
2013  int64_t mls_size,
2014  int32_t* mls_ls_sizes,
2015  int64_t mls_ls_num,
2016  int32_t ic1,
2017  int32_t isr1,
2018  int32_t ic2,
2019  int32_t isr2,
2020  int32_t osr,
2021  double threshold) {
2022  if (mls_ls_num <= 0) {
2023  return 0.0;
2024  }
2025  double min_distance = 0.0;
2026
2027  auto next_linestring_coords = mls;
2028
2029  for (auto l = 0; l < mls_ls_num; l++) {
2030  auto linestring_coords = next_linestring_coords;
2031  auto linestring_num_points = mls_ls_sizes[l];
2032  auto linestring_num_coords = 2 * linestring_num_points;
2033  auto linestring_coords_size = linestring_num_coords * compression_unit_size(ic2);
2034  next_linestring_coords += linestring_coords_size;
2035  double distance = distance_point_linestring(p,
2036  psize,
2037  linestring_coords,
2038  linestring_coords_size,
2039  ic1,
2040  isr1,
2041  ic2,
2042  isr2,
2043  osr,
2044  false,
2045  threshold);
2046  if (l == 0 || min_distance > distance) {
2047  min_distance = distance;
2048  if (tol_zero(min_distance)) {
2049  min_distance = 0.0;
2050  break;
2051  }
2052  if (min_distance <= threshold) {
2053  break;
2054  }
2055  }
2056  }
2057
2058  return min_distance;
2059 }
2060
2063  int64_t psize,
2064  int8_t* mls,
2065  int64_t mls_size,
2066  int32_t* mls_ls_sizes,
2067  int64_t mls_ls_num,
2068  int32_t ic1,
2069  int32_t isr1,
2070  int32_t ic2,
2071  int32_t isr2,
2072  int32_t osr,
2073  double threshold) {
2075  psize,
2076  mls,
2077  mls_size,
2078  mls_ls_sizes,
2079  mls_ls_num,
2080  ic1,
2081  isr1,
2082  ic2,
2083  isr2,
2084  osr,
2085  threshold);
2086 }
2087
2089  int64_t psize,
2090  int8_t* poly,
2091  int64_t polysize,
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 exterior_ring_num_coords = polysize / compression_unit_size(ic2);
2101  if (poly_num_rings > 0) {
2102  exterior_ring_num_coords = poly_ring_sizes[0] * 2;
2103  }
2104  auto exterior_ring_coords_size = exterior_ring_num_coords * compression_unit_size(ic2);
2105
2106  Point2D pt = get_point(p, 0, ic1, isr1, osr);
2108  poly, exterior_ring_num_coords, pt.x, pt.y, ic2, isr2, osr)) {
2109  // Outside the exterior ring
2111  p, psize, poly, exterior_ring_coords_size, ic1, isr1, ic2, isr2, osr, threshold);
2112  }
2113  // Inside exterior ring
2114  // Advance to first interior ring
2115  poly += exterior_ring_coords_size;
2116  // Check if one of the polygon's holes contains that point
2117  for (auto r = 1; r < poly_num_rings; r++) {
2118  auto interior_ring_num_coords = poly_ring_sizes[r] * 2;
2119  auto interior_ring_coords_size =
2120  interior_ring_num_coords * compression_unit_size(ic2);
2122  poly, interior_ring_num_coords, pt.x, pt.y, ic2, isr2, osr)) {
2123  // Inside an interior ring
2125  psize,
2126  poly,
2127  interior_ring_coords_size,
2128  ic1,
2129  isr1,
2130  ic2,
2131  isr2,
2132  osr,
2133  threshold);
2134  }
2135  poly += interior_ring_coords_size;
2136  }
2137  return 0.0;
2138 }
2139
2141  int64_t psize,
2142  int8_t* mpoly_coords,
2143  int64_t mpoly_coords_size,
2144  int32_t* mpoly_ring_sizes,
2145  int64_t mpoly_num_rings,
2146  int32_t* mpoly_poly_sizes,
2147  int64_t mpoly_num_polys,
2148  int32_t ic1,
2149  int32_t isr1,
2150  int32_t ic2,
2151  int32_t isr2,
2152  int32_t osr,
2153  double threshold) {
2154  if (mpoly_num_polys <= 0) {
2155  return 0.0;
2156  }
2157  double min_distance = 0.0;
2158
2159  // Set specific poly pointers as we move through the coords/ringsizes/polyrings arrays.
2160  auto next_poly_coords = mpoly_coords;
2161  auto next_poly_ring_sizes = mpoly_ring_sizes;
2162
2163  for (auto poly = 0; poly < mpoly_num_polys; poly++) {
2164  auto poly_coords = next_poly_coords;
2165  auto poly_ring_sizes = next_poly_ring_sizes;
2166  auto poly_num_rings = mpoly_poly_sizes[poly];
2167  // Count number of coords in all of poly's rings, advance ring size pointer.
2168  int32_t poly_num_coords = 0;
2169  for (auto ring = 0; ring < poly_num_rings; ring++) {
2170  poly_num_coords += 2 * *next_poly_ring_sizes++;
2171  }
2172  auto poly_coords_size = poly_num_coords * compression_unit_size(ic2);
2173  next_poly_coords += poly_coords_size;
2174  double distance = distance_point_polygon(p,
2175  psize,
2176  poly_coords,
2177  poly_coords_size,
2178  poly_ring_sizes,
2179  poly_num_rings,
2180  ic1,
2181  isr1,
2182  ic2,
2183  isr2,
2184  osr,
2185  threshold);
2186  if (poly == 0 || min_distance > distance) {
2187  min_distance = distance;
2188  if (tol_zero(min_distance)) {
2189  min_distance = 0.0;
2190  break;
2191  }
2192  if (min_distance <= threshold) {
2193  break;
2194  }
2195  }
2196  }
2197
2198  return min_distance;
2199 }
2200
2202 double ST_Distance_Point_Polygon(int8_t* p,
2203  int64_t psize,
2204  int8_t* poly,
2205  int64_t polysize,
2206  int32_t* poly_ring_sizes,
2207  int64_t poly_num_rings,
2208  int32_t ic1,
2209  int32_t isr1,
2210  int32_t ic2,
2211  int32_t isr2,
2212  int32_t osr,
2213  double threshold) {
2214  return distance_point_polygon(p,
2215  psize,
2216  poly,
2217  polysize,
2218  poly_ring_sizes,
2219  poly_num_rings,
2220  ic1,
2221  isr1,
2222  ic2,
2223  isr2,
2224  osr,
2225  threshold);
2226 }
2227
2230  int64_t psize,
2231  int8_t* mpoly_coords,
2232  int64_t mpoly_coords_size,
2233  int32_t* mpoly_ring_sizes,
2234  int64_t mpoly_num_rings,
2235  int32_t* mpoly_poly_sizes,
2236  int64_t mpoly_num_polys,
2237  int32_t ic1,
2238  int32_t isr1,
2239  int32_t ic2,
2240  int32_t isr2,
2241  int32_t osr,
2242  double threshold) {
2243  return distance_point_multipolygon(p,
2244  psize,
2245  mpoly_coords,
2246  mpoly_coords_size,
2247  mpoly_ring_sizes,
2248  mpoly_num_rings,
2249  mpoly_poly_sizes,
2250  mpoly_num_polys,
2251  ic1,
2252  isr1,
2253  ic2,
2254  isr2,
2255  osr,
2256  threshold);
2257 }
2258
2261  int64_t mpsize,
2262  int8_t* poly,
2263  int64_t polysize,
2264  int32_t* poly_ring_sizes,
2265  int64_t poly_num_rings,
2266  int32_t ic1,
2267  int32_t isr1,
2268  int32_t ic2,
2269  int32_t isr2,
2270  int32_t osr,
2271  double threshold) {
2272  auto mp_num_coords = mpsize / compression_unit_size(ic1);
2273  auto* p = mp;
2274  auto const psize = compression_unit_size(ic1) * 2;
2275  auto dist = distance_point_polygon(p,
2276  psize,
2277  poly,
2278  polysize,
2279  poly_ring_sizes,
2280  poly_num_rings,
2281  ic1,
2282  isr1,
2283  ic2,
2284  isr2,
2285  osr,
2286  threshold);
2287  p += psize;
2288  for (int32_t i = 2; i < mp_num_coords; i += 2) {
2289  auto ldist = distance_point_polygon(p,
2290  psize,
2291  poly,
2292  polysize,
2293  poly_ring_sizes,
2294  poly_num_rings,
2295  ic1,
2296  isr1,
2297  ic2,
2298  isr2,
2299  osr,
2300  threshold);
2301  p += psize;
2302  if (dist > ldist) {
2303  dist = ldist;
2304  }
2305  if (dist <= threshold) {
2306  return dist;
2307  }
2308  }
2309  return dist;
2310 }
2311
2314  int64_t mpsize,
2315  int8_t* mpoly_coords,
2316  int64_t mpoly_coords_size,
2317  int32_t* mpoly_ring_sizes,
2318  int64_t mpoly_num_rings,
2319  int32_t* mpoly_poly_sizes,
2320  int64_t mpoly_num_polys,
2321  int32_t ic1,
2322  int32_t isr1,
2323  int32_t ic2,
2324  int32_t isr2,
2325  int32_t osr,
2326  double threshold) {
2327  auto mp_num_coords = mpsize / compression_unit_size(ic1);
2328  auto* p = mp;
2329  auto const psize = compression_unit_size(ic1) * 2;
2330  auto dist = distance_point_multipolygon(p,
2331  psize,
2332  mpoly_coords,
2333  mpoly_coords_size,
2334  mpoly_ring_sizes,
2335  mpoly_num_rings,
2336  mpoly_poly_sizes,
2337  mpoly_num_polys,
2338  ic1,
2339  isr1,
2340  ic2,
2341  isr2,
2342  osr,
2343  threshold);
2344  p += psize;
2345  for (int32_t i = 2; i < mp_num_coords; i += 2) {
2346  auto ldist = distance_point_multipolygon(p,
2347  psize,
2348  mpoly_coords,
2349  mpoly_coords_size,
2350  mpoly_ring_sizes,
2351  mpoly_num_rings,
2352  mpoly_poly_sizes,
2353  mpoly_num_polys,
2354  ic1,
2355  isr1,
2356  ic2,
2357  isr2,
2358  osr,
2359  threshold);
2360  p += psize;
2361  if (dist > ldist) {
2362  dist = ldist;
2363  }
2364  if (dist <= threshold) {
2365  return dist;
2366  }
2367  }
2368  return dist;
2369 }
2370
2373  int64_t polysize,
2374  int32_t* poly_ring_sizes,
2375  int64_t poly_num_rings,
2376  int8_t* mp,
2377  int64_t mpsize,
2378  int32_t ic1,
2379  int32_t isr1,
2380  int32_t ic2,
2381  int32_t isr2,
2382  int32_t osr,
2383  double threshold) {
2385  mpsize,
2386  poly,
2387  polysize,
2388  poly_ring_sizes,
2389  poly_num_rings,
2390  ic2,
2391  isr2,
2392  ic1,
2393  isr1,
2394  osr,
2395  threshold);
2396 }
2397
2399 double ST_Distance_MultiPolygon_MultiPoint(int8_t* mpoly_coords,
2400  int64_t mpoly_coords_size,
2401  int32_t* mpoly_ring_sizes,
2402  int64_t mpoly_num_rings,
2403  int32_t* mpoly_poly_sizes,
2404  int64_t mpoly_num_polys,
2405  int8_t* mp,
2406  int64_t mpsize,
2407  int32_t ic1,
2408  int32_t isr1,
2409  int32_t ic2,
2410  int32_t isr2,
2411  int32_t osr,
2412  double threshold) {
2414  mpsize,
2415  mpoly_coords,
2416  mpoly_coords_size,
2417  mpoly_ring_sizes,
2418  mpoly_num_rings,
2419  mpoly_poly_sizes,
2420  mpoly_num_polys,
2421  ic2,
2422  isr2,
2423  ic1,
2424  isr1,
2425  osr,
2426  threshold);
2427 }
2428
2431  int64_t mpsize,
2432  int8_t* p,
2433  int64_t psize,
2434  int32_t ic1,
2435  int32_t isr1,
2436  int32_t ic2,
2437  int32_t isr2,
2438  int32_t osr,
2439  double threshold) {
2441  p, psize, mp, mpsize, ic2, isr2, ic1, isr1, osr, threshold);
2442 }
2443
2446  int64_t mpsize,
2447  int8_t* p,
2448  int64_t psize,
2449  int32_t ic1,
2450  int32_t isr1,
2451  int32_t ic2,
2452  int32_t isr2,
2453  int32_t osr,
2454  double threshold) {
2456  p, psize, mp, mpsize, ic2, isr2, ic1, isr1, osr, threshold);
2457 }
2458
2461  int64_t mp1size,
2462  int8_t* mp2,
2463  int64_t mp2size,
2464  int32_t ic1,
2465  int32_t isr1,
2466  int32_t ic2,
2467  int32_t isr2,
2468  int32_t osr,
2469  double threshold) {
2470  auto mp1_num_coords = mp1size / compression_unit_size(ic1);
2471  auto mp2_num_coords = mp2size / compression_unit_size(ic2);
2472  double dist = -1.0;
2473  for (int32_t i = 0; i < mp1_num_coords; i += 2) {
2474  Point2D mp1p = get_point(mp1, 0, ic1, isr1, osr);
2475  for (int32_t j = 0; j < mp2_num_coords; j += 2) {
2476  Point2D mp2p = get_point(mp2, 0, ic2, isr2, osr);
2477  double ldist = distance_point_point(mp1p.x, mp1p.y, mp2p.x, mp2p.y);
2478  if (dist > ldist || dist < 0.0) {
2479  dist = ldist;
2480  }
2481  if (dist <= threshold) {
2482  return dist;
2483  }
2484  }
2485  }
2486  return dist;
2487 }
2488
2491  int64_t mp1size,
2492  int8_t* mp2,
2493  int64_t mp2size,
2494  int32_t ic1,
2495  int32_t isr1,
2496  int32_t ic2,
2497  int32_t isr2,
2498  int32_t osr,
2499  double threshold) {
2500  auto mp1_num_coords = mp1size / compression_unit_size(ic1);
2501  auto mp2_num_coords = mp2size / compression_unit_size(ic2);
2502  double dist = -1.0;
2503  for (int32_t i = 0; i < mp1_num_coords; i += 2) {
2504  Point2D mp1p = get_point(mp1, i, ic1, isr1, osr);
2505  for (int32_t j = 0; j < mp2_num_coords; j += 2) {
2506  Point2D mp2p = get_point(mp2, j, ic2, isr2, osr);
2507  double ldist = distance_point_point_squared(mp1p.x, mp1p.y, mp2p.x, mp2p.y);
2508  if (dist > ldist || dist < 0.0) {
2509  dist = ldist;
2510  }
2511  if (dist <= threshold) {
2512  return dist;
2513  }
2514  }
2515  }
2516  return dist;
2517 }
2518
2521  int64_t mpsize,
2522  int8_t* l,
2523  int64_t lsize,
2524  int32_t ic1,
2525  int32_t isr1,
2526  int32_t ic2,
2527  int32_t isr2,
2528  int32_t osr,
2529  double threshold) {
2530  auto mp_num_coords = mpsize / compression_unit_size(ic1);
2531  auto* p = mp;
2532  auto const psize = compression_unit_size(ic1) * 2;
2533  auto dist = distance_point_linestring(
2534  p, psize, l, lsize, ic1, isr1, ic2, isr2, osr, true, threshold);
2535  p += psize;
2536  for (int32_t i = 2; i < mp_num_coords; i += 2) {
2537  auto ldist = distance_point_linestring(
2538  p, psize, l, lsize, ic1, isr1, ic2, isr2, osr, true, threshold);
2539  p += psize;
2540  if (dist > ldist) {
2541  dist = ldist;
2542  }
2543  if (dist <= threshold) {
2544  return dist;
2545  }
2546  }
2547  return dist;
2548 }
2549
2552  int64_t mpsize,
2553  int8_t* mls,
2554  int64_t mls_size,
2555  int32_t* mls_ls_sizes,
2556  int64_t mls_ls_num,
2557  int32_t ic1,
2558  int32_t isr1,
2559  int32_t ic2,
2560  int32_t isr2,
2561  int32_t osr,
2562  double threshold) {
2563  auto mp_num_coords = mpsize / compression_unit_size(ic1);
2564  auto* p = mp;
2565  auto const psize = compression_unit_size(ic1) * 2;
2566  auto dist = distance_point_multilinestring(p,
2567  psize,
2568  mls,
2569  mls_size,
2570  mls_ls_sizes,
2571  mls_ls_num,
2572  ic1,
2573  isr1,
2574  ic2,
2575  isr2,
2576  osr,
2577  threshold);
2578  p += psize;
2579  for (int32_t i = 2; i < mp_num_coords; i += 2) {
2580  auto ldist = distance_point_multilinestring(p,
2581  psize,
2582  mls,
2583  mls_size,
2584  mls_ls_sizes,
2585  mls_ls_num,
2586  ic1,
2587  isr1,
2588  ic2,
2589  isr2,
2590  osr,
2591  threshold);
2592  p += psize;
2593  if (dist > ldist) {
2594  dist = ldist;
2595  }
2596  if (dist <= threshold) {
2597  return dist;
2598  }
2599  }
2600  return dist;
2601 }
2602
2605  int64_t lsize,
2606  int8_t* p,
2607  int64_t psize,
2608  int32_t ic1,
2609  int32_t isr1,
2610  int32_t ic2,
2611  int32_t isr2,
2612  int32_t osr,
2613  double threshold) {
2615  p, psize, l, lsize, ic2, isr2, ic1, isr1, osr, threshold);
2616 }
2617
2620  int64_t lsize,
2621  int8_t* mp,
2622  int64_t mpsize,
2623  int32_t ic1,
2624  int32_t isr1,
2625  int32_t ic2,
2626  int32_t isr2,
2627  int32_t osr,
2628  double threshold) {
2630  mp, mpsize, l, lsize, ic2, isr2, ic1, isr1, osr, threshold);
2631 }
2632
2635  int64_t mls_size,
2636  int32_t* mls_ls_sizes,
2637  int64_t mls_ls_num,
2638  int8_t* mp,
2639  int64_t mpsize,
2640  int32_t ic1,
2641  int32_t isr1,
2642  int32_t ic2,
2643  int32_t isr2,
2644  int32_t osr,
2645  double threshold) {
2647  mpsize,
2648  mls,
2649  mls_size,
2650  mls_ls_sizes,
2651  mls_ls_num,
2652  ic2,
2653  isr2,
2654  ic1,
2655  isr1,
2656  osr,
2657  threshold);
2658 }
2659
2662  int64_t l1size,
2663  int8_t* l2,
2664  int64_t l2size,
2665  int32_t ic1,
2666  int32_t isr1,
2667  int32_t ic2,
2668  int32_t isr2,
2669  int32_t osr,
2670  double threshold) {
2671  auto l1_num_coords = l1size / compression_unit_size(ic1);
2672  auto l2_num_coords = l2size / compression_unit_size(ic2);
2673
2674  double threshold_squared = threshold * threshold;
2675  double dist_squared = 0.0;
2676  Point2D l11 = get_point(l1, 0, ic1, isr1, osr);
2677  for (int32_t i1 = 2; i1 < l1_num_coords; i1 += 2) {
2678  Point2D l12 = get_point(l1, i1, ic1, isr1, osr);
2679  Point2D l21 = get_point(l2, 0, ic2, isr2, osr);
2680  for (int32_t i2 = 2; i2 < l2_num_coords; i2 += 2) {
2681  Point2D l22 = get_point(l2, i2, ic2, isr2, osr);
2682
2683  // double ldist_squared =
2684  // distance_line_line_squared(l11x, l11y, l12x, l12y, l21x, l21y, l22x, l22y);
2685  // TODO: fix distance_line_line_squared
2686  double ldist =
2687  distance_line_line(l11.x, l11.y, l12.x, l12.y, l21.x, l21.y, l22.x, l22.y);
2688  double ldist_squared = ldist * ldist;
2689
2690  if (i1 == 2 && i2 == 2) {
2691  dist_squared = ldist_squared; // initialize dist with distance between the first
2692  // two segments
2693  } else if (dist_squared > ldist_squared) {
2694  dist_squared = ldist_squared;
2695  }
2696  if (tol_zero(dist_squared, TOLERANCE_DEFAULT_SQUARED)) {
2697  return 0.0; // Segments touch, short-circuit
2698  }
2699  if (dist_squared <= threshold_squared) {
2700  // If threashold_squared is defined and the calculated dist_squared dips under it,
2701  // short-circuit and return it
2702  return sqrt(dist_squared);
2703  }
2704
2705  l21 = l22; // advance to the next point on l2
2706  }
2707
2708  l11 = l12; // advance to the next point on l1
2709  }
2710  return sqrt(dist_squared);
2711 }
2712
2715  int64_t lsize,
2716  int8_t* poly_coords,
2717  int64_t poly_coords_size,
2718  int32_t* poly_ring_sizes,
2719  int64_t poly_num_rings,
2720  int32_t ic1,
2721  int32_t isr1,
2722  int32_t ic2,
2723  int32_t isr2,
2724  int32_t osr,
2725  double threshold) {
2726  auto lnum_coords = lsize / compression_unit_size(ic1);
2727  auto psize = 2 * compression_unit_size(ic1);
2728  auto min_distance =
2729  ST_Distance_Point_Polygon(l, // pointer to start of linestring for first point
2730  psize,
2731  poly_coords,
2732  poly_coords_size,
2733  poly_ring_sizes,
2734  poly_num_rings,
2735  ic1,
2736  isr1,
2737  ic2,
2738  isr2,
2739  osr,
2740  threshold);
2741  if (tol_zero(min_distance)) {
2742  // Linestring's first point is inside the poly
2743  return 0.0;
2744  }
2745  if (min_distance <= threshold) {
2746  return min_distance;
2747  }
2748
2749  // Otherwise, linestring's first point is outside the external ring or inside
2750  // an internal ring. Measure minimum distance between linestring segments and
2751  // poly rings. Crossing a ring zeroes the distance and causes an early return.
2752  auto poly_ring_coords = poly_coords;
2753  for (auto r = 0; r < poly_num_rings; r++) {
2754  int64_t poly_ring_num_coords = poly_ring_sizes[r] * 2;
2755
2756  auto distance = distance_ring_linestring(poly_ring_coords,
2757  poly_ring_num_coords,
2758  l,
2759  lnum_coords,
2760  ic2,
2761  isr2,
2762  ic1,
2763  isr1,
2764  osr,
2765  threshold);
2766  if (min_distance > distance) {
2767  min_distance = distance;
2768  if (tol_zero(min_distance)) {
2769  return 0.0;
2770  }
2771  if (min_distance <= threshold) {
2772  return min_distance;
2773  }
2774  }
2775
2776  poly_ring_coords += poly_ring_num_coords * compression_unit_size(ic2);
2777  }
2778
2779  return min_distance;
2780 }
2781
2784  int64_t lsize,
2785  int8_t* mpoly_coords,
2786  int64_t mpoly_coords_size,
2787  int32_t* mpoly_ring_sizes,
2788  int64_t mpoly_num_rings,
2789  int32_t* mpoly_poly_sizes,
2790  int64_t mpoly_num_polys,
2791  int32_t ic1,
2792  int32_t isr1,
2793  int32_t ic2,
2794  int32_t isr2,
2795  int32_t osr,
2796  double threshold) {
2797  // TODO: revisit implementation, cover all cases
2798  double min_distance = 0.0;
2799
2800  // Set specific poly pointers as we move through the coords/ringsizes/polyrings arrays.
2801  auto next_poly_coords = mpoly_coords;
2802  auto next_poly_ring_sizes = mpoly_ring_sizes;
2803
2804  for (auto poly = 0; poly < mpoly_num_polys; poly++) {
2805  auto poly_coords = next_poly_coords;
2806  auto poly_ring_sizes = next_poly_ring_sizes;
2807  auto poly_num_rings = mpoly_poly_sizes[poly];
2808  // Count number of coords in all of poly's rings, advance ring size pointer.
2809  int32_t poly_num_coords = 0;
2810  for (auto ring = 0; ring < poly_num_rings; ring++) {
2811  poly_num_coords += 2 * *next_poly_ring_sizes++;
2812  }
2813  auto poly_coords_size = poly_num_coords * compression_unit_size(ic2);
2814  next_poly_coords += poly_coords_size;
2815  double distance = ST_Distance_LineString_Polygon(l,
2816  lsize,
2817  poly_coords,
2818  poly_coords_size,
2819  poly_ring_sizes,
2820  poly_num_rings,
2821  ic1,
2822  isr1,
2823  ic2,
2824  isr2,
2825  osr,
2826  threshold);
2827  if (poly == 0 || min_distance > distance) {
2828  min_distance = distance;
2829  if (tol_zero(min_distance)) {
2830  min_distance = 0.0;
2831  break;
2832  }
2833  if (min_distance <= threshold) {
2834  return min_distance;
2835  }
2836  }
2837  }
2838
2839  return min_distance;
2840 }
2841
2844  int64_t mls_size,
2845  int32_t* mls_ls_sizes,
2846  int64_t mls_ls_num,
2847  int8_t* p,
2848  int64_t psize,
2849  int32_t ic1,
2850  int32_t isr1,
2851  int32_t ic2,
2852  int32_t isr2,
2853  int32_t osr,
2854  double threshold) {
2856  psize,
2857  mls,
2858  mls_size,
2859  mls_ls_sizes,
2860  mls_ls_num,
2861  ic2,
2862  isr2,
2863  ic1,
2864  isr1,
2865  osr,
2866  threshold);
2867 }
2868
2870 double ST_Distance_Polygon_Point(int8_t* poly_coords,
2871  int64_t poly_coords_size,
2872  int32_t* poly_ring_sizes,
2873  int64_t poly_num_rings,
2874  int8_t* p,
2875  int64_t psize,
2876  int32_t ic1,
2877  int32_t isr1,
2878  int32_t ic2,
2879  int32_t isr2,
2880  int32_t osr,
2881  double threshold) {
2882  return ST_Distance_Point_Polygon(p,
2883  psize,
2884  poly_coords,
2885  poly_coords_size,
2886  poly_ring_sizes,
2887  poly_num_rings,
2888  ic2,
2889  isr2,
2890  ic1,
2891  isr1,
2892  osr,
2893  threshold);
2894 }
2895
2897 double ST_Distance_Polygon_LineString(int8_t* poly_coords,
2898  int64_t poly_coords_size,
2899  int32_t* poly_ring_sizes,
2900  int64_t poly_num_rings,
2901  int8_t* l,
2902  int64_t lsize,
2903  int32_t ic1,
2904  int32_t isr1,
2905  int32_t ic2,
2906  int32_t isr2,
2907  int32_t osr,
2908  double threshold) {
2910  lsize,
2911  poly_coords,
2912  poly_coords_size,
2913  poly_ring_sizes,
2914  poly_num_rings,
2915  ic2,
2916  isr2,
2917  ic1,
2918  isr2,
2919  osr,
2920  threshold);
2921 }
2922
2924 double ST_Distance_Polygon_Polygon(int8_t* poly1_coords,
2925  int64_t poly1_coords_size,
2926  int32_t* poly1_ring_sizes,
2927  int64_t poly1_num_rings,
2928  int8_t* poly2_coords,
2929  int64_t poly2_coords_size,
2930  int32_t* poly2_ring_sizes,
2931  int64_t poly2_num_rings,
2932  int32_t ic1,
2933  int32_t isr1,
2934  int32_t ic2,
2935  int32_t isr2,
2936  int32_t osr,
2937  double threshold) {
2938  // Check if poly1 contains the first point of poly2's shape, i.e. the external ring
2939  auto poly2_first_point_coords = poly2_coords;
2940  auto poly2_first_point_coords_size = compression_unit_size(ic2) * 2;
2941  auto min_distance = ST_Distance_Polygon_Point(poly1_coords,
2942  poly1_coords_size,
2943  poly1_ring_sizes,
2944  poly1_num_rings,
2945  poly2_first_point_coords,
2946  poly2_first_point_coords_size,
2947  ic1,
2948  isr1,
2949  ic2,
2950  isr2,
2951  osr,
2952  threshold);
2953  if (tol_zero(min_distance)) {
2954  // Polygons overlap
2955  return 0.0;
2956  }
2957  if (min_distance <= threshold) {
2958  return min_distance;
2959  }
2960
2961  // Poly2's first point is either outside poly1's external ring or inside one of the
2962  // internal rings. Measure the smallest distance between a poly1 ring (external or
2963  // internal) and a poly2 ring (external or internal). If poly2 is completely outside
2964  // poly1, then the min distance would be between poly1's and poly2's external rings. If
2965  // poly2 is completely inside one of poly1 internal rings then the min distance would be
2966  // between that poly1 internal ring and poly2's external ring. If poly1 is completely
2967  // inside one of poly2 internal rings, min distance is between that internal ring and
2968  // poly1's external ring. In each case other rings don't get in the way. Any ring
2969  // intersection means zero distance - short-circuit and return.
2970
2971  auto poly1_ring_coords = poly1_coords;
2972  for (auto r1 = 0; r1 < poly1_num_rings; r1++) {
2973  int64_t poly1_ring_num_coords = poly1_ring_sizes[r1] * 2;
2974
2975  auto poly2_ring_coords = poly2_coords;
2976  for (auto r2 = 0; r2 < poly2_num_rings; r2++) {
2977  int64_t poly2_ring_num_coords = poly2_ring_sizes[r2] * 2;
2978
2979  auto distance = distance_ring_ring(poly1_ring_coords,
2980  poly1_ring_num_coords,
2981  poly2_ring_coords,
2982  poly2_ring_num_coords,
2983  ic1,
2984  isr1,
2985  ic2,
2986  isr2,
2987  osr,
2988  threshold);
2989  if (min_distance > distance) {
2990  min_distance = distance;
2991  if (tol_zero(min_distance)) {
2992  return 0.0;
2993  }
2994  if (min_distance <= threshold) {
2995  return min_distance;
2996  }
2997  }
2998
2999  poly2_ring_coords += poly2_ring_num_coords * compression_unit_size(ic2);
3000  }
3001
3002  poly1_ring_coords += poly1_ring_num_coords * compression_unit_size(ic1);
3003  }
3004
3005  return min_distance;
3006 }
3007
3009 double ST_Distance_Polygon_MultiPolygon(int8_t* poly1_coords,
3010  int64_t poly1_coords_size,
3011  int32_t* poly1_ring_sizes,
3012  int64_t poly1_num_rings,
3013  int8_t* mpoly_coords,
3014  int64_t mpoly_coords_size,
3015  int32_t* mpoly_ring_sizes,
3016  int64_t mpoly_num_rings,
3017  int32_t* mpoly_poly_sizes,
3018  int64_t mpoly_num_polys,
3019  int32_t ic1,
3020  int32_t isr1,
3021  int32_t ic2,
3022  int32_t isr2,
3023  int32_t osr,
3024  double threshold) {
3025  double min_distance = 0.0;
3026
3027  // Set specific poly pointers as we move through the coords/ringsizes/polyrings arrays.
3028  auto next_poly_coords = mpoly_coords;
3029  auto next_poly_ring_sizes = mpoly_ring_sizes;
3030
3031  for (auto poly = 0; poly < mpoly_num_polys; poly++) {
3032  auto poly_coords = next_poly_coords;
3033  auto poly_ring_sizes = next_poly_ring_sizes;
3034  auto poly_num_rings = mpoly_poly_sizes[poly];
3035  // Count number of coords in all of poly's rings, advance ring size pointer.
3036  int32_t poly_num_coords = 0;
3037  for (auto ring = 0; ring < poly_num_rings; ring++) {
3038  poly_num_coords += 2 * *next_poly_ring_sizes++;
3039  }
3040  auto poly_coords_size = poly_num_coords * compression_unit_size(ic2);
3041  next_poly_coords += poly_coords_size;
3042  double distance = ST_Distance_Polygon_Polygon(poly1_coords,
3043  poly1_coords_size,
3044  poly1_ring_sizes,
3045  poly1_num_rings,
3046  poly_coords,
3047  poly_coords_size,
3048  poly_ring_sizes,
3049  poly_num_rings,
3050  ic1,
3051  isr1,
3052  ic2,
3053  isr2,
3054  osr,
3055  threshold);
3056  if (poly == 0 || min_distance > distance) {
3057  min_distance = distance;
3058  if (tol_zero(min_distance)) {
3059  min_distance = 0.0;
3060  break;
3061  }
3062  if (min_distance <= threshold) {
3063  break;
3064  }
3065  }
3066  }
3067
3068  return min_distance;
3069 }
3070
3072 double ST_Distance_MultiPolygon_Point(int8_t* mpoly_coords,
3073  int64_t mpoly_coords_size,
3074  int32_t* mpoly_ring_sizes,
3075  int64_t mpoly_num_rings,
3076  int32_t* mpoly_poly_sizes,
3077  int64_t mpoly_num_polys,
3078  int8_t* p,
3079  int64_t psize,
3080  int32_t ic1,
3081  int32_t isr1,
3082  int32_t ic2,
3083  int32_t isr2,
3084  int32_t osr,
3085  double threshold) {
3087  psize,
3088  mpoly_coords,
3089  mpoly_coords_size,
3090  mpoly_ring_sizes,
3091  mpoly_num_rings,
3092  mpoly_poly_sizes,
3093  mpoly_num_polys,
3094  ic2,
3095  isr2,
3096  ic1,
3097  isr1,
3098  osr,
3099  threshold);
3100 }
3101
3103 double ST_Distance_MultiPolygon_LineString(int8_t* mpoly_coords,
3104  int64_t mpoly_coords_size,
3105  int32_t* mpoly_ring_sizes,
3106  int64_t mpoly_num_rings,
3107  int32_t* mpoly_poly_sizes,
3108  int64_t mpoly_num_polys,
3109  int8_t* l,
3110  int64_t lsize,
3111  int32_t ic1,
3112  int32_t isr1,
3113  int32_t ic2,
3114  int32_t isr2,
3115  int32_t osr,
3116  double threshold) {
3118  lsize,
3119  mpoly_coords,
3120  mpoly_coords_size,
3121  mpoly_ring_sizes,
3122  mpoly_num_rings,
3123  mpoly_poly_sizes,
3124  mpoly_num_polys,
3125  ic2,
3126  isr2,
3127  ic1,
3128  isr1,
3129  osr,
3130  threshold);
3131 }
3132
3134 double ST_Distance_MultiPolygon_Polygon(int8_t* mpoly_coords,
3135  int64_t mpoly_coords_size,
3136  int32_t* mpoly_ring_sizes,
3137  int64_t mpoly_num_rings,
3138  int32_t* mpoly_poly_sizes,
3139  int64_t mpoly_num_polys,
3140  int8_t* poly1_coords,
3141  int64_t poly1_coords_size,
3142  int32_t* poly1_ring_sizes,
3143  int64_t poly1_num_rings,
3144  int32_t ic1,
3145  int32_t isr1,
3146  int32_t ic2,
3147  int32_t isr2,
3148  int32_t osr,
3149  double threshold) {
3150  return ST_Distance_Polygon_MultiPolygon(poly1_coords,
3151  poly1_coords_size,
3152  poly1_ring_sizes,
3153  poly1_num_rings,
3154  mpoly_coords,
3155  mpoly_coords_size,
3156  mpoly_ring_sizes,
3157  mpoly_num_rings,
3158  mpoly_poly_sizes,
3159  mpoly_num_polys,
3160  ic2,
3161  isr2,
3162  ic1,
3163  isr1,
3164  osr,
3165  threshold);
3166 }
3167
3169 double ST_Distance_MultiPolygon_MultiPolygon(int8_t* mpoly1_coords,
3170  int64_t mpoly1_coords_size,
3171  int32_t* mpoly1_ring_sizes,
3172  int64_t mpoly1_num_rings,
3173  int32_t* mpoly1_poly_sizes,
3174  int64_t mpoly1_num_polys,
3175  int8_t* mpoly2_coords,
3176  int64_t mpoly2_coords_size,
3177  int32_t* mpoly2_ring_sizes,
3178  int64_t mpoly2_num_rings,
3179  int32_t* mpoly2_poly_sizes,
3180  int64_t mpoly2_num_polys,
3181  int32_t ic1,
3182  int32_t isr1,
3183  int32_t ic2,
3184  int32_t isr2,
3185  int32_t osr,
3186  double threshold) {
3187  double min_distance = 0.0;
3188
3189  // Set specific poly pointers as we move through mpoly1's coords/ringsizes/polyrings
3190  // arrays.
3191  auto next_poly_coords = mpoly1_coords;
3192  auto next_poly_ring_sizes = mpoly1_ring_sizes;
3193
3194  for (auto poly = 0; poly < mpoly1_num_polys; poly++) {
3195  auto poly_coords = next_poly_coords;
3196  auto poly_ring_sizes = next_poly_ring_sizes;
3197  auto poly_num_rings = mpoly1_poly_sizes[poly];
3198  // Count number of coords in all of poly's rings, advance ring size pointer.
3199  int32_t poly_num_coords = 0;
3200  for (auto ring = 0; ring < poly_num_rings; ring++) {
3201  poly_num_coords += 2 * *next_poly_ring_sizes++;
3202  }
3203  auto poly_coords_size = poly_num_coords * compression_unit_size(ic1);
3204  next_poly_coords += poly_coords_size;
3205  double distance = ST_Distance_Polygon_MultiPolygon(poly_coords,
3206  poly_coords_size,
3207  poly_ring_sizes,
3208  poly_num_rings,
3209  mpoly2_coords,
3210  mpoly2_coords_size,
3211  mpoly2_ring_sizes,
3212  mpoly2_num_rings,
3213  mpoly2_poly_sizes,
3214  mpoly2_num_polys,
3215  ic1,
3216  isr1,
3217  ic2,
3218  isr2,
3219  osr,
3220  threshold);
3221  if (poly == 0 || min_distance > distance) {
3222  min_distance = distance;
3223  if (tol_zero(min_distance)) {
3224  min_distance = 0.0;
3225  break;
3226  }
3227  if (min_distance <= threshold) {
3228  break;
3229  }
3230  }
3231  }
3232
3233  return min_distance;
3234 }
3235
3238  int64_t mls1size,
3239  int32_t* mls1_ls_sizes,
3240  int64_t mls1_ls_num,
3241  int8_t* mls2,
3242  int64_t mls2size,
3243  int32_t* mls2_ls_sizes,
3244  int64_t mls2_ls_num,
3245  int32_t ic1,
3246  int32_t isr1,
3247  int32_t ic2,
3248  int32_t isr2,
3249  int32_t osr,
3250  double threshold) {
3251  return 0.0;
3252 }
3253
3254 //
3255 // ST_DWithin
3256 //
3257
3259 bool ST_DWithin_Point_Point(int8_t* p1,
3260  int64_t p1size,
3261  int8_t* p2,
3262  int64_t p2size,
3263  int32_t ic1,
3264  int32_t isr1,
3265  int32_t ic2,
3266  int32_t isr2,
3267  int32_t osr,
3268  double distance_within) {
3270  p1, p1size, p2, p2size, ic1, isr1, ic2, isr2, osr) <=
3271  distance_within * distance_within;
3272 }
3273
3276  int64_t p1size,
3277  int8_t* p2,
3278  int64_t p2size,
3279  int32_t ic1,
3280  int32_t isr1,
3281  int32_t ic2,
3282  int32_t isr2,
3283  int32_t osr,
3284  double distance_within) {
3285  auto dist_meters =
3286  ST_Distance_Point_Point_Geodesic(p1, p1size, p2, p2size, ic1, isr1, ic2, isr2, osr);
3287  return (dist_meters <= distance_within);
3288 }
3289
3292  int64_t p1size,
3293  int8_t* l2,
3294  int64_t l2size,
3295  double* l2bounds,
3296  int64_t l2bounds_size,
3297  int32_t ic1,
3298  int32_t isr1,
3299  int32_t ic2,
3300  int32_t isr2,
3301  int32_t osr,
3302  double distance_within) {
3303  if (l2bounds) {
3304  // Bounding boxes need to be transformed to output SR before proximity check
3305  if (!point_dwithin_box(
3306  p1, p1size, ic1, isr1, l2bounds, l2bounds_size, isr2, osr, distance_within)) {
3307  return false;
3308  }
3309  }
3310
3311  // May need to adjust the threshold by TOLERANCE_DEFAULT
3312  const double threshold = distance_within;
3314  p1, p1size, l2, l2size, ic1, isr1, ic2, isr2, osr, threshold) <=
3315  distance_within;
3316 }
3317
3320  int64_t psize,
3321  int8_t* poly_coords,
3322  int64_t poly_coords_size,
3323  int32_t* poly_ring_sizes,
3324  int64_t poly_num_rings,
3325  double* poly_bounds,
3326  int64_t poly_bounds_size,
3327  int32_t ic1,
3328  int32_t isr1,
3329  int32_t ic2,
3330  int32_t isr2,
3331  int32_t osr,
3332  double distance_within) {
3333  if (poly_bounds) {
3334  if (!point_dwithin_box(p,
3335  psize,
3336  ic1,
3337  isr1,
3338  poly_bounds,
3339  poly_bounds_size,
3340  isr2,
3341  osr,
3342  distance_within)) {
3343  return false;
3344  }
3345  }
3346
3347  // May need to adjust the threshold by TOLERANCE_DEFAULT
3348  const double threshold = distance_within;
3349  return ST_Distance_Point_Polygon(p,
3350  psize,
3351  poly_coords,
3352  poly_coords_size,
3353  poly_ring_sizes,
3354  poly_num_rings,
3355  ic1,
3356  isr1,
3357  ic2,
3358  isr2,
3359  osr,
3360  threshold) <= distance_within;
3361 }
3362
3365  int64_t psize,
3366  int8_t* mpoly_coords,
3367  int64_t mpoly_coords_size,
3368  int32_t* mpoly_ring_sizes,
3369  int64_t mpoly_num_rings,
3370  int32_t* mpoly_poly_sizes,
3371  int64_t mpoly_num_polys,
3372  double* mpoly_bounds,
3373  int64_t mpoly_bounds_size,
3374  int32_t ic1,
3375  int32_t isr1,
3376  int32_t ic2,
3377  int32_t isr2,
3378  int32_t osr,
3379  double distance_within) {
3380  if (mpoly_bounds) {
3381  if (!point_dwithin_box(p,
3382  psize,
3383  ic1,
3384  isr1,
3385  mpoly_bounds,
3386  mpoly_bounds_size,
3387  isr2,
3388  osr,
3389  distance_within)) {
3390  return false;
3391  }
3392  }
3393
3394  // May need to adjust the threshold by TOLERANCE_DEFAULT
3395  const double threshold = distance_within;
3397  psize,
3398  mpoly_coords,
3399  mpoly_coords_size,
3400  mpoly_ring_sizes,
3401  mpoly_num_rings,
3402  mpoly_poly_sizes,
3403  mpoly_num_polys,
3404  ic1,
3405  isr1,
3406  ic2,
3407  isr2,
3408  osr,
3409  threshold) <= distance_within;
3410 }
3411
3414  int64_t l1size,
3415  double* l1bounds,
3416  int64_t l1bounds_size,
3417  int8_t* l2,
3418  int64_t l2size,
3419  double* l2bounds,
3420  int64_t l2bounds_size,
3421  int32_t ic1,
3422  int32_t isr1,
3423  int32_t ic2,
3424  int32_t isr2,
3425  int32_t osr,
3426  double distance_within) {
3427  if (l1bounds && l2bounds) {
3428  // Bounding boxes need to be transformed to output SR before proximity check
3429  if (!box_dwithin_box(l1bounds,
3430  l1bounds_size,
3431  isr1,
3432  l2bounds,
3433  l2bounds_size,
3434  isr2,
3435  osr,
3436  distance_within)) {
3437  return false;
3438  }
3439  }
3440
3441  // May need to adjust the threshold by TOLERANCE_DEFAULT
3442  const double threshold = distance_within;
3444  l1, l1size, l2, l2size, ic1, isr1, ic2, isr2, osr, threshold) <=
3445  distance_within;
3446 }
3447
3450  int64_t l1size,
3451  double* l1bounds,
3452  int64_t l1bounds_size,
3453  int8_t* poly_coords,
3454  int64_t poly_coords_size,
3455  int32_t* poly_ring_sizes,
3456  int64_t poly_num_rings,
3457  double* poly_bounds,
3458  int64_t poly_bounds_size,
3459  int32_t ic1,
3460  int32_t isr1,
3461  int32_t ic2,
3462  int32_t isr2,
3463  int32_t osr,
3464  double distance_within) {
3465  if (l1bounds && poly_bounds) {
3466  // Bounding boxes need to be transformed to output SR before proximity check
3467  if (!box_dwithin_box(l1bounds,
3468  l1bounds_size,
3469  isr1,
3470  poly_bounds,
3471  poly_bounds_size,
3472  isr2,
3473  osr,
3474  distance_within)) {
3475  return false;
3476  }
3477  }
3478
3479  // First, assume the entire linestring is relevant.
3480  CoordData l1_relevant_section{l1, l1size};
3481  // Mitigation for very long linestrings (5+ segments: think highways, powerlines)
3482  if (poly_bounds && l1size > 12 * compression_unit_size(ic1)) {
3483  // Before diving into distance calculations, try to trim the linestring just to
3484  // segments whose bounding boxes overlap the threshold-buffered poly bounding box.
3485  l1_relevant_section = trim_linestring_to_buffered_box(
3486  l1, l1size, ic1, isr1, poly_bounds, poly_bounds_size, isr2, osr, distance_within);
3487  if (!l1_relevant_section.ptr) {
3488  // The big short-circuit: if not a single segment's bounding box overlaped
3489  // with buffered poly bounding box, then it means that the entire linestring
3490  // doesn't come close enough to the polygon to be within distance.
3491  // The global bounding boxes for the input geometries may overlap but the linestring
3492  // never gets close enough to the poly to require an actual distance calculation.
3493  return false;
3494  }
3495  }
3496
3497  // May need to adjust the threshold by TOLERANCE_DEFAULT
3498  const double threshold = distance_within;
3499  return ST_Distance_LineString_Polygon(l1_relevant_section.ptr,
3500  l1_relevant_section.size,
3501  poly_coords,
3502  poly_coords_size,
3503  poly_ring_sizes,
3504  poly_num_rings,
3505  ic1,
3506  isr1,
3507  ic2,
3508  isr2,
3509  osr,
3510  threshold) <= distance_within;
3511 }
3512
3515  int64_t l1size,
3516  double* l1bounds,
3517  int64_t l1bounds_size,
3518  int8_t* mpoly_coords,
3519  int64_t mpoly_coords_size,
3520  int32_t* mpoly_ring_sizes,
3521  int64_t mpoly_num_rings,
3522  int32_t* mpoly_poly_sizes,
3523  int64_t mpoly_num_polys,
3524  double* mpoly_bounds,
3525  int64_t mpoly_bounds_size,
3526  int32_t ic1,
3527  int32_t isr1,
3528  int32_t ic2,
3529  int32_t isr2,
3530  int32_t osr,
3531  double distance_within) {
3532  if (l1bounds && mpoly_bounds) {
3533  // Bounding boxes need to be transformed to output SR before proximity check
3534  if (!box_dwithin_box(l1bounds,
3535  l1bounds_size,
3536  isr1,
3537  mpoly_bounds,
3538  mpoly_bounds_size,
3539  isr2,
3540  osr,
3541  distance_within)) {
3542  return false;
3543  }
3544  }
3545
3546  // First, assume the entire linestring is relevant.
3547  CoordData l1_relevant_section{l1, l1size};
3548  // Mitigation for very long linestrings (5+ segments: think highways, powerlines)
3549  if (mpoly_bounds && l1size > 12 * compression_unit_size(ic1)) {
3550  // Before diving into distance calculations, try to trim the linestring just to
3551  // segments whose bounding boxes overlap the threshold-buffered mpoly bounding box.
3552  l1_relevant_section = trim_linestring_to_buffered_box(l1,
3553  l1size,
3554  ic1,
3555  isr1,
3556  mpoly_bounds,
3557  mpoly_bounds_size,
3558  isr2,
3559  osr,
3560  distance_within);
3561  if (!l1_relevant_section.ptr) {
3562  // The big short-circuit: if not a single segment's bounding box overlaped
3563  // with buffered mpoly bounding box, then it means that the entire linestring
3564  // doesn't come close enough to the multipolygon to be within distance.
3565  // The global bounding boxes for the input geometries may overlap but the linestring
3566  // never gets close enough to the mpoly to require an actual distance calculation.
3567  return false;
3568  }
3569  }
3570
3571  // Run distance calculation but only on the linestring section that is relevant.
3572  // May need to adjust the threshold by TOLERANCE_DEFAULT
3573  const double threshold = distance_within;
3574  return ST_Distance_LineString_MultiPolygon(l1_relevant_section.ptr,
3575  l1_relevant_section.size,
3576  mpoly_coords,
3577  mpoly_coords_size,
3578  mpoly_ring_sizes,
3579  mpoly_num_rings,
3580  mpoly_poly_sizes,
3581  mpoly_num_polys,
3582  ic1,
3583  isr1,
3584  ic2,
3585  isr2,
3586  osr,
3587  threshold) <= distance_within;
3588 }
3589
3591 bool ST_DWithin_Polygon_Polygon(int8_t* poly1_coords,
3592  int64_t poly1_coords_size,
3593  int32_t* poly1_ring_sizes,
3594  int64_t poly1_num_rings,
3595  double* poly1_bounds,
3596  int64_t poly1_bounds_size,
3597  int8_t* poly2_coords,
3598  int64_t poly2_coords_size,
3599  int32_t* poly2_ring_sizes,
3600  int64_t poly2_num_rings,
3601  double* poly2_bounds,
3602  int64_t poly2_bounds_size,
3603  int32_t ic1,
3604  int32_t isr1,
3605  int32_t ic2,
3606  int32_t isr2,
3607  int32_t osr,
3608  double distance_within) {
3609  if (poly1_bounds && poly2_bounds) {
3610  // Bounding boxes need to be transformed to output SR before proximity check
3611  if (!box_dwithin_box(poly1_bounds,
3612  poly1_bounds_size,
3613  isr1,
3614  poly2_bounds,
3615  poly2_bounds_size,
3616  isr2,
3617  osr,
3618  distance_within)) {
3619  return false;
3620  }
3621  }
3622
3623  // May need to adjust the threshold by TOLERANCE_DEFAULT
3624  const double threshold = distance_within;
3625  return ST_Distance_Polygon_Polygon(poly1_coords,
3626  poly1_coords_size,
3627  poly1_ring_sizes,
3628  poly1_num_rings,
3629  poly2_coords,
3630  poly2_coords_size,
3631  poly2_ring_sizes,
3632  poly2_num_rings,
3633  ic1,
3634  isr1,
3635  ic2,
3636  isr2,
3637  osr,
3638  threshold) <= distance_within;
3639 }
3640
3642 bool ST_DWithin_Polygon_MultiPolygon(int8_t* poly_coords,
3643  int64_t poly_coords_size,
3644  int32_t* poly_ring_sizes,
3645  int64_t poly_num_rings,
3646  double* poly_bounds,
3647  int64_t poly_bounds_size,
3648  int8_t* mpoly_coords,
3649  int64_t mpoly_coords_size,
3650  int32_t* mpoly_ring_sizes,
3651  int64_t mpoly_num_rings,
3652  int32_t* mpoly_poly_sizes,
3653  int64_t mpoly_num_polys,
3654  double* mpoly_bounds,
3655  int64_t mpoly_bounds_size,
3656  int32_t ic1,
3657  int32_t isr1,
3658  int32_t ic2,
3659  int32_t isr2,
3660  int32_t osr,
3661  double distance_within) {
3662  if (poly_bounds && mpoly_bounds) {
3663  // Bounding boxes need to be transformed to output SR before proximity check
3664  if (!box_dwithin_box(poly_bounds,
3665  poly_bounds_size,
3666  isr1,
3667  mpoly_bounds,
3668  mpoly_bounds_size,
3669  isr2,
3670  osr,
3671  distance_within)) {
3672  return false;
3673  }
3674  }
3675
3676  // May need to adjust the threshold by TOLERANCE_DEFAULT
3677  const double threshold = distance_within;
3678  return ST_Distance_Polygon_MultiPolygon(poly_coords,
3679  poly_coords_size,
3680  poly_ring_sizes,
3681  poly_num_rings,
3682  mpoly_coords,
3683  mpoly_coords_size,
3684  mpoly_ring_sizes,
3685  mpoly_num_rings,
3686  mpoly_poly_sizes,
3687  mpoly_num_polys,
3688  ic1,
3689  isr1,
3690  ic2,
3691  isr2,
3692  osr,
3693  threshold) <= distance_within;
3694 }
3695
3697 bool ST_DWithin_MultiPolygon_MultiPolygon(int8_t* mpoly1_coords,
3698  int64_t mpoly1_coords_size,
3699  int32_t* mpoly1_ring_sizes,
3700  int64_t mpoly1_num_rings,
3701  int32_t* mpoly1_poly_sizes,
3702  int64_t mpoly1_num_polys,
3703  double* mpoly1_bounds,
3704  int64_t mpoly1_bounds_size,
3705  int8_t* mpoly2_coords,
3706  int64_t mpoly2_coords_size,
3707  int32_t* mpoly2_ring_sizes,
3708  int64_t mpoly2_num_rings,
3709  int32_t* mpoly2_poly_sizes,
3710  int64_t mpoly2_num_polys,
3711  double* mpoly2_bounds,
3712  int64_t mpoly2_bounds_size,
3713  int32_t ic1,
3714  int32_t isr1,
3715  int32_t ic2,
3716  int32_t isr2,
3717  int32_t osr,
3718  double distance_within) {
3719  if (mpoly1_bounds && mpoly2_bounds) {
3720  // Bounding boxes need to be transformed to output SR before proximity check
3721  if (!box_dwithin_box(mpoly1_bounds,
3722  mpoly1_bounds_size,
3723  isr1,
3724  mpoly2_bounds,
3725  mpoly2_bounds_size,
3726  isr2,
3727  osr,
3728  distance_within)) {
3729  return false;
3730  }
3731  }
3732
3733  // May need to adjust the threshold by TOLERANCE_DEFAULT
3734  const double threshold = distance_within;
3735  return ST_Distance_MultiPolygon_MultiPolygon(mpoly1_coords,
3736  mpoly1_coords_size,
3737  mpoly1_ring_sizes,
3738  mpoly1_num_rings,
3739  mpoly1_poly_sizes,
3740  mpoly1_num_polys,
3741  mpoly2_coords,
3742  mpoly2_coords_size,
3743  mpoly2_ring_sizes,
3744  mpoly2_num_rings,
3745  mpoly2_poly_sizes,
3746  mpoly2_num_polys,
3747  ic1,
3748  isr1,
3749  ic2,
3750  isr2,
3751  osr,
3752  threshold) <= distance_within;
3753 }
3754
3755 //
3756 // ST_MaxDistance
3757 //
3758
3759 // Max cartesian distance between a point and a line segment
3760 DEVICE
3761 double max_distance_point_line(double px,
3762  double py,
3763  double l1x,
3764  double l1y,
3765  double l2x,
3766  double l2y) {
3767  // TODO: switch to squared distances
3768  double length1 = distance_point_point(px, py, l1x, l1y);
3769  double length2 = distance_point_point(px, py, l2x, l2y);
3770  if (length1 > length2) {
3771  return length1;
3772  }
3773  return length2;
3774 }
3775
3777  int64_t psize,
3778  int8_t* l,
3779  int64_t lsize,
3780  int32_t ic1,
3781  int32_t isr1,
3782  int32_t ic2,
3783  int32_t isr2,
3784  int32_t osr,
3785  bool check_closed) {
3786  // TODO: switch to squared distances
3787  Point2D pt = get_point(p, 0, ic1, isr1, osr);
3788
3789  auto l_num_coords = lsize / compression_unit_size(ic2);
3790
3791  Point2D l1 = get_point(l, 0, ic2, isr2, osr);
3792  Point2D l2 = get_point(l, 2, ic2, isr2, osr);
3793
3794  double max_dist = max_distance_point_line(pt.x, pt.y, l1.x, l1.y, l2.x, l2.y);
3795  for (int32_t i = 4; i < l_num_coords; i += 2) {
3796  l1 = l2; // advance one point
3797  l2 = get_point(l, i, ic2, isr2, osr);
3798  double ldist = max_distance_point_line(pt.x, pt.y, l1.x, l1.y, l2.x, l2.y);
3799  if (max_dist < ldist) {
3800  max_dist = ldist;
3801  }
3802  }
3803  if (l_num_coords > 4 && check_closed) {
3804  // Also check distance to the closing edge between the first and the last points
3805  l1 = get_point(l, 0, ic2, isr2, osr);
3806  double ldist = max_distance_point_line(pt.x, pt.y, l1.x, l1.y, l2.x, l2.y);
3807  if (max_dist < ldist) {
3808  max_dist = ldist;
3809  }
3810  }
3811  return max_dist;
3812 }
3813
3816  int64_t psize,
3817  int8_t* l,
3818  int64_t lsize,
3819  int32_t ic1,
3820  int32_t isr1,
3821  int32_t ic2,
3822  int32_t isr2,
3823  int32_t osr) {
3825  p, psize, l, lsize, ic1, isr1, ic2, isr2, osr, false);
3826 }
3827
3830  int64_t lsize,
3831  int8_t* p,
3832  int64_t psize,
3833  int32_t ic1,
3834  int32_t isr1,
3835  int32_t ic2,
3836  int32_t isr2,
3837  int32_t osr) {
3839  p, psize, l, lsize, ic2, isr2, ic1, isr1, osr, false);
3840 }
3841
3842 // TODO: add ST_MaxDistance_LineString_LineString (with short-circuit threshold)
3843
3844 //
3845 // ST_Contains
3846 //
3847
3849 bool ST_Contains_Point_Point(int8_t* p1,
3850  int64_t p1size,
3851  int8_t* p2,
3852  int64_t p2size,
3853  int32_t ic1,
3854  int32_t isr1,
3855  int32_t ic2,
3856  int32_t isr2,
3857  int32_t osr) {
3858  Point2D const pt1 = get_point(p1, 0, ic1, isr1, osr);
3859  Point2D const pt2 = get_point(p2, 0, ic2, isr2, osr);
3860  double tolerance = tol(ic1, ic2);
3862 }
3863
3866  int64_t psize,
3867  int8_t* l,
3868  int64_t lsize,
3869  double* lbounds,
3870  int64_t lbounds_size,
3871  int32_t ic1,
3872  int32_t isr1,
3873  int32_t ic2,
3874  int32_t isr2,
3875  int32_t osr) {
3876  Point2D const pt = get_point(p, 0, ic1, isr1, osr);
3877
3878  if (lbounds) {
3879  if (tol_eq(pt.x, lbounds[0]) && tol_eq(pt.y, lbounds[1]) &&
3880  tol_eq(pt.x, lbounds[2]) && tol_eq(pt.y, lbounds[3])) {
3881  return true;
3882  }
3883  }
3884
3885  auto l_num_coords = lsize / compression_unit_size(ic2);
3886  for (int i = 0; i < l_num_coords; i += 2) {
3887  Point2D pl = get_point(l, i, ic2, isr2, osr);
3888  if (tol_eq(pt.x, pl.x) && tol_eq(pt.y, pl.y)) {
3889  continue;
3890  }
3891  return false;
3892  }
3893  return true;
3894 }
3895
3898  int64_t psize,
3899  int8_t* poly_coords,
3900  int64_t poly_coords_size,
3901  int32_t* poly_ring_sizes,
3902  int64_t poly_num_rings,
3903  double* poly_bounds,
3904  int64_t poly_bounds_size,
3905  int32_t ic1,
3906  int32_t isr1,
3907  int32_t ic2,
3908  int32_t isr2,
3909  int32_t osr) {
3910  auto exterior_ring_num_coords = poly_coords_size / compression_unit_size(ic2);
3911  if (poly_num_rings > 0) {
3912  exterior_ring_num_coords = poly_ring_sizes[0] * 2;
3913  }
3914  auto exterior_ring_coords_size = exterior_ring_num_coords * compression_unit_size(ic2);
3915
3917  psize,
3918  poly_coords,
3919  exterior_ring_coords_size,
3920  poly_bounds,
3921  poly_bounds_size,
3922  ic1,
3923  isr1,
3924  ic2,
3925  isr2,
3926  osr);
3927 }
3928
3931  int64_t lsize,
3932  double* lbounds,
3933  int64_t lbounds_size,
3934  int8_t* p,
3935  int64_t psize,
3936  int32_t ic1,
3937  int32_t isr1,
3938  int32_t ic2,
3939  int32_t isr2,
3940  int32_t osr) {
3942  ST_Distance_Point_LineString(p, psize, l, lsize, ic2, isr2, ic1, isr1, osr, 0.0));
3943 }
3944
3947  int64_t lsize,
3948  double* lbounds,
3949  int64_t lbounds_size,
3950  int8_t* poly_coords,
3951  int64_t poly_coords_size,
3952  int32_t* poly_ring_sizes,
3953  int64_t poly_num_rings,
3954  double* poly_bounds,
3955  int64_t poly_bounds_size,
3956  int32_t ic1,
3957  int32_t isr1,
3958  int32_t ic2,
3959  int32_t isr2,
3960  int32_t osr) {
3961  // TODO
3962  return false;
3963 }
3964
3965 template <typename T, EdgeBehavior TEdgeBehavior>
3966 DEVICE ALWAYS_INLINE bool Contains_Polygon_Point_Impl(const int8_t* poly_coords,
3967  const int64_t poly_coords_size,
3968  const int32_t* poly_ring_sizes,
3969  const int64_t poly_num_rings,
3970  const double* poly_bounds,
3971  const int64_t poly_bounds_size,
3972  const int8_t* p,
3973  const int64_t psize,
3974  const int32_t ic1,
3975  const int32_t isr1,
3976  const int32_t ic2,
3977  const int32_t isr2,
3978  const int32_t osr) {
3979  if (poly_bounds) {
3980  // TODO: codegen
3981  Point2D const pt = get_point(p, 0, ic2, isr2, osr);
3982  if (!box_contains_point(poly_bounds, poly_bounds_size, pt.x, pt.y)) {
3983  DEBUG_STMT(printf("Bounding box does not contain point, exiting.\n"));
3984  return false;
3985  }
3986  }
3987
3988  auto get_x_coord = [=]() -> T {
3989  if constexpr (std::is_floating_point<T>::value) {
3990  return get_point(p, 0, ic2, isr2, osr).x;
3991  } else {
3992  return compressed_coord(p, 0);
3993  }
3994  return T{}; // https://stackoverflow.com/a/64561686/2700898
3995  };
3996
3997  auto get_y_coord = [=]() -> T {
3998  if constexpr (std::is_floating_point<T>::value) {
3999  return get_point(p, 0, ic2, isr2, osr).y;
4000  } else {
4001  return compressed_coord(p, 1);
4002  }
4003  return T{}; // https://stackoverflow.com/a/64561686/2700898
4004  };
4005
4006  const T px = get_x_coord();
4007  const T py = get_y_coord();
4008  DEBUG_STMT(printf("pt: %f, %f\n", (double)px, (double)py));
4009
4010  const auto poly_num_coords = poly_coords_size / compression_unit_size(ic1);
4011  auto exterior_ring_num_coords =
4012  poly_num_rings > 0 ? poly_ring_sizes[0] * 2 : poly_num_coords;
4013
4014  auto poly = poly_coords;
4015  if (point_in_polygon_winding_number<T, TEdgeBehavior>(
4016  poly, exterior_ring_num_coords, px, py, ic1, isr1, osr)) {
4017  // Inside exterior ring
4018  poly += exterior_ring_num_coords * compression_unit_size(ic1);
4019  // Check that none of the polygon's holes contain that point
4020  for (auto r = 1; r < poly_num_rings; r++) {
4021  const int64_t interior_ring_num_coords = poly_ring_sizes[r] * 2;
4022  if (point_in_polygon_winding_number<T, TEdgeBehavior>(
4023  poly, interior_ring_num_coords, px, py, ic1, isr1, osr)) {
4024  return false;
4025  }
4026  poly += interior_ring_num_coords * compression_unit_size(ic1);
4027  }
4028  return true;
4029  }
4030  return false;
4031 }
4032
4033 EXTENSION_NOINLINE bool ST_Contains_Polygon_Point(const int8_t* poly_coords,
4034  const int64_t poly_coords_size,
4035  const int32_t* poly_ring_sizes,
4036  const int64_t poly_num_rings,
4037  const double* poly_bounds,
4038  const int64_t poly_bounds_size,
4039  const int8_t* p,
4040  const int64_t psize,
4041  const int32_t ic1,
4042  const int32_t isr1,
4043  const int32_t ic2,
4044  const int32_t isr2,
4045  const int32_t osr) {
4046  return Contains_Polygon_Point_Impl<double, EdgeBehavior::kExcludePointOnEdge>(
4047  poly_coords,
4048  poly_coords_size,
4049  poly_ring_sizes,
4050  poly_num_rings,
4051  poly_bounds,
4052  poly_bounds_size,
4053  p,
4054  psize,
4055  ic1,
4056  isr1,
4057  ic2,
4058  isr2,
4059  osr);
4060 }
4061
4062 EXTENSION_NOINLINE bool ST_cContains_Polygon_Point(const int8_t* poly_coords,
4063  const int64_t poly_coords_size,
4064  const int32_t* poly_ring_sizes,
4065  const int64_t poly_num_rings,
4066  const double* poly_bounds,
4067  const int64_t poly_bounds_size,
4068  const int8_t* p,
4069  const int64_t psize,
4070  const int32_t ic1,
4071  const int32_t isr1,
4072  const int32_t ic2,
4073  const int32_t isr2,
4074  const int32_t osr) {
4075  return Contains_Polygon_Point_Impl<int64_t, EdgeBehavior::kExcludePointOnEdge>(
4076  poly_coords,
4077  poly_coords_size,
4078  poly_ring_sizes,
4079  poly_num_rings,
4080  poly_bounds,
4081  poly_bounds_size,
4082  p,
4083  psize,
4084  ic1,
4085  isr1,
4086  ic2,
4087  isr2,
4088  osr);
4089 }
4090
4092 bool ST_Contains_Polygon_LineString(int8_t* poly_coords,
4093  int64_t poly_coords_size,
4094  int32_t* poly_ring_sizes,
4095  int64_t poly_num_rings,
4096  double* poly_bounds,
4097  int64_t poly_bounds_size,
4098  int8_t* l,
4099  int64_t lsize,
4100  double* lbounds,
4101  int64_t lbounds_size,
4102  int32_t ic1,
4103  int32_t isr1,
4104  int32_t ic2,
4105  int32_t isr2,
4106  int32_t osr) {
4107  if (poly_num_rings > 1) {
4108  return false; // TODO: support polygons with interior rings
4109  }
4110
4111  // Bail out if poly bounding box doesn't contain linestring bounding box
4112  if (poly_bounds && lbounds) {
4113  if (!box_contains_box(poly_bounds, poly_bounds_size, lbounds, lbounds_size)) {
4114  return false;
4115  }
4116  }
4117
4118  const auto poly_num_coords = poly_coords_size / compression_unit_size(ic1);
4119  const auto lnum_coords = lsize / compression_unit_size(ic2);
4120
4122  poly_coords, poly_num_coords, l, lnum_coords, ic1, isr1, ic2, isr2, osr);
4123 }
4124
4126 bool ST_Contains_Polygon_Polygon(int8_t* poly1_coords,
4127  int64_t poly1_coords_size,
4128  int32_t* poly1_ring_sizes,
4129  int64_t poly1_num_rings,
4130  double* poly1_bounds,
4131  int64_t poly1_bounds_size,
4132  int8_t* poly2_coords,
4133  int64_t poly2_coords_size,
4134  int32_t* poly2_ring_sizes,
4135  int64_t poly2_num_rings,
4136  double* poly2_bounds,
4137  int64_t poly2_bounds_size,
4138  int32_t ic1,
4139  int32_t isr1,
4140  int32_t ic2,
4141  int32_t isr2,
4142  int32_t osr) {
4143  // TODO: needs to be extended, cover more cases
4144  // Right now only checking if simple poly1 (no holes) contains poly2's exterior shape
4145  if (poly1_num_rings > 1) {
4146  return false; // TODO: support polygons with interior rings
4147  }
4148
4149  if (poly1_bounds && poly2_bounds) {
4150  if (!box_contains_box(
4151  poly1_bounds, poly1_bounds_size, poly2_bounds, poly2_bounds_size)) {
4152  return false;
4153  }
4154  }
4155
4156  int64_t poly2_exterior_ring_coords_size = poly2_coords_size;
4157  if (poly2_num_rings > 0) {
4158  poly2_exterior_ring_coords_size =
4159  2 * poly2_ring_sizes[0] * compression_unit_size(ic2);
4160  }
4161  return ST_Contains_Polygon_LineString(poly1_coords,
4162  poly1_coords_size,
4163  poly1_ring_sizes,
4164  poly1_num_rings,
4165  poly1_bounds,
4166  poly1_bounds_size,
4167  poly2_coords,
4168  poly2_exterior_ring_coords_size,
4169  poly2_bounds,
4170  poly2_bounds_size,
4171  ic1,
4172  isr1,
4173  ic2,
4174  isr2,
4175  osr);
4176 }
4177
4178 template <typename T, EdgeBehavior TEdgeBehavior>
4180  const int8_t* mpoly_coords,
4181  const int64_t mpoly_coords_size,
4182  const int32_t* mpoly_ring_sizes,
4183  const int64_t mpoly_num_rings,
4184  const int32_t* mpoly_poly_sizes,
4185  const int64_t mpoly_num_polys,
4186  const double* mpoly_bounds,
4187  const int64_t mpoly_bounds_size,
4188  const int8_t* p,
4189  const int64_t psize,
4190  const int32_t ic1,
4191  const int32_t isr1,
4192  const int32_t ic2,
4193  const int32_t isr2,
4194  const int32_t osr) {
4195  if (mpoly_num_polys <= 0) {
4196  return false;
4197  }
4198
4199  // TODO: mpoly_bounds could contain individual bounding boxes too:
4200  // first two points - box for the entire multipolygon, then a pair for each polygon
4201  if (mpoly_bounds) {
4202  Point2D const pt = get_point(p, 0, ic2, isr2, osr);
4203  if (!box_contains_point(mpoly_bounds, mpoly_bounds_size, pt.x, pt.y)) {
4204  return false;
4205  }
4206  }
4207
4208  // Set specific poly pointers as we move through the coords/ringsizes/polyrings
4209  // arrays.
4210  auto next_poly_coords = mpoly_coords;
4211  auto next_poly_ring_sizes = mpoly_ring_sizes;
4212
4213  for (auto poly = 0; poly < mpoly_num_polys; poly++) {
4214  const auto poly_coords = next_poly_coords;
4215  const auto poly_ring_sizes = next_poly_ring_sizes;
4216  const auto poly_num_rings = mpoly_poly_sizes[poly];
4217  // Count number of coords in all of poly's rings, advance ring size pointer.
4218  int32_t poly_num_coords = 0;
4219  for (auto ring = 0; ring < poly_num_rings; ring++) {
4220  poly_num_coords += 2 * *next_poly_ring_sizes++;
4221  }
4222  const auto poly_coords_size = poly_num_coords * compression_unit_size(ic1);
4223  next_poly_coords += poly_coords_size;
4224  // TODO: pass individual bounding boxes for each polygon
4225  if (Contains_Polygon_Point_Impl<T, TEdgeBehavior>(poly_coords,
4226  poly_coords_size,
4227  poly_ring_sizes,
4228  poly_num_rings,
4229  nullptr,
4230  0,
4231  p,
4232  psize,
4233  ic1,
4234  isr1,
4235  ic2,
4236  isr2,
4237  osr)) {
4238  return true;
4239  }
4240  }
4241
4242  return false;
4243 }
4244
4246 bool ST_Contains_MultiPolygon_Point(int8_t* mpoly_coords,
4247  int64_t mpoly_coords_size,
4248  int32_t* mpoly_ring_sizes,
4249  int64_t mpoly_num_rings,
4250  int32_t* mpoly_poly_sizes,
4251  int64_t mpoly_num_polys,
4252  double* mpoly_bounds,
4253  int64_t mpoly_bounds_size,
4254  int8_t* p,
4255  int64_t psize,
4256  int32_t ic1,
4257  int32_t isr1,
4258  int32_t ic2,
4259  int32_t isr2,
4260  int32_t osr) {
4261  return Contains_MultiPolygon_Point_Impl<double, EdgeBehavior::kExcludePointOnEdge>(
4262  mpoly_coords,
4263  mpoly_coords_size,
4264  mpoly_ring_sizes,
4265  mpoly_num_rings,
4266  mpoly_poly_sizes,
4267  mpoly_num_polys,
4268  mpoly_bounds,
4269  mpoly_bounds_size,
4270  p,
4271  psize,
4272  ic1,
4273  isr1,
4274  ic2,
4275  isr2,
4276  osr);
4277 }
4278
4280  int64_t mpoly_coords_size,
4281  int32_t* mpoly_ring_sizes,
4282  int64_t mpoly_num_rings,
4283  int32_t* mpoly_poly_sizes,
4284  int64_t mpoly_num_polys,
4285  double* mpoly_bounds,
4286  int64_t mpoly_bounds_size,
4287  int8_t* p,
4288  int64_t psize,
4289  int32_t ic1,
4290  int32_t isr1,
4291  int32_t ic2,
4292  int32_t isr2,
4293  int32_t osr) {
4294  return Contains_MultiPolygon_Point_Impl<int64_t, EdgeBehavior::kExcludePointOnEdge>(
4295  mpoly_coords,
4296  mpoly_coords_size,
4297  mpoly_ring_sizes,
4298  mpoly_num_rings,
4299  mpoly_poly_sizes,
4300  mpoly_num_polys,
4301  mpoly_bounds,
4302  mpoly_bounds_size,
4303  p,
4304  psize,
4305  ic1,
4306  isr1,
4307  ic2,
4308  isr2,
4309  osr);
4310 }
4311
4313 bool ST_Contains_MultiPolygon_LineString(int8_t* mpoly_coords,
4314  int64_t mpoly_coords_size,
4315  int32_t* mpoly_ring_sizes,
4316  int64_t mpoly_num_rings,
4317  int32_t* mpoly_poly_sizes,
4318  int64_t mpoly_num_polys,
4319  double* mpoly_bounds,
4320  int64_t mpoly_bounds_size,
4321  int8_t* l,
4322  int64_t lsize,
4323  double* lbounds,
4324  int64_t lbounds_size,
4325  int32_t ic1,
4326  int32_t isr1,
4327  int32_t ic2,
4328  int32_t isr2,
4329  int32_t osr) {
4330  if (mpoly_num_polys <= 0) {
4331  return false;
4332  }
4333
4334  if (mpoly_bounds && lbounds) {
4335  if (!box_contains_box(mpoly_bounds, mpoly_bounds_size, lbounds, lbounds_size)) {
4336  return false;
4337  }
4338  }
4339
4340  // Set specific poly pointers as we move through the coords/ringsizes/polyrings arrays.
4341  auto next_poly_coords = mpoly_coords;
4342  auto next_poly_ring_sizes = mpoly_ring_sizes;
4343
4344  for (auto poly = 0; poly < mpoly_num_polys; poly++) {
4345  auto poly_coords = next_poly_coords;
4346  auto poly_ring_sizes = next_poly_ring_sizes;
4347  auto poly_num_rings = mpoly_poly_sizes[poly];
4348  // Count number of coords in all of poly's rings, advance ring size pointer.
4349  int32_t poly_num_coords = 0;
4350  for (auto ring = 0; ring < poly_num_rings; ring++) {
4351  poly_num_coords += 2 * *next_poly_ring_sizes++;
4352  }
4353  auto poly_coords_size = poly_num_coords * compression_unit_size(ic1);
4354  next_poly_coords += poly_coords_size;
4355
4356  if (ST_Contains_Polygon_LineString(poly_coords,
4357  poly_coords_size,
4358  poly_ring_sizes,
4359  poly_num_rings,
4360  nullptr,
4361  0,
4362  l,
4363  lsize,
4364  nullptr,
4365  0,
4366  ic1,
4367  isr1,
4368  ic2,
4369  isr2,
4370  osr)) {
4371  return true;
4372  }
4373  }
4374
4375  return false;
4376 }
4377
4378 //
4379 // ST_Intersects
4380 //
4381
4384  int64_t p1size,
4385  int8_t* p2,
4386  int64_t p2size,
4387  int32_t ic1,
4388  int32_t isr1,
4389  int32_t ic2,
4390  int32_t isr2,
4391  int32_t osr) {
4393  ST_Distance_Point_Point(p1, p1size, p2, p2size, ic1, isr1, ic2, isr2, osr));
4394 }
4395
4398  int64_t psize,
4399  int8_t* l,
4400  int64_t lsize,
4401  double* lbounds,
4402  int64_t lbounds_size,
4403  int32_t ic1,
4404  int32_t isr1,
4405  int32_t ic2,
4406  int32_t isr2,
4407  int32_t osr) {
4408  if (lbounds) {
4409  Point2D const pt = get_point(p, 0, ic1, isr1, osr);
4410  if (!box_contains_point(lbounds, lbounds_size, pt.x, pt.y)) {
4411  return false;
4412  }
4413  }
4415  ST_Distance_Point_LineString(p, psize, l, lsize, ic1, isr1, ic2, isr2, osr, 0.0));
4416 }
4417
4420  int64_t psize,
4421  int8_t* poly,
4422  int64_t polysize,
4423  int32_t* poly_ring_sizes,
4424  int64_t poly_num_rings,
4425  double* poly_bounds,
4426  int64_t poly_bounds_size,
4427  int32_t ic1,
4428  int32_t isr1,
4429  int32_t ic2,
4430  int32_t isr2,
4431  int32_t osr) {
4432  return Contains_Polygon_Point_Impl<double, EdgeBehavior::kIncludePointOnEdge>(
4433  poly,
4434  polysize,
4435  poly_ring_sizes,
4436  poly_num_rings,
4437  poly_bounds,
4438  poly_bounds_size,
4439  p,
4440  psize,
4441  ic2,
4442  isr2,
4443  ic1,
4444  isr1,
4445  osr);
4446 }
4447
4450  int64_t psize,
4451  int8_t* mpoly_coords,
4452  int64_t mpoly_coords_size,
4453  int32_t* mpoly_ring_sizes,
4454  int64_t mpoly_num_rings,
4455  int32_t* mpoly_poly_sizes,
4456  int64_t mpoly_num_polys,
4457  double* mpoly_bounds,
4458  int64_t mpoly_bounds_size,
4459  int32_t ic1,
4460  int32_t isr1,
4461  int32_t ic2,
4462  int32_t isr2,
4463  int32_t osr) {
4464  return Contains_MultiPolygon_Point_Impl<double, EdgeBehavior::kIncludePointOnEdge>(
4465  mpoly_coords,
4466  mpoly_coords_size,
4467  mpoly_ring_sizes,
4468  mpoly_num_rings,
4469  mpoly_poly_sizes,
4470  mpoly_num_polys,
4471  mpoly_bounds,
4472  mpoly_bounds_size,
4473  p,
4474  psize,
4475  ic1,
4476  isr1,
4477  ic2,
4478  isr2,
4479  osr);
4480 }
4481
4484  int64_t lsize,
4485  double* lbounds,
4486  int64_t lbounds_size,
4487  int8_t* p,
4488  int64_t psize,
4489  int32_t ic1,
4490  int32_t isr1,
4491  int32_t ic2,
4492  int32_t isr2,
4493  int32_t osr) {
4495  p, psize, l, lsize, lbounds, lbounds_size, ic2, isr2, ic1, isr1, osr);
4496 }
4497
4500  int64_t l1size,
4501  double* l1bounds,
4502  int64_t l1bounds_size,
4503  int8_t* l2,
4504  int64_t l2size,
4505  double* l2bounds,
4506  int64_t l2bounds_size,
4507  int32_t ic1,
4508  int32_t isr1,
4509  int32_t ic2,
4510  int32_t isr2,
4511  int32_t osr) {
4512  if (l1bounds && l2bounds) {
4513  if (!box_overlaps_box(l1bounds, l1bounds_size, l2bounds, l2bounds_size)) {
4514  return false;
4515  }
4516  }
4517
4519  l1, l1size, l2, l2size, ic1, isr1, ic2, isr2, osr, 0.0));
4520 }
4521
4524  int64_t lsize,
4525  double* lbounds,
4526  int64_t lbounds_size,
4527  int8_t* poly,
4528  int64_t polysize,
4529  int32_t* poly_ring_sizes,
4530  int64_t poly_num_rings,
4531  double* poly_bounds,
4532  int64_t poly_bounds_size,
4533  int32_t ic1,
4534  int32_t isr1,
4535  int32_t ic2,
4536  int32_t isr2,
4537  int32_t osr) {
4538  if (lbounds && poly_bounds) {
4539  if (!box_overlaps_box(lbounds, lbounds_size, poly_bounds, poly_bounds_size)) {
4540  return false;
4541  }
4542  }
4543
4544  // Check for spatial intersection.
4545  // One way to do that would be to start with linestring's first point, if it's inside
4546  // the polygon - it means intersection. Otherwise follow the linestring, segment by
4547  // segment, checking for intersections with polygon rings, bail as soon as we cross into
4548  // the polygon.
4549
4550  // Or, alternatively, just measure the distance:
4552  lsize,
4553  poly,
4554  polysize,
4555  poly_ring_sizes,
4556  poly_num_rings,
4557  ic1,
4558  isr1,
4559  ic2,
4560  isr2,
4561  osr,
4562  0.0));
4563 }
4564
4567  int64_t lsize,
4568  double* lbounds,
4569  int64_t lbounds_size,
4570  int8_t* mpoly_coords,
4571  int64_t mpoly_coords_size,
4572  int32_t* mpoly_ring_sizes,
4573  int64_t mpoly_num_rings,
4574  int32_t* mpoly_poly_sizes,
4575  int64_t mpoly_num_polys,
4576  double* mpoly_bounds,
4577  int64_t mpoly_bounds_size,
4578  int32_t ic1,
4579  int32_t isr1,
4580  int32_t ic2,
4581  int32_t isr2,
4582  int32_t osr) {
4583  if (lbounds && mpoly_bounds) {
4584  if (!box_overlaps_box(lbounds, lbounds_size, mpoly_bounds, mpoly_bounds_size)) {
4585  return false;
4586  }
4587  }
4588
4589  // Check for spatial intersection.
4590  // One way to do that would be to start with linestring's first point, if it's inside
4591  // any of the polygons - it means intersection. Otherwise follow the linestring, segment
4592  // by segment, checking for intersections with polygon shapes/holes, bail as soon as we
4593  // cross into a polygon.
4594
4595  // Or, alternatively, just measure the distance:
4597  lsize,
4598  mpoly_coords,
4599  mpoly_coords_size,
4600  mpoly_ring_sizes,
4601  mpoly_num_rings,
4602  mpoly_poly_sizes,
4603  mpoly_num_polys,
4604  ic1,
4605  isr1,
4606  ic2,
4607  isr2,
4608  osr,
4609  0.0));
4610 }
4611
4614  int64_t polysize,
4615  int32_t* poly_ring_sizes,
4616  int64_t poly_num_rings,
4617  double* poly_bounds,
4618  int64_t poly_bounds_size,
4619  int8_t* p,
4620  int64_t psize,
4621  int32_t ic1,
4622  int32_t isr1,
4623  int32_t ic2,
4624  int32_t isr2,
4625  int32_t osr) {
4626  return Contains_Polygon_Point_Impl<double, EdgeBehavior::kIncludePointOnEdge>(
4627  poly,
4628  polysize,
4629  poly_ring_sizes,
4630  poly_num_rings,
4631  poly_bounds,
4632  poly_bounds_size,
4633  p,
4634  psize,
4635  ic1,
4636  isr1,
4637  ic2,
4638  isr2,
4639  osr);
4640 }
4641
4644  int64_t polysize,
4645  int32_t* poly_ring_sizes,
4646  int64_t poly_num_rings,
4647  double* poly_bounds,
4648  int64_t poly_bounds_size,
4649  int8_t* p,
4650  int64_t psize,
4651  int32_t ic1,
4652  int32_t isr1,
4653  int32_t ic2,
4654  int32_t isr2,
4655  int32_t osr) {
4656  return Contains_Polygon_Point_Impl<int64_t, EdgeBehavior::kIncludePointOnEdge>(
4657  poly,
4658  polysize,
4659  poly_ring_sizes,
4660  poly_num_rings,
4661  poly_bounds,
4662  poly_bounds_size,
4663  p,
4664  psize,
4665  ic1,
4666  isr1,
4667  ic2,
4668  isr2,
4669  osr);
4670 }
4671
4674  int64_t polysize,
4675  int32_t* poly_ring_sizes,
4676  int64_t poly_num_rings,
4677  double* poly_bounds,
4678  int64_t poly_bounds_size,
4679  int8_t* l,
4680  int64_t lsize,
4681  double* lbounds,
4682  int64_t lbounds_size,
4683  int32_t ic1,
4684  int32_t isr1,
4685  int32_t ic2,
4686  int32_t isr2,
4687  int32_t osr) {
4689  lsize,
4690  lbounds,
4691  lbounds_size,
4692  poly,
4693  polysize,
4694  poly_ring_sizes,
4695  poly_num_rings,
4696  poly_bounds,
4697  poly_bounds_size,
4698  ic2,
4699  isr2,
4700  ic1,
4701  isr1,
4702  osr);
4703 }
4704
4706 bool ST_Intersects_Polygon_Polygon(int8_t* poly1_coords,
4707  int64_t poly1_coords_size,
4708  int32_t* poly1_ring_sizes,
4709  int64_t poly1_num_rings,
4710  double* poly1_bounds,
4711  int64_t poly1_bounds_size,
4712  int8_t* poly2_coords,
4713  int64_t poly2_coords_size,
4714  int32_t* poly2_ring_sizes,
4715  int64_t poly2_num_rings,
4716  double* poly2_bounds,
4717  int64_t poly2_bounds_size,
4718  int32_t ic1,
4719  int32_t isr1,
4720  int32_t ic2,
4721  int32_t isr2,
4722  int32_t osr) {
4723  if (poly1_bounds && poly2_bounds) {
4724  if (!box_overlaps_box(
4725  poly1_bounds, poly1_bounds_size, poly2_bounds, poly2_bounds_size)) {
4726  return false;
4727  }
4728  }
4729
4731  poly1_coords_size,
4732  poly1_ring_sizes,
4733  poly1_num_rings,
4734  poly2_coords,
4735  poly2_coords_size,
4736  poly2_ring_sizes,
4737  poly2_num_rings,
4738  ic1,
4739  isr1,
4740  ic2,
4741  isr2,
4742  osr,
4743  0.0));
4744 }
4745
4747 bool ST_Intersects_Polygon_MultiPolygon(int8_t* poly_coords,
4748  int64_t poly_coords_size,
4749  int32_t* poly_ring_sizes,
4750  int64_t poly_num_rings,
4751  double* poly_bounds,
4752  int64_t poly_bounds_size,
4753  int8_t* mpoly_coords,
4754  int64_t mpoly_coords_size,
4755  int32_t* mpoly_ring_sizes,
4756  int64_t mpoly_num_rings,
4757  int32_t* mpoly_poly_sizes,
4758  int64_t mpoly_num_polys,
4759  double* mpoly_bounds,
4760  int64_t mpoly_bounds_size,
4761  int32_t ic1,
4762  int32_t isr1,
4763  int32_t ic2,
4764  int32_t isr2,
4765  int32_t osr) {
4766  if (poly_bounds && mpoly_bounds) {
4767  if (!box_overlaps_box(
4768  poly_bounds, poly_bounds_size, mpoly_bounds, mpoly_bounds_size)) {
4769  return false;
4770  }
4771  }
4772
4774  poly_coords_size,
4775  poly_ring_sizes,
4776  poly_num_rings,
4777  mpoly_coords,
4778  mpoly_coords_size,
4779  mpoly_ring_sizes,
4780  mpoly_num_rings,
4781  mpoly_poly_sizes,
4782  mpoly_num_polys,
4783  ic1,
4784  isr1,
4785  ic2,
4786  isr2,
4787  osr,
4788  0.0));
4789 }
4790
4792 bool ST_Intersects_MultiPolygon_Point(int8_t* mpoly_coords,
4793  int64_t mpoly_coords_size,
4794  int32_t* mpoly_ring_sizes,
4795  int64_t mpoly_num_rings,
4796  int32_t* mpoly_poly_sizes,
4797  int64_t mpoly_num_polys,
4798  double* mpoly_bounds,
4799  int64_t mpoly_bounds_size,
4800  int8_t* p,
4801  int64_t psize,
4802  int32_t ic1,
4803  int32_t isr1,
4804  int32_t ic2,
4805  int32_t isr2,
4806  int32_t osr) {
4807  return Contains_MultiPolygon_Point_Impl<double, EdgeBehavior::kIncludePointOnEdge>(
4808  mpoly_coords,
4809  mpoly_coords_size,
4810  mpoly_ring_sizes,
4811  mpoly_num_rings,
4812  mpoly_poly_sizes,
4813  mpoly_num_polys,
4814  mpoly_bounds,
4815  mpoly_bounds_size,
4816  p,
4817  psize,
4818  ic1,
4819  isr1,
4820  ic2,
4821  isr2,
4822  osr);
4823 }
4824
4826 bool ST_cIntersects_MultiPolygon_Point(int8_t* mpoly_coords,
4827  int64_t mpoly_coords_size,
4828  int32_t* mpoly_ring_sizes,
4829  int64_t mpoly_num_rings,
4830  int32_t* mpoly_poly_sizes,
4831  int64_t mpoly_num_polys,
4832  double* mpoly_bounds,
4833  int64_t mpoly_bounds_size,
4834  int8_t* p,
4835  int64_t psize,
4836  int32_t ic1,
4837  int32_t isr1,
4838  int32_t ic2,
4839  int32_t isr2,
4840  int32_t osr) {
4841  return Contains_MultiPolygon_Point_Impl<int64_t, EdgeBehavior::kIncludePointOnEdge>(
4842  mpoly_coords,
4843  mpoly_coords_size,
4844  mpoly_ring_sizes,
4845  mpoly_num_rings,
4846  mpoly_poly_sizes,
4847  mpoly_num_polys,
4848  mpoly_bounds,
4849  mpoly_bounds_size,
4850  p,
4851  psize,
4852  ic1,
4853  isr1,
4854  ic2,
4855  isr2,
4856  osr);
4857 }
4858
4860 bool ST_Intersects_MultiPolygon_LineString(int8_t* mpoly_coords,
4861  int64_t mpoly_coords_size,
4862  int32_t* mpoly_ring_sizes,
4863  int64_t mpoly_num_rings,
4864  int32_t* mpoly_poly_sizes,
4865  int64_t mpoly_num_polys,
4866  double* mpoly_bounds,
4867  int64_t mpoly_bounds_size,
4868  int8_t* l,
4869  int64_t lsize,
4870  double* lbounds,
4871  int64_t lbounds_size,
4872  int32_t ic1,
4873  int32_t isr1,
4874  int32_t ic2,
4875  int32_t isr2,
4876  int32_t osr) {
4878  lsize,
4879  lbounds,
4880  lbounds_size,
4881  mpoly_coords,
4882  mpoly_coords_size,
4883  mpoly_ring_sizes,
4884  mpoly_num_rings,
4885  mpoly_poly_sizes,
4886  mpoly_num_polys,
4887  mpoly_bounds,
4888  mpoly_bounds_size,
4889  ic2,
4890  isr2,
4891  ic1,
4892  isr1,
4893  osr);
4894 }
4895
4897 bool ST_Intersects_MultiPolygon_Polygon(int8_t* mpoly_coords,
4898  int64_t mpoly_coords_size,
4899  int32_t* mpoly_ring_sizes,
4900  int64_t mpoly_num_rings,
4901  int32_t* mpoly_poly_sizes,
4902  int64_t mpoly_num_polys,
4903  double* mpoly_bounds,
4904  int64_t mpoly_bounds_size,
4905  int8_t* poly_coords,
4906  int64_t poly_coords_size,
4907  int32_t* poly_ring_sizes,
4908  int64_t poly_num_rings,
4909  double* poly_bounds,
4910  int64_t poly_bounds_size,
4911  int32_t ic1,
4912  int32_t isr1,
4913  int32_t ic2,
4914  int32_t isr2,
4915  int32_t osr) {
4916  return ST_Intersects_Polygon_MultiPolygon(poly_coords,
4917  poly_coords_size,
4918  poly_ring_sizes,
4919  poly_num_rings,
4920  poly_bounds,
4921  poly_bounds_size,
4922  mpoly_coords,
4923  mpoly_coords_size,
4924  mpoly_ring_sizes,
4925  mpoly_num_rings,
4926  mpoly_poly_sizes,
4927  mpoly_num_polys,
4928  mpoly_bounds,
4929  mpoly_bounds_size,
4930  ic2,
4931  isr2,
4932  ic1,
4933  isr1,
4934  osr);
4935 }
4936
4938 bool ST_Intersects_MultiPolygon_MultiPolygon(int8_t* mpoly1_coords,
4939  int64_t mpoly1_coords_size,
4940  int32_t* mpoly1_ring_sizes,
4941  int64_t mpoly1_num_rings,
4942  int32_t* mpoly1_poly_sizes,
4943  int64_t mpoly1_num_polys,
4944  double* mpoly1_bounds,
4945  int64_t mpoly1_bounds_size,
4946  int8_t* mpoly2_coords,
4947  int64_t mpoly2_coords_size,
4948  int32_t* mpoly2_ring_sizes,
4949  int64_t mpoly2_num_rings,
4950  int32_t* mpoly2_poly_sizes,
4951  int64_t mpoly2_num_polys,
4952  double* mpoly2_bounds,
4953  int64_t mpoly2_bounds_size,
4954  int32_t ic1,
4955  int32_t isr1,
4956  int32_t ic2,
4957  int32_t isr2,
4958  int32_t osr) {
4959  if (mpoly1_bounds && mpoly2_bounds) {
4960  if (!box_overlaps_box(
4961  mpoly1_bounds, mpoly1_bounds_size, mpoly2_bounds, mpoly2_bounds_size)) {
4962  return false;
4963  }
4964  }
4965
4967  mpoly1_coords_size,
4968  mpoly1_ring_sizes,
4969  mpoly1_num_rings,
4970  mpoly1_poly_sizes,
4971  mpoly1_num_polys,
4972  mpoly2_coords,
4973  mpoly2_coords_size,
4974  mpoly2_ring_sizes,
4975  mpoly2_num_rings,
4976  mpoly2_poly_sizes,
4977  mpoly2_num_polys,
4978  ic1,
4979  isr1,
4980  ic2,
4981  isr2,
4982  osr,
4983  0.0));
4984 }
4985
4986 //
4987 // Accessors for poly bounds and render group for in-situ poly render queries
4988 //
4989
4991 int64_t HeavyDB_Geo_PolyBoundsPtr(double* bounds, int64_t size) {
4992  return reinterpret_cast<int64_t>(bounds);
4993 }
4994
4996 int32_t HeavyDB_Geo_PolyRenderGroup(int32_t render_group) {
4997  return render_group;
4998 }
4999
5000 // TODO Update for UTM. This assumes x and y are independent, which is not true for UTM.
5002 double convert_meters_to_pixel_width(const double meters,
5003  int8_t* p,
5004  const int64_t psize,
5005  const int32_t ic,
5006  const int32_t isr,
5007  const int32_t osr,
5008  const double min_lon,
5009  const double max_lon,
5010  const int32_t img_width,
5011  const double min_width) {
5012  const double const1 = 0.017453292519943295769236907684886;
5013  const double const2 = 6372797.560856;
5014  const auto lon = decompress_coord<X>(p, 0, ic);
5015  const auto lat = decompress_coord<Y>(p, 0, ic);
5016  double t1 = sinf(meters / (2.0 * const2));
5017  double t2 = cosf(const1 * lat);
5018  const double newlon = lon - (2.0 * asinf(t1 / t2)) / const1;
5019  t1 = transform_point<X>({lon, {}}, isr, osr).x;
5020  t2 = transform_point<X>({newlon, {}}, isr, osr).x;
5021  const double min_domain_x = transform_point<X>({min_lon, {}}, isr, osr).x;
5022  const double max_domain_x = transform_point<X>({max_lon, {}}, isr, osr).x;
5023  const double domain_diff = max_domain_x - min_domain_x;
5024  t1 = ((t1 - min_domain_x) / domain_diff) * static_cast<double>(img_width);
5025  t2 = ((t2 - min_domain_x) / domain_diff) * static_cast<double>(img_width);
5026
5027  // TODO(croot): need to account for edge cases, such as getting close to the poles.
5028  const double sz = fabs(t1 - t2);
5029  return (sz < min_width ? min_width : sz);
5030 }
5031
5032 // TODO Update for UTM. This assumes x and y are independent, which is not true for UTM.
5034 double convert_meters_to_pixel_height(const double meters,
5035  int8_t* p,
5036  const int64_t psize,
5037  const int32_t ic,
5038  const int32_t isr,
5039  const int32_t osr,
5040  const double min_lat,
5041  const double max_lat,
5042  const int32_t img_height,
5043  const double min_height) {
5044  const double const1 = 0.017453292519943295769236907684886;
5045  const double const2 = 6372797.560856;
5046  const auto lat = decompress_coord<Y>(p, 0, ic);
5047  const double latdiff = meters / (const1 * const2);
5048  const double newlat =
5049  (lat < 0) ? lat + latdiff : lat - latdiff; // assumes a lat range of [-90, 90]
5050  double t1 = transform_point<Y>({{}, lat}, isr, osr).y;
5051  double t2 = transform_point<Y>({{}, newlat}, isr, osr).y;
5052  const double min_domain_y = transform_point<Y>({{}, min_lat}, isr, osr).y;
5053  const double max_domain_y = transform_point<Y>({{}, max_lat}, isr, osr).y;
5054  const double domain_diff = max_domain_y - min_domain_y;
5055  t1 = ((t1 - min_domain_y) / domain_diff) * static_cast<double>(img_height);
5056  t2 = ((t2 - min_domain_y) / domain_diff) * static_cast<double>(img_height);
5057
5058  // TODO(croot): need to account for edge cases, such as getting close to the poles.
5059  const double sz = fabs(t1 - t2);
5060  return (sz < min_height ? min_height : sz);
5061 }
5062
5064  const int64_t psize,
5065  const int32_t ic,
5066  const double min_lon,
5067  const double max_lon,
5068  const double min_lat,
5069  const double max_lat) {
5070  const auto lon = decompress_coord<X>(p, 0, ic);
5071  const auto lat = decompress_coord<Y>(p, 0, ic);
5072  return !(lon < min_lon || lon > max_lon || lat < min_lat || lat > max_lat);
5073 }
5074
5076  const int64_t psize,
5077  const int32_t ic,
5078  const double meters,
5079  const double min_lon,
5080  const double max_lon,
5081  const double min_lat,
5082  const double max_lat) {
5083  const double const1 = 0.017453292519943295769236907684886;
5084  const double const2 = 6372797.560856;
5085  const auto lon = decompress_coord<X>(p, 0, ic);
5086  const auto lat = decompress_coord<Y>(p, 0, ic);
5087  const double latdiff = meters / (const1 * const2);
5088  const double t1 = sinf(meters / (2.0 * const2));
5089  const double t2 = cosf(const1 * lat);
5090  const double londiff = (2.0 * asinf(t1 / t2)) / const1;
5091  return !(lon + londiff < min_lon || lon - londiff > max_lon ||
5092  lat + latdiff < min_lat || lat - latdiff > max_lat);
5093 }
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)
EXTENSION_NOINLINE double ST_Distance_MultiPoint_MultiPoint_Squared(int8_t *mp1, int64_t mp1size, int8_t *mp2, int64_t mp2size, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr, double threshold)
EXTENSION_INLINE double ST_Distance_MultiPolygon_MultiPoint(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 *mp, int64_t mpsize, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr, double threshold)
DEVICE ALWAYS_INLINE double distance_point_multilinestring(int8_t *p, int64_t psize, int8_t *mls, int64_t mls_size, int32_t *mls_ls_sizes, int64_t mls_ls_num, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr, double threshold)
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_MultiPoint_Polygon(int8_t *mp, int64_t mpsize, 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_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 double ST_Distance_MultiPoint_MultiPolygon(int8_t *mp, int64_t mpsize, 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_INLINE double ST_Distance_MultiPoint_Point_Squared(int8_t *mp, int64_t mpsize, 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_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:52
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_NOINLINE double ST_Distance_MultiLineString_MultiLineString(int8_t *mls1, int64_t mls1size, int32_t *mls1_ls_sizes, int64_t mls1_ls_num, int8_t *mls2, int64_t mls2size, int32_t *mls2_ls_sizes, int64_t mls2_ls_num, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr, double threshold)
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 double ST_Distance_MultiPoint_MultiLineString(int8_t *mp, int64_t mpsize, int8_t *mls, int64_t mls_size, int32_t *mls_ls_sizes, int64_t mls_ls_num, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr, double threshold)
EXTENSION_INLINE bool ST_Intersects_Polygon_Point(int8_t *poly, int64_t polysize, int32_t *poly_ring_sizes, int64_t poly_num_rings, double *poly_bounds, int64_t poly_bounds_size, int8_t *p, int64_t psize, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr)
EXTENSION_NOINLINE double ST_Distance_Point_Polygon(int8_t *p, int64_t psize, int8_t *poly, int64_t polysize, int32_t *poly_ring_sizes, int64_t poly_num_rings, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr, double threshold)
EXTENSION_NOINLINE double ST_Distance_Point_Point_Squared(int8_t *p1, int64_t p1size, int8_t *p2, int64_t p2size, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr)
EXTENSION_INLINE bool ST_DWithin_Point_Point(int8_t *p1, int64_t p1size, int8_t *p2, int64_t p2size, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr, double distance_within)
EXTENSION_NOINLINE double ST_Y_Point(int8_t *p, int64_t psize, int32_t ic, int32_t isr, int32_t osr)
EXTENSION_NOINLINE double ST_Distance_Point_MultiPoint_Squared(int8_t *p, int64_t psize, int8_t *mp, int64_t mpsize, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr, double threshold)
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)
EXTENSION_INLINE double ST_Distance_Polygon_MultiPoint(int8_t *poly, int64_t polysize, int32_t *poly_ring_sizes, int64_t poly_num_rings, int8_t *mp, int64_t mpsize, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr, double threshold)
EXTENSION_NOINLINE void ST_Centroid_MultiPoint(int8_t *mp, int32_t mpsize, int32_t ic, int32_t isr, int32_t osr, double *multipoint_centroid)
#define DEVICE
EXTENSION_INLINE double ST_Distance_MultiLineString_Point(int8_t *mls, int64_t mls_size, int32_t *mls_ls_sizes, int64_t mls_ls_num, 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_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)
DEVICE ALWAYS_INLINE double 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_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:51
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)
EXTENSION_NOINLINE double ST_Distance_Point_MultiPoint(int8_t *p, int64_t psize, int8_t *mp, int64_t mpsize, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr, double threshold)
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 double ST_Length_MultiLineString(int8_t *coords, int64_t coords_sz, int8_t *linestring_sizes_in, int64_t linestring_sizes_sz, int32_t ic, int32_t isr, int32_t osr)
EXTENSION_NOINLINE bool ST_Intersects_Polygon_MultiPolygon(int8_t *poly_coords, int64_t poly_coords_size, int32_t *poly_ring_sizes, int64_t poly_num_rings, double *poly_bounds, int64_t poly_bounds_size, int8_t *mpoly_coords, int64_t mpoly_coords_size, int32_t *mpoly_ring_sizes, int64_t mpoly_num_rings, int32_t *mpoly_poly_sizes, int64_t mpoly_num_polys, double *mpoly_bounds, int64_t mpoly_bounds_size, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr)
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 double ST_Distance_LineString_MultiPoint(int8_t *l, int64_t lsize, int8_t *mp, int64_t mpsize, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr, double threshold)
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 double 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)
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_INLINE double ST_Distance_MultiPoint_LineString(int8_t *mp, int64_t mpsize, 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 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_NOINLINE double ST_Distance_Point_MultiLineString(int8_t *p, int64_t psize, int8_t *mls, int64_t mls_size, int32_t *mls_ls_sizes, int64_t mls_ls_num, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr, double threshold)
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)
EXTENSION_NOINLINE double ST_Distance_MultiPoint_MultiPoint(int8_t *mp1, int64_t mp1size, int8_t *mp2, int64_t mp2size, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr, double threshold)
EXTENSION_INLINE double ST_Distance_MultiLineString_MultiPoint(int8_t *mls, int64_t mls_size, int32_t *mls_ls_sizes, int64_t mls_ls_num, int8_t *mp, int64_t mpsize, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr, double threshold)
EXTENSION_INLINE double ST_Distance_MultiPoint_Point(int8_t *mp, int64_t mpsize, int8_t *p, int64_t psize, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr, double threshold)
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)