OmniSciDB  95562058bd
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
ExtensionFunctionsGeo.hpp
Go to the documentation of this file.
2 
3 // Adjustable tolerance, determined by compression mode.
4 // The criteria is to still recognize a compressed+decompressed number.
5 // For example 1.0 longitude compressed with GEOINT32 and then decompressed
6 // comes back as 0.99999994086101651, which is still within GEOINT32
7 // tolerance val 0.0000001
8 DEVICE ALWAYS_INLINE double tol(int32_t ic) {
9  if (ic == COMPRESSION_GEOINT32) {
10  return TOLERANCE_GEOINT32;
11  }
12  return TOLERANCE_DEFAULT;
13 }
14 
15 // Combine tolerances for two-component calculations
16 DEVICE ALWAYS_INLINE double tol(int32_t ic1, int32_t ic2) {
17  return fmax(tol(ic1), tol(ic2));
18 }
19 
20 DEVICE ALWAYS_INLINE bool tol_zero(double x, double tolerance = TOLERANCE_DEFAULT) {
21  return (-tolerance <= x) && (x <= tolerance);
22 }
23 
24 DEVICE ALWAYS_INLINE bool tol_eq(double x,
25  double y,
26  double tolerance = TOLERANCE_DEFAULT) {
27  auto diff = x - y;
28  return (-tolerance <= diff) && (diff <= tolerance);
29 }
30 
31 DEVICE ALWAYS_INLINE bool tol_le(double x,
32  double y,
33  double tolerance = TOLERANCE_DEFAULT) {
34  return x <= (y + tolerance);
35 }
36 
37 DEVICE ALWAYS_INLINE bool tol_ge(double x,
38  double y,
39  double tolerance = TOLERANCE_DEFAULT) {
40  return (x + tolerance) >= y;
41 }
42 
44  int32_t index,
45  int32_t ic,
46  bool x) {
47  if (ic == COMPRESSION_GEOINT32) {
48  auto compressed_coords = reinterpret_cast<int32_t*>(data);
49  auto compressed_coord = compressed_coords[index];
50  if (x) {
51  return Geospatial::decompress_longitude_coord_geoint32(compressed_coord);
52  } else {
53  return Geospatial::decompress_lattitude_coord_geoint32(compressed_coord);
54  }
55  }
56  auto double_coords = reinterpret_cast<double*>(data);
57  return double_coords[index];
58 }
59 
61  if (ic == COMPRESSION_GEOINT32) {
62  return 4;
63  }
64  return 8;
65 }
66 
67 DEVICE ALWAYS_INLINE double transform_coord(double coord,
68  int32_t isr,
69  int32_t osr,
70  bool x) {
71  if (isr == 4326) {
72  if (osr == 900913) {
73  // WGS 84 --> Web Mercator
74  if (x) {
75  return conv_4326_900913_x(coord);
76  } else {
77  return conv_4326_900913_y(coord);
78  }
79  }
80  }
81  return coord;
82 }
83 
84 // X coord accessor handling on-the-fly decommpression and transforms
85 DEVICE ALWAYS_INLINE double coord_x(int8_t* data,
86  int32_t index,
87  int32_t ic,
88  int32_t isr,
89  int32_t osr) {
90  auto decompressed_coord_x = decompress_coord(data, index, ic, true);
91  auto decompressed_transformed_coord_x =
92  transform_coord(decompressed_coord_x, isr, osr, true);
93  return decompressed_transformed_coord_x;
94 }
95 
96 // Y coord accessor handling on-the-fly decommpression and transforms
97 DEVICE ALWAYS_INLINE double coord_y(int8_t* data,
98  int32_t index,
99  int32_t ic,
100  int32_t isr,
101  int32_t osr) {
102  auto decompressed_coord_y = decompress_coord(data, index, ic, false);
103  auto decompressed_transformed_coord_y =
104  transform_coord(decompressed_coord_y, isr, osr, false);
105  return decompressed_transformed_coord_y;
106 }
107 
108 // Cartesian distance between points, squared
110  double p1y,
111  double p2x,
112  double p2y) {
113  auto x = p1x - p2x;
114  auto y = p1y - p2y;
115  auto x2 = x * x;
116  auto y2 = y * y;
117  auto d2 = x2 + y2;
119  return 0.0;
120  }
121  return d2;
122 }
123 
124 // Cartesian distance between points
126  double p1y,
127  double p2x,
128  double p2y) {
129  auto d2 = distance_point_point_squared(p1x, p1y, p2x, p2y);
130  return sqrt(d2);
131 }
132 
133 // Cartesian distance between a point and a line segment
134 DEVICE
135 double distance_point_line(double px,
136  double py,
137  double l1x,
138  double l1y,
139  double l2x,
140  double l2y) {
141  double length = distance_point_point(l1x, l1y, l2x, l2y);
142  if (tol_zero(length)) {
143  return distance_point_point(px, py, l1x, l1y);
144  }
145 
146  // Find projection of point P onto the line segment AB:
147  // Line containing that segment: A + k * (B - A)
148  // Projection of point P onto the line touches it at
149  // k = dot(P-A,B-A) / length^2
150  // AB segment is represented by k = [0,1]
151  // Clamping k to [0,1] will give the shortest distance from P to AB segment
152  double dotprod = (px - l1x) * (l2x - l1x) + (py - l1y) * (l2y - l1y);
153  double k = dotprod / (length * length);
154  k = fmax(0.0, fmin(1.0, k));
155  double projx = l1x + k * (l2x - l1x);
156  double projy = l1y + k * (l2y - l1y);
157  return distance_point_point(px, py, projx, projy);
158 }
159 
160 // Given three colinear points p, q, r, the function checks if
161 // point q lies on line segment 'pr'
163  double py,
164  double qx,
165  double qy,
166  double rx,
167  double ry) {
168  return (tol_le(qx, fmax(px, rx)) && tol_ge(qx, fmin(px, rx)) &&
169  tol_le(qy, fmax(py, ry)) && tol_ge(qy, fmin(py, ry)));
170 }
171 
172 DEVICE ALWAYS_INLINE int16_t
173 orientation(double px, double py, double qx, double qy, double rx, double ry) {
174  auto val = ((qy - py) * (rx - qx) - (qx - px) * (ry - qy));
175  if (tol_zero(val)) {
176  return 0; // Points p, q and r are colinear
177  }
178  if (val > 0.0) {
179  return 1; // Clockwise point orientation
180  }
181  return 2; // Counterclockwise point orientation
182 }
183 
184 // Cartesian intersection of two line segments l11-l12 and l21-l22
185 DEVICE
186 bool line_intersects_line(double l11x,
187  double l11y,
188  double l12x,
189  double l12y,
190  double l21x,
191  double l21y,
192  double l22x,
193  double l22y) {
194  auto o1 = orientation(l11x, l11y, l12x, l12y, l21x, l21y);
195  auto o2 = orientation(l11x, l11y, l12x, l12y, l22x, l22y);
196  auto o3 = orientation(l21x, l21y, l22x, l22y, l11x, l11y);
197  auto o4 = orientation(l21x, l21y, l22x, l22y, l12x, l12y);
198 
199  // General case
200  if (o1 != o2 && o3 != o4) {
201  return true;
202  }
203 
204  // Special Cases
205  // l11, l12 and l21 are colinear and l21 lies on segment l11-l12
206  if (o1 == 0 && on_segment(l11x, l11y, l21x, l21y, l12x, l12y)) {
207  return true;
208  }
209 
210  // l11, l12 and l21 are colinear and l22 lies on segment l11-l12
211  if (o2 == 0 && on_segment(l11x, l11y, l22x, l22y, l12x, l12y)) {
212  return true;
213  }
214 
215  // l21, l22 and l11 are colinear and l11 lies on segment l21-l22
216  if (o3 == 0 && on_segment(l21x, l21y, l11x, l11y, l22x, l22y)) {
217  return true;
218  }
219 
220  // l21, l22 and l12 are colinear and l12 lies on segment l21-l22
221  if (o4 == 0 && on_segment(l21x, l21y, l12x, l12y, l22x, l22y)) {
222  return true;
223  }
224 
225  return false;
226 }
227 
228 DEVICE
230  int32_t lnum_coords,
231  double l1x,
232  double l1y,
233  double l2x,
234  double l2y,
235  int32_t ic1,
236  int32_t isr1,
237  int32_t osr) {
238  double e1x = coord_x(l, 0, ic1, isr1, osr);
239  double e1y = coord_y(l, 1, ic1, isr1, osr);
240  for (int64_t i = 2; i < lnum_coords; i += 2) {
241  double e2x = coord_x(l, i, ic1, isr1, osr);
242  double e2y = coord_y(l, i + 1, ic1, isr1, osr);
243  if (line_intersects_line(e1x, e1y, e2x, e2y, l1x, l1y, l2x, l2y)) {
244  return true;
245  }
246  e1x = e2x;
247  e1y = e2y;
248  }
249  return false;
250 }
251 
252 DEVICE
253 bool ring_intersects_line(int8_t* ring,
254  int32_t ring_num_coords,
255  double l1x,
256  double l1y,
257  double l2x,
258  double l2y,
259  int32_t ic1,
260  int32_t isr1,
261  int32_t osr) {
262  double e1x = coord_x(ring, ring_num_coords - 2, ic1, isr1, osr);
263  double e1y = coord_y(ring, ring_num_coords - 1, ic1, isr1, osr);
264  double e2x = coord_x(ring, 0, ic1, isr1, osr);
265  double e2y = coord_y(ring, 1, ic1, isr1, osr);
266  if (line_intersects_line(e1x, e1y, e2x, e2y, l1x, l1y, l2x, l2y)) {
267  return true;
268  }
270  ring, ring_num_coords, l1x, l1y, l2x, l2y, ic1, isr1, osr);
271 }
272 
273 DEVICE
275  int32_t lnum_coords,
276  double l1x,
277  double l1y,
278  double l2x,
279  double l2y,
280  int32_t ic1,
281  int32_t isr1,
282  int32_t osr) {
283  double e1x = coord_x(l, 0, ic1, isr1, osr);
284  double e1y = coord_y(l, 1, ic1, isr1, osr);
285  for (int64_t i = 2; i < lnum_coords; i += 2) {
286  double e2x = coord_x(l, i, ic1, isr1, osr);
287  double e2y = coord_y(l, i + 1, ic1, isr1, osr);
288  if (line_intersects_line(e1x, e1y, e2x, e2y, l1x, l1y, l2x, l2y)) {
289  return true;
290  }
291  e1x = e2x;
292  e1y = e2y;
293  }
294  return false;
295 }
296 
297 // Cartesian distance between two line segments l11-l12 and l21-l22
298 DEVICE
299 double distance_line_line(double l11x,
300  double l11y,
301  double l12x,
302  double l12y,
303  double l21x,
304  double l21y,
305  double l22x,
306  double l22y) {
307  if (line_intersects_line(l11x, l11y, l12x, l12y, l21x, l21y, l22x, l22y)) {
308  return 0.0;
309  }
310  double dist12 = fmin(distance_point_line(l11x, l11y, l21x, l21y, l22x, l22y),
311  distance_point_line(l12x, l12y, l21x, l21y, l22x, l22y));
312  double dist21 = fmin(distance_point_line(l21x, l21y, l11x, l11y, l12x, l12y),
313  distance_point_line(l22x, l22y, l11x, l11y, l12x, l12y));
314  return fmin(dist12, dist21);
315 }
316 
317 DEVICE
318 double distance_ring_linestring(int8_t* ring,
319  int32_t ring_num_coords,
320  int8_t* l,
321  int32_t lnum_coords,
322  int32_t ic1,
323  int32_t isr1,
324  int32_t ic2,
325  int32_t isr2,
326  int32_t osr) {
327  double min_distance = 0.0;
328 
329  double re1x = coord_x(ring, ring_num_coords - 2, ic1, isr1, osr);
330  double re1y = coord_y(ring, ring_num_coords - 1, ic1, isr1, osr);
331  for (auto i = 0; i < ring_num_coords; i += 2) {
332  double re2x = coord_x(ring, i, ic1, isr1, osr);
333  double re2y = coord_y(ring, i + 1, ic1, isr1, osr);
334 
335  double le1x = coord_x(l, 0, ic2, isr2, osr);
336  double le1y = coord_y(l, 1, ic2, isr2, osr);
337  for (auto j = 2; j < lnum_coords; j += 2) {
338  double le2x = coord_x(l, j, ic2, isr2, osr);
339  double le2y = coord_y(l, j + 1, ic2, isr2, osr);
340 
341  auto distance = distance_line_line(re1x, re1y, re2x, re2y, le1x, le1y, le2x, le2y);
342  if ((i == 0 && j == 2) || min_distance > distance) {
343  min_distance = distance;
344  if (tol_zero(min_distance)) {
345  return 0.0;
346  }
347  }
348  le1x = le2x;
349  le1y = le2y;
350  }
351  re1x = re2x;
352  re1y = re2y;
353  }
354 
355  return min_distance;
356 }
357 
358 DEVICE
359 double distance_ring_ring(int8_t* ring1,
360  int32_t ring1_num_coords,
361  int8_t* ring2,
362  int32_t ring2_num_coords,
363  int32_t ic1,
364  int32_t isr1,
365  int32_t ic2,
366  int32_t isr2,
367  int32_t osr) {
368  double min_distance = 0.0;
369 
370  double e11x = coord_x(ring1, ring1_num_coords - 2, ic1, isr1, osr);
371  double e11y = coord_y(ring1, ring1_num_coords - 1, ic1, isr1, osr);
372  for (auto i = 0; i < ring1_num_coords; i += 2) {
373  double e12x = coord_x(ring1, i, ic1, isr1, osr);
374  double e12y = coord_y(ring1, i + 1, ic1, isr1, osr);
375 
376  double e21x = coord_x(ring2, ring2_num_coords - 2, ic2, isr2, osr);
377  double e21y = coord_y(ring2, ring2_num_coords - 1, ic2, isr2, osr);
378  for (auto j = 0; j < ring2_num_coords; j += 2) {
379  double e22x = coord_x(ring2, j, ic2, isr2, osr);
380  double e22y = coord_y(ring2, j + 1, ic2, isr2, osr);
381 
382  auto distance = distance_line_line(e11x, e11y, e12x, e12y, e21x, e21y, e22x, e22y);
383  if ((i == 0 && j == 0) || min_distance > distance) {
384  min_distance = distance;
385  if (tol_zero(min_distance)) {
386  return 0.0;
387  }
388  }
389  e21x = e22x;
390  e21y = e22y;
391  }
392  e11x = e12x;
393  e11y = e12y;
394  }
395 
396  return min_distance;
397 }
398 
399 // Checks if a simple polygon (no holes) contains a point.
400 //
401 // Poly coords are extracted from raw data, based on compression (ic1) and input/output
402 // SRIDs (isr1/osr).
403 //
404 // Shoot a ray from point P to the right, register intersections with any of polygon's
405 // edges. Each intersection means entrance into or exit from the polygon. Odd number of
406 // intersections means the polygon does contain P. Account for special cases: touch+cross,
407 // touch+leave, touch+overlay+cross, touch+overlay+leave, on edge, etc.
408 //
409 // Secondary ray is shot from point P down for simple redundancy, to reduce main probe's
410 // chance of error. No intersections means P is outside, irrespective of main probe's
411 // result.
412 //
413 DEVICE
414 bool polygon_contains_point(int8_t* poly,
415  int32_t poly_num_coords,
416  double px,
417  double py,
418  int32_t ic1,
419  int32_t isr1,
420  int32_t osr) {
421  bool result = false;
422  int xray_touch = 0;
423  bool horizontal_edge = false;
424  bool yray_intersects = false;
425 
426  double e1x = coord_x(poly, poly_num_coords - 2, ic1, isr1, osr);
427  double e1y = coord_y(poly, poly_num_coords - 1, ic1, isr1, osr);
428  for (int64_t i = 0; i < poly_num_coords; i += 2) {
429  double e2x = coord_x(poly, i, ic1, isr1, osr);
430  double e2y = coord_y(poly, i + 1, ic1, isr1, osr);
431 
432  // Check if point sits on an edge.
433  if (tol_zero(distance_point_line(px, py, e1x, e1y, e2x, e2y))) {
434  return true;
435  }
436 
437  // Before flipping the switch, check if xray hit a horizontal edge
438  // - If an edge lays on the xray, one of the previous edges touched it
439  // so while moving horizontally we're in 'xray_touch' state
440  // - Last edge that touched xray at (e2x,e2y) didn't register intersection
441  // - Next edge that diverges from xray at (e1,e1y) will register intersection
442  // - Can have several horizontal edges, one after the other, keep moving though
443  // in 'xray_touch' state without flipping the switch
444  horizontal_edge = (xray_touch != 0) && tol_eq(py, e1y) && tol_eq(py, e2y);
445 
446  // Main probe: xray
447  // Overshoot the xray to detect an intersection if there is one.
448  double xray = fmax(e2x, e1x) + 1.0;
449  if (px <= xray && // Only check for intersection if the edge is on the right
450  !horizontal_edge && // Keep moving through horizontal edges
451  line_intersects_line(px, // xray shooting from point p to the right
452  py,
453  xray,
454  py,
455  e1x, // polygon edge
456  e1y,
457  e2x,
458  e2y)) {
459  // Register intersection
460  result = !result;
461 
462  // Adjust for special cases
463  if (xray_touch == 0) {
464  if (tol_zero(distance_point_line(e2x, e2y, px, py, xray + 1.0, py))) {
465  // Xray goes through the edge's second vertex, unregister intersection -
466  // that vertex will be crossed again when we look at the following edge(s)
467  result = !result;
468  // Enter the xray-touch state:
469  // (1) - xray was touched by the edge from above, (-1) from below
470  xray_touch = (e1y > py) ? 1 : -1;
471  }
472  } else {
473  // Previous edge touched the xray, intersection hasn't been registered,
474  // it has to be registered now if this edge continues across the xray.
475  if (xray_touch > 0) {
476  // Previous edge touched the xray from above
477  if (e2y <= py) {
478  // Current edge crosses under xray: intersection is already registered
479  } else {
480  // Current edge just touched the xray and pulled up: unregister intersection
481  result = !result;
482  }
483  } else {
484  // Previous edge touched the xray from below
485  if (e2y > py) {
486  // Current edge crosses over xray: intersection is already registered
487  } else {
488  // Current edge just touched the xray and pulled down: unregister intersection
489  result = !result;
490  }
491  }
492  // Exit the xray-touch state
493  xray_touch = 0;
494  }
495  }
496 
497  // Redundancy: vertical yray down
498  // Main probe xray may hit multiple complex fragments which increases a chance of
499  // error. Perform a simple secondary check for edge intersections to see if point is
500  // outside.
501  if (!yray_intersects) { // Continue checking on yray until intersection is found
502  double yray = fmin(e2y, e1y) - 1.0;
503  if (yray <= py) { // Only check for yray intersection if point P is above the edge
504  yray_intersects = line_intersects_line(px, // yray shooting from point P down
505  py,
506  px,
507  yray,
508  e1x, // polygon edge
509  e1y,
510  e2x,
511  e2y);
512  }
513  }
514 
515  // Advance to the next vertex
516  e1x = e2x;
517  e1y = e2y;
518  }
519  if (!yray_intersects) {
520  // yray has zero intersections - point is outside the polygon
521  return false;
522  }
523  // Otherwise rely on the main probe
524  return result;
525 }
526 
527 // Returns true if simple polygon (no holes) contains a linestring
528 DEVICE
529 bool polygon_contains_linestring(int8_t* poly,
530  int32_t poly_num_coords,
531  int8_t* l,
532  int64_t lnum_coords,
533  int32_t ic1,
534  int32_t isr1,
535  int32_t ic2,
536  int32_t isr2,
537  int32_t osr) {
538  // Check that the first point is in the polygon
539  double l1x = coord_x(l, 0, ic2, isr2, osr);
540  double l1y = coord_y(l, 1, ic2, isr2, osr);
541  if (!polygon_contains_point(poly, poly_num_coords, l1x, l1y, ic1, isr1, osr)) {
542  return false;
543  }
544 
545  // Go through line segments and check if there are no intersections with poly edges,
546  // i.e. linestring doesn't escape
547  for (int32_t i = 2; i < lnum_coords; i += 2) {
548  double l2x = coord_x(l, i, ic2, isr2, osr);
549  double l2y = coord_y(l, i + 1, ic2, isr2, osr);
550  if (ring_intersects_line(poly, poly_num_coords, l1x, l1y, l2x, l2y, ic1, isr1, osr)) {
551  return false;
552  }
553  l1x = l2x;
554  l1y = l2y;
555  }
556  return true;
557 }
558 
560  int64_t bounds_size,
561  double px,
562  double py) {
563  return (tol_ge(px, bounds[0]) && tol_ge(py, bounds[1]) && tol_le(px, bounds[2]) &&
564  tol_le(py, bounds[3]));
565 }
566 
568  int64_t bounds_size,
569  double px,
570  double py) {
571  return box_contains_point(bounds, bounds_size, px, py);
572 }
573 
574 DEVICE ALWAYS_INLINE bool box_contains_box(double* bounds1,
575  int64_t bounds1_size,
576  double* bounds2,
577  int64_t bounds2_size) {
578  return (
580  bounds1, bounds1_size, bounds2[0], bounds2[1]) && // box1 <- box2: xmin, ymin
582  bounds1, bounds1_size, bounds2[2], bounds2[3])); // box1 <- box2: xmax, ymax
583 }
584 
586  int64_t bounds1_size,
587  double* bounds2,
588  int64_t bounds2_size) {
589  return (
591  bounds1, bounds1_size, bounds2[0], bounds2[1]) || // box1 <- box2: xmin, ymin
593  bounds1, bounds1_size, bounds2[2], bounds2[3]) || // box1 <- box2: xmax, ymax
595  bounds1, bounds1_size, bounds2[0], bounds2[3]) || // box1 <- box2: xmin, ymax
597  bounds1, bounds1_size, bounds2[2], bounds2[1])); // box1 <- box2: xmax, ymin
598 }
599 
600 DEVICE ALWAYS_INLINE bool box_overlaps_box(double* bounds1,
601  int64_t bounds1_size,
602  double* bounds2,
603  int64_t bounds2_size) {
604  // TODO: tolerance
605  if (bounds1[2] < bounds2[0] || // box1 is left of box2: box1.xmax < box2.xmin
606  bounds1[0] > bounds2[2] || // box1 is right of box2: box1.xmin > box2.xmax
607  bounds1[3] < bounds2[1] || // box1 is below box2: box1.ymax < box2.miny
608  bounds1[1] > bounds2[3]) { // box1 is above box2: box1.ymin > box1.ymax
609  return false;
610  }
611  return true;
612 }
613 
615 double ST_X_Point(int8_t* p, int64_t psize, int32_t ic, int32_t isr, int32_t osr) {
616  return coord_x(p, 0, ic, isr, osr);
617 }
618 
620 double ST_Y_Point(int8_t* p, int64_t psize, int32_t ic, int32_t isr, int32_t osr) {
621  return coord_y(p, 1, ic, isr, osr);
622 }
623 
625 double ST_X_LineString(int8_t* l,
626  int64_t lsize,
627  int32_t lindex,
628  int32_t ic,
629  int32_t isr,
630  int32_t osr) {
631  auto l_num_points = lsize / (2 * compression_unit_size(ic));
632  if (lindex < 0 || lindex > l_num_points) {
633  lindex = l_num_points; // Endpoint
634  }
635  return coord_x(l, 2 * (lindex - 1), ic, isr, osr);
636 }
637 
639 double ST_Y_LineString(int8_t* l,
640  int64_t lsize,
641  int32_t lindex,
642  int32_t ic,
643  int32_t isr,
644  int32_t osr) {
645  auto l_num_points = lsize / (2 * compression_unit_size(ic));
646  if (lindex < 0 || lindex > l_num_points) {
647  lindex = l_num_points; // Endpoint
648  }
649  return coord_y(l, 2 * (lindex - 1) + 1, ic, isr, osr);
650 }
651 
653 double ST_XMin(int8_t* coords, int64_t size, int32_t ic, int32_t isr, int32_t osr) {
654  auto num_coords = size / compression_unit_size(ic);
655  double xmin = 0.0;
656  for (int32_t i = 0; i < num_coords; i += 2) {
657  double x = coord_x(coords, i, ic, isr, osr);
658  if (i == 0 || x < xmin) {
659  xmin = x;
660  }
661  }
662  return xmin;
663 }
664 
666 double ST_YMin(int8_t* coords, int64_t size, int32_t ic, int32_t isr, int32_t osr) {
667  auto num_coords = size / compression_unit_size(ic);
668  double ymin = 0.0;
669  for (int32_t i = 1; i < num_coords; i += 2) {
670  double y = coord_y(coords, i, ic, isr, osr);
671  if (i == 1 || y < ymin) {
672  ymin = y;
673  }
674  }
675  return ymin;
676 }
677 
679 double ST_XMax(int8_t* coords, int64_t size, int32_t ic, int32_t isr, int32_t osr) {
680  auto num_coords = size / compression_unit_size(ic);
681  double xmax = 0.0;
682  for (int32_t i = 0; i < num_coords; i += 2) {
683  double x = coord_x(coords, i, ic, isr, osr);
684  if (i == 0 || x > xmax) {
685  xmax = x;
686  }
687  }
688  return xmax;
689 }
690 
692 double ST_YMax(int8_t* coords, int64_t size, int32_t ic, int32_t isr, int32_t osr) {
693  auto num_coords = size / compression_unit_size(ic);
694  double ymax = 0.0;
695  for (int32_t i = 1; i < num_coords; i += 2) {
696  double y = coord_y(coords, i, ic, isr, osr);
697  if (i == 1 || y > ymax) {
698  ymax = y;
699  }
700  }
701  return ymax;
702 }
703 
705 double ST_XMin_Bounds(double* bounds, int64_t size, int32_t isr, int32_t osr) {
706  return transform_coord(bounds[0], isr, osr, true);
707 }
708 
710 double ST_YMin_Bounds(double* bounds, int64_t size, int32_t isr, int32_t osr) {
711  return transform_coord(bounds[1], isr, osr, false);
712 }
713 
715 double ST_XMax_Bounds(double* bounds, int64_t size, int32_t isr, int32_t osr) {
716  return transform_coord(bounds[2], isr, osr, true);
717 }
718 
720 double ST_YMax_Bounds(double* bounds, int64_t size, int32_t isr, int32_t osr) {
721  return transform_coord(bounds[3], isr, osr, false);
722 }
723 
724 //
725 // ST_Length
726 //
727 
729  int64_t lsize,
730  int32_t ic,
731  int32_t isr,
732  int32_t osr,
733  bool geodesic,
734  bool check_closed) {
735  auto l_num_coords = lsize / compression_unit_size(ic);
736 
737  double length = 0.0;
738 
739  double l0x = coord_x(l, 0, ic, isr, osr);
740  double l0y = coord_y(l, 1, ic, isr, osr);
741  double l2x = l0x;
742  double l2y = l0y;
743  for (int32_t i = 2; i < l_num_coords; i += 2) {
744  double l1x = l2x;
745  double l1y = l2y;
746  l2x = coord_x(l, i, ic, isr, osr);
747  l2y = coord_y(l, i + 1, ic, isr, osr);
748  double ldist = geodesic ? distance_in_meters(l1x, l1y, l2x, l2y)
749  : distance_point_point(l1x, l1y, l2x, l2y);
750  length += ldist;
751  }
752  if (check_closed) {
753  double ldist = geodesic ? distance_in_meters(l2x, l2y, l0x, l0y)
754  : distance_point_point(l2x, l2y, l0x, l0y);
755  length += ldist;
756  }
757  return length;
758 }
759 
761 double ST_Length_LineString(int8_t* coords,
762  int64_t coords_sz,
763  int32_t ic,
764  int32_t isr,
765  int32_t osr) {
766  return length_linestring(coords, coords_sz, ic, isr, osr, false, false);
767 }
768 
770 double ST_Length_LineString_Geodesic(int8_t* coords,
771  int64_t coords_sz,
772  int32_t ic,
773  int32_t isr,
774  int32_t osr) {
775  return length_linestring(coords, coords_sz, ic, isr, osr, true, false);
776 }
777 
778 //
779 // ST_Perimeter
780 //
781 
783 double ST_Perimeter_Polygon(int8_t* poly,
784  int64_t polysize,
785  int32_t* poly_ring_sizes,
786  int64_t poly_num_rings,
787  int32_t ic,
788  int32_t isr,
789  int32_t osr) {
790  if (poly_num_rings <= 0) {
791  return 0.0;
792  }
793 
794  auto exterior_ring_num_coords = poly_ring_sizes[0] * 2;
795  auto exterior_ring_coords_size = exterior_ring_num_coords * compression_unit_size(ic);
796 
797  return length_linestring(poly, exterior_ring_coords_size, ic, isr, osr, false, true);
798 }
799 
801 double ST_Perimeter_Polygon_Geodesic(int8_t* poly,
802  int64_t polysize,
803  int32_t* poly_ring_sizes,
804  int64_t poly_num_rings,
805  int32_t ic,
806  int32_t isr,
807  int32_t osr) {
808  if (poly_num_rings <= 0) {
809  return 0.0;
810  }
811 
812  auto exterior_ring_num_coords = poly_ring_sizes[0] * 2;
813  auto exterior_ring_coords_size = exterior_ring_num_coords * compression_unit_size(ic);
814 
815  return length_linestring(poly, exterior_ring_coords_size, ic, isr, osr, true, true);
816 }
817 
818 DEVICE ALWAYS_INLINE double perimeter_multipolygon(int8_t* mpoly_coords,
819  int64_t mpoly_coords_size,
820  int32_t* mpoly_ring_sizes,
821  int64_t mpoly_num_rings,
822  int32_t* mpoly_poly_sizes,
823  int64_t mpoly_num_polys,
824  int32_t ic,
825  int32_t isr,
826  int32_t osr,
827  bool geodesic) {
828  if (mpoly_num_polys <= 0 || mpoly_num_rings <= 0) {
829  return 0.0;
830  }
831 
832  double perimeter = 0.0;
833 
834  // Set specific poly pointers as we move through the coords/ringsizes/polyrings arrays.
835  auto next_poly_coords = mpoly_coords;
836  auto next_poly_ring_sizes = mpoly_ring_sizes;
837 
838  for (auto poly = 0; poly < mpoly_num_polys; poly++) {
839  auto poly_coords = next_poly_coords;
840  auto poly_ring_sizes = next_poly_ring_sizes;
841  auto poly_num_rings = mpoly_poly_sizes[poly];
842  // Count number of coords in all of poly's rings, advance ring size pointer.
843  int32_t poly_num_coords = 0;
844  for (auto ring = 0; ring < poly_num_rings; ring++) {
845  poly_num_coords += 2 * *next_poly_ring_sizes++;
846  }
847  auto poly_coords_size = poly_num_coords * compression_unit_size(ic);
848  next_poly_coords += poly_coords_size;
849 
850  auto exterior_ring_num_coords = poly_ring_sizes[0] * 2;
851  auto exterior_ring_coords_size = exterior_ring_num_coords * compression_unit_size(ic);
852 
853  perimeter += length_linestring(
854  poly_coords, exterior_ring_coords_size, ic, isr, osr, geodesic, true);
855  }
856 
857  return perimeter;
858 }
859 
861 double ST_Perimeter_MultiPolygon(int8_t* mpoly_coords,
862  int64_t mpoly_coords_size,
863  int32_t* mpoly_ring_sizes,
864  int64_t mpoly_num_rings,
865  int32_t* mpoly_poly_sizes,
866  int64_t mpoly_num_polys,
867  int32_t ic,
868  int32_t isr,
869  int32_t osr) {
870  return perimeter_multipolygon(mpoly_coords,
871  mpoly_coords_size,
872  mpoly_ring_sizes,
873  mpoly_num_rings,
874  mpoly_poly_sizes,
875  mpoly_num_polys,
876  ic,
877  isr,
878  osr,
879  false);
880 }
881 
883 double ST_Perimeter_MultiPolygon_Geodesic(int8_t* mpoly_coords,
884  int64_t mpoly_coords_size,
885  int32_t* mpoly_ring_sizes,
886  int64_t mpoly_num_rings,
887  int32_t* mpoly_poly_sizes,
888  int64_t mpoly_num_polys,
889  int32_t ic,
890  int32_t isr,
891  int32_t osr) {
892  return perimeter_multipolygon(mpoly_coords,
893  mpoly_coords_size,
894  mpoly_ring_sizes,
895  mpoly_num_rings,
896  mpoly_poly_sizes,
897  mpoly_num_polys,
898  ic,
899  isr,
900  osr,
901  true);
902 }
903 
904 //
905 // ST_Area
906 //
907 
909  double y1,
910  double x2,
911  double y2,
912  double x3,
913  double y3) {
914  return (x1 * y2 - x2 * y1 + x3 * y1 - x1 * y3 + x2 * y3 - x3 * y2) / 2.0;
915 }
916 
917 DEVICE ALWAYS_INLINE double area_ring(int8_t* ring,
918  int64_t ringsize,
919  int32_t ic,
920  int32_t isr,
921  int32_t osr) {
922  auto ring_num_coords = ringsize / compression_unit_size(ic);
923 
924  if (ring_num_coords < 6) {
925  return 0.0;
926  }
927 
928  double area = 0.0;
929 
930  double x1 = coord_x(ring, 0, ic, isr, osr);
931  double y1 = coord_y(ring, 1, ic, isr, osr);
932  double x2 = coord_x(ring, 2, ic, isr, osr);
933  double y2 = coord_y(ring, 3, ic, isr, osr);
934  for (int32_t i = 4; i < ring_num_coords; i += 2) {
935  double x3 = coord_x(ring, i, ic, isr, osr);
936  double y3 = coord_y(ring, i + 1, ic, isr, osr);
937  area += area_triangle(x1, y1, x2, y2, x3, y3);
938  x2 = x3;
939  y2 = y3;
940  }
941  return area;
942 }
943 
944 DEVICE ALWAYS_INLINE double area_polygon(int8_t* poly_coords,
945  int64_t poly_coords_size,
946  int32_t* poly_ring_sizes,
947  int64_t poly_num_rings,
948  int32_t ic,
949  int32_t isr,
950  int32_t osr) {
951  if (poly_num_rings <= 0) {
952  return 0.0;
953  }
954 
955  double area = 0.0;
956  auto ring_coords = poly_coords;
957 
958  // Add up the areas of all rings.
959  // External ring is CCW, open - positive area.
960  // Internal rings (holes) are CW, open - negative areas.
961  for (auto r = 0; r < poly_num_rings; r++) {
962  auto ring_coords_size = poly_ring_sizes[r] * 2 * compression_unit_size(ic);
963  area += area_ring(ring_coords, ring_coords_size, ic, isr, osr);
964  // Advance to the next ring.
965  ring_coords += ring_coords_size;
966  }
967  return area;
968 }
969 
971 double ST_Area_Polygon(int8_t* poly_coords,
972  int64_t poly_coords_size,
973  int32_t* poly_ring_sizes,
974  int64_t poly_num_rings,
975  int32_t ic,
976  int32_t isr,
977  int32_t osr) {
978  return area_polygon(
979  poly_coords, poly_coords_size, poly_ring_sizes, poly_num_rings, ic, isr, osr);
980 }
981 
983 double ST_Area_Polygon_Geodesic(int8_t* poly_coords,
984  int64_t poly_coords_size,
985  int32_t* poly_ring_sizes,
986  int64_t poly_num_rings,
987  int32_t ic,
988  int32_t isr,
989  int32_t osr) {
990  return ST_Area_Polygon(
991  poly_coords, poly_coords_size, poly_ring_sizes, poly_num_rings, ic, isr, osr);
992 }
993 
995 double ST_Area_MultiPolygon(int8_t* mpoly_coords,
996  int64_t mpoly_coords_size,
997  int32_t* mpoly_ring_sizes,
998  int64_t mpoly_num_rings,
999  int32_t* mpoly_poly_sizes,
1000  int64_t mpoly_num_polys,
1001  int32_t ic,
1002  int32_t isr,
1003  int32_t osr) {
1004  if (mpoly_num_rings <= 0 || mpoly_num_polys <= 0) {
1005  return 0.0;
1006  }
1007 
1008  double area = 0.0;
1009 
1010  // Set specific poly pointers as we move through the coords/ringsizes/polyrings arrays.
1011  auto next_poly_coords = mpoly_coords;
1012  auto next_poly_ring_sizes = mpoly_ring_sizes;
1013 
1014  for (auto poly = 0; poly < mpoly_num_polys; poly++) {
1015  auto poly_coords = next_poly_coords;
1016  auto poly_ring_sizes = next_poly_ring_sizes;
1017  auto poly_num_rings = mpoly_poly_sizes[poly];
1018  // Count number of coords in all of poly's rings, advance ring size pointer.
1019  int32_t poly_num_coords = 0;
1020  for (auto ring = 0; ring < poly_num_rings; ring++) {
1021  poly_num_coords += 2 * *next_poly_ring_sizes++;
1022  }
1023  auto poly_coords_size = poly_num_coords * compression_unit_size(ic);
1024  next_poly_coords += poly_coords_size;
1025 
1026  area += area_polygon(
1027  poly_coords, poly_coords_size, poly_ring_sizes, poly_num_rings, ic, isr, osr);
1028  }
1029  return area;
1030 }
1031 
1033 double ST_Area_MultiPolygon_Geodesic(int8_t* mpoly_coords,
1034  int64_t mpoly_coords_size,
1035  int32_t* mpoly_ring_sizes,
1036  int64_t mpoly_num_rings,
1037  int32_t* mpoly_poly_sizes,
1038  int64_t mpoly_num_polys,
1039  int32_t ic,
1040  int32_t isr,
1041  int32_t osr) {
1042  return ST_Area_MultiPolygon(mpoly_coords,
1043  mpoly_coords_size,
1044  mpoly_ring_sizes,
1045  mpoly_num_rings,
1046  mpoly_poly_sizes,
1047  mpoly_num_polys,
1048  ic,
1049  isr,
1050  osr);
1051 }
1052 
1053 //
1054 // ST_Centroid
1055 //
1056 
1057 // Point centroid
1058 // - calculate average coordinates for all points
1059 // LineString centroid
1060 // - take midpoint of each line segment
1061 // - calculate average coordinates of all midpoints weighted by segment length
1062 // Polygon, MultiPolygon centroid, holes..
1063 // - weighted sum of centroids of triangles coming out of area decomposition
1064 //
1065 // On zero area fall back to linestring centroid
1066 // On zero length fall back to point centroid
1067 
1069 double ST_Centroid_Point(int8_t* p,
1070  int64_t psize,
1071  int32_t ic,
1072  int32_t isr,
1073  int32_t osr,
1074  bool ycoord) {
1075  if (ycoord) {
1076  return coord_y(p, 1, ic, isr, osr);
1077  }
1078  return coord_x(p, 0, ic, isr, osr);
1079 }
1080 
1082  double y1,
1083  double x2,
1084  double y2,
1085  double* length,
1086  double* linestring_centroid_sum) {
1087  double ldist = distance_point_point(x1, y1, x2, y2);
1088  *length += ldist;
1089  double segment_midpoint_x = (x1 + x2) / 2.0;
1090  double segment_midpoint_y = (y1 + y2) / 2.0;
1091  linestring_centroid_sum[0] += ldist * segment_midpoint_x;
1092  linestring_centroid_sum[1] += ldist * segment_midpoint_y;
1093  return true;
1094 }
1095 
1097  int64_t lsize,
1098  int32_t ic,
1099  int32_t isr,
1100  int32_t osr,
1101  bool closed,
1102  double* total_length,
1103  double* linestring_centroid_sum,
1104  int64_t* num_points,
1105  double* point_centroid_sum) {
1106  auto l_num_coords = lsize / compression_unit_size(ic);
1107  double length = 0.0;
1108  double l0x = coord_x(l, 0, ic, isr, osr);
1109  double l0y = coord_y(l, 1, ic, isr, osr);
1110  double l2x = l0x;
1111  double l2y = l0y;
1112  for (int32_t i = 2; i < l_num_coords; i += 2) {
1113  double l1x = l2x;
1114  double l1y = l2y;
1115  l2x = coord_x(l, i, ic, isr, osr);
1116  l2y = coord_y(l, i + 1, ic, isr, osr);
1117  centroid_add_segment(l1x, l1y, l2x, l2y, &length, linestring_centroid_sum);
1118  }
1119  if (l_num_coords > 4 && closed) {
1120  // Also add the closing segment between the last and the first points
1121  centroid_add_segment(l2x, l2y, l0x, l0y, &length, linestring_centroid_sum);
1122  }
1123  *total_length += length;
1124  if (length == 0.0 && l_num_coords > 0) {
1125  *num_points += 1;
1126  point_centroid_sum[0] += l0x;
1127  point_centroid_sum[1] += l0y;
1128  }
1129  return true;
1130 }
1131 
1133 double ST_Centroid_LineString(int8_t* coords,
1134  int64_t coords_sz,
1135  int32_t ic,
1136  int32_t isr,
1137  int32_t osr,
1138  bool ycoord) {
1139  double length = 0.0;
1140  double linestring_centroid_sum[2] = {0.0, 0.0};
1141  int64_t num_points = 0;
1142  double point_centroid_sum[2] = {0.0, 0.0};
1143  centroid_add_linestring(coords,
1144  coords_sz,
1145  ic,
1146  isr,
1147  osr,
1148  false, // not closed
1149  &length,
1150  &linestring_centroid_sum[0],
1151  &num_points,
1152  &point_centroid_sum[0]);
1153  double linestring_centroid[2] = {0.0, 0.0};
1154  if (length > 0) {
1155  linestring_centroid[0] = linestring_centroid_sum[0] / length;
1156  linestring_centroid[1] = linestring_centroid_sum[1] / length;
1157  } else if (num_points > 0) {
1158  linestring_centroid[0] = point_centroid_sum[0] / num_points;
1159  linestring_centroid[1] = point_centroid_sum[1] / num_points;
1160  }
1161  if (ycoord) {
1162  return linestring_centroid[1];
1163  }
1164  return linestring_centroid[0];
1165 }
1166 
1168  double y1,
1169  double x2,
1170  double y2,
1171  double x3,
1172  double y3,
1173  double sign,
1174  double* total_area2,
1175  double* cg3) {
1176  double cx = x1 + x2 + x3;
1177  double cy = y1 + y2 + y3;
1178  double area2 = x1 * y2 - x2 * y1 + x3 * y1 - x1 * y3 + x2 * y3 - x3 * y2;
1179  cg3[0] += sign * area2 * cx;
1180  cg3[1] += sign * area2 * cy;
1181  *total_area2 += sign * area2;
1182  return true;
1183 }
1184 
1186  int64_t ringsize,
1187  int32_t ic,
1188  int32_t isr,
1189  int32_t osr,
1190  double sign,
1191  double* total_area2,
1192  double* cg3,
1193  double* total_length,
1194  double* linestring_centroid_sum,
1195  int64_t* num_points,
1196  double* point_centroid_sum) {
1197  auto ring_num_coords = ringsize / compression_unit_size(ic);
1198 
1199  if (ring_num_coords < 6) {
1200  return false;
1201  }
1202 
1203  double x1 = coord_x(ring, 0, ic, isr, osr);
1204  double y1 = coord_y(ring, 1, ic, isr, osr);
1205  double x2 = coord_x(ring, 2, ic, isr, osr);
1206  double y2 = coord_y(ring, 3, ic, isr, osr);
1207  for (int32_t i = 4; i < ring_num_coords; i += 2) {
1208  double x3 = coord_x(ring, i, ic, isr, osr);
1209  double y3 = coord_y(ring, i + 1, ic, isr, osr);
1210  centroid_add_triangle(x1, y1, x2, y2, x3, y3, sign, total_area2, cg3);
1211  x2 = x3;
1212  y2 = y3;
1213  }
1214 
1216  ringsize,
1217  ic,
1218  isr,
1219  osr,
1220  true, // closed
1221  total_length,
1222  linestring_centroid_sum,
1223  num_points,
1224  point_centroid_sum);
1225  return true;
1226 }
1227 
1228 DEVICE ALWAYS_INLINE bool centroid_add_polygon(int8_t* poly_coords,
1229  int64_t poly_coords_size,
1230  int32_t* poly_ring_sizes,
1231  int64_t poly_num_rings,
1232  int32_t ic,
1233  int32_t isr,
1234  int32_t osr,
1235  double* total_area2,
1236  double* cg3,
1237  double* total_length,
1238  double* linestring_centroid_sum,
1239  int64_t* num_points,
1240  double* point_centroid_sum) {
1241  if (poly_num_rings <= 0) {
1242  return false;
1243  }
1244 
1245  auto ring_coords = poly_coords;
1246 
1247  for (auto r = 0; r < poly_num_rings; r++) {
1248  auto ring_coords_size = poly_ring_sizes[r] * 2 * compression_unit_size(ic);
1249  // Shape - positive area, holes - negative areas
1250  double sign = (r == 0) ? 1.0 : -1.0;
1251  centroid_add_ring(ring_coords,
1252  ring_coords_size,
1253  ic,
1254  isr,
1255  osr,
1256  sign,
1257  total_area2,
1258  cg3,
1259  total_length,
1260  linestring_centroid_sum,
1261  num_points,
1262  point_centroid_sum);
1263  // Advance to the next ring.
1264  ring_coords += ring_coords_size;
1265  }
1266  return true;
1267 }
1268 
1270 double ST_Centroid_Polygon(int8_t* poly_coords,
1271  int64_t poly_coords_size,
1272  int32_t* poly_ring_sizes,
1273  int64_t poly_num_rings,
1274  int32_t ic,
1275  int32_t isr,
1276  int32_t osr,
1277  bool ycoord) {
1278  if (poly_num_rings <= 0) {
1279  return 0.0;
1280  }
1281  double total_area2 = 0.0;
1282  double cg3[2] = {0.0, 0.0};
1283  double total_length = 0.0;
1284  double linestring_centroid_sum[2] = {0.0, 0.0};
1285  int64_t num_points = 0;
1286  double point_centroid_sum[2] = {0.0, 0.0};
1287  centroid_add_polygon(poly_coords,
1288  poly_coords_size,
1289  poly_ring_sizes,
1290  poly_num_rings,
1291  ic,
1292  isr,
1293  osr,
1294  &total_area2,
1295  &cg3[0],
1296  &total_length,
1297  &linestring_centroid_sum[0],
1298  &num_points,
1299  &point_centroid_sum[0]);
1300 
1301  double x1 = coord_x(poly_coords, 0, ic, isr, osr);
1302  double y1 = coord_y(poly_coords, 1, ic, isr, osr);
1303  double poly_centroid[2] = {x1, y1};
1304  if (total_area2 != 0.0) {
1305  poly_centroid[0] = cg3[0] / 3 / total_area2;
1306  poly_centroid[1] = cg3[1] / 3 / total_area2;
1307  } else if (total_length > 0.0) {
1308  // zero-area polygon, fall back to linear centroid
1309  poly_centroid[0] = linestring_centroid_sum[0] / total_length;
1310  poly_centroid[1] = linestring_centroid_sum[1] / total_length;
1311  } else if (num_points > 0) {
1312  // zero-area zero-length polygon, fall further back to point centroid
1313  poly_centroid[0] = point_centroid_sum[0] / num_points;
1314  poly_centroid[1] = point_centroid_sum[1] / num_points;
1315  }
1316 
1317  if (ycoord) {
1318  return poly_centroid[1];
1319  }
1320  return poly_centroid[0];
1321 }
1322 
1324 double ST_Centroid_MultiPolygon(int8_t* mpoly_coords,
1325  int64_t mpoly_coords_size,
1326  int32_t* mpoly_ring_sizes,
1327  int64_t mpoly_num_rings,
1328  int32_t* mpoly_poly_sizes,
1329  int64_t mpoly_num_polys,
1330  int32_t ic,
1331  int32_t isr,
1332  int32_t osr,
1333  bool ycoord) {
1334  if (mpoly_num_rings <= 0 || mpoly_num_polys <= 0) {
1335  return 0.0;
1336  }
1337 
1338  double total_area2 = 0.0;
1339  double cg3[2] = {0.0, 0.0};
1340  double total_length = 0.0;
1341  double linestring_centroid_sum[2] = {0.0, 0.0};
1342  int64_t num_points = 0;
1343  double point_centroid_sum[2] = {0.0, 0.0};
1344 
1345  // Set specific poly pointers as we move through the coords/ringsizes/polyrings arrays.
1346  auto next_poly_coords = mpoly_coords;
1347  auto next_poly_ring_sizes = mpoly_ring_sizes;
1348 
1349  for (auto poly = 0; poly < mpoly_num_polys; poly++) {
1350  auto poly_coords = next_poly_coords;
1351  auto poly_ring_sizes = next_poly_ring_sizes;
1352  auto poly_num_rings = mpoly_poly_sizes[poly];
1353  // Count number of coords in all of poly's rings, advance ring size pointer.
1354  int32_t poly_num_coords = 0;
1355  for (auto ring = 0; ring < poly_num_rings; ring++) {
1356  poly_num_coords += 2 * *next_poly_ring_sizes++;
1357  }
1358  auto poly_coords_size = poly_num_coords * compression_unit_size(ic);
1359  next_poly_coords += poly_coords_size;
1360 
1361  centroid_add_polygon(poly_coords,
1362  poly_coords_size,
1363  poly_ring_sizes,
1364  poly_num_rings,
1365  ic,
1366  isr,
1367  osr,
1368  &total_area2,
1369  &cg3[0],
1370  &total_length,
1371  &linestring_centroid_sum[0],
1372  &num_points,
1373  &point_centroid_sum[0]);
1374  }
1375 
1376  double x1 = coord_x(mpoly_coords, 0, ic, isr, osr);
1377  double y1 = coord_y(mpoly_coords, 1, ic, isr, osr);
1378  double mpoly_centroid[2] = {x1, y1};
1379  if (total_area2 != 0.0) {
1380  mpoly_centroid[0] = cg3[0] / 3 / total_area2;
1381  mpoly_centroid[1] = cg3[1] / 3 / total_area2;
1382  } else if (total_length > 0.0) {
1383  // zero-area multipolygon, fall back to linear centroid
1384  mpoly_centroid[0] = linestring_centroid_sum[0] / total_length;
1385  mpoly_centroid[1] = linestring_centroid_sum[1] / total_length;
1386  } else if (num_points > 0) {
1387  // zero-area zero-length multipolygon, fall further back to point centroid
1388  mpoly_centroid[0] = point_centroid_sum[0] / num_points;
1389  mpoly_centroid[1] = point_centroid_sum[1] / num_points;
1390  }
1391 
1392  if (ycoord) {
1393  return mpoly_centroid[1];
1394  }
1395  return mpoly_centroid[0];
1396 }
1397 
1399 int32_t ST_NPoints(int8_t* coords, int64_t coords_sz, int32_t ic) {
1400  auto num_pts = coords_sz / compression_unit_size(ic);
1401  return static_cast<int32_t>(num_pts / 2);
1402 }
1403 
1405 int32_t ST_NRings(int32_t* poly_ring_sizes, int64_t poly_num_rings) {
1406  return static_cast<int32_t>(poly_num_rings);
1407 }
1408 
1409 //
1410 // ST_Distance
1411 //
1412 
1414 double ST_Distance_Point_Point(int8_t* p1,
1415  int64_t p1size,
1416  int8_t* p2,
1417  int64_t p2size,
1418  int32_t ic1,
1419  int32_t isr1,
1420  int32_t ic2,
1421  int32_t isr2,
1422  int32_t osr) {
1423  double p1x = coord_x(p1, 0, ic1, isr1, osr);
1424  double p1y = coord_y(p1, 1, ic1, isr1, osr);
1425  double p2x = coord_x(p2, 0, ic2, isr2, osr);
1426  double p2y = coord_y(p2, 1, ic2, isr2, osr);
1427  return distance_point_point(p1x, p1y, p2x, p2y);
1428 }
1429 
1432  int64_t p1size,
1433  int8_t* p2,
1434  int64_t p2size,
1435  int32_t ic1,
1436  int32_t isr1,
1437  int32_t ic2,
1438  int32_t isr2,
1439  int32_t osr) {
1440  double p1x = coord_x(p1, 0, ic1, isr1, osr);
1441  double p1y = coord_y(p1, 1, ic1, isr1, osr);
1442  double p2x = coord_x(p2, 0, ic2, isr2, osr);
1443  double p2y = coord_y(p2, 1, ic2, isr2, osr);
1444  return distance_point_point_squared(p1x, p1y, p2x, p2y);
1445 }
1446 
1449  int64_t p1size,
1450  int8_t* p2,
1451  int64_t p2size,
1452  int32_t ic1,
1453  int32_t isr1,
1454  int32_t ic2,
1455  int32_t isr2,
1456  int32_t osr) {
1457  double p1x = coord_x(p1, 0, ic1, 4326, 4326);
1458  double p1y = coord_y(p1, 1, ic1, 4326, 4326);
1459  double p2x = coord_x(p2, 0, ic2, 4326, 4326);
1460  double p2y = coord_y(p2, 1, ic2, 4326, 4326);
1461  return distance_in_meters(p1x, p1y, p2x, p2y);
1462 }
1463 
1466  int64_t psize,
1467  int8_t* l,
1468  int64_t lsize,
1469  int32_t lindex,
1470  int32_t ic1,
1471  int32_t isr1,
1472  int32_t ic2,
1473  int32_t isr2,
1474  int32_t osr) {
1475  // Currently only statically indexed LineString is supported
1476  double px = coord_x(p, 0, ic1, 4326, 4326);
1477  double py = coord_y(p, 1, ic1, 4326, 4326);
1478  auto lpoints = lsize / (2 * compression_unit_size(ic2));
1479  if (lindex < 0 || lindex > lpoints) {
1480  lindex = lpoints; // Endpoint
1481  }
1482  double lx = coord_x(l, 2 * (lindex - 1), ic2, 4326, 4326);
1483  double ly = coord_y(l, 2 * (lindex - 1) + 1, ic2, 4326, 4326);
1484  return distance_in_meters(px, py, lx, ly);
1485 }
1486 
1489  int64_t lsize,
1490  int32_t lindex,
1491  int8_t* p,
1492  int64_t psize,
1493  int32_t ic1,
1494  int32_t isr1,
1495  int32_t ic2,
1496  int32_t isr2,
1497  int32_t osr) {
1498  // Currently only statically indexed LineString is supported
1500  p, psize, l, lsize, lindex, ic2, isr2, ic1, isr1, osr);
1501 }
1502 
1505  int64_t l1size,
1506  int32_t l1index,
1507  int8_t* l2,
1508  int64_t l2size,
1509  int32_t l2index,
1510  int32_t ic1,
1511  int32_t isr1,
1512  int32_t ic2,
1513  int32_t isr2,
1514  int32_t osr) {
1515  // Currently only statically indexed LineStrings are supported
1516  auto l1points = l1size / (2 * compression_unit_size(ic1));
1517  if (l1index < 0 || l1index > l1points) {
1518  l1index = l1points; // Endpoint
1519  }
1520  double l1x = coord_x(l1, 2 * (l1index - 1), ic1, 4326, 4326);
1521  double l1y = coord_y(l1, 2 * (l1index - 1) + 1, ic1, 4326, 4326);
1522  auto l2points = l2size / (2 * compression_unit_size(ic2));
1523  if (l2index < 0 || l2index > l2points) {
1524  l2index = l2points; // Endpoint
1525  }
1526  double l2x = coord_x(l2, 2 * (l2index - 1), ic2, 4326, 4326);
1527  double l2y = coord_y(l2, 2 * (l2index - 1) + 1, ic2, 4326, 4326);
1528  return distance_in_meters(l1x, l1y, l2x, l2y);
1529 }
1530 
1532  int64_t psize,
1533  int8_t* l,
1534  int64_t lsize,
1535  int32_t lindex,
1536  int32_t ic1,
1537  int32_t isr1,
1538  int32_t ic2,
1539  int32_t isr2,
1540  int32_t osr,
1541  bool check_closed) {
1542  double px = coord_x(p, 0, ic1, isr1, osr);
1543  double py = coord_y(p, 1, ic1, isr1, osr);
1544 
1545  auto l_num_coords = lsize / compression_unit_size(ic2);
1546  auto l_num_points = l_num_coords / 2;
1547  if (lindex != 0) { // Statically indexed linestring
1548  if (lindex < 0 || lindex > l_num_points) {
1549  lindex = l_num_points; // Endpoint
1550  }
1551  double lx = coord_x(l, 2 * (lindex - 1), ic2, isr2, osr);
1552  double ly = coord_y(l, 2 * (lindex - 1) + 1, ic2, isr2, osr);
1553  return distance_point_point(px, py, lx, ly);
1554  }
1555 
1556  double l1x = coord_x(l, 0, ic2, isr2, osr);
1557  double l1y = coord_y(l, 1, ic2, isr2, osr);
1558  double l2x = coord_x(l, 2, ic2, isr2, osr);
1559  double l2y = coord_y(l, 3, ic2, isr2, osr);
1560 
1561  double dist = distance_point_line(px, py, l1x, l1y, l2x, l2y);
1562  for (int32_t i = 4; i < l_num_coords; i += 2) {
1563  l1x = l2x; // advance one point
1564  l1y = l2y;
1565  l2x = coord_x(l, i, ic2, isr2, osr);
1566  l2y = coord_y(l, i + 1, ic2, isr2, osr);
1567  double ldist = distance_point_line(px, py, l1x, l1y, l2x, l2y);
1568  if (dist > ldist) {
1569  dist = ldist;
1570  }
1571  }
1572  if (l_num_coords > 4 && check_closed) {
1573  // Also check distance to the closing edge between the first and the last points
1574  l1x = coord_x(l, 0, ic2, isr2, osr);
1575  l1y = coord_y(l, 1, ic2, isr2, osr);
1576  double ldist = distance_point_line(px, py, l1x, l1y, l2x, l2y);
1577  if (dist > ldist) {
1578  dist = ldist;
1579  }
1580  }
1581  return dist;
1582 }
1583 
1586  int64_t psize,
1587  int8_t* l,
1588  int64_t lsize,
1589  int32_t lindex,
1590  int32_t ic1,
1591  int32_t isr1,
1592  int32_t ic2,
1593  int32_t isr2,
1594  int32_t osr) {
1596  p, psize, l, lsize, lindex, ic1, isr1, ic2, isr2, osr, true);
1597 }
1598 
1601  int64_t psize,
1602  int8_t* l,
1603  int64_t lsize,
1604  int32_t lindex,
1605  int32_t ic1,
1606  int32_t isr1,
1607  int32_t ic2,
1608  int32_t isr2,
1609  int32_t osr) {
1610  if (lindex != 0) { // Statically indexed linestring
1611  auto l_num_coords = lsize / compression_unit_size(ic2);
1612  auto l_num_points = l_num_coords / 2;
1613  if (lindex < 0 || lindex > l_num_points) {
1614  lindex = l_num_points; // Endpoint
1615  }
1616  double px = coord_x(p, 0, ic1, isr1, osr);
1617  double py = coord_y(p, 1, ic1, isr1, osr);
1618  double lx = coord_x(l, 2 * (lindex - 1), ic2, isr2, osr);
1619  double ly = coord_y(l, 2 * (lindex - 1) + 1, ic2, isr2, osr);
1620  return distance_point_point(px, py, lx, ly);
1621  }
1622 
1624  p, psize, l, lsize, lindex, ic1, isr1, ic2, isr2, osr, false);
1625 }
1626 
1628 double ST_Distance_Point_Polygon(int8_t* p,
1629  int64_t psize,
1630  int8_t* poly,
1631  int64_t polysize,
1632  int32_t* poly_ring_sizes,
1633  int64_t poly_num_rings,
1634  int32_t ic1,
1635  int32_t isr1,
1636  int32_t ic2,
1637  int32_t isr2,
1638  int32_t osr) {
1639  auto exterior_ring_num_coords = polysize / compression_unit_size(ic2);
1640  if (poly_num_rings > 0) {
1641  exterior_ring_num_coords = poly_ring_sizes[0] * 2;
1642  }
1643  auto exterior_ring_coords_size = exterior_ring_num_coords * compression_unit_size(ic2);
1644 
1645  double px = coord_x(p, 0, ic1, isr1, osr);
1646  double py = coord_y(p, 1, ic1, isr1, osr);
1647  if (!polygon_contains_point(poly, exterior_ring_num_coords, px, py, ic2, isr2, osr)) {
1648  // Outside the exterior ring
1650  p, psize, poly, exterior_ring_coords_size, 0, ic1, isr1, ic2, isr2, osr);
1651  }
1652  // Inside exterior ring
1653  // Advance to first interior ring
1654  poly += exterior_ring_coords_size;
1655  // Check if one of the polygon's holes contains that point
1656  for (auto r = 1; r < poly_num_rings; r++) {
1657  auto interior_ring_num_coords = poly_ring_sizes[r] * 2;
1658  auto interior_ring_coords_size =
1659  interior_ring_num_coords * compression_unit_size(ic2);
1660  if (polygon_contains_point(poly, interior_ring_num_coords, px, py, ic2, isr2, osr)) {
1661  // Inside an interior ring
1663  p, psize, poly, interior_ring_coords_size, 0, ic1, isr1, ic2, isr2, osr);
1664  }
1665  poly += interior_ring_coords_size;
1666  }
1667  return 0.0;
1668 }
1669 
1672  int64_t psize,
1673  int8_t* mpoly_coords,
1674  int64_t mpoly_coords_size,
1675  int32_t* mpoly_ring_sizes,
1676  int64_t mpoly_num_rings,
1677  int32_t* mpoly_poly_sizes,
1678  int64_t mpoly_num_polys,
1679  int32_t ic1,
1680  int32_t isr1,
1681  int32_t ic2,
1682  int32_t isr2,
1683  int32_t osr) {
1684  if (mpoly_num_polys <= 0) {
1685  return 0.0;
1686  }
1687  double min_distance = 0.0;
1688 
1689  // Set specific poly pointers as we move through the coords/ringsizes/polyrings arrays.
1690  auto next_poly_coords = mpoly_coords;
1691  auto next_poly_ring_sizes = mpoly_ring_sizes;
1692 
1693  for (auto poly = 0; poly < mpoly_num_polys; poly++) {
1694  auto poly_coords = next_poly_coords;
1695  auto poly_ring_sizes = next_poly_ring_sizes;
1696  auto poly_num_rings = mpoly_poly_sizes[poly];
1697  // Count number of coords in all of poly's rings, advance ring size pointer.
1698  int32_t poly_num_coords = 0;
1699  for (auto ring = 0; ring < poly_num_rings; ring++) {
1700  poly_num_coords += 2 * *next_poly_ring_sizes++;
1701  }
1702  auto poly_coords_size = poly_num_coords * compression_unit_size(ic2);
1703  next_poly_coords += poly_coords_size;
1704  double distance = ST_Distance_Point_Polygon(p,
1705  psize,
1706  poly_coords,
1707  poly_coords_size,
1708  poly_ring_sizes,
1709  poly_num_rings,
1710  ic1,
1711  isr1,
1712  ic2,
1713  isr2,
1714  osr);
1715  if (poly == 0 || min_distance > distance) {
1716  min_distance = distance;
1717  if (tol_zero(min_distance)) {
1718  min_distance = 0.0;
1719  break;
1720  }
1721  }
1722  }
1723 
1724  return min_distance;
1725 }
1726 
1729  int64_t lsize,
1730  int32_t lindex,
1731  int8_t* p,
1732  int64_t psize,
1733  int32_t ic1,
1734  int32_t isr1,
1735  int32_t ic2,
1736  int32_t isr2,
1737  int32_t osr) {
1739  p, psize, l, lsize, lindex, ic2, isr2, ic1, isr1, osr);
1740 }
1741 
1744  int64_t l1size,
1745  int32_t l1index,
1746  int8_t* l2,
1747  int64_t l2size,
1748  int32_t l2index,
1749  int32_t ic1,
1750  int32_t isr1,
1751  int32_t ic2,
1752  int32_t isr2,
1753  int32_t osr) {
1754  auto l1_num_coords = l1size / compression_unit_size(ic1);
1755  if (l1index != 0) {
1756  // l1 is a statically indexed linestring
1757  auto l1_num_points = l1_num_coords / 2;
1758  if (l1index < 0 || l1index > l1_num_points) {
1759  l1index = l1_num_points;
1760  }
1761  int8_t* p = l1 + 2 * (l1index - 1) * compression_unit_size(ic1);
1762  int64_t psize = 2 * compression_unit_size(ic1);
1764  p, psize, l2, l2size, l2index, ic1, isr1, ic2, isr2, osr);
1765  }
1766 
1767  auto l2_num_coords = l2size / compression_unit_size(ic2);
1768  if (l2index != 0) {
1769  // l2 is a statically indexed linestring
1770  auto l2_num_points = l2_num_coords / 2;
1771  if (l2index < 0 || l2index > l2_num_points) {
1772  l2index = l2_num_points;
1773  }
1774  int8_t* p = l2 + 2 * (l2index - 1) * compression_unit_size(ic2);
1775  int64_t psize = 2 * compression_unit_size(ic2);
1777  p, psize, l1, l1size, l1index, ic2, isr2, ic1, isr1, osr);
1778  }
1779 
1780  double dist = 0.0;
1781  double l11x = coord_x(l1, 0, ic1, isr1, osr);
1782  double l11y = coord_y(l1, 1, ic1, isr1, osr);
1783  for (int32_t i1 = 2; i1 < l1_num_coords; i1 += 2) {
1784  double l12x = coord_x(l1, i1, ic1, isr1, osr);
1785  double l12y = coord_y(l1, i1 + 1, ic1, isr1, osr);
1786 
1787  double l21x = coord_x(l2, 0, ic2, isr2, osr);
1788  double l21y = coord_y(l2, 1, ic2, isr2, osr);
1789  for (int32_t i2 = 2; i2 < l2_num_coords; i2 += 2) {
1790  double l22x = coord_x(l2, i2, ic2, isr2, osr);
1791  double l22y = coord_y(l2, i2 + 1, ic2, isr2, osr);
1792 
1793  double ldist = distance_line_line(l11x, l11y, l12x, l12y, l21x, l21y, l22x, l22y);
1794  if (i1 == 2 && i2 == 2) {
1795  dist = ldist; // initialize dist with distance between the first two segments
1796  } else if (dist > ldist) {
1797  dist = ldist;
1798  }
1799  if (tol_zero(dist)) {
1800  return 0.0; // segments touch
1801  }
1802 
1803  l21x = l22x; // advance to the next point on l2
1804  l21y = l22y;
1805  }
1806 
1807  l11x = l12x; // advance to the next point on l1
1808  l11y = l12y;
1809  }
1810  return dist;
1811 }
1812 
1815  int64_t lsize,
1816  int32_t lindex,
1817  int8_t* poly_coords,
1818  int64_t poly_coords_size,
1819  int32_t* poly_ring_sizes,
1820  int64_t poly_num_rings,
1821  int32_t ic1,
1822  int32_t isr1,
1823  int32_t ic2,
1824  int32_t isr2,
1825  int32_t osr) {
1826  auto lnum_coords = lsize / compression_unit_size(ic1);
1827  auto lnum_points = lnum_coords / 2;
1828  if (lindex < 0 || lindex > lnum_points) {
1829  lindex = lnum_points;
1830  }
1831  auto p = l + lindex * compression_unit_size(ic1);
1832  auto psize = 2 * compression_unit_size(ic1);
1833  auto min_distance = ST_Distance_Point_Polygon(p,
1834  psize,
1835  poly_coords,
1836  poly_coords_size,
1837  poly_ring_sizes,
1838  poly_num_rings,
1839  ic1,
1840  isr1,
1841  ic2,
1842  isr2,
1843  osr);
1844  if (lindex != 0) {
1845  // Statically indexed linestring: return distance from the indexed point to poly
1846  return min_distance;
1847  }
1848  if (tol_zero(min_distance)) {
1849  // Linestring's first point is inside the poly
1850  return 0.0;
1851  }
1852 
1853  // Otherwise, linestring's first point is outside the external ring or inside
1854  // an internal ring. Measure minimum distance between linestring segments and
1855  // poly rings. Crossing a ring zeroes the distance and causes an early return.
1856  auto poly_ring_coords = poly_coords;
1857  for (auto r = 0; r < poly_num_rings; r++) {
1858  int64_t poly_ring_num_coords = poly_ring_sizes[r] * 2;
1859 
1860  auto distance = distance_ring_linestring(poly_ring_coords,
1861  poly_ring_num_coords,
1862  l,
1863  lnum_coords,
1864  ic2,
1865  isr2,
1866  ic1,
1867  isr1,
1868  osr);
1869  if (min_distance > distance) {
1870  min_distance = distance;
1871  if (tol_zero(min_distance)) {
1872  return 0.0;
1873  }
1874  }
1875 
1876  poly_ring_coords += poly_ring_num_coords * compression_unit_size(ic2);
1877  }
1878 
1879  return min_distance;
1880 }
1881 
1884  int64_t lsize,
1885  int32_t lindex,
1886  int8_t* mpoly_coords,
1887  int64_t mpoly_coords_size,
1888  int32_t* mpoly_ring_sizes,
1889  int64_t mpoly_num_rings,
1890  int32_t* mpoly_poly_sizes,
1891  int64_t mpoly_num_polys,
1892  int32_t ic1,
1893  int32_t isr1,
1894  int32_t ic2,
1895  int32_t isr2,
1896  int32_t osr) {
1897  // TODO: revisit implementation, cover all cases
1898 
1899  auto lnum_coords = lsize / compression_unit_size(ic1);
1900  auto lnum_points = lnum_coords / 2;
1901  if (lindex != 0) {
1902  // Statically indexed linestring
1903  if (lindex < 0 || lindex > lnum_points) {
1904  lindex = lnum_points;
1905  }
1906  auto p = l + lindex * compression_unit_size(ic1);
1907  auto psize = 2 * compression_unit_size(ic1);
1909  psize,
1910  mpoly_coords,
1911  mpoly_coords_size,
1912  mpoly_ring_sizes,
1913  mpoly_num_rings,
1914  mpoly_poly_sizes,
1915  mpoly_num_polys,
1916  ic1,
1917  isr1,
1918  ic2,
1919  isr2,
1920  osr);
1921  }
1922 
1923  double min_distance = 0.0;
1924 
1925  // Set specific poly pointers as we move through the coords/ringsizes/polyrings arrays.
1926  auto next_poly_coords = mpoly_coords;
1927  auto next_poly_ring_sizes = mpoly_ring_sizes;
1928 
1929  for (auto poly = 0; poly < mpoly_num_polys; poly++) {
1930  auto poly_coords = next_poly_coords;
1931  auto poly_ring_sizes = next_poly_ring_sizes;
1932  auto poly_num_rings = mpoly_poly_sizes[poly];
1933  // Count number of coords in all of poly's rings, advance ring size pointer.
1934  int32_t poly_num_coords = 0;
1935  for (auto ring = 0; ring < poly_num_rings; ring++) {
1936  poly_num_coords += 2 * *next_poly_ring_sizes++;
1937  }
1938  auto poly_coords_size = poly_num_coords * compression_unit_size(ic2);
1939  next_poly_coords += poly_coords_size;
1940  double distance = ST_Distance_LineString_Polygon(l,
1941  lsize,
1942  lindex,
1943  poly_coords,
1944  poly_coords_size,
1945  poly_ring_sizes,
1946  poly_num_rings,
1947  ic1,
1948  isr1,
1949  ic2,
1950  isr2,
1951  osr);
1952  if (poly == 0 || min_distance > distance) {
1953  min_distance = distance;
1954  if (tol_zero(min_distance)) {
1955  min_distance = 0.0;
1956  break;
1957  }
1958  }
1959  }
1960 
1961  return min_distance;
1962 }
1963 
1965 double ST_Distance_Polygon_Point(int8_t* poly_coords,
1966  int64_t poly_coords_size,
1967  int32_t* poly_ring_sizes,
1968  int64_t poly_num_rings,
1969  int8_t* p,
1970  int64_t psize,
1971  int32_t ic1,
1972  int32_t isr1,
1973  int32_t ic2,
1974  int32_t isr2,
1975  int32_t osr) {
1976  return ST_Distance_Point_Polygon(p,
1977  psize,
1978  poly_coords,
1979  poly_coords_size,
1980  poly_ring_sizes,
1981  poly_num_rings,
1982  ic2,
1983  isr2,
1984  ic1,
1985  isr1,
1986  osr);
1987 }
1988 
1990 double ST_Distance_Polygon_LineString(int8_t* poly_coords,
1991  int64_t poly_coords_size,
1992  int32_t* poly_ring_sizes,
1993  int64_t poly_num_rings,
1994  int8_t* l,
1995  int64_t lsize,
1996  int32_t li,
1997  int32_t ic1,
1998  int32_t isr1,
1999  int32_t ic2,
2000  int32_t isr2,
2001  int32_t osr) {
2003  lsize,
2004  li,
2005  poly_coords,
2006  poly_coords_size,
2007  poly_ring_sizes,
2008  poly_num_rings,
2009  ic2,
2010  isr2,
2011  ic1,
2012  isr2,
2013  osr);
2014 }
2015 
2017 double ST_Distance_Polygon_Polygon(int8_t* poly1_coords,
2018  int64_t poly1_coords_size,
2019  int32_t* poly1_ring_sizes,
2020  int64_t poly1_num_rings,
2021  int8_t* poly2_coords,
2022  int64_t poly2_coords_size,
2023  int32_t* poly2_ring_sizes,
2024  int64_t poly2_num_rings,
2025  int32_t ic1,
2026  int32_t isr1,
2027  int32_t ic2,
2028  int32_t isr2,
2029  int32_t osr) {
2030  // Check if poly1 contains the first point of poly2's shape, i.e. the external ring
2031  auto poly2_first_point_coords = poly2_coords;
2032  auto poly2_first_point_coords_size = compression_unit_size(ic2) * 2;
2033  auto min_distance = ST_Distance_Polygon_Point(poly1_coords,
2034  poly1_coords_size,
2035  poly1_ring_sizes,
2036  poly1_num_rings,
2037  poly2_first_point_coords,
2038  poly2_first_point_coords_size,
2039  ic1,
2040  isr1,
2041  ic2,
2042  isr2,
2043  osr);
2044  if (tol_zero(min_distance)) {
2045  // Polygons overlap
2046  return 0.0;
2047  }
2048 
2049  // Poly2's first point is either outside poly1's external ring or inside one of the
2050  // internal rings. Measure the smallest distance between a poly1 ring (external or
2051  // internal) and a poly2 ring (external or internal). If poly2 is completely outside
2052  // poly1, then the min distance would be between poly1's and poly2's external rings. If
2053  // poly2 is completely inside one of poly1 internal rings then the min distance would be
2054  // between that poly1 internal ring and poly2's external ring. If poly1 is completely
2055  // inside one of poly2 internal rings, min distance is between that internal ring and
2056  // poly1's external ring. In each case other rings don't get in the way. Any ring
2057  // intersection means zero distance - short-circuit and return.
2058 
2059  auto poly1_ring_coords = poly1_coords;
2060  for (auto r1 = 0; r1 < poly1_num_rings; r1++) {
2061  int64_t poly1_ring_num_coords = poly1_ring_sizes[r1] * 2;
2062 
2063  auto poly2_ring_coords = poly2_coords;
2064  for (auto r2 = 0; r2 < poly2_num_rings; r2++) {
2065  int64_t poly2_ring_num_coords = poly2_ring_sizes[r2] * 2;
2066 
2067  auto distance = distance_ring_ring(poly1_ring_coords,
2068  poly1_ring_num_coords,
2069  poly2_ring_coords,
2070  poly2_ring_num_coords,
2071  ic1,
2072  isr1,
2073  ic2,
2074  isr2,
2075  osr);
2076  if (min_distance > distance) {
2077  min_distance = distance;
2078  if (tol_zero(min_distance)) {
2079  return 0.0;
2080  }
2081  }
2082 
2083  poly2_ring_coords += poly2_ring_num_coords * compression_unit_size(ic2);
2084  }
2085 
2086  poly1_ring_coords += poly1_ring_num_coords * compression_unit_size(ic1);
2087  }
2088 
2089  return min_distance;
2090 }
2091 
2093 double ST_Distance_Polygon_MultiPolygon(int8_t* poly1_coords,
2094  int64_t poly1_coords_size,
2095  int32_t* poly1_ring_sizes,
2096  int64_t poly1_num_rings,
2097  int8_t* mpoly_coords,
2098  int64_t mpoly_coords_size,
2099  int32_t* mpoly_ring_sizes,
2100  int64_t mpoly_num_rings,
2101  int32_t* mpoly_poly_sizes,
2102  int64_t mpoly_num_polys,
2103  int32_t ic1,
2104  int32_t isr1,
2105  int32_t ic2,
2106  int32_t isr2,
2107  int32_t osr) {
2108  double min_distance = 0.0;
2109 
2110  // Set specific poly pointers as we move through the coords/ringsizes/polyrings arrays.
2111  auto next_poly_coords = mpoly_coords;
2112  auto next_poly_ring_sizes = mpoly_ring_sizes;
2113 
2114  for (auto poly = 0; poly < mpoly_num_polys; poly++) {
2115  auto poly_coords = next_poly_coords;
2116  auto poly_ring_sizes = next_poly_ring_sizes;
2117  auto poly_num_rings = mpoly_poly_sizes[poly];
2118  // Count number of coords in all of poly's rings, advance ring size pointer.
2119  int32_t poly_num_coords = 0;
2120  for (auto ring = 0; ring < poly_num_rings; ring++) {
2121  poly_num_coords += 2 * *next_poly_ring_sizes++;
2122  }
2123  auto poly_coords_size = poly_num_coords * compression_unit_size(ic2);
2124  next_poly_coords += poly_coords_size;
2125  double distance = ST_Distance_Polygon_Polygon(poly1_coords,
2126  poly1_coords_size,
2127  poly1_ring_sizes,
2128  poly1_num_rings,
2129  poly_coords,
2130  poly_coords_size,
2131  poly_ring_sizes,
2132  poly_num_rings,
2133  ic1,
2134  isr1,
2135  ic2,
2136  isr2,
2137  osr);
2138  if (poly == 0 || min_distance > distance) {
2139  min_distance = distance;
2140  if (tol_zero(min_distance)) {
2141  min_distance = 0.0;
2142  break;
2143  }
2144  }
2145  }
2146 
2147  return min_distance;
2148 }
2149 
2151 double ST_Distance_MultiPolygon_Point(int8_t* mpoly_coords,
2152  int64_t mpoly_coords_size,
2153  int32_t* mpoly_ring_sizes,
2154  int64_t mpoly_num_rings,
2155  int32_t* mpoly_poly_sizes,
2156  int64_t mpoly_num_polys,
2157  int8_t* p,
2158  int64_t psize,
2159  int32_t ic1,
2160  int32_t isr1,
2161  int32_t ic2,
2162  int32_t isr2,
2163  int32_t osr) {
2165  psize,
2166  mpoly_coords,
2167  mpoly_coords_size,
2168  mpoly_ring_sizes,
2169  mpoly_num_rings,
2170  mpoly_poly_sizes,
2171  mpoly_num_polys,
2172  ic2,
2173  isr2,
2174  ic1,
2175  isr1,
2176  osr);
2177 }
2178 
2180 double ST_Distance_MultiPolygon_LineString(int8_t* mpoly_coords,
2181  int64_t mpoly_coords_size,
2182  int32_t* mpoly_ring_sizes,
2183  int64_t mpoly_num_rings,
2184  int32_t* mpoly_poly_sizes,
2185  int64_t mpoly_num_polys,
2186  int8_t* l,
2187  int64_t lsize,
2188  int32_t lindex,
2189  int32_t ic1,
2190  int32_t isr1,
2191  int32_t ic2,
2192  int32_t isr2,
2193  int32_t osr) {
2195  lsize,
2196  lindex,
2197  mpoly_coords,
2198  mpoly_coords_size,
2199  mpoly_ring_sizes,
2200  mpoly_num_rings,
2201  mpoly_poly_sizes,
2202  mpoly_num_polys,
2203  ic2,
2204  isr2,
2205  ic1,
2206  isr1,
2207  osr);
2208 }
2209 
2211 double ST_Distance_MultiPolygon_Polygon(int8_t* mpoly_coords,
2212  int64_t mpoly_coords_size,
2213  int32_t* mpoly_ring_sizes,
2214  int64_t mpoly_num_rings,
2215  int32_t* mpoly_poly_sizes,
2216  int64_t mpoly_num_polys,
2217  int8_t* poly1_coords,
2218  int64_t poly1_coords_size,
2219  int32_t* poly1_ring_sizes,
2220  int64_t poly1_num_rings,
2221  int32_t ic1,
2222  int32_t isr1,
2223  int32_t ic2,
2224  int32_t isr2,
2225  int32_t osr) {
2226  return ST_Distance_Polygon_MultiPolygon(poly1_coords,
2227  poly1_coords_size,
2228  poly1_ring_sizes,
2229  poly1_num_rings,
2230  mpoly_coords,
2231  mpoly_coords_size,
2232  mpoly_ring_sizes,
2233  mpoly_num_rings,
2234  mpoly_poly_sizes,
2235  mpoly_num_polys,
2236  ic2,
2237  isr2,
2238  ic1,
2239  isr1,
2240  osr);
2241 }
2242 
2244 double ST_Distance_MultiPolygon_MultiPolygon(int8_t* mpoly1_coords,
2245  int64_t mpoly1_coords_size,
2246  int32_t* mpoly1_ring_sizes,
2247  int64_t mpoly1_num_rings,
2248  int32_t* mpoly1_poly_sizes,
2249  int64_t mpoly1_num_polys,
2250  int8_t* mpoly2_coords,
2251  int64_t mpoly2_coords_size,
2252  int32_t* mpoly2_ring_sizes,
2253  int64_t mpoly2_num_rings,
2254  int32_t* mpoly2_poly_sizes,
2255  int64_t mpoly2_num_polys,
2256  int32_t ic1,
2257  int32_t isr1,
2258  int32_t ic2,
2259  int32_t isr2,
2260  int32_t osr) {
2261  double min_distance = 0.0;
2262 
2263  // Set specific poly pointers as we move through mpoly1's coords/ringsizes/polyrings
2264  // arrays.
2265  auto next_poly_coords = mpoly1_coords;
2266  auto next_poly_ring_sizes = mpoly1_ring_sizes;
2267 
2268  for (auto poly = 0; poly < mpoly1_num_polys; poly++) {
2269  auto poly_coords = next_poly_coords;
2270  auto poly_ring_sizes = next_poly_ring_sizes;
2271  auto poly_num_rings = mpoly1_poly_sizes[poly];
2272  // Count number of coords in all of poly's rings, advance ring size pointer.
2273  int32_t poly_num_coords = 0;
2274  for (auto ring = 0; ring < poly_num_rings; ring++) {
2275  poly_num_coords += 2 * *next_poly_ring_sizes++;
2276  }
2277  auto poly_coords_size = poly_num_coords * compression_unit_size(ic1);
2278  next_poly_coords += poly_coords_size;
2279  double distance = ST_Distance_Polygon_MultiPolygon(poly_coords,
2280  poly_coords_size,
2281  poly_ring_sizes,
2282  poly_num_rings,
2283  mpoly2_coords,
2284  mpoly2_coords_size,
2285  mpoly2_ring_sizes,
2286  mpoly2_num_rings,
2287  mpoly2_poly_sizes,
2288  mpoly2_num_polys,
2289  ic1,
2290  isr1,
2291  ic2,
2292  isr2,
2293  osr);
2294  if (poly == 0 || min_distance > distance) {
2295  min_distance = distance;
2296  if (tol_zero(min_distance)) {
2297  min_distance = 0.0;
2298  break;
2299  }
2300  }
2301  }
2302 
2303  return min_distance;
2304 }
2305 
2306 //
2307 // ST_MaxDistance
2308 //
2309 
2310 // Max cartesian distance between a point and a line segment
2311 DEVICE
2312 double max_distance_point_line(double px,
2313  double py,
2314  double l1x,
2315  double l1y,
2316  double l2x,
2317  double l2y) {
2318  double length1 = distance_point_point(px, py, l1x, l1y);
2319  double length2 = distance_point_point(px, py, l2x, l2y);
2320  if (length1 > length2) {
2321  return length1;
2322  }
2323  return length2;
2324 }
2325 
2327  int64_t psize,
2328  int8_t* l,
2329  int64_t lsize,
2330  int32_t lindex,
2331  int32_t ic1,
2332  int32_t isr1,
2333  int32_t ic2,
2334  int32_t isr2,
2335  int32_t osr,
2336  bool check_closed) {
2337  double px = coord_x(p, 0, ic1, isr1, osr);
2338  double py = coord_y(p, 1, ic1, isr1, osr);
2339 
2340  auto l_num_coords = lsize / compression_unit_size(ic2);
2341  auto l_num_points = l_num_coords / 2;
2342  if (lindex != 0) { // Statically indexed linestring
2343  if (lindex < 0 || lindex > l_num_points) {
2344  lindex = l_num_points; // Endpoint
2345  }
2346  double lx = coord_x(l, 2 * (lindex - 1), ic2, isr2, osr);
2347  double ly = coord_y(l, 2 * (lindex - 1) + 1, ic2, isr2, osr);
2348  return distance_point_point(px, py, lx, ly);
2349  }
2350 
2351  double l1x = coord_x(l, 0, ic2, isr2, osr);
2352  double l1y = coord_y(l, 1, ic2, isr2, osr);
2353  double l2x = coord_x(l, 2, ic2, isr2, osr);
2354  double l2y = coord_y(l, 3, ic2, isr2, osr);
2355 
2356  double max_dist = max_distance_point_line(px, py, l1x, l1y, l2x, l2y);
2357  for (int32_t i = 4; i < l_num_coords; i += 2) {
2358  l1x = l2x; // advance one point
2359  l1y = l2y;
2360  l2x = coord_x(l, i, ic2, isr2, osr);
2361  l2y = coord_y(l, i + 1, ic2, isr2, osr);
2362  double ldist = max_distance_point_line(px, py, l1x, l1y, l2x, l2y);
2363  if (max_dist < ldist) {
2364  max_dist = ldist;
2365  }
2366  }
2367  if (l_num_coords > 4 && check_closed) {
2368  // Also check distance to the closing edge between the first and the last points
2369  l1x = coord_x(l, 0, ic2, isr2, osr);
2370  l1y = coord_y(l, 1, ic2, isr2, osr);
2371  double ldist = max_distance_point_line(px, py, l1x, l1y, l2x, l2y);
2372  if (max_dist < ldist) {
2373  max_dist = ldist;
2374  }
2375  }
2376  return max_dist;
2377 }
2378 
2381  int64_t psize,
2382  int8_t* l,
2383  int64_t lsize,
2384  int32_t lindex,
2385  int32_t ic1,
2386  int32_t isr1,
2387  int32_t ic2,
2388  int32_t isr2,
2389  int32_t osr) {
2391  p, psize, l, lsize, lindex, ic1, isr1, ic2, isr2, osr, false);
2392 }
2393 
2396  int64_t lsize,
2397  int32_t lindex,
2398  int8_t* p,
2399  int64_t psize,
2400  int32_t ic1,
2401  int32_t isr1,
2402  int32_t ic2,
2403  int32_t isr2,
2404  int32_t osr) {
2406  p, psize, l, lsize, lindex, ic2, isr2, ic1, isr1, osr, false);
2407 }
2408 
2409 //
2410 // ST_Contains
2411 //
2412 
2414 bool ST_Contains_Point_Point(int8_t* p1,
2415  int64_t p1size,
2416  int8_t* p2,
2417  int64_t p2size,
2418  int32_t ic1,
2419  int32_t isr1,
2420  int32_t ic2,
2421  int32_t isr2,
2422  int32_t osr) {
2423  double p1x = coord_x(p1, 0, ic1, isr1, osr);
2424  double p1y = coord_y(p1, 1, ic1, isr1, osr);
2425  double p2x = coord_x(p2, 0, ic2, isr2, osr);
2426  double p2y = coord_y(p2, 1, ic2, isr2, osr);
2427  double tolerance = tol(ic1, ic2);
2428  return tol_eq(p1x, p2x, tolerance) && tol_eq(p1y, p2y, tolerance);
2429 }
2430 
2433  int64_t psize,
2434  int8_t* l,
2435  int64_t lsize,
2436  double* lbounds,
2437  int64_t lbounds_size,
2438  int32_t li,
2439  int32_t ic1,
2440  int32_t isr1,
2441  int32_t ic2,
2442  int32_t isr2,
2443  int32_t osr) {
2444  double px = coord_x(p, 0, ic1, isr1, osr);
2445  double py = coord_y(p, 1, ic1, isr1, osr);
2446 
2447  if (lbounds) {
2448  if (tol_eq(px, lbounds[0]) && tol_eq(py, lbounds[1]) && tol_eq(px, lbounds[2]) &&
2449  tol_eq(py, lbounds[3])) {
2450  return true;
2451  }
2452  }
2453 
2454  auto l_num_coords = lsize / compression_unit_size(ic2);
2455  for (int i = 0; i < l_num_coords; i += 2) {
2456  double lx = coord_x(l, i, ic2, isr2, osr);
2457  double ly = coord_y(l, i + 1, ic2, isr2, osr);
2458  if (tol_eq(px, lx) && tol_eq(py, ly)) {
2459  continue;
2460  }
2461  return false;
2462  }
2463  return true;
2464 }
2465 
2468  int64_t psize,
2469  int8_t* poly_coords,
2470  int64_t poly_coords_size,
2471  int32_t* poly_ring_sizes,
2472  int64_t poly_num_rings,
2473  double* poly_bounds,
2474  int64_t poly_bounds_size,
2475  int32_t ic1,
2476  int32_t isr1,
2477  int32_t ic2,
2478  int32_t isr2,
2479  int32_t osr) {
2480  auto exterior_ring_num_coords = poly_coords_size / compression_unit_size(ic2);
2481  if (poly_num_rings > 0) {
2482  exterior_ring_num_coords = poly_ring_sizes[0] * 2;
2483  }
2484  auto exterior_ring_coords_size = exterior_ring_num_coords * compression_unit_size(ic2);
2485 
2487  psize,
2488  poly_coords,
2489  exterior_ring_coords_size,
2490  poly_bounds,
2491  poly_bounds_size,
2492  0,
2493  ic1,
2494  isr1,
2495  ic2,
2496  isr2,
2497  osr);
2498 }
2499 
2502  int64_t lsize,
2503  double* lbounds,
2504  int64_t lbounds_size,
2505  int32_t li,
2506  int8_t* p,
2507  int64_t psize,
2508  int32_t ic1,
2509  int32_t isr1,
2510  int32_t ic2,
2511  int32_t isr2,
2512  int32_t osr) {
2513  return tol_zero(
2514  ST_Distance_Point_LineString(p, psize, l, lsize, li, ic2, isr2, ic1, isr1, osr));
2515 }
2516 
2519  int64_t l1size,
2520  double* l1bounds,
2521  int64_t l1bounds_size,
2522  int32_t l1i,
2523  int8_t* l2,
2524  int64_t l2size,
2525  double* l2bounds,
2526  int64_t l2bounds_size,
2527  int32_t l2i,
2528  int32_t ic1,
2529  int32_t isr1,
2530  int32_t ic2,
2531  int32_t isr2,
2532  int32_t osr) {
2533  if (l1i != 0 || l2i != 0) {
2534  // At least one linestring is indexed, can rely on distance
2536  l1, l1size, l1i, l2, l2size, l2i, ic1, isr1, ic2, isr2, osr));
2537  }
2538 
2539  // TODO: sublinestring
2540  // For each line segment in l2 check if there is a segment in l1
2541  // that it's colinear with and both l2 vertices are on l1 segment.
2542  // Bail if any line segment deviates from the path.
2543  return false;
2544 }
2545 
2548  int64_t lsize,
2549  double* lbounds,
2550  int64_t lbounds_size,
2551  int32_t li,
2552  int8_t* poly_coords,
2553  int64_t poly_coords_size,
2554  int32_t* poly_ring_sizes,
2555  int64_t poly_num_rings,
2556  double* poly_bounds,
2557  int64_t poly_bounds_size,
2558  int32_t ic1,
2559  int32_t isr1,
2560  int32_t ic2,
2561  int32_t isr2,
2562  int32_t osr) {
2563  // TODO
2564  return false;
2565 }
2566 
2568 bool ST_Contains_Polygon_Point(int8_t* poly_coords,
2569  int64_t poly_coords_size,
2570  int32_t* poly_ring_sizes,
2571  int64_t poly_num_rings,
2572  double* poly_bounds,
2573  int64_t poly_bounds_size,
2574  int8_t* p,
2575  int64_t psize,
2576  int32_t ic1,
2577  int32_t isr1,
2578  int32_t ic2,
2579  int32_t isr2,
2580  int32_t osr) {
2581  double px = coord_x(p, 0, ic2, isr2, osr);
2582  double py = coord_y(p, 1, ic2, isr2, osr);
2583 
2584  if (poly_bounds) {
2585  if (!box_contains_point(poly_bounds, poly_bounds_size, px, py)) {
2586  return false;
2587  }
2588  }
2589 
2590  auto poly_num_coords = poly_coords_size / compression_unit_size(ic1);
2591  auto exterior_ring_num_coords = poly_num_coords;
2592  if (poly_num_rings > 0) {
2593  exterior_ring_num_coords = poly_ring_sizes[0] * 2;
2594  }
2595 
2596  auto poly = poly_coords;
2597  if (polygon_contains_point(poly, exterior_ring_num_coords, px, py, ic1, isr1, osr)) {
2598  // Inside exterior ring
2599  poly += exterior_ring_num_coords * compression_unit_size(ic1);
2600  // Check that none of the polygon's holes contain that point
2601  for (auto r = 1; r < poly_num_rings; r++) {
2602  int64_t interior_ring_num_coords = poly_ring_sizes[r] * 2;
2604  poly, interior_ring_num_coords, px, py, ic1, isr1, osr)) {
2605  return false;
2606  }
2607  poly += interior_ring_num_coords * compression_unit_size(ic1);
2608  }
2609  return true;
2610  }
2611  return false;
2612 }
2613 
2615 bool ST_Contains_Polygon_LineString(int8_t* poly_coords,
2616  int64_t poly_coords_size,
2617  int32_t* poly_ring_sizes,
2618  int64_t poly_num_rings,
2619  double* poly_bounds,
2620  int64_t poly_bounds_size,
2621  int8_t* l,
2622  int64_t lsize,
2623  double* lbounds,
2624  int64_t lbounds_size,
2625  int32_t li,
2626  int32_t ic1,
2627  int32_t isr1,
2628  int32_t ic2,
2629  int32_t isr2,
2630  int32_t osr) {
2631  if (poly_num_rings > 1) {
2632  return false; // TODO: support polygons with interior rings
2633  }
2634 
2635  auto poly_num_coords = poly_coords_size / compression_unit_size(ic1);
2636  auto lnum_coords = lsize / compression_unit_size(ic2);
2637  auto lnum_points = lnum_coords / 2;
2638  if (li != 0) {
2639  // Statically indexed linestring
2640  if (li < 0 || li > lnum_points) {
2641  li = lnum_points;
2642  }
2643  double lx = coord_x(l, 2 * (li - 1), ic2, isr2, osr);
2644  double ly = coord_y(l, 2 * (li - 1) + 1, ic2, isr2, osr);
2645 
2646  if (poly_bounds) {
2647  if (!box_contains_point(poly_bounds, poly_bounds_size, lx, ly)) {
2648  return false;
2649  }
2650  }
2651  return polygon_contains_point(poly_coords, poly_num_coords, lx, ly, ic1, isr1, osr);
2652  }
2653 
2654  // Bail out if poly bounding box doesn't contain linestring bounding box
2655  if (poly_bounds && lbounds) {
2656  if (!box_contains_box(poly_bounds, poly_bounds_size, lbounds, lbounds_size)) {
2657  return false;
2658  }
2659  }
2660 
2662  poly_coords, poly_num_coords, l, lnum_coords, ic1, isr1, ic2, isr2, osr);
2663 }
2664 
2666 bool ST_Contains_Polygon_Polygon(int8_t* poly1_coords,
2667  int64_t poly1_coords_size,
2668  int32_t* poly1_ring_sizes,
2669  int64_t poly1_num_rings,
2670  double* poly1_bounds,
2671  int64_t poly1_bounds_size,
2672  int8_t* poly2_coords,
2673  int64_t poly2_coords_size,
2674  int32_t* poly2_ring_sizes,
2675  int64_t poly2_num_rings,
2676  double* poly2_bounds,
2677  int64_t poly2_bounds_size,
2678  int32_t ic1,
2679  int32_t isr1,
2680  int32_t ic2,
2681  int32_t isr2,
2682  int32_t osr) {
2683  // TODO: needs to be extended, cover more cases
2684  // Right now only checking if simple poly1 (no holes) contains poly2's exterior shape
2685  if (poly1_num_rings > 1) {
2686  return false; // TODO: support polygons with interior rings
2687  }
2688 
2689  if (poly1_bounds && poly2_bounds) {
2690  if (!box_contains_box(
2691  poly1_bounds, poly1_bounds_size, poly2_bounds, poly2_bounds_size)) {
2692  return false;
2693  }
2694  }
2695 
2696  int64_t poly2_exterior_ring_coords_size = poly2_coords_size;
2697  if (poly2_num_rings > 0) {
2698  poly2_exterior_ring_coords_size =
2699  2 * poly2_ring_sizes[0] * compression_unit_size(ic2);
2700  }
2701  return ST_Contains_Polygon_LineString(poly1_coords,
2702  poly1_coords_size,
2703  poly1_ring_sizes,
2704  poly1_num_rings,
2705  poly1_bounds,
2706  poly1_bounds_size,
2707  poly2_coords,
2708  poly2_exterior_ring_coords_size,
2709  poly2_bounds,
2710  poly2_bounds_size,
2711  0,
2712  ic1,
2713  isr1,
2714  ic2,
2715  isr2,
2716  osr);
2717 }
2718 
2720 bool ST_Contains_MultiPolygon_Point(int8_t* mpoly_coords,
2721  int64_t mpoly_coords_size,
2722  int32_t* mpoly_ring_sizes,
2723  int64_t mpoly_num_rings,
2724  int32_t* mpoly_poly_sizes,
2725  int64_t mpoly_num_polys,
2726  double* mpoly_bounds,
2727  int64_t mpoly_bounds_size,
2728  int8_t* p,
2729  int64_t psize,
2730  int32_t ic1,
2731  int32_t isr1,
2732  int32_t ic2,
2733  int32_t isr2,
2734  int32_t osr) {
2735  if (mpoly_num_polys <= 0) {
2736  return false;
2737  }
2738 
2739  double px = coord_x(p, 0, ic2, isr2, osr);
2740  double py = coord_y(p, 1, ic2, isr2, osr);
2741 
2742  // TODO: mpoly_bounds could contain individual bounding boxes too:
2743  // first two points - box for the entire multipolygon, then a pair for each polygon
2744  if (mpoly_bounds) {
2745  if (!box_contains_point(mpoly_bounds, mpoly_bounds_size, px, py)) {
2746  return false;
2747  }
2748  }
2749 
2750  // Set specific poly pointers as we move through the coords/ringsizes/polyrings arrays.
2751  auto next_poly_coords = mpoly_coords;
2752  auto next_poly_ring_sizes = mpoly_ring_sizes;
2753 
2754  for (auto poly = 0; poly < mpoly_num_polys; poly++) {
2755  auto poly_coords = next_poly_coords;
2756  auto poly_ring_sizes = next_poly_ring_sizes;
2757  auto poly_num_rings = mpoly_poly_sizes[poly];
2758  // Count number of coords in all of poly's rings, advance ring size pointer.
2759  int32_t poly_num_coords = 0;
2760  for (auto ring = 0; ring < poly_num_rings; ring++) {
2761  poly_num_coords += 2 * *next_poly_ring_sizes++;
2762  }
2763  auto poly_coords_size = poly_num_coords * compression_unit_size(ic1);
2764  next_poly_coords += poly_coords_size;
2765  // TODO: pass individual bounding boxes for each polygon
2766  if (ST_Contains_Polygon_Point(poly_coords,
2767  poly_coords_size,
2768  poly_ring_sizes,
2769  poly_num_rings,
2770  nullptr,
2771  0,
2772  p,
2773  psize,
2774  ic1,
2775  isr1,
2776  ic2,
2777  isr2,
2778  osr)) {
2779  return true;
2780  }
2781  }
2782 
2783  return false;
2784 }
2785 
2787 bool ST_Contains_MultiPolygon_LineString(int8_t* mpoly_coords,
2788  int64_t mpoly_coords_size,
2789  int32_t* mpoly_ring_sizes,
2790  int64_t mpoly_num_rings,
2791  int32_t* mpoly_poly_sizes,
2792  int64_t mpoly_num_polys,
2793  double* mpoly_bounds,
2794  int64_t mpoly_bounds_size,
2795  int8_t* l,
2796  int64_t lsize,
2797  double* lbounds,
2798  int64_t lbounds_size,
2799  int32_t li,
2800  int32_t ic1,
2801  int32_t isr1,
2802  int32_t ic2,
2803  int32_t isr2,
2804  int32_t osr) {
2805  if (mpoly_num_polys <= 0) {
2806  return false;
2807  }
2808 
2809  auto lnum_coords = lsize / compression_unit_size(ic2);
2810  auto lnum_points = lnum_coords / 2;
2811  if (li != 0) {
2812  // Statically indexed linestring
2813  if (li < 0 || li > lnum_points) {
2814  li = lnum_points;
2815  }
2816  double lx = coord_x(l, 2 * (li - 1), ic2, isr2, osr);
2817  double ly = coord_y(l, 2 * (li - 1) + 1, ic2, isr2, osr);
2818 
2819  if (mpoly_bounds) {
2820  if (!box_contains_point(mpoly_bounds, mpoly_bounds_size, lx, ly)) {
2821  return false;
2822  }
2823  }
2824  auto p = l + li * compression_unit_size(ic2);
2825  auto psize = 2 * compression_unit_size(ic2);
2826  return ST_Contains_MultiPolygon_Point(mpoly_coords,
2827  mpoly_coords_size,
2828  mpoly_ring_sizes,
2829  mpoly_num_rings,
2830  mpoly_poly_sizes,
2831  mpoly_num_polys,
2832  mpoly_bounds,
2833  mpoly_bounds_size,
2834  p,
2835  psize,
2836  ic1,
2837  isr1,
2838  ic2,
2839  isr2,
2840  osr);
2841  }
2842 
2843  if (mpoly_bounds && lbounds) {
2844  if (!box_contains_box(mpoly_bounds, mpoly_bounds_size, lbounds, lbounds_size)) {
2845  return false;
2846  }
2847  }
2848 
2849  // Set specific poly pointers as we move through the coords/ringsizes/polyrings arrays.
2850  auto next_poly_coords = mpoly_coords;
2851  auto next_poly_ring_sizes = mpoly_ring_sizes;
2852 
2853  for (auto poly = 0; poly < mpoly_num_polys; poly++) {
2854  auto poly_coords = next_poly_coords;
2855  auto poly_ring_sizes = next_poly_ring_sizes;
2856  auto poly_num_rings = mpoly_poly_sizes[poly];
2857  // Count number of coords in all of poly's rings, advance ring size pointer.
2858  int32_t poly_num_coords = 0;
2859  for (auto ring = 0; ring < poly_num_rings; ring++) {
2860  poly_num_coords += 2 * *next_poly_ring_sizes++;
2861  }
2862  auto poly_coords_size = poly_num_coords * compression_unit_size(ic1);
2863  next_poly_coords += poly_coords_size;
2864 
2865  if (ST_Contains_Polygon_LineString(poly_coords,
2866  poly_coords_size,
2867  poly_ring_sizes,
2868  poly_num_rings,
2869  nullptr,
2870  0,
2871  l,
2872  lsize,
2873  nullptr,
2874  0,
2875  li,
2876  ic1,
2877  isr1,
2878  ic2,
2879  isr2,
2880  osr)) {
2881  return true;
2882  }
2883  }
2884 
2885  return false;
2886 }
2887 
2888 //
2889 // ST_Intersects
2890 //
2891 
2894  int64_t p1size,
2895  int8_t* p2,
2896  int64_t p2size,
2897  int32_t ic1,
2898  int32_t isr1,
2899  int32_t ic2,
2900  int32_t isr2,
2901  int32_t osr) {
2902  return tol_zero(
2903  ST_Distance_Point_Point(p1, p1size, p2, p2size, ic1, isr1, ic2, isr2, osr));
2904 }
2905 
2908  int64_t psize,
2909  int8_t* l,
2910  int64_t lsize,
2911  double* lbounds,
2912  int64_t lbounds_size,
2913  int32_t li,
2914  int32_t ic1,
2915  int32_t isr1,
2916  int32_t ic2,
2917  int32_t isr2,
2918  int32_t osr) {
2919  double px = coord_x(p, 0, ic1, isr1, osr);
2920  double py = coord_y(p, 1, ic1, isr1, osr);
2921 
2922  auto lnum_coords = lsize / compression_unit_size(ic2);
2923  auto lnum_points = lnum_coords / 2;
2924  if (li != 0) {
2925  // Statically indexed linestring
2926  if (li < 0 || li > lnum_points) {
2927  li = lnum_points;
2928  }
2929  auto p2 = l + li * compression_unit_size(ic2);
2930  auto p2size = 2 * compression_unit_size(ic2);
2931  return tol_zero(
2932  ST_Distance_Point_Point(p2, p2size, p, psize, ic2, isr2, ic1, isr1, osr));
2933  }
2934 
2935  if (lbounds) {
2936  if (!box_contains_point(lbounds, lbounds_size, px, py)) {
2937  return false;
2938  }
2939  }
2940  return tol_zero(
2941  ST_Distance_Point_LineString(p, psize, l, lsize, li, ic1, isr1, ic2, isr2, osr));
2942 }
2943 
2946  int64_t psize,
2947  int8_t* poly,
2948  int64_t polysize,
2949  int32_t* poly_ring_sizes,
2950  int64_t poly_num_rings,
2951  double* poly_bounds,
2952  int64_t poly_bounds_size,
2953  int32_t ic1,
2954  int32_t isr1,
2955  int32_t ic2,
2956  int32_t isr2,
2957  int32_t osr) {
2958  return ST_Contains_Polygon_Point(poly,
2959  polysize,
2960  poly_ring_sizes,
2961  poly_num_rings,
2962  poly_bounds,
2963  poly_bounds_size,
2964  p,
2965  psize,
2966  ic2,
2967  isr2,
2968  ic1,
2969  isr1,
2970  osr);
2971 }
2972 
2975  int64_t psize,
2976  int8_t* mpoly_coords,
2977  int64_t mpoly_coords_size,
2978  int32_t* mpoly_ring_sizes,
2979  int64_t mpoly_num_rings,
2980  int32_t* mpoly_poly_sizes,
2981  int64_t mpoly_num_polys,
2982  double* mpoly_bounds,
2983  int64_t mpoly_bounds_size,
2984  int32_t ic1,
2985  int32_t isr1,
2986  int32_t ic2,
2987  int32_t isr2,
2988  int32_t osr) {
2989  return ST_Contains_MultiPolygon_Point(mpoly_coords,
2990  mpoly_coords_size,
2991  mpoly_ring_sizes,
2992  mpoly_num_rings,
2993  mpoly_poly_sizes,
2994  mpoly_num_polys,
2995  mpoly_bounds,
2996  mpoly_bounds_size,
2997  p,
2998  psize,
2999  ic2,
3000  isr2,
3001  ic1,
3002  isr1,
3003  osr);
3004 }
3005 
3008  int64_t lsize,
3009  double* lbounds,
3010  int64_t lbounds_size,
3011  int32_t li,
3012  int8_t* p,
3013  int64_t psize,
3014  int32_t ic1,
3015  int32_t isr1,
3016  int32_t ic2,
3017  int32_t isr2,
3018  int32_t osr) {
3020  p, psize, l, lsize, lbounds, lbounds_size, li, ic2, isr2, ic1, isr1, osr);
3021 }
3022 
3025  int64_t l1size,
3026  double* l1bounds,
3027  int64_t l1bounds_size,
3028  int32_t l1i,
3029  int8_t* l2,
3030  int64_t l2size,
3031  double* l2bounds,
3032  int64_t l2bounds_size,
3033  int32_t l2i,
3034  int32_t ic1,
3035  int32_t isr1,
3036  int32_t ic2,
3037  int32_t isr2,
3038  int32_t osr) {
3039  auto l2num_coords = l2size / compression_unit_size(ic2);
3040  auto l2num_points = l2num_coords / 2;
3041  if (l2i != 0) {
3042  // Statically indexed linestring
3043  if (l2i < 0 || l2i > l2num_points) {
3044  l2i = l2num_points;
3045  }
3046  auto p2 = l2 + l2i * compression_unit_size(ic2);
3047  auto p2size = 2 * compression_unit_size(ic2);
3049  l1, l1size, l1bounds, l1bounds_size, l1i, p2, p2size, ic1, isr1, ic2, isr2, osr);
3050  }
3051  auto l1num_coords = l1size / compression_unit_size(ic1);
3052  auto l1num_points = l1num_coords / 2;
3053  if (l1i != 0) {
3054  // Statically indexed linestring
3055  if (l1i < 0 || l1i > l1num_points) {
3056  l1i = l1num_points;
3057  }
3058  auto p1 = l1 + l1i * compression_unit_size(ic1);
3059  auto p1size = 2 * compression_unit_size(ic1);
3061  l2, l2size, l2bounds, l2bounds_size, l2i, p1, p1size, ic2, isr2, ic1, isr1, osr);
3062  }
3063 
3064  if (l1bounds && l2bounds) {
3065  if (!box_overlaps_box(l1bounds, l1bounds_size, l2bounds, l2bounds_size)) {
3066  return false;
3067  }
3068  }
3069 
3071  l1, l1size, l1i, l2, l2size, l2i, ic1, isr1, ic2, isr2, osr));
3072 }
3073 
3076  int64_t lsize,
3077  double* lbounds,
3078  int64_t lbounds_size,
3079  int32_t li,
3080  int8_t* poly,
3081  int64_t polysize,
3082  int32_t* poly_ring_sizes,
3083  int64_t poly_num_rings,
3084  double* poly_bounds,
3085  int64_t poly_bounds_size,
3086  int32_t ic1,
3087  int32_t isr1,
3088  int32_t ic2,
3089  int32_t isr2,
3090  int32_t osr) {
3091  auto lnum_coords = lsize / compression_unit_size(ic1);
3092  auto lnum_points = lnum_coords / 2;
3093  if (li != 0) {
3094  // Statically indexed linestring
3095  if (li < 0 || li > lnum_points) {
3096  li = lnum_points;
3097  }
3098  auto p = l + li * compression_unit_size(ic1);
3099  auto psize = 2 * compression_unit_size(ic1);
3100  return ST_Contains_Polygon_Point(poly,
3101  polysize,
3102  poly_ring_sizes,
3103  poly_num_rings,
3104  poly_bounds,
3105  poly_bounds_size,
3106  p,
3107  psize,
3108  ic2,
3109  isr2,
3110  ic1,
3111  isr1,
3112  osr);
3113  }
3114 
3115  if (lbounds && poly_bounds) {
3116  if (!box_overlaps_box(lbounds, lbounds_size, poly_bounds, poly_bounds_size)) {
3117  return false;
3118  }
3119  }
3120 
3121  // Check for spatial intersection.
3122  // One way to do that would be to start with linestring's first point, if it's inside
3123  // the polygon - it means intersection. Otherwise follow the linestring, segment by
3124  // segment, checking for intersections with polygon rings, bail as soon as we cross into
3125  // the polygon.
3126 
3127  // Or, alternatively, just measure the distance:
3129  lsize,
3130  li,
3131  poly,
3132  polysize,
3133  poly_ring_sizes,
3134  poly_num_rings,
3135  ic1,
3136  isr1,
3137  ic2,
3138  isr2,
3139  osr));
3140 }
3141 
3144  int64_t lsize,
3145  double* lbounds,
3146  int64_t lbounds_size,
3147  int32_t li,
3148  int8_t* mpoly_coords,
3149  int64_t mpoly_coords_size,
3150  int32_t* mpoly_ring_sizes,
3151  int64_t mpoly_num_rings,
3152  int32_t* mpoly_poly_sizes,
3153  int64_t mpoly_num_polys,
3154  double* mpoly_bounds,
3155  int64_t mpoly_bounds_size,
3156  int32_t ic1,
3157  int32_t isr1,
3158  int32_t ic2,
3159  int32_t isr2,
3160  int32_t osr) {
3161  auto lnum_coords = lsize / compression_unit_size(ic1);
3162  auto lnum_points = lnum_coords / 2;
3163  if (li != 0) {
3164  // Statically indexed linestring
3165  if (li < 0 || li > lnum_points) {
3166  li = lnum_points;
3167  }
3168  auto p = l + li * compression_unit_size(ic1);
3169  auto psize = 2 * compression_unit_size(ic1);
3170  return ST_Contains_MultiPolygon_Point(mpoly_coords,
3171  mpoly_coords_size,
3172  mpoly_ring_sizes,
3173  mpoly_num_rings,
3174  mpoly_poly_sizes,
3175  mpoly_num_polys,
3176  mpoly_bounds,
3177  mpoly_bounds_size,
3178  p,
3179  psize,
3180  ic2,
3181  isr2,
3182  ic1,
3183  isr1,
3184  osr);
3185  }
3186 
3187  if (lbounds && mpoly_bounds) {
3188  if (!box_overlaps_box(lbounds, lbounds_size, mpoly_bounds, mpoly_bounds_size)) {
3189  return false;
3190  }
3191  }
3192 
3193  // Check for spatial intersection.
3194  // One way to do that would be to start with linestring's first point, if it's inside
3195  // any of the polygons - it means intersection. Otherwise follow the linestring, segment
3196  // by segment, checking for intersections with polygon shapes/holes, bail as soon as we
3197  // cross into a polygon.
3198 
3199  // Or, alternatively, just measure the distance:
3201  lsize,
3202  li,
3203  mpoly_coords,
3204  mpoly_coords_size,
3205  mpoly_ring_sizes,
3206  mpoly_num_rings,
3207  mpoly_poly_sizes,
3208  mpoly_num_polys,
3209  ic1,
3210  isr1,
3211  ic2,
3212  isr2,
3213  osr));
3214 }
3215 
3218  int64_t polysize,
3219  int32_t* poly_ring_sizes,
3220  int64_t poly_num_rings,
3221  double* poly_bounds,
3222  int64_t poly_bounds_size,
3223  int8_t* p,
3224  int64_t psize,
3225  int32_t ic1,
3226  int32_t isr1,
3227  int32_t ic2,
3228  int32_t isr2,
3229  int32_t osr) {
3230  return ST_Contains_Polygon_Point(poly,
3231  polysize,
3232  poly_ring_sizes,
3233  poly_num_rings,
3234  poly_bounds,
3235  poly_bounds_size,
3236  p,
3237  psize,
3238  ic1,
3239  isr1,
3240  ic2,
3241  isr2,
3242  osr);
3243 }
3244 
3247  int64_t polysize,
3248  int32_t* poly_ring_sizes,
3249  int64_t poly_num_rings,
3250  double* poly_bounds,
3251  int64_t poly_bounds_size,
3252  int8_t* l,
3253  int64_t lsize,
3254  double* lbounds,
3255  int64_t lbounds_size,
3256  int32_t li,
3257  int32_t ic1,
3258  int32_t isr1,
3259  int32_t ic2,
3260  int32_t isr2,
3261  int32_t osr) {
3263  lsize,
3264  lbounds,
3265  lbounds_size,
3266  li,
3267  poly,
3268  polysize,
3269  poly_ring_sizes,
3270  poly_num_rings,
3271  poly_bounds,
3272  poly_bounds_size,
3273  ic2,
3274  isr2,
3275  ic1,
3276  isr1,
3277  osr);
3278 }
3279 
3281 bool ST_Intersects_Polygon_Polygon(int8_t* poly1_coords,
3282  int64_t poly1_coords_size,
3283  int32_t* poly1_ring_sizes,
3284  int64_t poly1_num_rings,
3285  double* poly1_bounds,
3286  int64_t poly1_bounds_size,
3287  int8_t* poly2_coords,
3288  int64_t poly2_coords_size,
3289  int32_t* poly2_ring_sizes,
3290  int64_t poly2_num_rings,
3291  double* poly2_bounds,
3292  int64_t poly2_bounds_size,
3293  int32_t ic1,
3294  int32_t isr1,
3295  int32_t ic2,
3296  int32_t isr2,
3297  int32_t osr) {
3298  if (poly1_bounds && poly2_bounds) {
3299  if (!box_overlaps_box(
3300  poly1_bounds, poly1_bounds_size, poly2_bounds, poly2_bounds_size)) {
3301  return false;
3302  }
3303  }
3304 
3305  return tol_zero(ST_Distance_Polygon_Polygon(poly1_coords,
3306  poly1_coords_size,
3307  poly1_ring_sizes,
3308  poly1_num_rings,
3309  poly2_coords,
3310  poly2_coords_size,
3311  poly2_ring_sizes,
3312  poly2_num_rings,
3313  ic1,
3314  isr1,
3315  ic2,
3316  isr2,
3317  osr));
3318 }
3319 
3321 bool ST_Intersects_Polygon_MultiPolygon(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  int8_t* mpoly_coords,
3328  int64_t mpoly_coords_size,
3329  int32_t* mpoly_ring_sizes,
3330  int64_t mpoly_num_rings,
3331  int32_t* mpoly_poly_sizes,
3332  int64_t mpoly_num_polys,
3333  double* mpoly_bounds,
3334  int64_t mpoly_bounds_size,
3335  int32_t ic1,
3336  int32_t isr1,
3337  int32_t ic2,
3338  int32_t isr2,
3339  int32_t osr) {
3340  if (poly_bounds && mpoly_bounds) {
3341  if (!box_overlaps_box(
3342  poly_bounds, poly_bounds_size, mpoly_bounds, mpoly_bounds_size)) {
3343  return false;
3344  }
3345  }
3346 
3347  return tol_zero(ST_Distance_Polygon_MultiPolygon(poly_coords,
3348  poly_coords_size,
3349  poly_ring_sizes,
3350  poly_num_rings,
3351  mpoly_coords,
3352  mpoly_coords_size,
3353  mpoly_ring_sizes,
3354  mpoly_num_rings,
3355  mpoly_poly_sizes,
3356  mpoly_num_polys,
3357  ic1,
3358  isr1,
3359  ic2,
3360  isr2,
3361  osr));
3362 }
3363 
3365 bool ST_Intersects_MultiPolygon_Point(int8_t* mpoly_coords,
3366  int64_t mpoly_coords_size,
3367  int32_t* mpoly_ring_sizes,
3368  int64_t mpoly_num_rings,
3369  int32_t* mpoly_poly_sizes,
3370  int64_t mpoly_num_polys,
3371  double* mpoly_bounds,
3372  int64_t mpoly_bounds_size,
3373  int8_t* p,
3374  int64_t psize,
3375  int32_t ic1,
3376  int32_t isr1,
3377  int32_t ic2,
3378  int32_t isr2,
3379  int32_t osr) {
3380  return ST_Contains_MultiPolygon_Point(mpoly_coords,
3381  mpoly_coords_size,
3382  mpoly_ring_sizes,
3383  mpoly_num_rings,
3384  mpoly_poly_sizes,
3385  mpoly_num_polys,
3386  mpoly_bounds,
3387  mpoly_bounds_size,
3388  p,
3389  psize,
3390  ic1,
3391  isr1,
3392  ic2,
3393  isr2,
3394  osr);
3395 }
3396 
3398 bool ST_Intersects_MultiPolygon_LineString(int8_t* mpoly_coords,
3399  int64_t mpoly_coords_size,
3400  int32_t* mpoly_ring_sizes,
3401  int64_t mpoly_num_rings,
3402  int32_t* mpoly_poly_sizes,
3403  int64_t mpoly_num_polys,
3404  double* mpoly_bounds,
3405  int64_t mpoly_bounds_size,
3406  int8_t* l,
3407  int64_t lsize,
3408  double* lbounds,
3409  int64_t lbounds_size,
3410  int32_t li,
3411  int32_t ic1,
3412  int32_t isr1,
3413  int32_t ic2,
3414  int32_t isr2,
3415  int32_t osr) {
3417  lsize,
3418  lbounds,
3419  lbounds_size,
3420  li,
3421  mpoly_coords,
3422  mpoly_coords_size,
3423  mpoly_ring_sizes,
3424  mpoly_num_rings,
3425  mpoly_poly_sizes,
3426  mpoly_num_polys,
3427  mpoly_bounds,
3428  mpoly_bounds_size,
3429  ic2,
3430  isr2,
3431  ic1,
3432  isr1,
3433  osr);
3434 }
3435 
3437 bool ST_Intersects_MultiPolygon_Polygon(int8_t* mpoly_coords,
3438  int64_t mpoly_coords_size,
3439  int32_t* mpoly_ring_sizes,
3440  int64_t mpoly_num_rings,
3441  int32_t* mpoly_poly_sizes,
3442  int64_t mpoly_num_polys,
3443  double* mpoly_bounds,
3444  int64_t mpoly_bounds_size,
3445  int8_t* poly_coords,
3446  int64_t poly_coords_size,
3447  int32_t* poly_ring_sizes,
3448  int64_t poly_num_rings,
3449  double* poly_bounds,
3450  int64_t poly_bounds_size,
3451  int32_t ic1,
3452  int32_t isr1,
3453  int32_t ic2,
3454  int32_t isr2,
3455  int32_t osr) {
3456  return ST_Intersects_Polygon_MultiPolygon(poly_coords,
3457  poly_coords_size,
3458  poly_ring_sizes,
3459  poly_num_rings,
3460  poly_bounds,
3461  poly_bounds_size,
3462  mpoly_coords,
3463  mpoly_coords_size,
3464  mpoly_ring_sizes,
3465  mpoly_num_rings,
3466  mpoly_poly_sizes,
3467  mpoly_num_polys,
3468  mpoly_bounds,
3469  mpoly_bounds_size,
3470  ic2,
3471  isr2,
3472  ic1,
3473  isr1,
3474  osr);
3475 }
3476 
3478 bool ST_Intersects_MultiPolygon_MultiPolygon(int8_t* mpoly1_coords,
3479  int64_t mpoly1_coords_size,
3480  int32_t* mpoly1_ring_sizes,
3481  int64_t mpoly1_num_rings,
3482  int32_t* mpoly1_poly_sizes,
3483  int64_t mpoly1_num_polys,
3484  double* mpoly1_bounds,
3485  int64_t mpoly1_bounds_size,
3486  int8_t* mpoly2_coords,
3487  int64_t mpoly2_coords_size,
3488  int32_t* mpoly2_ring_sizes,
3489  int64_t mpoly2_num_rings,
3490  int32_t* mpoly2_poly_sizes,
3491  int64_t mpoly2_num_polys,
3492  double* mpoly2_bounds,
3493  int64_t mpoly2_bounds_size,
3494  int32_t ic1,
3495  int32_t isr1,
3496  int32_t ic2,
3497  int32_t isr2,
3498  int32_t osr) {
3499  if (mpoly1_bounds && mpoly2_bounds) {
3500  if (!box_overlaps_box(
3501  mpoly1_bounds, mpoly1_bounds_size, mpoly2_bounds, mpoly2_bounds_size)) {
3502  return false;
3503  }
3504  }
3505 
3506  return tol_zero(ST_Distance_MultiPolygon_MultiPolygon(mpoly1_coords,
3507  mpoly1_coords_size,
3508  mpoly1_ring_sizes,
3509  mpoly1_num_rings,
3510  mpoly1_poly_sizes,
3511  mpoly1_num_polys,
3512  mpoly2_coords,
3513  mpoly2_coords_size,
3514  mpoly2_ring_sizes,
3515  mpoly2_num_rings,
3516  mpoly2_poly_sizes,
3517  mpoly2_num_polys,
3518  ic1,
3519  isr1,
3520  ic2,
3521  isr2,
3522  osr));
3523 }
3524 
3525 //
3526 // Accessors for poly bounds and render group for in-situ poly render queries
3527 //
3528 // The MapD_* varieties are deprecated and renamed to "OmniSci_Geo*"
3529 // There may be some clients out there who are playing with the MapD_* so leaving
3530 // them for backwards compatibility.
3531 //
3532 
3534 int64_t OmniSci_Geo_PolyBoundsPtr(double* bounds, int64_t size) {
3535  return reinterpret_cast<int64_t>(bounds);
3536 }
3537 
3539 int32_t OmniSci_Geo_PolyRenderGroup(int32_t render_group) {
3540  return render_group;
3541 }
3542 
3544 int64_t MapD_GeoPolyBoundsPtr(double* bounds, int64_t size) {
3545  return OmniSci_Geo_PolyBoundsPtr(bounds, size);
3546 }
3547 
3549 int32_t MapD_GeoPolyRenderGroup(int32_t render_group) {
3550  return OmniSci_Geo_PolyRenderGroup(render_group);
3551 }
3552 
3554 double convert_meters_to_pixel_width(const double meters,
3555  int8_t* p,
3556  const int64_t psize,
3557  const int32_t ic,
3558  const int32_t isr,
3559  const int32_t osr,
3560  const double min_lon,
3561  const double max_lon,
3562  const int32_t img_width,
3563  const double min_width) {
3564  const double const1 = 0.017453292519943295769236907684886;
3565  const double const2 = 6372797.560856;
3566  const auto lon = decompress_coord(p, 0, ic, true);
3567  const auto lat = decompress_coord(p, 1, ic, false);
3568  double t1 = sinf(meters / (2.0 * const2));
3569  double t2 = cosf(const1 * lat);
3570  const double newlon = lon - (2.0 * asinf(t1 / t2)) / const1;
3571  t1 = transform_coord(lon, isr, osr, true);
3572  t2 = transform_coord(newlon, isr, osr, true);
3573  const double min_domain_x = transform_coord(min_lon, isr, osr, true);
3574  const double max_domain_x = transform_coord(max_lon, isr, osr, true);
3575  const double domain_diff = max_domain_x - min_domain_x;
3576  t1 = ((t1 - min_domain_x) / domain_diff) * static_cast<double>(img_width);
3577  t2 = ((t2 - min_domain_x) / domain_diff) * static_cast<double>(img_width);
3578 
3579  // TODO(croot): need to account for edge cases, such as getting close to the poles.
3580  const double sz = fabs(t1 - t2);
3581  return (sz < min_width ? min_width : sz);
3582 }
3583 
3585 double convert_meters_to_pixel_height(const double meters,
3586  int8_t* p,
3587  const int64_t psize,
3588  const int32_t ic,
3589  const int32_t isr,
3590  const int32_t osr,
3591  const double min_lat,
3592  const double max_lat,
3593  const int32_t img_height,
3594  const double min_height) {
3595  const double const1 = 0.017453292519943295769236907684886;
3596  const double const2 = 6372797.560856;
3597  const auto lat = decompress_coord(p, 1, ic, false);
3598  const double latdiff = meters / (const1 * const2);
3599  const double newlat =
3600  (lat < 0) ? lat + latdiff : lat - latdiff; // assumes a lat range of [-90, 90]
3601  double t1 = transform_coord(lat, isr, osr, false);
3602  double t2 = transform_coord(newlat, isr, osr, false);
3603  const double min_domain_y = transform_coord(min_lat, isr, osr, false);
3604  const double max_domain_y = transform_coord(max_lat, isr, osr, false);
3605  const double domain_diff = max_domain_y - min_domain_y;
3606  t1 = ((t1 - min_domain_y) / domain_diff) * static_cast<double>(img_height);
3607  t2 = ((t2 - min_domain_y) / domain_diff) * static_cast<double>(img_height);
3608 
3609  // TODO(croot): need to account for edge cases, such as getting close to the poles.
3610  const double sz = fabs(t1 - t2);
3611  return (sz < min_height ? min_height : sz);
3612 }
3613 
3615  const int64_t psize,
3616  const int32_t ic,
3617  const double min_lon,
3618  const double max_lon,
3619  const double min_lat,
3620  const double max_lat) {
3621  const auto lon = decompress_coord(p, 0, ic, true);
3622  const auto lat = decompress_coord(p, 1, ic, false);
3623  return !(lon < min_lon || lon > max_lon || lat < min_lat || lat > max_lat);
3624 }
3625 
3627  const int64_t psize,
3628  const int32_t ic,
3629  const double meters,
3630  const double min_lon,
3631  const double max_lon,
3632  const double min_lat,
3633  const double max_lat) {
3634  const double const1 = 0.017453292519943295769236907684886;
3635  const double const2 = 6372797.560856;
3636  const auto lon = decompress_coord(p, 0, ic, true);
3637  const auto lat = decompress_coord(p, 1, ic, false);
3638  const double latdiff = meters / (const1 * const2);
3639  const double t1 = sinf(meters / (2.0 * const2));
3640  const double t2 = cosf(const1 * lat);
3641  const double londiff = (2.0 * asinf(t1 / t2)) / const1;
3642  return !(lon + londiff < min_lon || lon - londiff > max_lon ||
3643  lat + latdiff < min_lat || lat - latdiff > max_lat);
3644 }
EXTENSION_NOINLINE double ST_Distance_LineString_LineString_Geodesic(int8_t *l1, int64_t l1size, int32_t l1index, int8_t *l2, int64_t l2size, int32_t l2index, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr)
EXTENSION_NOINLINE double convert_meters_to_pixel_width(const double meters, int8_t *p, const int64_t psize, const int32_t ic, const int32_t isr, const int32_t osr, const double min_lon, const double max_lon, const int32_t img_width, const double min_width)
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)
EXTENSION_INLINE bool ST_Intersects_Point_Polygon(int8_t *p, int64_t psize, int8_t *poly, int64_t polysize, int32_t *poly_ring_sizes, int64_t poly_num_rings, double *poly_bounds, int64_t poly_bounds_size, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr)
DEVICE ALWAYS_INLINE double coord_x(int8_t *data, int32_t index, int32_t ic, int32_t isr, int32_t osr)
DEVICE ALWAYS_INLINE double coord_y(int8_t *data, int32_t index, int32_t ic, int32_t isr, int32_t osr)
DEVICE ALWAYS_INLINE double length_linestring(int8_t *l, int64_t lsize, int32_t ic, int32_t isr, int32_t osr, bool geodesic, bool check_closed)
DEVICE ALWAYS_INLINE double transform_coord(double coord, int32_t isr, int32_t osr, bool x)
EXTENSION_NOINLINE double ST_Distance_Point_Point(int8_t *p1, int64_t p1size, int8_t *p2, int64_t p2size, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr)
EXTENSION_NOINLINE double ST_X_Point(int8_t *p, int64_t psize, int32_t ic, int32_t isr, int32_t osr)
DEVICE ALWAYS_INLINE bool on_segment(double px, double py, double qx, double qy, double rx, double ry)
EXTENSION_NOINLINE double ST_Centroid_LineString(int8_t *coords, int64_t coords_sz, int32_t ic, int32_t isr, int32_t osr, bool ycoord)
EXTENSION_INLINE 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)
EXTENSION_INLINE bool ST_Intersects_MultiPolygon_Polygon(int8_t *mpoly_coords, int64_t mpoly_coords_size, int32_t *mpoly_ring_sizes, int64_t mpoly_num_rings, int32_t *mpoly_poly_sizes, int64_t mpoly_num_polys, double *mpoly_bounds, int64_t mpoly_bounds_size, int8_t *poly_coords, int64_t poly_coords_size, int32_t *poly_ring_sizes, int64_t poly_num_rings, double *poly_bounds, int64_t poly_bounds_size, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr)
DEVICE ALWAYS_INLINE double area_polygon(int8_t *poly_coords, int64_t poly_coords_size, int32_t *poly_ring_sizes, int64_t poly_num_rings, int32_t ic, int32_t isr, int32_t osr)
EXTENSION_INLINE double ST_Distance_Polygon_LineString(int8_t *poly_coords, int64_t poly_coords_size, int32_t *poly_ring_sizes, int64_t poly_num_rings, int8_t *l, int64_t lsize, int32_t li, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr)
EXTENSION_NOINLINE bool ST_Contains_Point_Point(int8_t *p1, int64_t p1size, int8_t *p2, int64_t p2size, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr)
EXTENSION_NOINLINE double ST_Perimeter_Polygon(int8_t *poly, int64_t polysize, int32_t *poly_ring_sizes, int64_t poly_num_rings, int32_t ic, int32_t isr, int32_t osr)
EXTENSION_NOINLINE bool ST_Contains_LineString_Polygon(int8_t *l, int64_t lsize, double *lbounds, int64_t lbounds_size, int32_t li, int8_t *poly_coords, int64_t poly_coords_size, int32_t *poly_ring_sizes, int64_t poly_num_rings, double *poly_bounds, int64_t poly_bounds_size, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr)
EXTENSION_INLINE int32_t OmniSci_Geo_PolyRenderGroup(int32_t render_group)
DEVICE ALWAYS_INLINE bool centroid_add_triangle(double x1, double y1, double x2, double y2, double x3, double y3, double sign, double *total_area2, double *cg3)
EXTENSION_NOINLINE double ST_XMax(int8_t *coords, int64_t size, int32_t ic, int32_t isr, int32_t osr)
EXTENSION_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)
EXTENSION_NOINLINE bool ST_Contains_Polygon_Point(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 *p, int64_t psize, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr)
EXTENSION_NOINLINE bool ST_Intersects_LineString_MultiPolygon(int8_t *l, int64_t lsize, double *lbounds, int64_t lbounds_size, int32_t li, int8_t *mpoly_coords, int64_t mpoly_coords_size, int32_t *mpoly_ring_sizes, int64_t mpoly_num_rings, int32_t *mpoly_poly_sizes, int64_t mpoly_num_polys, double *mpoly_bounds, int64_t mpoly_bounds_size, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr)
EXTENSION_NOINLINE bool ST_Intersects_Polygon_Polygon(int8_t *poly1_coords, int64_t poly1_coords_size, int32_t *poly1_ring_sizes, int64_t poly1_num_rings, double *poly1_bounds, int64_t poly1_bounds_size, int8_t *poly2_coords, int64_t poly2_coords_size, int32_t *poly2_ring_sizes, int64_t poly2_num_rings, double *poly2_bounds, int64_t poly2_bounds_size, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr)
EXTENSION_NOINLINE bool is_point_size_in_view(int8_t *p, const int64_t psize, const int32_t ic, const double meters, const double min_lon, const double max_lon, const double min_lat, const double max_lat)
DEVICE double distance_point_line(double px, double py, double l1x, double l1y, double l2x, double l2y)
EXTENSION_NOINLINE double ST_X_LineString(int8_t *l, int64_t lsize, int32_t lindex, int32_t ic, int32_t isr, int32_t osr)
EXTENSION_INLINE double ST_YMin_Bounds(double *bounds, int64_t size, int32_t isr, int32_t osr)
EXTENSION_NOINLINE double ST_YMax(int8_t *coords, int64_t size, int32_t ic, int32_t isr, int32_t osr)
DEVICE ALWAYS_INLINE bool tol_eq(double x, double y, double tolerance=TOLERANCE_DEFAULT)
EXTENSION_INLINE double ST_XMin_Bounds(double *bounds, int64_t size, int32_t isr, int32_t osr)
DEVICE ALWAYS_INLINE double perimeter_multipolygon(int8_t *mpoly_coords, int64_t mpoly_coords_size, int32_t *mpoly_ring_sizes, int64_t mpoly_num_rings, int32_t *mpoly_poly_sizes, int64_t mpoly_num_polys, int32_t ic, int32_t isr, int32_t osr, bool geodesic)
DEVICE ALWAYS_INLINE double max_distance_point_linestring(int8_t *p, int64_t psize, int8_t *l, int64_t lsize, int32_t lindex, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr, bool check_closed)
DEVICE double decompress_lattitude_coord_geoint32(const int32_t compressed)
EXTENSION_NOINLINE double ST_Area_Polygon(int8_t *poly_coords, int64_t poly_coords_size, int32_t *poly_ring_sizes, int64_t poly_num_rings, int32_t ic, int32_t isr, int32_t osr)
EXTENSION_NOINLINE bool ST_Contains_MultiPolygon_LineString(int8_t *mpoly_coords, int64_t mpoly_coords_size, int32_t *mpoly_ring_sizes, int64_t mpoly_num_rings, int32_t *mpoly_poly_sizes, int64_t mpoly_num_polys, double *mpoly_bounds, int64_t mpoly_bounds_size, int8_t *l, int64_t lsize, double *lbounds, int64_t lbounds_size, int32_t li, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr)
DEVICE ALWAYS_INLINE double area_triangle(double x1, double y1, double x2, double y2, double x3, double y3)
EXTENSION_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)
EXTENSION_NOINLINE double ST_XMin(int8_t *coords, int64_t size, int32_t ic, int32_t isr, int32_t osr)
EXTENSION_NOINLINE bool ST_Contains_Polygon_LineString(int8_t *poly_coords, int64_t poly_coords_size, int32_t *poly_ring_sizes, int64_t poly_num_rings, double *poly_bounds, int64_t poly_bounds_size, int8_t *l, int64_t lsize, double *lbounds, int64_t lbounds_size, int32_t li, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr)
EXTENSION_INLINE bool ST_Intersects_Polygon_Point(int8_t *poly, int64_t polysize, int32_t *poly_ring_sizes, int64_t poly_num_rings, double *poly_bounds, int64_t poly_bounds_size, int8_t *p, int64_t psize, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr)
EXTENSION_NOINLINE double ST_Distance_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)
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_NOINLINE double ST_Distance_Point_LineString(int8_t *p, int64_t psize, int8_t *l, int64_t lsize, int32_t lindex, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr)
EXTENSION_NOINLINE double ST_Y_Point(int8_t *p, int64_t psize, int32_t ic, int32_t isr, int32_t osr)
DEVICE ALWAYS_INLINE bool tol_le(double x, double y, double tolerance=TOLERANCE_DEFAULT)
EXTENSION_NOINLINE double ST_Perimeter_MultiPolygon_Geodesic(int8_t *mpoly_coords, int64_t mpoly_coords_size, int32_t *mpoly_ring_sizes, int64_t mpoly_num_rings, int32_t *mpoly_poly_sizes, int64_t mpoly_num_polys, int32_t ic, int32_t isr, int32_t osr)
EXTENSION_NOINLINE double convert_meters_to_pixel_height(const double meters, int8_t *p, const int64_t psize, const int32_t ic, const int32_t isr, const int32_t osr, const double min_lat, const double max_lat, const int32_t img_height, const double min_height)
#define TOLERANCE_GEOINT32
EXTENSION_NOINLINE double ST_Y_LineString(int8_t *l, int64_t lsize, int32_t lindex, int32_t ic, int32_t isr, int32_t osr)
DEVICE double max_distance_point_line(double px, double py, double l1x, double l1y, double l2x, double l2y)
EXTENSION_INLINE int32_t MapD_GeoPolyRenderGroup(int32_t render_group)
EXTENSION_NOINLINE bool Point_Overlaps_Box(double *bounds, int64_t bounds_size, double px, double py)
EXTENSION_INLINE int64_t OmniSci_Geo_PolyBoundsPtr(double *bounds, int64_t size)
DEVICE ALWAYS_INLINE int32_t compression_unit_size(int32_t ic)
#define DEVICE
EXTENSION_NOINLINE bool ST_Intersects_LineString_Linestring(int8_t *l1, int64_t l1size, double *l1bounds, int64_t l1bounds_size, int32_t l1i, int8_t *l2, int64_t l2size, double *l2bounds, int64_t l2bounds_size, int32_t l2i, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr)
#define EXTENSION_NOINLINE
Definition: OmniSciTypes.h:25
EXTENSION_INLINE double ST_Distance_LineString_Point(int8_t *l, int64_t lsize, int32_t lindex, int8_t *p, int64_t psize, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr)
EXTENSION_INLINE bool ST_Intersects_Polygon_LineString(int8_t *poly, int64_t polysize, int32_t *poly_ring_sizes, int64_t poly_num_rings, double *poly_bounds, int64_t poly_bounds_size, int8_t *l, int64_t lsize, double *lbounds, int64_t lbounds_size, int32_t li, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr)
DEVICE ALWAYS_INLINE int16_t orientation(double px, double py, double qx, double qy, double rx, double ry)
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)
EXTENSION_NOINLINE double ST_Centroid_MultiPolygon(int8_t *mpoly_coords, int64_t mpoly_coords_size, int32_t *mpoly_ring_sizes, int64_t mpoly_num_rings, int32_t *mpoly_poly_sizes, int64_t mpoly_num_polys, int32_t ic, int32_t isr, int32_t osr, bool ycoord)
EXTENSION_NOINLINE double ST_YMin(int8_t *coords, int64_t size, int32_t ic, int32_t isr, int32_t osr)
EXTENSION_INLINE bool ST_Intersects_LineString_Point(int8_t *l, int64_t lsize, double *lbounds, int64_t lbounds_size, int32_t li, int8_t *p, int64_t psize, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr)
DEVICE ALWAYS_INLINE double area_ring(int8_t *ring, int64_t ringsize, int32_t ic, int32_t isr, int32_t osr)
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)
EXTENSION_INLINE double ST_Area_Polygon_Geodesic(int8_t *poly_coords, int64_t poly_coords_size, int32_t *poly_ring_sizes, int64_t poly_num_rings, int32_t ic, int32_t isr, int32_t osr)
EXTENSION_NOINLINE bool ST_Intersects_LineString_Polygon(int8_t *l, int64_t lsize, double *lbounds, int64_t lbounds_size, int32_t li, int8_t *poly, int64_t polysize, int32_t *poly_ring_sizes, int64_t poly_num_rings, double *poly_bounds, int64_t poly_bounds_size, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr)
EXTENSION_INLINE double ST_Distance_MultiPolygon_LineString(int8_t *mpoly_coords, int64_t mpoly_coords_size, int32_t *mpoly_ring_sizes, int64_t mpoly_num_rings, int32_t *mpoly_poly_sizes, int64_t mpoly_num_polys, int8_t *l, int64_t lsize, int32_t lindex, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr)
DEVICE ALWAYS_INLINE bool box_overlaps_box(double *bounds1, int64_t bounds1_size, double *bounds2, int64_t bounds2_size)
EXTENSION_INLINE int64_t MapD_GeoPolyBoundsPtr(double *bounds, int64_t size)
EXTENSION_INLINE bool ST_Intersects_Point_MultiPolygon(int8_t *p, int64_t psize, int8_t *mpoly_coords, int64_t mpoly_coords_size, int32_t *mpoly_ring_sizes, int64_t mpoly_num_rings, int32_t *mpoly_poly_sizes, int64_t mpoly_num_polys, double *mpoly_bounds, int64_t mpoly_bounds_size, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr)
DEVICE bool linestring_intersects_line(int8_t *l, int32_t lnum_coords, double l1x, double l1y, double l2x, double l2y, int32_t ic1, int32_t isr1, int32_t osr)
EXTENSION_NOINLINE double ST_Distance_LineString_MultiPolygon(int8_t *l, int64_t lsize, int32_t lindex, int8_t *mpoly_coords, int64_t mpoly_coords_size, int32_t *mpoly_ring_sizes, int64_t mpoly_num_rings, int32_t *mpoly_poly_sizes, int64_t mpoly_num_polys, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr)
EXTENSION_NOINLINE double ST_Centroid_Polygon(int8_t *poly_coords, int64_t poly_coords_size, int32_t *poly_ring_sizes, int64_t poly_num_rings, int32_t ic, int32_t isr, int32_t osr, bool ycoord)
#define TOLERANCE_DEFAULT
EXTENSION_INLINE double ST_Distance_LineString_Point_Geodesic(int8_t *l, int64_t lsize, int32_t lindex, int8_t *p, int64_t psize, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr)
DEVICE bool line_intersects_line(double l11x, double l11y, double l12x, double l12y, double l21x, double l21y, double l22x, double l22y)
EXTENSION_NOINLINE double ST_Area_MultiPolygon(int8_t *mpoly_coords, int64_t mpoly_coords_size, int32_t *mpoly_ring_sizes, int64_t mpoly_num_rings, int32_t *mpoly_poly_sizes, int64_t mpoly_num_polys, int32_t ic, int32_t isr, int32_t osr)
EXTENSION_NOINLINE bool ST_Intersects_Polygon_MultiPolygon(int8_t *poly_coords, int64_t poly_coords_size, int32_t *poly_ring_sizes, int64_t poly_num_rings, double *poly_bounds, int64_t poly_bounds_size, int8_t *mpoly_coords, int64_t mpoly_coords_size, int32_t *mpoly_ring_sizes, int64_t mpoly_num_rings, int32_t *mpoly_poly_sizes, int64_t mpoly_num_polys, double *mpoly_bounds, int64_t mpoly_bounds_size, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr)
DEVICE ALWAYS_INLINE bool box_contains_point(double *bounds, int64_t bounds_size, double px, double py)
DEVICE bool polygon_contains_point(int8_t *poly, int32_t poly_num_coords, double px, double py, int32_t ic1, int32_t isr1, int32_t osr)
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_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)
EXTENSION_INLINE bool ST_Contains_LineString_Point(int8_t *l, int64_t lsize, double *lbounds, int64_t lbounds_size, int32_t li, int8_t *p, int64_t psize, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr)
EXTENSION_NOINLINE bool ST_Contains_LineString_LineString(int8_t *l1, int64_t l1size, double *l1bounds, int64_t l1bounds_size, int32_t l1i, int8_t *l2, int64_t l2size, double *l2bounds, int64_t l2bounds_size, int32_t l2i, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr)
DEVICE ALWAYS_INLINE bool tol_zero(double x, double tolerance=TOLERANCE_DEFAULT)
DEVICE ALWAYS_INLINE bool centroid_add_polygon(int8_t *poly_coords, int64_t poly_coords_size, int32_t *poly_ring_sizes, int64_t poly_num_rings, int32_t ic, int32_t isr, int32_t osr, double *total_area2, double *cg3, double *total_length, double *linestring_centroid_sum, int64_t *num_points, double *point_centroid_sum)
EXTENSION_NOINLINE double ST_Distance_Point_ClosedLineString(int8_t *p, int64_t psize, int8_t *l, int64_t lsize, int32_t lindex, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr)
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)
EXTENSION_NOINLINE bool ST_Contains_Point_Polygon(int8_t *p, int64_t psize, int8_t *poly_coords, int64_t poly_coords_size, int32_t *poly_ring_sizes, int64_t poly_num_rings, double *poly_bounds, int64_t poly_bounds_size, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr)
EXTENSION_NOINLINE double ST_Distance_Point_Point_Geodesic(int8_t *p1, int64_t p1size, int8_t *p2, int64_t p2size, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr)
EXTENSION_NOINLINE double ST_MaxDistance_Point_LineString(int8_t *p, int64_t psize, int8_t *l, int64_t lsize, int32_t lindex, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr)
EXTENSION_NOINLINE double ST_Length_LineString(int8_t *coords, int64_t coords_sz, int32_t ic, int32_t isr, int32_t osr)
DEVICE double decompress_longitude_coord_geoint32(const int32_t compressed)
EXTENSION_NOINLINE double ST_Perimeter_MultiPolygon(int8_t *mpoly_coords, int64_t mpoly_coords_size, int32_t *mpoly_ring_sizes, int64_t mpoly_num_rings, int32_t *mpoly_poly_sizes, int64_t mpoly_num_polys, int32_t ic, int32_t isr, int32_t osr)
EXTENSION_INLINE bool ST_Intersects_Point_Point(int8_t *p1, int64_t p1size, int8_t *p2, int64_t p2size, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr)
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)
DEVICE ALWAYS_INLINE bool centroid_add_segment(double x1, double y1, double x2, double y2, double *length, double *linestring_centroid_sum)
EXTENSION_INLINE double ST_Area_MultiPolygon_Geodesic(int8_t *mpoly_coords, int64_t mpoly_coords_size, int32_t *mpoly_ring_sizes, int64_t mpoly_num_rings, int32_t *mpoly_poly_sizes, int64_t mpoly_num_polys, int32_t ic, int32_t isr, int32_t osr)
EXTENSION_INLINE int32_t ST_NRings(int32_t *poly_ring_sizes, int64_t poly_num_rings)
DEVICE ALWAYS_INLINE bool centroid_add_linestring(int8_t *l, int64_t lsize, int32_t ic, int32_t isr, int32_t osr, bool closed, double *total_length, double *linestring_centroid_sum, int64_t *num_points, double *point_centroid_sum)
EXTENSION_NOINLINE double ST_Perimeter_Polygon_Geodesic(int8_t *poly, int64_t polysize, int32_t *poly_ring_sizes, int64_t poly_num_rings, int32_t ic, int32_t isr, int32_t osr)
DEVICE ALWAYS_INLINE double distance_point_point_squared(double p1x, double p1y, double p2x, double p2y)
DEVICE bool linestring_intersects_linestring(int8_t *l, int32_t lnum_coords, double l1x, double l1y, double l2x, double l2y, int32_t ic1, int32_t isr1, int32_t osr)
EXTENSION_INLINE bool ST_Intersects_MultiPolygon_MultiPolygon(int8_t *mpoly1_coords, int64_t mpoly1_coords_size, int32_t *mpoly1_ring_sizes, int64_t mpoly1_num_rings, int32_t *mpoly1_poly_sizes, int64_t mpoly1_num_polys, double *mpoly1_bounds, int64_t mpoly1_bounds_size, int8_t *mpoly2_coords, int64_t mpoly2_coords_size, int32_t *mpoly2_ring_sizes, int64_t mpoly2_num_rings, int32_t *mpoly2_poly_sizes, int64_t mpoly2_num_polys, double *mpoly2_bounds, int64_t mpoly2_bounds_size, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr)
EXTENSION_NOINLINE double distance_in_meters(const double fromlon, const double fromlat, const double tolon, const double tolat)
Computes the distance, in meters, between two WGS-84 positions.
EXTENSION_NOINLINE double ST_Length_LineString_Geodesic(int8_t *coords, int64_t coords_sz, int32_t ic, int32_t isr, int32_t osr)
EXTENSION_NOINLINE double ST_Distance_Point_LineString_Geodesic(int8_t *p, int64_t psize, int8_t *l, int64_t lsize, int32_t lindex, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr)
#define EXTENSION_INLINE
Definition: OmniSciTypes.h:24
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 tol_ge(double x, double y, double tolerance=TOLERANCE_DEFAULT)
#define COMPRESSION_GEOINT32
EXTENSION_NOINLINE double ST_Distance_LineString_Polygon(int8_t *l, int64_t lsize, int32_t lindex, int8_t *poly_coords, int64_t poly_coords_size, int32_t *poly_ring_sizes, int64_t poly_num_rings, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr)
EXTENSION_NOINLINE double ST_MaxDistance_LineString_Point(int8_t *l, int64_t lsize, int32_t lindex, int8_t *p, int64_t psize, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr)
DEVICE ALWAYS_INLINE double tol(int32_t ic)
DEVICE double distance_line_line(double l11x, double l11y, double l12x, double l12y, double l21x, double l21y, double l22x, double l22y)
DEVICE bool polygon_contains_linestring(int8_t *poly, int32_t poly_num_coords, int8_t *l, int64_t lnum_coords, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr)
EXTENSION_NOINLINE double ST_Centroid_Point(int8_t *p, int64_t psize, int32_t ic, int32_t isr, int32_t osr, bool ycoord)
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)
#define ALWAYS_INLINE
EXTENSION_INLINE bool ST_Intersects_MultiPolygon_LineString(int8_t *mpoly_coords, int64_t mpoly_coords_size, int32_t *mpoly_ring_sizes, int64_t mpoly_num_rings, int32_t *mpoly_poly_sizes, int64_t mpoly_num_polys, double *mpoly_bounds, int64_t mpoly_bounds_size, int8_t *l, int64_t lsize, double *lbounds, int64_t lbounds_size, int32_t li, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr)
EXTENSION_INLINE int32_t ST_NPoints(int8_t *coords, int64_t coords_sz, int32_t ic)
EXTENSION_NOINLINE double ST_Distance_LineString_LineString(int8_t *l1, int64_t l1size, int32_t l1index, int8_t *l2, int64_t l2size, int32_t l2index, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr)
EXTENSION_NOINLINE double conv_4326_900913_y(const double y)
DEVICE ALWAYS_INLINE double decompress_coord(int8_t *data, int32_t index, int32_t ic, bool x)
EXTENSION_NOINLINE double conv_4326_900913_x(const double x)
DEVICE ALWAYS_INLINE bool box_contains_box_vertex(double *bounds1, int64_t bounds1_size, double *bounds2, int64_t bounds2_size)
DEVICE ALWAYS_INLINE bool box_contains_box(double *bounds1, int64_t bounds1_size, double *bounds2, int64_t bounds2_size)
EXTENSION_INLINE double ST_YMax_Bounds(double *bounds, int64_t size, int32_t isr, int32_t osr)
EXTENSION_INLINE double ST_XMax_Bounds(double *bounds, int64_t size, int32_t isr, int32_t osr)
EXTENSION_NOINLINE bool ST_Contains_Point_LineString(int8_t *p, int64_t psize, int8_t *l, int64_t lsize, double *lbounds, int64_t lbounds_size, int32_t li, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr)
DEVICE ALWAYS_INLINE double distance_point_linestring(int8_t *p, int64_t psize, int8_t *l, int64_t lsize, int32_t lindex, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr, bool check_closed)
EXTENSION_NOINLINE bool ST_Intersects_Point_LineString(int8_t *p, int64_t psize, int8_t *l, int64_t lsize, double *lbounds, int64_t lbounds_size, int32_t li, int32_t ic1, int32_t isr1, int32_t ic2, int32_t isr2, int32_t osr)
DEVICE ALWAYS_INLINE double distance_point_point(double p1x, double p1y, double p2x, double p2y)
DEVICE ALWAYS_INLINE bool centroid_add_ring(int8_t *ring, int64_t ringsize, int32_t ic, int32_t isr, int32_t osr, double sign, double *total_area2, double *cg3, double *total_length, double *linestring_centroid_sum, int64_t *num_points, double *point_centroid_sum)