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