Line data Source code
1 : /*
2 : * SPDX-License-Identifier: MPL-2.0
3 : *
4 : * This Source Code Form is subject to the terms of the Mozilla Public
5 : * License, v. 2.0. If a copy of the MPL was not distributed with this
6 : * file, You can obtain one at http://mozilla.org/MPL/2.0/.
7 : *
8 : * Copyright 2024 MonetDB Foundation;
9 : * Copyright August 2008 - 2023 MonetDB B.V.;
10 : * Copyright 1997 - July 2008 CWI.
11 : */
12 :
13 : #include "geod.h"
14 :
15 : /**
16 : * Conversions
17 : *
18 : **/
19 :
20 : const double earth_radius = 6371.009;
21 : const double earth_radius_meters = 6371009;
22 : #ifndef M_PI
23 : #define M_PI ((double) 3.14159265358979323846) /* pi */
24 : #endif
25 : #ifndef M_PI_2
26 : #define M_PI_2 1.57079632679489661923132169163975144 /* pi/2 */
27 : #endif
28 :
29 : /* Converts a longitude value in degrees to radians */
30 : static double
31 0 : deg2RadLongitude(double lon_degrees)
32 : {
33 : //Convert
34 0 : double lon = M_PI * lon_degrees / 180.0;
35 : //Normalize
36 0 : if (lon == -1.0 * M_PI)
37 : return M_PI;
38 0 : if (lon == -2.0 * M_PI)
39 : return 0.0;
40 0 : if (lon > 2.0 * M_PI)
41 0 : lon = remainder(lon, 2.0 * M_PI);
42 :
43 0 : if (lon < -2.0 * M_PI)
44 0 : lon = remainder(lon, -2.0 * M_PI);
45 :
46 0 : if (lon > M_PI)
47 0 : lon = -2.0 * M_PI + lon;
48 :
49 0 : if (lon < -1.0 * M_PI)
50 0 : lon = 2.0 * M_PI + lon;
51 :
52 0 : if (lon == -2.0 * M_PI)
53 0 : lon *= -1.0;
54 :
55 : return lon;
56 : }
57 :
58 : /* Converts a latitude value in degrees to radians */
59 : static double
60 0 : deg2RadLatitude(double lat_degrees)
61 : {
62 : //Convert
63 0 : double lat = M_PI * lat_degrees / 180.0;
64 : //Normalize
65 0 : if (lat > 2.0 * M_PI)
66 0 : lat = remainder(lat, 2.0 * M_PI);
67 :
68 0 : if (lat < -2.0 * M_PI)
69 0 : lat = remainder(lat, -2.0 * M_PI);
70 :
71 0 : if (lat > M_PI)
72 0 : lat = M_PI - lat;
73 :
74 0 : if (lat < -1.0 * M_PI)
75 0 : lat = -1.0 * M_PI - lat;
76 :
77 0 : if (lat > M_PI_2)
78 0 : lat = M_PI - lat;
79 :
80 0 : if (lat < -1.0 * M_PI_2)
81 0 : lat = -1.0 * M_PI - lat;
82 :
83 0 : return lat;
84 : }
85 :
86 : /* Converts the GeoPoint from degrees to radians latitude and longitude*/
87 : static GeoPoint
88 0 : deg2RadPoint(GeoPoint geo)
89 : {
90 0 : geo.lon = deg2RadLongitude(geo.lon);
91 0 : geo.lat = deg2RadLatitude(geo.lat);
92 0 : return geo;
93 : }
94 :
95 : /**
96 : * Converts a longitude value in radians to degrees
97 : */
98 : static double
99 0 : rad2DegLongitude(double lon_radians)
100 : {
101 : //Convert
102 0 : double lon = lon_radians * 180.0 / M_PI;
103 : //Normalize
104 0 : if (lon > 360.0)
105 0 : lon = remainder(lon, 360.0);
106 :
107 0 : if (lon < -360.0)
108 0 : lon = remainder(lon, -360.0);
109 :
110 0 : if (lon > 180.0)
111 0 : lon = -360.0 + lon;
112 :
113 0 : if (lon < -180.0)
114 0 : lon = 360 + lon;
115 :
116 0 : if (lon == -180.0)
117 : return 180.0;
118 :
119 0 : if (lon == -360.0)
120 0 : return 0.0;
121 :
122 : return lon;
123 : }
124 :
125 : /**
126 : * Converts a latitude value in radians to degrees
127 : */
128 : static double
129 0 : rad2DegLatitude(double lat_radians)
130 : {
131 : //Convert
132 0 : double lat = lat_radians * 180.0 / M_PI;
133 : //Normalize
134 0 : if (lat > 360.0)
135 0 : lat = remainder(lat, 360.0);
136 :
137 0 : if (lat < -360.0)
138 0 : lat = remainder(lat, -360.0);
139 :
140 0 : if (lat > 180.0)
141 0 : lat = 180.0 - lat;
142 :
143 0 : if (lat < -180.0)
144 0 : lat = -180.0 - lat;
145 :
146 0 : if (lat > 90.0)
147 0 : lat = 180.0 - lat;
148 :
149 0 : if (lat < -90.0)
150 0 : lat = -180.0 - lat;
151 :
152 0 : return lat;
153 : }
154 :
155 : /* Converts the GeoPoint from degrees to radians latitude and longitude*/
156 : static GeoPoint
157 0 : rad2DegPoint(GeoPoint geo)
158 : {
159 0 : geo.lon = rad2DegLongitude(geo.lon);
160 0 : geo.lat = rad2DegLatitude(geo.lat);
161 0 : return geo;
162 : }
163 :
164 : /* Converts the a GEOSGeom Point into a GeoPoint */
165 : static GeoPoint
166 0 : geoPointFromGeom(GEOSGeom geom)
167 : {
168 0 : GeoPoint geo;
169 0 : GEOSGeomGetX_r(geoshandle, geom, &(geo.lon));
170 0 : GEOSGeomGetY_r(geoshandle, geom, &(geo.lat));
171 0 : return geo;
172 : }
173 :
174 : /* Converts the a GEOSGeom Line into GeoLines
175 : Argument must be a Line geometry. */
176 : static GeoLines
177 0 : geoLinesFromGeom(GEOSGeom geom)
178 : {
179 0 : const GEOSCoordSequence *gcs = GEOSGeom_getCoordSeq_r(geoshandle, geom);
180 0 : GeoLines geo;
181 0 : geo.pointCount = GEOSGeomGetNumPoints_r(geoshandle, geom);
182 0 : geo.points = GDKmalloc(sizeof(GeoPoint) * geo.pointCount);
183 0 : for (int i = 0; i < geo.pointCount; i++)
184 0 : GEOSCoordSeq_getXY_r(geoshandle, gcs, i, &geo.points[i].lon, &geo.points[i].lat);
185 0 : geo.bbox = NULL;
186 0 : return geo;
187 : }
188 :
189 : static BoundingBox * boundingBoxLines(GeoLines lines);
190 :
191 : /* Converts the a GEOSGeom Line into GeoPolygon (with exterior ring and zero-to-multiple interior rings)
192 : Argument must be a Polygon geometry. */
193 : static GeoPolygon
194 0 : geoPolygonFromGeom(GEOSGeom geom)
195 : {
196 0 : GeoPolygon geo;
197 : //Get exterior ring GeoLines
198 0 : geo.exteriorRing = geoLinesFromGeom((GEOSGeom)GEOSGetExteriorRing_r(geoshandle, geom));
199 0 : geo.interiorRingsCount = GEOSGetNumInteriorRings_r(geoshandle, geom);
200 : //If there are interior rings, allocate space to their GeoLines representation
201 0 : if (geo.interiorRingsCount > 0)
202 : //TODO Malloc fail exception?
203 0 : geo.interiorRings = GDKmalloc(sizeof(GeoLines) * geo.interiorRingsCount);
204 : else
205 : geo.interiorRings = NULL;
206 : //Get interior rings GeoLines
207 0 : for (int i = 0; i < geo.interiorRingsCount; i++)
208 0 : geo.interiorRings[i] = geoLinesFromGeom((GEOSGeom)GEOSGetInteriorRingN_r(geoshandle, geom, i));
209 : // If the geometry doesn't have BoundingBoxe, calculate it
210 0 : geo.bbox = boundingBoxLines(geo.exteriorRing);
211 0 : return geo;
212 : }
213 :
214 : static GeoPoint
215 0 : geoPointFromLatLon(double lon, double lat)
216 : {
217 0 : GeoPoint geo;
218 0 : geo.lon = lon;
219 0 : geo.lat = lat;
220 0 : return geo;
221 : }
222 :
223 : static str
224 0 : freeGeoLines(GeoLines lines) {
225 0 : str msg = MAL_SUCCEED;
226 0 : GDKfree(lines.points);
227 0 : if (lines.bbox)
228 0 : GDKfree(lines.bbox);
229 0 : return msg;
230 : }
231 :
232 : static str
233 0 : freeGeoPolygon(GeoPolygon polygon) {
234 0 : str msg = MAL_SUCCEED;
235 0 : msg = freeGeoLines(polygon.exteriorRing);
236 0 : if (polygon.bbox)
237 0 : GDKfree(polygon.bbox);
238 0 : for (int i = 0; i < polygon.interiorRingsCount; i++)
239 0 : msg = freeGeoLines(polygon.interiorRings[i]);
240 0 : if (polygon.interiorRings)
241 0 : GDKfree(polygon.interiorRings);
242 0 : return msg;
243 : }
244 :
245 : static CartPoint3D
246 0 : cartPointFromXYZ(double x, double y, double z)
247 : {
248 0 : CartPoint3D cart;
249 0 : cart.x = x;
250 0 : cart.y = y;
251 0 : cart.z = z;
252 0 : return cart;
253 : }
254 :
255 : /* Converts Well-Known Bytes into Geos Geometries, if they are not NULL and have the same SRID (used for geographic functions) */
256 : str
257 13 : wkbGetCompatibleGeometries(wkb * const *a, wkb * const *b, GEOSGeom *ga, GEOSGeom *gb)
258 : {
259 13 : str err = MAL_SUCCEED;
260 :
261 13 : if (is_wkb_nil(*a) || is_wkb_nil(*b)) {
262 0 : (*ga) = NULL;
263 0 : (*gb) = NULL;
264 0 : return MAL_SUCCEED;
265 : }
266 13 : (*ga) = wkb2geos(*a);
267 13 : (*gb) = wkb2geos(*b);
268 13 : if ((*ga) == NULL || (*gb) == NULL)
269 0 : err = createException(MAL, "geom.wkbGetComplatibleGeometries", SQLSTATE(38000) "Geos operation wkb2geos failed");
270 13 : else if (GEOSGetSRID_r(geoshandle, (*ga)) != GEOSGetSRID_r(geoshandle, *gb)) {
271 0 : GEOSGeom_destroy_r(geoshandle, *ga);
272 0 : GEOSGeom_destroy_r(geoshandle, *gb);
273 0 : err = createException(MAL, "geom.wkbGetComplatibleGeometries", SQLSTATE(38000) "Geometries of different SRID");
274 : }
275 : return err;
276 : }
277 :
278 : /**
279 : * Convert spherical coordinates to cartesian coordinates on unit sphere.
280 : * The inputs have to be in radians.
281 : */
282 : static CartPoint3D
283 0 : geo2cart(GeoPoint geo)
284 : {
285 0 : CartPoint3D cart;
286 0 : cart.x = cos(geo.lat) * cos(geo.lon);
287 0 : cart.y = cos(geo.lat) * sin(geo.lon);
288 0 : cart.z = sin(geo.lat);
289 0 : return cart;
290 : }
291 :
292 : /**
293 : * Convert spherical coordinates to cartesian coordinates on unit sphere.
294 : * The inputs have to be in degrees.
295 : */
296 : static CartPoint3D
297 0 : geo2cartFromDegrees(GeoPoint geo)
298 : {
299 0 : return geo2cart(deg2RadPoint(geo));
300 : }
301 :
302 : /* Convert cartesian coordinates to spherical coordinates on unit sphere */
303 : static GeoPoint
304 0 : cart2geo(CartPoint3D cart)
305 : {
306 0 : GeoPoint geo;
307 0 : geo.lon = atan2(cart.y, cart.x);
308 0 : geo.lat = asin(cart.z);
309 0 : return geo;
310 : }
311 :
312 : /* Converts two lat/lon points into cartesian coordinates and creates a Line geometry */
313 : static GEOSGeom
314 0 : cartesianLineFromGeoPoints(GeoPoint p1, GeoPoint p2)
315 : {
316 0 : CartPoint3D p1_cart, p2_cart;
317 0 : p1_cart = geo2cartFromDegrees(p1);
318 0 : p2_cart = geo2cartFromDegrees(p2);
319 0 : GEOSCoordSequence *lineSeq = GEOSCoordSeq_create_r(geoshandle, 2, 3);
320 0 : GEOSCoordSeq_setXYZ_r(geoshandle, lineSeq, 0, p1_cart.x, p1_cart.y, p1_cart.z);
321 0 : GEOSCoordSeq_setXYZ_r(geoshandle, lineSeq, 1, p2_cart.x, p2_cart.y, p2_cart.z);
322 0 : return GEOSGeom_createLineString_r(geoshandle, lineSeq);
323 : }
324 :
325 : /**
326 : * Bounding Box functions
327 : *
328 : **/
329 : /* Adds a Cartesian Point to the BoundingBox */
330 : static void
331 0 : boundingBoxAddPoint(BoundingBox *bb, CartPoint3D p)
332 : {
333 0 : if (bb->xmin > p.x)
334 0 : bb->xmin = p.x;
335 0 : if (bb->xmax < p.x)
336 0 : bb->xmax = p.x;
337 0 : if (bb->ymin > p.y)
338 0 : bb->ymin = p.y;
339 0 : if (bb->ymax < p.y)
340 0 : bb->ymax = p.y;
341 0 : if (bb->zmin > p.z)
342 0 : bb->zmin = p.z;
343 0 : if (bb->zmax < p.z)
344 0 : bb->zmax = p.z;
345 0 : }
346 :
347 : /* Builds the BoundingBox for a GeoLines geometry */
348 : static BoundingBox *
349 0 : boundingBoxLines(GeoLines lines)
350 : {
351 0 : CartPoint3D c;
352 0 : BoundingBox *bb;
353 :
354 : //If there are no segments, return NULL
355 0 : if (lines.pointCount == 0)
356 : return NULL;
357 :
358 0 : bb = GDKmalloc(sizeof(BoundingBox));
359 0 : if (bb == NULL)
360 : return NULL;
361 :
362 0 : c = geo2cartFromDegrees(lines.points[0]);
363 :
364 : //Initialize the bounding box with the first point
365 0 : *bb = (BoundingBox) {
366 : .xmin = c.x,
367 : .xmax = c.x,
368 : .ymin = c.y,
369 : .ymax = c.y,
370 : .zmin = c.z,
371 : .zmax = c.z,
372 : };
373 :
374 0 : for (int i = 1; i < lines.pointCount; i++) {
375 0 : c = geo2cartFromDegrees(lines.points[i]);
376 0 : boundingBoxAddPoint(bb, c);
377 : }
378 : return bb;
379 : }
380 :
381 : static int
382 0 : boundingBoxContainsPoint(BoundingBox bb, CartPoint3D pt)
383 : {
384 0 : return bb.xmin <= pt.x && bb.xmax >= pt.x && bb.ymin <= pt.y && bb.ymax >= pt.y && bb.zmin <= pt.z && bb.zmax >= pt.z;
385 : }
386 :
387 : /* Returns a point outside of the polygon's bounding box, for Point-In-Polygon calculation */
388 : static GeoPoint
389 0 : pointOutsidePolygon(GeoPolygon polygon)
390 : {
391 0 : BoundingBox bb = *polygon.bbox;
392 0 : BoundingBox bb2 = *polygon.bbox;
393 :
394 : //TODO: From POSTGIS -> CHANGE
395 0 : double grow = M_PI / 180.0 / 60.0;
396 0 : CartPoint3D corners[8];
397 0 : while (grow < M_PI) {
398 0 : if (bb.xmin > -1)
399 0 : bb.xmin -= grow;
400 0 : if (bb.ymin > -1)
401 0 : bb.ymin -= grow;
402 0 : if (bb.zmin > -1)
403 0 : bb.zmin -= grow;
404 0 : if (bb.xmax < 1)
405 0 : bb.xmax += grow;
406 0 : if (bb.ymax < 1)
407 0 : bb.ymax += grow;
408 0 : if (bb.zmax < 1)
409 0 : bb.zmax += grow;
410 :
411 0 : corners[0].x = bb.xmin;
412 0 : corners[0].y = bb.ymin;
413 0 : corners[0].z = bb.zmin;
414 :
415 0 : corners[1].x = bb.xmin;
416 0 : corners[1].y = bb.ymax;
417 0 : corners[1].z = bb.zmin;
418 :
419 0 : corners[2].x = bb.xmin;
420 0 : corners[2].y = bb.ymin;
421 0 : corners[2].z = bb.zmax;
422 :
423 0 : corners[3].x = bb.xmax;
424 0 : corners[3].y = bb.ymin;
425 0 : corners[3].z = bb.zmin;
426 :
427 0 : corners[4].x = bb.xmax;
428 0 : corners[4].y = bb.ymax;
429 0 : corners[4].z = bb.zmin;
430 :
431 0 : corners[5].x = bb.xmax;
432 0 : corners[5].y = bb.ymin;
433 0 : corners[5].z = bb.zmax;
434 :
435 0 : corners[6].x = bb.xmin;
436 0 : corners[6].y = bb.ymax;
437 0 : corners[6].z = bb.zmax;
438 :
439 0 : corners[7].x = bb.xmax;
440 0 : corners[7].y = bb.ymax;
441 0 : corners[7].z = bb.zmax;
442 :
443 0 : for (int i = 0; i < 8; i++)
444 0 : if (!boundingBoxContainsPoint(bb2, corners[i])) {
445 0 : CartPoint3D pt_cart = corners[i];
446 0 : return rad2DegPoint(cart2geo(pt_cart));
447 : }
448 0 : grow *= 2.0;
449 : }
450 : //TODO: Should this be the return value in case no point is found?
451 0 : return geoPointFromLatLon(0, 0);
452 : }
453 :
454 : /**
455 : * Distance functions
456 : *
457 : **/
458 : /**
459 : * The haversine formula calculate the distance in meters between two lat/lon points
460 : * The points must be measured in radians.
461 : * This formula assumes a spherical model of the earth, which can lead to an error of about 0.3% compared to a ellipsoidal model.
462 : */
463 : static double
464 0 : haversine(GeoPoint a, GeoPoint b)
465 : {
466 0 : double d_lon = b.lon - a.lon;
467 0 : double d_lat = b.lat - a.lat;
468 0 : double d = sin(d_lat / 2) * sin(d_lat / 2) + sin(d_lon / 2) * sin(d_lon / 2) * cos(b.lat) * cos(a.lat);
469 0 : double c = 2 * atan2(sqrt(d), sqrt(1 - d));
470 : //TODO: Same as the previous line (which one is best?) -> 2 * asin(sqrt(d));
471 0 : return earth_radius_meters * c;
472 : }
473 :
474 : /* Distance between two Points */
475 : static double
476 0 : geoDistancePointPoint(GeoPoint a, GeoPoint b)
477 : {
478 0 : return haversine(deg2RadPoint(a), deg2RadPoint(b));
479 : }
480 :
481 : /* Calculates the distance between the perpendicular projection of a point in the Line */
482 : static double
483 0 : calculatePerpendicularDistance(GeoPoint p_geo, GeoPoint l1_geo, GeoPoint l2_geo)
484 : {
485 0 : CartPoint3D l1, l2, p, projection;
486 0 : GeoPoint projection_geo;
487 :
488 : //First, convert the points to 3D cartesian coordinates
489 0 : l1 = geo2cartFromDegrees(l1_geo);
490 0 : l2 = geo2cartFromDegrees(l2_geo);
491 0 : p = geo2cartFromDegrees(p_geo);
492 :
493 : //Calculate the projection of point into the line
494 0 : double d_ab = (l2.z - l1.z) * (l2.z - l1.z) + (l2.y - l1.y) * (l2.y - l1.y) + (l2.x - l1.x) * (l2.x - l1.x);
495 0 : double t = (((p.x - l1.x) * (l2.x - l1.x)) + ((p.y - l1.y) * (l2.y - l1.y)) + ((p.z - l1.z) * (l2.z - l1.z))) / d_ab;
496 :
497 : //If t is not between 0 and 1, the projected point is not in the line, so there is no perpendicular, return huge number
498 0 : if (t < 0 || t > 1)
499 : return INT_MAX;
500 :
501 : //If the projection is in the line segment, build the point -> projection = l1 + t * (l2-l1)
502 0 : projection = cartPointFromXYZ(l1.x + t * (l2.x - l1.x), l1.y + t * (l2.y - l1.y), l1.z + t * (l2.z - l1.z));
503 :
504 : //Convert into geographic coordinates (radians)
505 0 : projection_geo = cart2geo(projection);
506 :
507 : //Calculate distance from original point to the projection
508 0 : return haversine(deg2RadPoint(p_geo), projection_geo);
509 : }
510 :
511 : /* Distance between Point and Line
512 : The returned distance is the minimum distance between the point and the line vertices
513 : and the perpendicular projection of the point in each line segment. */
514 : static double
515 0 : geoDistancePointLine(GeoPoint point, GeoLines lines, double distance_min_limit)
516 : {
517 0 : double distancePoint, distancePerpendicular, min_distance = INT_MAX;
518 0 : for (int i = 0; i < lines.pointCount-1; i++) {
519 0 : distancePoint = geoDistancePointPoint(point, lines.points[i]);
520 0 : distancePerpendicular = calculatePerpendicularDistance(point,lines.points[i],lines.points[i+1]);
521 0 : if (distancePoint < min_distance)
522 0 : min_distance = distancePoint;
523 0 : if (distancePerpendicular < min_distance)
524 0 : min_distance = distancePerpendicular;
525 : //Shortcut, if the geometries are already at their minimum distance
526 0 : if (min_distance <= distance_min_limit)
527 0 : return min_distance;
528 : }
529 0 : distancePoint = geoDistancePointPoint(point, lines.points[lines.pointCount-1]);
530 0 : return distancePoint < min_distance ? distancePoint : min_distance;
531 : }
532 :
533 : /* Distance between two Lines. */
534 : static double
535 0 : geoDistanceLineLine(GeoLines line1, GeoLines line2, double distance_min_limit)
536 : {
537 0 : double distance, min_distance = INT_MAX;
538 0 : for (int i = 0; i < line1.pointCount; i++) {
539 0 : distance = geoDistancePointLine(line1.points[i], line2, distance_min_limit);
540 0 : if (distance < min_distance)
541 0 : min_distance = distance;
542 : //Shortcut, if the geometries are already at their minimum distance
543 0 : if (min_distance <= distance_min_limit)
544 0 : return min_distance;
545 : }
546 0 : for (int i = 0; i < line2.pointCount; i++) {
547 0 : for (int j = 0; j < line1.pointCount - 1; j++) {
548 0 : distance = calculatePerpendicularDistance(line2.points[i],line1.points[j],line1.points[j+1]);
549 0 : if (distance < min_distance)
550 0 : min_distance = distance;
551 : //Shortcut, if the geometries are already at their minimum distance
552 0 : if (min_distance <= distance_min_limit)
553 0 : return min_distance;
554 : }
555 : }
556 : return min_distance;
557 : }
558 :
559 : /* Checks if a Point is within a Polygon */
560 : static bool
561 0 : pointWithinPolygon(GeoPolygon polygon, GeoPoint point)
562 : {
563 0 : int intersectionNum = 0;
564 0 : GEOSGeometry *segmentPolygon, *intersectionPoints;
565 0 : GeoLines polygonRing;
566 :
567 : //Get an point that's outside the polygon
568 0 : GeoPoint outsidePoint = pointOutsidePolygon(polygon);
569 :
570 : //No outside point was found, return false
571 0 : if (outsidePoint.lat == 0 && outsidePoint.lon == 0)
572 : return false;
573 :
574 : /*printf("Outside point: (%f %f)\n",outsidePoint.lon, outsidePoint.lat);
575 : fflush(stdout);*/
576 :
577 : //Construct a line between the outside point and the input point
578 0 : GEOSGeometry *outInLine = cartesianLineFromGeoPoints(point, outsidePoint);
579 :
580 : //TODO This is producing wrong results, review the intersection conditional
581 : //Count the number of intersections between the polygon exterior ring and the constructed line
582 0 : polygonRing = polygon.exteriorRing;
583 0 : for (int i = 0; i < polygonRing.pointCount-1; i++) {
584 0 : segmentPolygon = cartesianLineFromGeoPoints(polygonRing.points[i], polygonRing.points[i+1]);
585 0 : intersectionPoints = GEOSIntersection_r(geoshandle, segmentPolygon, outInLine);
586 :
587 : //If there is an intersection, a point will be returned (line when there is none)
588 0 : if (GEOSGeomTypeId_r(geoshandle, intersectionPoints) == GEOS_POINT)
589 0 : intersectionNum++;
590 :
591 0 : if (intersectionPoints != NULL)
592 0 : GEOSGeom_destroy_r(geoshandle, intersectionPoints);
593 0 : if (segmentPolygon != NULL)
594 0 : GEOSGeom_destroy_r(geoshandle, segmentPolygon);
595 : }
596 :
597 : //Count the number of intersections between the polygon interior rings and the constructed line
598 0 : for (int j = 0; j < polygon.interiorRingsCount; j++) {
599 0 : polygonRing = polygon.interiorRings[j];
600 0 : for (int i = 0; i < polygonRing.pointCount-1; i++) {
601 0 : segmentPolygon = cartesianLineFromGeoPoints(polygonRing.points[i], polygonRing.points[i+1]);
602 0 : intersectionPoints = GEOSIntersection_r(geoshandle, segmentPolygon, outInLine);
603 :
604 : //If there is an intersection, a point will be returned (line when there is none)
605 0 : if (GEOSGeomTypeId_r(geoshandle, intersectionPoints) == GEOS_POINT)
606 0 : intersectionNum++;
607 :
608 0 : if (intersectionPoints != NULL)
609 0 : GEOSGeom_destroy_r(geoshandle, intersectionPoints);
610 0 : if (segmentPolygon != NULL)
611 0 : GEOSGeom_destroy_r(geoshandle, segmentPolygon);
612 : }
613 : }
614 :
615 0 : if (outInLine != NULL)
616 0 : GEOSGeom_destroy_r(geoshandle, outInLine);
617 :
618 : //If even, the point is not within the polygon. If odd, it is within
619 0 : return intersectionNum % 2 == 1;
620 : }
621 :
622 : /* Distance between Point and Polygon.*/
623 : static double
624 0 : geoDistancePointPolygon(GeoPoint point, GeoPolygon polygon, double distance_min_limit)
625 : {
626 : //Check if point is in polygon
627 0 : if (pointWithinPolygon(polygon, point))
628 : return 0;
629 :
630 : //Calculate distance from Point to the exterior and interior rings of the polygon
631 0 : double distance, min_distance = INT_MAX;
632 : //First, calculate distance to the exterior ring
633 0 : min_distance = geoDistancePointLine(point, polygon.exteriorRing, distance_min_limit);
634 : //Then, calculate distance to the interior rings
635 0 : for (int i = 0; i < polygon.interiorRingsCount; i++) {
636 : //Shortcut, if the geometries are already at their minimum distance
637 0 : if (min_distance <= distance_min_limit)
638 0 : return min_distance;
639 0 : distance = geoDistancePointLine(point, polygon.interiorRings[i], distance_min_limit);
640 0 : if (distance < min_distance)
641 0 : min_distance = distance;
642 : }
643 : return min_distance;
644 : }
645 :
646 : /* Distance between Line and Polygon. */
647 : static double
648 0 : geoDistanceLinePolygon(GeoLines line, GeoPolygon polygon, double distance_min_limit)
649 : {
650 0 : double distance, min_distance = INT_MAX;
651 : //Calculate distance to all start vertices of the line
652 0 : for (int i = 0; i < line.pointCount; i++) {
653 0 : distance = geoDistancePointPolygon(line.points[i], polygon, distance_min_limit);
654 :
655 : //Short-cut in case the point is within the polygon
656 0 : if (distance <= distance_min_limit)
657 0 : return distance;
658 :
659 0 : if (distance < min_distance)
660 0 : min_distance = distance;
661 : }
662 : return min_distance;
663 : }
664 :
665 : /* Distance between two Polygons. */
666 : static double
667 0 : geoDistancePolygonPolygon(GeoPolygon polygon1, GeoPolygon polygon2, double distance_min_limit)
668 : {
669 0 : double distance1, distance2;
670 : //Calculate the distance between the exterior ring of polygon1 and all segments of polygon2 (including the interior rings)
671 0 : distance1 = geoDistanceLinePolygon(polygon1.exteriorRing, polygon2, distance_min_limit);
672 : //Shortcut, if the geometries are already at their minimum distance
673 0 : if (distance1 <= distance_min_limit)
674 : return distance1;
675 0 : distance2 = geoDistanceLinePolygon(polygon2.exteriorRing, polygon1, distance_min_limit);
676 0 : return distance1 < distance2 ? distance1 : distance2;
677 : }
678 :
679 : /* Distance between two (non-collection) geometries. */
680 : static double
681 0 : geoDistanceSingle(GEOSGeom aGeom, GEOSGeom bGeom, double distance_min_limit)
682 : {
683 0 : int dimA, dimB;
684 0 : double distance = INT_MAX;
685 0 : dimA = GEOSGeom_getDimensions_r(geoshandle, aGeom);
686 0 : dimB = GEOSGeom_getDimensions_r(geoshandle, bGeom);
687 0 : if (dimA == 0 && dimB == 0) {
688 : /* Point and Point */
689 0 : GeoPoint a = geoPointFromGeom(aGeom);
690 0 : GeoPoint b = geoPointFromGeom(bGeom);
691 0 : distance = geoDistancePointPoint(a, b);
692 0 : } else if (dimA == 0 && dimB == 1) {
693 : /* Point and Line/LinearRing */
694 0 : GeoPoint a = geoPointFromGeom(aGeom);
695 0 : GeoLines b = geoLinesFromGeom(bGeom);
696 0 : distance = geoDistancePointLine(a, b, distance_min_limit);
697 0 : freeGeoLines(b);
698 0 : } else if (dimA == 1 && dimB == 0) {
699 : /* Line/LinearRing and Point */
700 0 : GeoLines a = geoLinesFromGeom(aGeom);
701 0 : GeoPoint b = geoPointFromGeom(bGeom);
702 0 : distance = geoDistancePointLine(b, a, distance_min_limit);
703 0 : freeGeoLines(a);
704 0 : } else if (dimA == 1 && dimB == 1) {
705 : /* Line/LinearRing and Line/LinearRing */
706 0 : GeoLines a = geoLinesFromGeom(aGeom);
707 0 : GeoLines b = geoLinesFromGeom(bGeom);
708 0 : distance = geoDistanceLineLine(a, b, distance_min_limit);
709 0 : freeGeoLines(a);
710 0 : freeGeoLines(b);
711 0 : } else if (dimA == 0 && dimB == 2) {
712 : /* Point and Polygon */
713 0 : GeoPoint a = geoPointFromGeom(aGeom);
714 0 : GeoPolygon b = geoPolygonFromGeom(bGeom);
715 0 : distance = geoDistancePointPolygon(a, b, distance_min_limit);
716 0 : freeGeoPolygon(b);
717 0 : } else if (dimA == 2 && dimB == 0) {
718 : /* Polygon and Point */
719 0 : GeoPolygon a = geoPolygonFromGeom(aGeom);
720 0 : GeoPoint b = geoPointFromGeom(bGeom);
721 0 : distance = geoDistancePointPolygon(b, a, distance_min_limit);
722 0 : freeGeoPolygon(a);
723 0 : } else if (dimA == 1 && dimB == 2) {
724 : /* Line/LinearRing and Polygon */
725 0 : GeoLines a = geoLinesFromGeom(aGeom);
726 0 : GeoPolygon b = geoPolygonFromGeom(bGeom);
727 0 : distance = geoDistanceLinePolygon(a, b, distance_min_limit);
728 0 : freeGeoLines(a);
729 0 : freeGeoPolygon(b);
730 0 : } else if (dimA == 2 && dimB == 1) {
731 : /* Polygon and Line/LinearRing */
732 0 : GeoPolygon a = geoPolygonFromGeom(aGeom);
733 0 : GeoLines b = geoLinesFromGeom(bGeom);
734 0 : distance = geoDistanceLinePolygon(b, a, distance_min_limit);
735 0 : freeGeoPolygon(a);
736 0 : freeGeoLines(b);
737 0 : } else if (dimA == 2 && dimB == 2) {
738 : /* Polygon and Polygon */
739 0 : GeoPolygon a = geoPolygonFromGeom(aGeom);
740 0 : GeoPolygon b = geoPolygonFromGeom(bGeom);
741 0 : distance = geoDistancePolygonPolygon(a, b, distance_min_limit);
742 0 : freeGeoPolygon(a);
743 0 : freeGeoPolygon(b);
744 : }
745 0 : return distance;
746 : }
747 :
748 : //The distance_min_limit argument is used for DWithin and Intersects.
749 : //If we get to the minimum distance for the predicate, return immediately
750 : //It is equal to 0 if the operation is Distance
751 : static double
752 0 : geoDistanceInternal(GEOSGeom a, GEOSGeom b, double distance_min_limit)
753 : {
754 0 : int numGeomsA = GEOSGetNumGeometries_r(geoshandle, a), numGeomsB = GEOSGetNumGeometries_r(geoshandle, b);
755 0 : double distance, min_distance = INT_MAX;
756 0 : GEOSGeometry *geo1, *geo2;
757 0 : for (int i = 0; i < numGeomsA; i++) {
758 0 : geo1 = (GEOSGeometry *)GEOSGetGeometryN_r(geoshandle, (const GEOSGeometry *)a, i);
759 0 : for (int j = 0; j < numGeomsB; j++) {
760 0 : geo2 = (GEOSGeometry *)GEOSGetGeometryN_r(geoshandle, (const GEOSGeometry *)b, j);
761 0 : distance = geoDistanceSingle(geo1, geo2, distance_min_limit);
762 : //Shortcut, if the geometries are already at their minimum distance (0 in the case of normal Distance)
763 0 : if (distance <= distance_min_limit)
764 0 : return distance;
765 0 : if (distance < min_distance)
766 0 : min_distance = distance;
767 : }
768 : }
769 : return min_distance;
770 : }
771 :
772 : /**
773 : * Distance
774 : *
775 : **/
776 : /* Calculates the distance, in meters, between two geographic geometries with latitude/longitude coordinates */
777 : str
778 0 : wkbDistanceGeographic(dbl *out, wkb * const *a, wkb * const *b)
779 : {
780 0 : str err = MAL_SUCCEED;
781 0 : GEOSGeom ga, gb;
782 0 : err = wkbGetCompatibleGeometries(a, b, &ga, &gb);
783 0 : if (ga && gb) {
784 0 : (*out) = geoDistanceInternal(ga, gb, 0);
785 : }
786 0 : GEOSGeom_destroy_r(geoshandle, ga);
787 0 : GEOSGeom_destroy_r(geoshandle, gb);
788 0 : return err;
789 : }
790 :
791 : /**
792 : * Distance Within
793 : *
794 : **/
795 : /* Checks if two geographic geometries are within d meters of one another */
796 : str
797 0 : wkbDWithinGeographic(bit *out, wkb * const *a, wkb * const *b, const dbl *d)
798 : {
799 0 : str err = MAL_SUCCEED;
800 0 : GEOSGeom ga, gb;
801 0 : double distance;
802 0 : err = wkbGetCompatibleGeometries(a, b, &ga, &gb);
803 0 : if (ga && gb) {
804 0 : distance = geoDistanceInternal(ga, gb, *d);
805 0 : (*out) = (distance <= (*d));
806 : }
807 0 : GEOSGeom_destroy_r(geoshandle, ga);
808 0 : GEOSGeom_destroy_r(geoshandle, gb);
809 0 : return err;
810 : }
811 :
812 : /**
813 : * Intersects
814 : *
815 : **/
816 : /* Checks if two geographic geometries intersect at any point */
817 : str
818 0 : wkbIntersectsGeographic(bit *out, wkb * const *a, wkb * const *b)
819 : {
820 0 : str err = MAL_SUCCEED;
821 0 : GEOSGeom ga, gb;
822 0 : double distance;
823 0 : err = wkbGetCompatibleGeometries(a, b, &ga, &gb);
824 0 : if (ga && gb) {
825 0 : distance = geoDistanceInternal(ga, gb, 0);
826 0 : (*out) = (distance == 0);
827 : }
828 0 : GEOSGeom_destroy_r(geoshandle, ga);
829 0 : GEOSGeom_destroy_r(geoshandle, gb);
830 0 : return err;
831 : }
832 :
833 : /* Checks if a Polygon covers a Line geometry */
834 : static bool
835 0 : geoPolygonCoversLine(GeoPolygon polygon, GeoLines lines)
836 : {
837 0 : for (int i = 0; i < lines.pointCount; i++) {
838 0 : if (pointWithinPolygon(polygon, lines.points[i]) == false)
839 : return false;
840 : }
841 : return true;
842 : }
843 :
844 : /* Compares two GeoPoints, returns true if they're equal */
845 : static bool
846 0 : geoPointEquals(GeoPoint pointA, GeoPoint pointB)
847 : {
848 0 : return (pointA.lat == pointB.lat) && (pointA.lon = pointB.lon);
849 : }
850 :
851 : //TODO Check if this works correctly
852 : static bool
853 0 : geoCoversSingle(GEOSGeom a, GEOSGeom b)
854 : {
855 0 : int dimA = GEOSGeom_getDimensions_r(geoshandle, a), dimB = GEOSGeom_getDimensions_r(geoshandle, b);
856 0 : if (dimA < dimB)
857 : //If the dimension of A is smaller than B, then it must not cover it
858 : return false;
859 :
860 0 : if (dimA == 0){
861 : //A and B are Points
862 0 : GeoPoint pointA = geoPointFromGeom(a);
863 0 : GeoPoint pointB = geoPointFromGeom(b);
864 0 : return geoPointEquals(pointA, pointB);
865 0 : } else if (dimA == 1) {
866 : //A is Line
867 : //GeoLines lineA = geoLinesFromGeom(a);
868 : if (dimB == 0) {
869 : //B is Point
870 : //GeoPoint pointB = geoPointFromGeom(b);
871 : //return geoLineCoversPoint(lineA,pointB);
872 : return false;
873 : } else {
874 : //B is Line
875 : //GeoLines lineB = geoLinesFromGeom(b);
876 : //return geoLineCoversLine(lineA, lineB);
877 : return false;
878 : }
879 0 : } else if (dimA == 2) {
880 : //A is Polygon
881 0 : GeoPolygon polygonA = geoPolygonFromGeom(a);
882 0 : if (dimB == 0){
883 : //B is Point
884 0 : GeoPoint pointB = geoPointFromGeom(b);
885 0 : return pointWithinPolygon(polygonA, pointB);
886 0 : } else if (dimB == 1) {
887 : //B is Line
888 0 : GeoLines lineB = geoLinesFromGeom(b);
889 0 : return geoPolygonCoversLine(polygonA, lineB);
890 : } else {
891 : //B is Polygon
892 0 : GeoPolygon polygonB = geoPolygonFromGeom(b);
893 : //If every point in the exterior ring of B is covered, polygon B is covered by polygon A
894 0 : return geoPolygonCoversLine(polygonA, polygonB.exteriorRing);
895 : }
896 : } else
897 : return false;
898 : }
899 :
900 : static bool
901 0 : geoCoversInternal(GEOSGeom a, GEOSGeom b)
902 : {
903 0 : int numGeomsA = GEOSGetNumGeometries_r(geoshandle, a), numGeomsB = GEOSGetNumGeometries_r(geoshandle, b);
904 0 : GEOSGeometry *geo1, *geo2;
905 0 : for (int i = 0; i < numGeomsA; i++) {
906 0 : geo1 = (GEOSGeometry *)GEOSGetGeometryN_r(geoshandle, (const GEOSGeometry *)a, i);
907 0 : for (int j = 0; j < numGeomsB; j++) {
908 0 : geo2 = (GEOSGeometry *)GEOSGetGeometryN_r(geoshandle, (const GEOSGeometry *)b, j);
909 0 : if (geoCoversSingle(geo1, geo2) == 0)
910 : return 0;
911 : }
912 : }
913 : return 1;
914 : }
915 :
916 : /**
917 : * Covers
918 : *
919 : **/
920 : /* Checks if no point of Geometry B is outside Geometry A */
921 : str
922 0 : wkbCoversGeographic(bit *out, wkb * const *a, wkb * const *b)
923 : {
924 0 : str err = MAL_SUCCEED;
925 0 : GEOSGeom ga, gb;
926 0 : err = wkbGetCompatibleGeometries(a, b, &ga, &gb);
927 0 : if (ga && gb)
928 0 : (*out) = geoCoversInternal(ga, gb);
929 :
930 0 : GEOSGeom_destroy_r(geoshandle, ga);
931 0 : GEOSGeom_destroy_r(geoshandle, gb);
932 :
933 0 : return err;
934 : }
935 :
936 : /**
937 : * FILTER FUNCTIONS
938 : **/
939 :
940 : static inline bit
941 0 : geosDistanceWithin (GEOSGeom geom1, GEOSGeom geom2, dbl distance_within) {
942 0 : dbl distance = geoDistanceInternal(geom1, geom2, distance_within);
943 0 : return distance <= distance_within;
944 : }
945 :
946 : //TODO Change BUNappend with manual insertion into the result BAT
947 : static str
948 0 : filterSelectGeomGeomDoubleToBit(bat* outid, const bat *bid , const bat *sid, const wkb *wkb_const, dbl double_flag, bit anti, bit (*func) (GEOSGeom, GEOSGeom, dbl), const char *name)
949 : {
950 0 : BAT *out = NULL, *b = NULL, *s = NULL;
951 0 : BATiter b_iter;
952 0 : struct canditer ci;
953 0 : GEOSGeom col_geom, const_geom;
954 :
955 : //Check if the geometry is null and convert to GEOS
956 0 : if ((const_geom = wkb2geos(wkb_const)) == NULL) {
957 0 : if ((out = BATdense(0, 0, 0)) == NULL)
958 0 : throw(MAL, name, GDK_EXCEPTION);
959 0 : *outid = out->batCacheid;
960 0 : BBPkeepref(out);
961 0 : return MAL_SUCCEED;
962 : }
963 : //get the BATs
964 0 : if ((b = BATdescriptor(*bid)) == NULL)
965 0 : throw(MAL, name, SQLSTATE(HY002) RUNTIME_OBJECT_MISSING);
966 : //get the candidate lists
967 0 : if (sid && !is_bat_nil(*sid) && !(s = BATdescriptor(*sid))) {
968 0 : BBPunfix(b->batCacheid);
969 0 : throw(MAL, name, SQLSTATE(HY002) RUNTIME_OBJECT_MISSING);
970 : }
971 0 : canditer_init(&ci, b, s);
972 : //create a new BAT for the output
973 0 : if ((out = COLnew(0, ATOMindex("oid"), ci.ncand, TRANSIENT)) == NULL) {
974 0 : BBPunfix(b->batCacheid);
975 0 : if (s)
976 0 : BBPunfix(s->batCacheid);
977 0 : throw(MAL, name, SQLSTATE(HY013) MAL_MALLOC_FAIL);
978 : }
979 0 : b_iter = bat_iterator(b);
980 : //Loop through column and compare with constant
981 0 : for (BUN i = 0; i < ci.ncand; i++) {
982 0 : oid c_oid = canditer_next(&ci);
983 0 : const wkb *col_wkb = BUNtvar(b_iter, c_oid - b->hseqbase);
984 0 : if ((col_geom = wkb2geos(col_wkb)) == NULL)
985 0 : continue;
986 0 : if (GEOSGetSRID_r(geoshandle, col_geom) != GEOSGetSRID_r(geoshandle, const_geom)) {
987 0 : GEOSGeom_destroy_r(geoshandle, col_geom);
988 0 : GEOSGeom_destroy_r(geoshandle, const_geom);
989 0 : bat_iterator_end(&b_iter);
990 0 : BBPunfix(b->batCacheid);
991 0 : if (s)
992 0 : BBPunfix(s->batCacheid);
993 0 : BBPreclaim(out);
994 0 : throw(MAL, name, SQLSTATE(38000) "Geometries of different SRID");
995 : }
996 : //Apply the (Geom, Geom, double) -> bit function
997 0 : bit cond = (*func)(col_geom, const_geom, double_flag);
998 0 : if (cond != anti) {
999 0 : if (BUNappend(out, &c_oid, false) != GDK_SUCCEED) {
1000 0 : if (col_geom)
1001 0 : GEOSGeom_destroy_r(geoshandle, col_geom);
1002 0 : if (const_geom)
1003 0 : GEOSGeom_destroy_r(geoshandle, const_geom);
1004 0 : bat_iterator_end(&b_iter);
1005 0 : BBPunfix(b->batCacheid);
1006 0 : if (s)
1007 0 : BBPunfix(s->batCacheid);
1008 0 : BBPreclaim(out);
1009 0 : throw(MAL, name, SQLSTATE(HY013) MAL_MALLOC_FAIL);
1010 : }
1011 : }
1012 0 : GEOSGeom_destroy_r(geoshandle, col_geom);
1013 : }
1014 0 : GEOSGeom_destroy_r(geoshandle, const_geom);
1015 0 : bat_iterator_end(&b_iter);
1016 0 : BBPunfix(b->batCacheid);
1017 0 : if (s)
1018 0 : BBPunfix(s->batCacheid);
1019 0 : *outid = out->batCacheid;
1020 0 : BBPkeepref(out);
1021 0 : return MAL_SUCCEED;
1022 : }
1023 :
1024 : static str
1025 0 : filterJoinGeomGeomDoubleToBit(bat *lres_id, bat *rres_id, const bat *l_id, const bat *r_id, double double_flag, const bat *ls_id, const bat *rs_id, bit nil_matches, const lng *estimate, bit anti, bit (*func) (GEOSGeom, GEOSGeom, dbl), const char *name)
1026 : {
1027 0 : BAT *lres = NULL, *rres = NULL, *l = NULL, *r = NULL, *ls = NULL, *rs = NULL;
1028 0 : BATiter l_iter, r_iter;
1029 0 : str msg = MAL_SUCCEED;
1030 0 : struct canditer l_ci, r_ci;
1031 0 : GEOSGeom l_geom, r_geom;
1032 0 : GEOSGeom *l_geoms = NULL, *r_geoms = NULL;
1033 0 : BUN est;
1034 :
1035 : //get the input BATs
1036 0 : l = BATdescriptor(*l_id);
1037 0 : r = BATdescriptor(*r_id);
1038 0 : if (l == NULL || r == NULL) {
1039 0 : if (l)
1040 0 : BBPunfix(l->batCacheid);
1041 0 : if (r)
1042 0 : BBPunfix(r->batCacheid);
1043 0 : throw(MAL, name, SQLSTATE(HY002) RUNTIME_OBJECT_MISSING);
1044 : }
1045 : //get the candidate lists
1046 0 : if (ls_id && !is_bat_nil(*ls_id) && !(ls = BATdescriptor(*ls_id)) && rs_id && !is_bat_nil(*rs_id) && !(rs = BATdescriptor(*rs_id))) {
1047 0 : msg = createException(MAL, name, SQLSTATE(HY002) RUNTIME_OBJECT_MISSING);
1048 0 : goto free;
1049 : }
1050 0 : canditer_init(&l_ci, l, ls);
1051 0 : canditer_init(&r_ci, r, rs);
1052 : //create new BATs for the output
1053 0 : est = is_lng_nil(*estimate) || *estimate == 0 ? l_ci.ncand : (BUN) *estimate;
1054 0 : if ((lres = COLnew(0, ATOMindex("oid"), est, TRANSIENT)) == NULL) {
1055 0 : msg = createException(MAL, name, SQLSTATE(HY013) MAL_MALLOC_FAIL);
1056 0 : goto free;
1057 : }
1058 0 : if ((rres = COLnew(0, ATOMindex("oid"), est, TRANSIENT)) == NULL) {
1059 0 : msg = createException(MAL, name, SQLSTATE(HY013) MAL_MALLOC_FAIL);
1060 0 : goto free;
1061 : }
1062 :
1063 : //Allocate arrays for reutilizing GEOS type conversion
1064 0 : if ((l_geoms = GDKmalloc(l_ci.ncand * sizeof(GEOSGeometry *))) == NULL || (r_geoms = GDKmalloc(r_ci.ncand * sizeof(GEOSGeometry *))) == NULL) {
1065 0 : msg = createException(MAL, name, SQLSTATE(HY013) MAL_MALLOC_FAIL);
1066 0 : goto free;
1067 : }
1068 :
1069 0 : l_iter = bat_iterator(l);
1070 0 : r_iter = bat_iterator(r);
1071 :
1072 : //Convert wkb to GEOS only once
1073 0 : for (BUN i = 0; i < l_ci.ncand; i++) {
1074 0 : oid l_oid = canditer_next(&l_ci);
1075 0 : l_geoms[i] = wkb2geos((const wkb*) BUNtvar(l_iter, l_oid - l->hseqbase));
1076 : }
1077 0 : for (BUN j = 0; j < r_ci.ncand; j++) {
1078 0 : oid r_oid = canditer_next(&r_ci);
1079 0 : r_geoms[j] = wkb2geos((const wkb*)BUNtvar(r_iter, r_oid - r->hseqbase));
1080 : }
1081 :
1082 0 : canditer_reset(&l_ci);
1083 0 : for (BUN i = 0; i < l_ci.ncand; i++) {
1084 0 : oid l_oid = canditer_next(&l_ci);
1085 0 : l_geom = l_geoms[i];
1086 0 : if (!nil_matches && l_geom == NULL)
1087 0 : continue;
1088 0 : canditer_reset(&r_ci);
1089 0 : for (BUN j = 0; j < r_ci.ncand; j++) {
1090 0 : oid r_oid = canditer_next(&r_ci);
1091 0 : r_geom = r_geoms[j];
1092 : //Null handling
1093 0 : if (r_geom == NULL) {
1094 0 : if (nil_matches && l_geom == NULL) {
1095 0 : if (BUNappend(lres, &l_oid, false) != GDK_SUCCEED || BUNappend(rres, &r_oid, false) != GDK_SUCCEED) {
1096 0 : msg = createException(MAL, name, SQLSTATE(HY013) MAL_MALLOC_FAIL);
1097 0 : bat_iterator_end(&l_iter);
1098 0 : bat_iterator_end(&r_iter);
1099 0 : goto free;
1100 : }
1101 : }
1102 : else
1103 0 : continue;
1104 : }
1105 : //TODO Do we need to do this check for every element?
1106 0 : if (GEOSGetSRID_r(geoshandle, l_geom) != GEOSGetSRID_r(geoshandle, r_geom)) {
1107 0 : msg = createException(MAL, name, SQLSTATE(38000) "Geometries of different SRID");
1108 0 : bat_iterator_end(&l_iter);
1109 0 : bat_iterator_end(&r_iter);
1110 0 : goto free;
1111 : }
1112 : //Apply the (Geom, Geom) -> bit function
1113 0 : bit cond = (*func)(l_geom, r_geom, double_flag);
1114 0 : if (cond != anti) {
1115 0 : if (BUNappend(lres, &l_oid, false) != GDK_SUCCEED || BUNappend(rres, &r_oid, false) != GDK_SUCCEED) {
1116 0 : msg = createException(MAL, name, SQLSTATE(HY013) MAL_MALLOC_FAIL);
1117 0 : bat_iterator_end(&l_iter);
1118 0 : bat_iterator_end(&r_iter);
1119 0 : goto free;
1120 : }
1121 : }
1122 : }
1123 : }
1124 : if (l_geoms) {
1125 0 : for (BUN i = 0; i < l_ci.ncand; i++) {
1126 0 : GEOSGeom_destroy_r(geoshandle, l_geoms[i]);
1127 : }
1128 0 : GDKfree(l_geoms);
1129 : }
1130 0 : if (r_geoms) {
1131 0 : for (BUN i = 0; i < r_ci.ncand; i++) {
1132 0 : GEOSGeom_destroy_r(geoshandle, r_geoms[i]);
1133 : }
1134 0 : GDKfree(r_geoms);
1135 : }
1136 0 : bat_iterator_end(&l_iter);
1137 0 : bat_iterator_end(&r_iter);
1138 0 : BBPunfix(l->batCacheid);
1139 0 : BBPunfix(r->batCacheid);
1140 0 : if (ls)
1141 0 : BBPunfix(ls->batCacheid);
1142 0 : if (rs)
1143 0 : BBPunfix(rs->batCacheid);
1144 0 : *lres_id = lres->batCacheid;
1145 0 : BBPkeepref(lres);
1146 0 : *rres_id = rres->batCacheid;
1147 0 : BBPkeepref(rres);
1148 0 : return MAL_SUCCEED;
1149 0 : free:
1150 0 : if (l_geoms) {
1151 0 : for (BUN i = 0; i < l_ci.ncand; i++) {
1152 0 : GEOSGeom_destroy_r(geoshandle, l_geoms[i]);
1153 : }
1154 0 : GDKfree(l_geoms);
1155 : }
1156 0 : if (r_geoms) {
1157 0 : for (BUN i = 0; i < r_ci.ncand; i++) {
1158 0 : GEOSGeom_destroy_r(geoshandle, r_geoms[i]);
1159 : }
1160 0 : GDKfree(r_geoms);
1161 : }
1162 0 : BBPunfix(l->batCacheid);
1163 0 : BBPunfix(r->batCacheid);
1164 0 : if (ls)
1165 0 : BBPunfix(ls->batCacheid);
1166 0 : if (rs)
1167 0 : BBPunfix(rs->batCacheid);
1168 0 : if (lres)
1169 0 : BBPreclaim(lres);
1170 0 : if (rres)
1171 0 : BBPreclaim(rres);
1172 : return msg;
1173 : }
1174 :
1175 : str
1176 0 : wkbDWithinGeographicJoin(bat *lres_id, bat *rres_id, const bat *l_id, const bat *r_id, const bat *d_id, const bat *ls_id, const bat *rs_id, const bit *nil_matches, const lng *estimate, const bit *anti) {
1177 0 : double distance_within = 0;
1178 0 : BAT *d = NULL;
1179 : //Get the distance BAT and get the double value
1180 0 : if ((d = BATdescriptor(*d_id)) == NULL) {
1181 0 : throw(MAL, "geom.wkbDWithinGeographicJoin", SQLSTATE(HY002) RUNTIME_OBJECT_MISSING);
1182 : }
1183 0 : if (BATcount(d) == 1) {
1184 0 : distance_within = *((double*) Tloc(d, 0));
1185 : }
1186 0 : BBPunfix(d->batCacheid);
1187 0 : return filterJoinGeomGeomDoubleToBit(lres_id,rres_id,l_id,r_id,distance_within,ls_id,rs_id,*nil_matches,estimate,*anti,geosDistanceWithin,"geom.wkbDWithinGeographicJoin");
1188 : }
1189 :
1190 : str
1191 0 : wkbDWithinGeographicSelect(bat* outid, const bat *bid , const bat *sid, wkb * const *wkb_const, const dbl *distance_within, const bit *anti) {
1192 0 : return filterSelectGeomGeomDoubleToBit(outid,bid,sid,*wkb_const,*distance_within,*anti,geosDistanceWithin,"geom.wkbDWithinGeographicSelect");
1193 : }
1194 :
1195 : str
1196 0 : wkbIntersectsGeographicJoin(bat *lres_id, bat *rres_id, const bat *l_id, const bat *r_id, const bat *ls_id, const bat *rs_id, const bit *nil_matches, const lng *estimate, const bit *anti) {
1197 0 : return filterJoinGeomGeomDoubleToBit(lres_id,rres_id,l_id,r_id,0,ls_id,rs_id,*nil_matches,estimate,*anti,geosDistanceWithin,"geom.wkbIntersectsGeographicJoin");
1198 : }
1199 :
1200 : str
1201 0 : wkbIntersectsGeographicSelect(bat* outid, const bat *bid , const bat *sid, wkb * const *wkb_const, const bit *anti) {
1202 0 : return filterSelectGeomGeomDoubleToBit(outid,bid,sid,*wkb_const,0,*anti,geosDistanceWithin,"geom.wkbIntersectsGeographicSelect");
1203 : }
1204 :
1205 : static inline CartPoint3D
1206 0 : crossProduct (const CartPoint3D *p1, const CartPoint3D *p2)
1207 : {
1208 0 : CartPoint3D p3;
1209 0 : p3.x = p1->y * p2->z - p1->z * p2->y;
1210 0 : p3.y = p1->z * p2->x - p1->x * p2->z;
1211 0 : p3.z = p1->x * p2->y - p1->y * p2->x;
1212 0 : return p3;
1213 : }
1214 :
1215 : static inline double
1216 0 : dotProduct (const CartPoint3D *p1, const CartPoint3D *p2)
1217 : {
1218 0 : return (p1->x * p2->x) + (p1->y * p2->y) + (p1->z * p2->z);
1219 : }
1220 :
1221 : //
1222 : static inline int
1223 0 : angleRotation (const CartPoint2D p1, const CartPoint2D p2, const CartPoint2D p3)
1224 : {
1225 : // vector P = P1P2
1226 : // vector Q = P1P3
1227 : // side = || Q x P ||
1228 : // Q and P are 2D vectors, so z component is 0
1229 0 : double side = ((p3.x - p1.x) * (p2.y - p1.y) - (p2.x - p1.x) * (p3.y - p1.y));
1230 : //
1231 0 : if (side < 0)
1232 : return -1;
1233 0 : else if (side > 0)
1234 : return 1;
1235 : else
1236 0 : return 0;
1237 : }
1238 :
1239 : static void
1240 0 : normalize2D (CartPoint2D *p) {
1241 0 : double c = sqrt(p->x * p->x + p->y * p->y);
1242 0 : if (c == 0) {
1243 0 : p->x = p->y = 0;
1244 0 : return;
1245 : }
1246 0 : p->x = p->x / c;
1247 0 : p->y = p->y / c;
1248 : }
1249 :
1250 : static void
1251 0 : normalize3D (CartPoint3D *p) {
1252 0 : double c = sqrt(p->x * p->x + p->y * p->y + p->z * p->z);
1253 0 : if (c == 0) {
1254 0 : p->x = p->y = p->z = 0;
1255 0 : return;
1256 : }
1257 0 : p->x = p->x / c;
1258 0 : p->y = p->y / c;
1259 0 : p->z = p->z / c;
1260 : }
1261 :
1262 : static inline bool
1263 0 : FP_EQUALS (double x, double y)
1264 : {
1265 0 : return fabs(x-y) < 1e-12;
1266 : }
1267 :
1268 : str
1269 0 : geodeticEdgeBoundingBox(const CartPoint3D* p1, const CartPoint3D* p2, BoundingBox* mbox)
1270 : {
1271 0 : CartPoint3D p3, pn, ep3d;
1272 0 : CartPoint3D e[6];
1273 0 : CartPoint2D s1, s2, ep, o;
1274 0 : int rotation_to_origin, rotation_to_ep;
1275 :
1276 : // check coinciding points
1277 0 : if (FP_EQUALS(p1->x,p2->x) && FP_EQUALS(p1->y,p2->y) && FP_EQUALS(p1->z,p2->z))
1278 : return MAL_SUCCEED;
1279 : // check antipodal points
1280 0 : if (FP_EQUALS(p1->x,-p2->x) && FP_EQUALS(p1->y,-p2->y) && FP_EQUALS(p1->z,-p2->z))
1281 0 : throw(MAL, "geom.geodeticEdgeBoundingBox", SQLSTATE(38000) "Antipodal edge");
1282 :
1283 : // create the great circle plane coord system (p1, p3)
1284 : // pn = p1 x p2
1285 : // p3 = pn x p1
1286 : // TODO handle the narrow and wide angle cases
1287 0 : pn = crossProduct(p1,p2);
1288 0 : normalize3D(&pn);
1289 0 : p3 = crossProduct(&pn,p1);
1290 0 : normalize3D(&p3);
1291 :
1292 : // represent p1, p2 with (s1, s2) 2-D space
1293 : // s1.x = 1, s1.y = 0
1294 : // s2.x = p2 * p1, s2.y = p2 * p3
1295 0 : s1.x = 1;
1296 0 : s1.y = 0;
1297 0 : s2.x = dotProduct(p2, p1);
1298 0 : s2.y = dotProduct(p2, &p3);
1299 : // 2-D space origin
1300 : // O.x = 0, O.y = 0
1301 0 : o.x = 0;
1302 0 : o.y = 0;
1303 :
1304 : // create 3D endpoints E.x, E.-x, ...
1305 : // E.x = (1, 0, 0), E.-x = (-1, 0, 0) ...
1306 0 : memset(e, 0, sizeof(CartPoint2D) * 6);
1307 0 : e[0].x = e[1].y = e[2].z = 1;
1308 0 : e[3].x = e[4].y = e[5].z = -1;
1309 :
1310 : // find the rotation between s1->s2 and s1->O
1311 : // rot = norm( vec(s1,s2) x vec(s1,0))
1312 0 : rotation_to_origin = angleRotation(s1,s2,o);
1313 :
1314 : // for every endpoint E
1315 0 : for (int i = 0; i < 6; i++) {
1316 : // project the endpoint in the 2-D space
1317 0 : ep.x = dotProduct(&e[i],p1);
1318 0 : ep.y = dotProduct(&e[i],&p3);
1319 : // re-normalize it e.g. EP (for endpoint_projection)
1320 0 : normalize2D(&ep);
1321 : // ep_rot = norm( vec(s1,s2) x vec(s1,EP_end) )
1322 0 : rotation_to_ep = angleRotation(s1,s2,ep);
1323 0 : if (rotation_to_origin != rotation_to_ep) {
1324 : // convert the 2-D EP into 3-D space
1325 0 : ep3d.x = ep.x * p1->x + ep.y * p3.x;
1326 0 : ep3d.y = ep.x * p1->y + ep.y * p3.y;
1327 0 : ep3d.z = ep.x * p1->z + ep.y * p3.z;
1328 : // expand the mbox in order to include 3-D representation of EP
1329 0 : boundingBoxAddPoint(mbox,ep3d);
1330 : }
1331 : }
1332 : return MAL_SUCCEED;
1333 : }
|
