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 = GDKzalloc(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->xmin = bb->xmax = c.x;
366 0 : bb->ymin = bb->ymax = c.y;
367 0 : bb->zmin = bb->zmax = c.z;
368 :
369 0 : for (int i = 1; i < lines.pointCount; i++) {
370 0 : c = geo2cartFromDegrees(lines.points[i]);
371 0 : boundingBoxAddPoint(bb, c);
372 : }
373 : return bb;
374 : }
375 :
376 : static int
377 0 : boundingBoxContainsPoint(BoundingBox bb, CartPoint3D pt)
378 : {
379 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;
380 : }
381 :
382 : /* Returns a point outside of the polygon's bounding box, for Point-In-Polygon calculation */
383 : static GeoPoint
384 0 : pointOutsidePolygon(GeoPolygon polygon)
385 : {
386 0 : BoundingBox bb = *polygon.bbox;
387 0 : BoundingBox bb2 = *polygon.bbox;
388 :
389 : //TODO: From POSTGIS -> CHANGE
390 0 : double grow = M_PI / 180.0 / 60.0;
391 0 : CartPoint3D corners[8];
392 0 : while (grow < M_PI) {
393 0 : if (bb.xmin > -1)
394 0 : bb.xmin -= grow;
395 0 : if (bb.ymin > -1)
396 0 : bb.ymin -= grow;
397 0 : if (bb.zmin > -1)
398 0 : bb.zmin -= grow;
399 0 : if (bb.xmax < 1)
400 0 : bb.xmax += grow;
401 0 : if (bb.ymax < 1)
402 0 : bb.ymax += grow;
403 0 : if (bb.zmax < 1)
404 0 : bb.zmax += grow;
405 :
406 0 : corners[0].x = bb.xmin;
407 0 : corners[0].y = bb.ymin;
408 0 : corners[0].z = bb.zmin;
409 :
410 0 : corners[1].x = bb.xmin;
411 0 : corners[1].y = bb.ymax;
412 0 : corners[1].z = bb.zmin;
413 :
414 0 : corners[2].x = bb.xmin;
415 0 : corners[2].y = bb.ymin;
416 0 : corners[2].z = bb.zmax;
417 :
418 0 : corners[3].x = bb.xmax;
419 0 : corners[3].y = bb.ymin;
420 0 : corners[3].z = bb.zmin;
421 :
422 0 : corners[4].x = bb.xmax;
423 0 : corners[4].y = bb.ymax;
424 0 : corners[4].z = bb.zmin;
425 :
426 0 : corners[5].x = bb.xmax;
427 0 : corners[5].y = bb.ymin;
428 0 : corners[5].z = bb.zmax;
429 :
430 0 : corners[6].x = bb.xmin;
431 0 : corners[6].y = bb.ymax;
432 0 : corners[6].z = bb.zmax;
433 :
434 0 : corners[7].x = bb.xmax;
435 0 : corners[7].y = bb.ymax;
436 0 : corners[7].z = bb.zmax;
437 :
438 0 : for (int i = 0; i < 8; i++)
439 0 : if (!boundingBoxContainsPoint(bb2, corners[i])) {
440 0 : CartPoint3D pt_cart = corners[i];
441 0 : return rad2DegPoint(cart2geo(pt_cart));
442 : }
443 0 : grow *= 2.0;
444 : }
445 : //TODO: Should this be the return value in case no point is found?
446 0 : return geoPointFromLatLon(0, 0);
447 : }
448 :
449 : /**
450 : * Distance functions
451 : *
452 : **/
453 : /**
454 : * The haversine formula calculate the distance in meters between two lat/lon points
455 : * The points must be measured in radians.
456 : * This formula assumes a spherical model of the earth, which can lead to an error of about 0.3% compared to a ellipsoidal model.
457 : */
458 : static double
459 0 : haversine(GeoPoint a, GeoPoint b)
460 : {
461 0 : double d_lon = b.lon - a.lon;
462 0 : double d_lat = b.lat - a.lat;
463 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);
464 0 : double c = 2 * atan2(sqrt(d), sqrt(1 - d));
465 : //TODO: Same as the previous line (which one is best?) -> 2 * asin(sqrt(d));
466 0 : return earth_radius_meters * c;
467 : }
468 :
469 : /* Distance between two Points */
470 : static double
471 0 : geoDistancePointPoint(GeoPoint a, GeoPoint b)
472 : {
473 0 : return haversine(deg2RadPoint(a), deg2RadPoint(b));
474 : }
475 :
476 : /* Calculates the distance between the perpendicular projection of a point in the Line */
477 : static double
478 0 : calculatePerpendicularDistance(GeoPoint p_geo, GeoPoint l1_geo, GeoPoint l2_geo)
479 : {
480 0 : CartPoint3D l1, l2, p, projection;
481 0 : GeoPoint projection_geo;
482 :
483 : //First, convert the points to 3D cartesian coordinates
484 0 : l1 = geo2cartFromDegrees(l1_geo);
485 0 : l2 = geo2cartFromDegrees(l2_geo);
486 0 : p = geo2cartFromDegrees(p_geo);
487 :
488 : //Calculate the projection of point into the line
489 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);
490 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;
491 :
492 : //If t is not between 0 and 1, the projected point is not in the line, so there is no perpendicular, return huge number
493 0 : if (t < 0 || t > 1)
494 : return INT_MAX;
495 :
496 : //If the projection is in the line segment, build the point -> projection = l1 + t * (l2-l1)
497 0 : projection = cartPointFromXYZ(l1.x + t * (l2.x - l1.x), l1.y + t * (l2.y - l1.y), l1.z + t * (l2.z - l1.z));
498 :
499 : //Convert into geographic coordinates (radians)
500 0 : projection_geo = cart2geo(projection);
501 :
502 : //Calculate distance from original point to the projection
503 0 : return haversine(deg2RadPoint(p_geo), projection_geo);
504 : }
505 :
506 : /* Distance between Point and Line
507 : The returned distance is the minimum distance between the point and the line vertices
508 : and the perpendicular projection of the point in each line segment. */
509 : static double
510 0 : geoDistancePointLine(GeoPoint point, GeoLines lines, double distance_min_limit)
511 : {
512 0 : double distancePoint, distancePerpendicular, min_distance = INT_MAX;
513 0 : for (int i = 0; i < lines.pointCount-1; i++) {
514 0 : distancePoint = geoDistancePointPoint(point, lines.points[i]);
515 0 : distancePerpendicular = calculatePerpendicularDistance(point,lines.points[i],lines.points[i+1]);
516 0 : if (distancePoint < min_distance)
517 0 : min_distance = distancePoint;
518 0 : if (distancePerpendicular < min_distance)
519 0 : min_distance = distancePerpendicular;
520 : //Shortcut, if the geometries are already at their minimum distance
521 0 : if (min_distance <= distance_min_limit)
522 0 : return min_distance;
523 : }
524 0 : distancePoint = geoDistancePointPoint(point, lines.points[lines.pointCount-1]);
525 0 : return distancePoint < min_distance ? distancePoint : min_distance;
526 : }
527 :
528 : /* Distance between two Lines. */
529 : static double
530 0 : geoDistanceLineLine(GeoLines line1, GeoLines line2, double distance_min_limit)
531 : {
532 0 : double distance, min_distance = INT_MAX;
533 0 : for (int i = 0; i < line1.pointCount; i++) {
534 0 : distance = geoDistancePointLine(line1.points[i], line2, distance_min_limit);
535 0 : if (distance < min_distance)
536 0 : min_distance = distance;
537 : //Shortcut, if the geometries are already at their minimum distance
538 0 : if (min_distance <= distance_min_limit)
539 0 : return min_distance;
540 : }
541 0 : for (int i = 0; i < line2.pointCount; i++) {
542 0 : for (int j = 0; j < line1.pointCount - 1; j++) {
543 0 : distance = calculatePerpendicularDistance(line2.points[i],line1.points[j],line1.points[j+1]);
544 0 : if (distance < min_distance)
545 0 : min_distance = distance;
546 : //Shortcut, if the geometries are already at their minimum distance
547 0 : if (min_distance <= distance_min_limit)
548 0 : return min_distance;
549 : }
550 : }
551 : return min_distance;
552 : }
553 :
554 : /* Checks if a Point is within a Polygon */
555 : static bool
556 0 : pointWithinPolygon(GeoPolygon polygon, GeoPoint point)
557 : {
558 0 : int intersectionNum = 0;
559 0 : GEOSGeometry *segmentPolygon, *intersectionPoints;
560 0 : GeoLines polygonRing;
561 :
562 : //Get an point that's outside the polygon
563 0 : GeoPoint outsidePoint = pointOutsidePolygon(polygon);
564 :
565 : //No outside point was found, return false
566 0 : if (outsidePoint.lat == 0 && outsidePoint.lon == 0)
567 : return false;
568 :
569 : /*printf("Outside point: (%f %f)\n",outsidePoint.lon, outsidePoint.lat);
570 : fflush(stdout);*/
571 :
572 : //Construct a line between the outside point and the input point
573 0 : GEOSGeometry *outInLine = cartesianLineFromGeoPoints(point, outsidePoint);
574 :
575 : //TODO This is producing wrong results, review the intersection conditional
576 : //Count the number of intersections between the polygon exterior ring and the constructed line
577 0 : polygonRing = polygon.exteriorRing;
578 0 : for (int i = 0; i < polygonRing.pointCount-1; i++) {
579 0 : segmentPolygon = cartesianLineFromGeoPoints(polygonRing.points[i], polygonRing.points[i+1]);
580 0 : intersectionPoints = GEOSIntersection_r(geoshandle, segmentPolygon, outInLine);
581 :
582 : //If there is an intersection, a point will be returned (line when there is none)
583 0 : if (GEOSGeomTypeId_r(geoshandle, intersectionPoints) == GEOS_POINT)
584 0 : intersectionNum++;
585 :
586 0 : if (intersectionPoints != NULL)
587 0 : GEOSGeom_destroy_r(geoshandle, intersectionPoints);
588 0 : if (segmentPolygon != NULL)
589 0 : GEOSGeom_destroy_r(geoshandle, segmentPolygon);
590 : }
591 :
592 : //Count the number of intersections between the polygon interior rings and the constructed line
593 0 : for (int j = 0; j < polygon.interiorRingsCount; j++) {
594 0 : polygonRing = polygon.interiorRings[j];
595 0 : for (int i = 0; i < polygonRing.pointCount-1; i++) {
596 0 : segmentPolygon = cartesianLineFromGeoPoints(polygonRing.points[i], polygonRing.points[i+1]);
597 0 : intersectionPoints = GEOSIntersection_r(geoshandle, segmentPolygon, outInLine);
598 :
599 : //If there is an intersection, a point will be returned (line when there is none)
600 0 : if (GEOSGeomTypeId_r(geoshandle, intersectionPoints) == GEOS_POINT)
601 0 : intersectionNum++;
602 :
603 0 : if (intersectionPoints != NULL)
604 0 : GEOSGeom_destroy_r(geoshandle, intersectionPoints);
605 0 : if (segmentPolygon != NULL)
606 0 : GEOSGeom_destroy_r(geoshandle, segmentPolygon);
607 : }
608 : }
609 :
610 0 : if (outInLine != NULL)
611 0 : GEOSGeom_destroy_r(geoshandle, outInLine);
612 :
613 : //If even, the point is not within the polygon. If odd, it is within
614 0 : return intersectionNum % 2 == 1;
615 : }
616 :
617 : /* Distance between Point and Polygon.*/
618 : static double
619 0 : geoDistancePointPolygon(GeoPoint point, GeoPolygon polygon, double distance_min_limit)
620 : {
621 : //Check if point is in polygon
622 0 : if (pointWithinPolygon(polygon, point))
623 : return 0;
624 :
625 : //Calculate distance from Point to the exterior and interior rings of the polygon
626 0 : double distance, min_distance = INT_MAX;
627 : //First, calculate distance to the exterior ring
628 0 : min_distance = geoDistancePointLine(point, polygon.exteriorRing, distance_min_limit);
629 : //Then, calculate distance to the interior rings
630 0 : for (int i = 0; i < polygon.interiorRingsCount; i++) {
631 : //Shortcut, if the geometries are already at their minimum distance
632 0 : if (min_distance <= distance_min_limit)
633 0 : return min_distance;
634 0 : distance = geoDistancePointLine(point, polygon.interiorRings[i], distance_min_limit);
635 0 : if (distance < min_distance)
636 0 : min_distance = distance;
637 : }
638 : return min_distance;
639 : }
640 :
641 : /* Distance between Line and Polygon. */
642 : static double
643 0 : geoDistanceLinePolygon(GeoLines line, GeoPolygon polygon, double distance_min_limit)
644 : {
645 0 : double distance, min_distance = INT_MAX;
646 : //Calculate distance to all start vertices of the line
647 0 : for (int i = 0; i < line.pointCount; i++) {
648 0 : distance = geoDistancePointPolygon(line.points[i], polygon, distance_min_limit);
649 :
650 : //Short-cut in case the point is within the polygon
651 0 : if (distance <= distance_min_limit)
652 0 : return distance;
653 :
654 0 : if (distance < min_distance)
655 0 : min_distance = distance;
656 : }
657 : return min_distance;
658 : }
659 :
660 : /* Distance between two Polygons. */
661 : static double
662 0 : geoDistancePolygonPolygon(GeoPolygon polygon1, GeoPolygon polygon2, double distance_min_limit)
663 : {
664 0 : double distance1, distance2;
665 : //Calculate the distance between the exterior ring of polygon1 and all segments of polygon2 (including the interior rings)
666 0 : distance1 = geoDistanceLinePolygon(polygon1.exteriorRing, polygon2, distance_min_limit);
667 : //Shortcut, if the geometries are already at their minimum distance
668 0 : if (distance1 <= distance_min_limit)
669 : return distance1;
670 0 : distance2 = geoDistanceLinePolygon(polygon2.exteriorRing, polygon1, distance_min_limit);
671 0 : return distance1 < distance2 ? distance1 : distance2;
672 : }
673 :
674 : /* Distance between two (non-collection) geometries. */
675 : static double
676 0 : geoDistanceSingle(GEOSGeom aGeom, GEOSGeom bGeom, double distance_min_limit)
677 : {
678 0 : int dimA, dimB;
679 0 : double distance = INT_MAX;
680 0 : dimA = GEOSGeom_getDimensions_r(geoshandle, aGeom);
681 0 : dimB = GEOSGeom_getDimensions_r(geoshandle, bGeom);
682 0 : if (dimA == 0 && dimB == 0) {
683 : /* Point and Point */
684 0 : GeoPoint a = geoPointFromGeom(aGeom);
685 0 : GeoPoint b = geoPointFromGeom(bGeom);
686 0 : distance = geoDistancePointPoint(a, b);
687 0 : } else if (dimA == 0 && dimB == 1) {
688 : /* Point and Line/LinearRing */
689 0 : GeoPoint a = geoPointFromGeom(aGeom);
690 0 : GeoLines b = geoLinesFromGeom(bGeom);
691 0 : distance = geoDistancePointLine(a, b, distance_min_limit);
692 0 : freeGeoLines(b);
693 0 : } else if (dimA == 1 && dimB == 0) {
694 : /* Line/LinearRing and Point */
695 0 : GeoLines a = geoLinesFromGeom(aGeom);
696 0 : GeoPoint b = geoPointFromGeom(bGeom);
697 0 : distance = geoDistancePointLine(b, a, distance_min_limit);
698 0 : freeGeoLines(a);
699 0 : } else if (dimA == 1 && dimB == 1) {
700 : /* Line/LinearRing and Line/LinearRing */
701 0 : GeoLines a = geoLinesFromGeom(aGeom);
702 0 : GeoLines b = geoLinesFromGeom(bGeom);
703 0 : distance = geoDistanceLineLine(a, b, distance_min_limit);
704 0 : freeGeoLines(a);
705 0 : freeGeoLines(b);
706 0 : } else if (dimA == 0 && dimB == 2) {
707 : /* Point and Polygon */
708 0 : GeoPoint a = geoPointFromGeom(aGeom);
709 0 : GeoPolygon b = geoPolygonFromGeom(bGeom);
710 0 : distance = geoDistancePointPolygon(a, b, distance_min_limit);
711 0 : freeGeoPolygon(b);
712 0 : } else if (dimA == 2 && dimB == 0) {
713 : /* Polygon and Point */
714 0 : GeoPolygon a = geoPolygonFromGeom(aGeom);
715 0 : GeoPoint b = geoPointFromGeom(bGeom);
716 0 : distance = geoDistancePointPolygon(b, a, distance_min_limit);
717 0 : freeGeoPolygon(a);
718 0 : } else if (dimA == 1 && dimB == 2) {
719 : /* Line/LinearRing and Polygon */
720 0 : GeoLines a = geoLinesFromGeom(aGeom);
721 0 : GeoPolygon b = geoPolygonFromGeom(bGeom);
722 0 : distance = geoDistanceLinePolygon(a, b, distance_min_limit);
723 0 : freeGeoLines(a);
724 0 : freeGeoPolygon(b);
725 0 : } else if (dimA == 2 && dimB == 1) {
726 : /* Polygon and Line/LinearRing */
727 0 : GeoPolygon a = geoPolygonFromGeom(aGeom);
728 0 : GeoLines b = geoLinesFromGeom(bGeom);
729 0 : distance = geoDistanceLinePolygon(b, a, distance_min_limit);
730 0 : freeGeoPolygon(a);
731 0 : freeGeoLines(b);
732 0 : } else if (dimA == 2 && dimB == 2) {
733 : /* Polygon and Polygon */
734 0 : GeoPolygon a = geoPolygonFromGeom(aGeom);
735 0 : GeoPolygon b = geoPolygonFromGeom(bGeom);
736 0 : distance = geoDistancePolygonPolygon(a, b, distance_min_limit);
737 0 : freeGeoPolygon(a);
738 0 : freeGeoPolygon(b);
739 : }
740 0 : return distance;
741 : }
742 :
743 : //The distance_min_limit argument is used for DWithin and Intersects.
744 : //If we get to the minimum distance for the predicate, return immediately
745 : //It is equal to 0 if the operation is Distance
746 : static double
747 0 : geoDistanceInternal(GEOSGeom a, GEOSGeom b, double distance_min_limit)
748 : {
749 0 : int numGeomsA = GEOSGetNumGeometries_r(geoshandle, a), numGeomsB = GEOSGetNumGeometries_r(geoshandle, b);
750 0 : double distance, min_distance = INT_MAX;
751 0 : GEOSGeometry *geo1, *geo2;
752 0 : for (int i = 0; i < numGeomsA; i++) {
753 0 : geo1 = (GEOSGeometry *)GEOSGetGeometryN_r(geoshandle, (const GEOSGeometry *)a, i);
754 0 : for (int j = 0; j < numGeomsB; j++) {
755 0 : geo2 = (GEOSGeometry *)GEOSGetGeometryN_r(geoshandle, (const GEOSGeometry *)b, j);
756 0 : distance = geoDistanceSingle(geo1, geo2, distance_min_limit);
757 : //Shortcut, if the geometries are already at their minimum distance (0 in the case of normal Distance)
758 0 : if (distance <= distance_min_limit)
759 0 : return distance;
760 0 : if (distance < min_distance)
761 0 : min_distance = distance;
762 : }
763 : }
764 : return min_distance;
765 : }
766 :
767 : /**
768 : * Distance
769 : *
770 : **/
771 : /* Calculates the distance, in meters, between two geographic geometries with latitude/longitude coordinates */
772 : str
773 0 : wkbDistanceGeographic(dbl *out, wkb * const *a, wkb * const *b)
774 : {
775 0 : str err = MAL_SUCCEED;
776 0 : GEOSGeom ga, gb;
777 0 : err = wkbGetCompatibleGeometries(a, b, &ga, &gb);
778 0 : if (ga && gb) {
779 0 : (*out) = geoDistanceInternal(ga, gb, 0);
780 : }
781 0 : GEOSGeom_destroy_r(geoshandle, ga);
782 0 : GEOSGeom_destroy_r(geoshandle, gb);
783 0 : return err;
784 : }
785 :
786 : /**
787 : * Distance Within
788 : *
789 : **/
790 : /* Checks if two geographic geometries are within d meters of one another */
791 : str
792 0 : wkbDWithinGeographic(bit *out, wkb * const *a, wkb * const *b, const dbl *d)
793 : {
794 0 : str err = MAL_SUCCEED;
795 0 : GEOSGeom ga, gb;
796 0 : double distance;
797 0 : err = wkbGetCompatibleGeometries(a, b, &ga, &gb);
798 0 : if (ga && gb) {
799 0 : distance = geoDistanceInternal(ga, gb, *d);
800 0 : (*out) = (distance <= (*d));
801 : }
802 0 : GEOSGeom_destroy_r(geoshandle, ga);
803 0 : GEOSGeom_destroy_r(geoshandle, gb);
804 0 : return err;
805 : }
806 :
807 : /**
808 : * Intersects
809 : *
810 : **/
811 : /* Checks if two geographic geometries intersect at any point */
812 : str
813 0 : wkbIntersectsGeographic(bit *out, wkb * const *a, wkb * const *b)
814 : {
815 0 : str err = MAL_SUCCEED;
816 0 : GEOSGeom ga, gb;
817 0 : double distance;
818 0 : err = wkbGetCompatibleGeometries(a, b, &ga, &gb);
819 0 : if (ga && gb) {
820 0 : distance = geoDistanceInternal(ga, gb, 0);
821 0 : (*out) = (distance == 0);
822 : }
823 0 : GEOSGeom_destroy_r(geoshandle, ga);
824 0 : GEOSGeom_destroy_r(geoshandle, gb);
825 0 : return err;
826 : }
827 :
828 : /* Checks if a Polygon covers a Line geometry */
829 : static bool
830 0 : geoPolygonCoversLine(GeoPolygon polygon, GeoLines lines)
831 : {
832 0 : for (int i = 0; i < lines.pointCount; i++) {
833 0 : if (pointWithinPolygon(polygon, lines.points[i]) == false)
834 : return false;
835 : }
836 : return true;
837 : }
838 :
839 : /* Compares two GeoPoints, returns true if they're equal */
840 : static bool
841 0 : geoPointEquals(GeoPoint pointA, GeoPoint pointB)
842 : {
843 0 : return (pointA.lat == pointB.lat) && (pointA.lon = pointB.lon);
844 : }
845 :
846 : //TODO Check if this works correctly
847 : static bool
848 0 : geoCoversSingle(GEOSGeom a, GEOSGeom b)
849 : {
850 0 : int dimA = GEOSGeom_getDimensions_r(geoshandle, a), dimB = GEOSGeom_getDimensions_r(geoshandle, b);
851 0 : if (dimA < dimB)
852 : //If the dimension of A is smaller than B, then it must not cover it
853 : return false;
854 :
855 0 : if (dimA == 0){
856 : //A and B are Points
857 0 : GeoPoint pointA = geoPointFromGeom(a);
858 0 : GeoPoint pointB = geoPointFromGeom(b);
859 0 : return geoPointEquals(pointA, pointB);
860 0 : } else if (dimA == 1) {
861 : //A is Line
862 : //GeoLines lineA = geoLinesFromGeom(a);
863 : if (dimB == 0) {
864 : //B is Point
865 : //GeoPoint pointB = geoPointFromGeom(b);
866 : //return geoLineCoversPoint(lineA,pointB);
867 : return false;
868 : } else {
869 : //B is Line
870 : //GeoLines lineB = geoLinesFromGeom(b);
871 : //return geoLineCoversLine(lineA, lineB);
872 : return false;
873 : }
874 0 : } else if (dimA == 2) {
875 : //A is Polygon
876 0 : GeoPolygon polygonA = geoPolygonFromGeom(a);
877 0 : if (dimB == 0){
878 : //B is Point
879 0 : GeoPoint pointB = geoPointFromGeom(b);
880 0 : return pointWithinPolygon(polygonA, pointB);
881 0 : } else if (dimB == 1) {
882 : //B is Line
883 0 : GeoLines lineB = geoLinesFromGeom(b);
884 0 : return geoPolygonCoversLine(polygonA, lineB);
885 : } else {
886 : //B is Polygon
887 0 : GeoPolygon polygonB = geoPolygonFromGeom(b);
888 : //If every point in the exterior ring of B is covered, polygon B is covered by polygon A
889 0 : return geoPolygonCoversLine(polygonA, polygonB.exteriorRing);
890 : }
891 : } else
892 : return false;
893 : }
894 :
895 : static bool
896 0 : geoCoversInternal(GEOSGeom a, GEOSGeom b)
897 : {
898 0 : int numGeomsA = GEOSGetNumGeometries_r(geoshandle, a), numGeomsB = GEOSGetNumGeometries_r(geoshandle, b);
899 0 : GEOSGeometry *geo1, *geo2;
900 0 : for (int i = 0; i < numGeomsA; i++) {
901 0 : geo1 = (GEOSGeometry *)GEOSGetGeometryN_r(geoshandle, (const GEOSGeometry *)a, i);
902 0 : for (int j = 0; j < numGeomsB; j++) {
903 0 : geo2 = (GEOSGeometry *)GEOSGetGeometryN_r(geoshandle, (const GEOSGeometry *)b, j);
904 0 : if (geoCoversSingle(geo1, geo2) == 0)
905 : return 0;
906 : }
907 : }
908 : return 1;
909 : }
910 :
911 : /**
912 : * Covers
913 : *
914 : **/
915 : /* Checks if no point of Geometry B is outside Geometry A */
916 : str
917 0 : wkbCoversGeographic(bit *out, wkb * const *a, wkb * const *b)
918 : {
919 0 : str err = MAL_SUCCEED;
920 0 : GEOSGeom ga, gb;
921 0 : err = wkbGetCompatibleGeometries(a, b, &ga, &gb);
922 0 : if (ga && gb)
923 0 : (*out) = geoCoversInternal(ga, gb);
924 :
925 0 : GEOSGeom_destroy_r(geoshandle, ga);
926 0 : GEOSGeom_destroy_r(geoshandle, gb);
927 :
928 0 : return err;
929 : }
930 :
931 : /**
932 : * FILTER FUNCTIONS
933 : **/
934 :
935 : static inline bit
936 0 : geosDistanceWithin (GEOSGeom geom1, GEOSGeom geom2, dbl distance_within) {
937 0 : dbl distance = geoDistanceInternal(geom1, geom2, distance_within);
938 0 : return distance <= distance_within;
939 : }
940 :
941 : //TODO Change BUNappend with manual insertion into the result BAT
942 : static str
943 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)
944 : {
945 0 : BAT *out = NULL, *b = NULL, *s = NULL;
946 0 : BATiter b_iter;
947 0 : struct canditer ci;
948 0 : GEOSGeom col_geom, const_geom;
949 :
950 : //Check if the geometry is null and convert to GEOS
951 0 : if ((const_geom = wkb2geos(wkb_const)) == NULL) {
952 0 : if ((out = BATdense(0, 0, 0)) == NULL)
953 0 : throw(MAL, name, GDK_EXCEPTION);
954 0 : *outid = out->batCacheid;
955 0 : BBPkeepref(out);
956 0 : return MAL_SUCCEED;
957 : }
958 : //get the BATs
959 0 : if ((b = BATdescriptor(*bid)) == NULL)
960 0 : throw(MAL, name, SQLSTATE(HY002) RUNTIME_OBJECT_MISSING);
961 : //get the candidate lists
962 0 : if (sid && !is_bat_nil(*sid) && !(s = BATdescriptor(*sid))) {
963 0 : BBPunfix(b->batCacheid);
964 0 : throw(MAL, name, SQLSTATE(HY002) RUNTIME_OBJECT_MISSING);
965 : }
966 0 : canditer_init(&ci, b, s);
967 : //create a new BAT for the output
968 0 : if ((out = COLnew(0, ATOMindex("oid"), ci.ncand, TRANSIENT)) == NULL) {
969 0 : BBPunfix(b->batCacheid);
970 0 : if (s)
971 0 : BBPunfix(s->batCacheid);
972 0 : throw(MAL, name, SQLSTATE(HY013) MAL_MALLOC_FAIL);
973 : }
974 0 : b_iter = bat_iterator(b);
975 : //Loop through column and compare with constant
976 0 : for (BUN i = 0; i < ci.ncand; i++) {
977 0 : oid c_oid = canditer_next(&ci);
978 0 : const wkb *col_wkb = BUNtvar(b_iter, c_oid - b->hseqbase);
979 0 : if ((col_geom = wkb2geos(col_wkb)) == NULL)
980 0 : continue;
981 0 : if (GEOSGetSRID_r(geoshandle, col_geom) != GEOSGetSRID_r(geoshandle, const_geom)) {
982 0 : GEOSGeom_destroy_r(geoshandle, col_geom);
983 0 : GEOSGeom_destroy_r(geoshandle, const_geom);
984 0 : bat_iterator_end(&b_iter);
985 0 : BBPunfix(b->batCacheid);
986 0 : if (s)
987 0 : BBPunfix(s->batCacheid);
988 0 : BBPreclaim(out);
989 0 : throw(MAL, name, SQLSTATE(38000) "Geometries of different SRID");
990 : }
991 : //Apply the (Geom, Geom, double) -> bit function
992 0 : bit cond = (*func)(col_geom, const_geom, double_flag);
993 0 : if (cond != anti) {
994 0 : if (BUNappend(out, &c_oid, false) != GDK_SUCCEED) {
995 0 : if (col_geom)
996 0 : GEOSGeom_destroy_r(geoshandle, col_geom);
997 0 : if (const_geom)
998 0 : GEOSGeom_destroy_r(geoshandle, const_geom);
999 0 : bat_iterator_end(&b_iter);
1000 0 : BBPunfix(b->batCacheid);
1001 0 : if (s)
1002 0 : BBPunfix(s->batCacheid);
1003 0 : BBPreclaim(out);
1004 0 : throw(MAL, name, SQLSTATE(HY013) MAL_MALLOC_FAIL);
1005 : }
1006 : }
1007 0 : GEOSGeom_destroy_r(geoshandle, col_geom);
1008 : }
1009 0 : GEOSGeom_destroy_r(geoshandle, const_geom);
1010 0 : bat_iterator_end(&b_iter);
1011 0 : BBPunfix(b->batCacheid);
1012 0 : if (s)
1013 0 : BBPunfix(s->batCacheid);
1014 0 : *outid = out->batCacheid;
1015 0 : BBPkeepref(out);
1016 0 : return MAL_SUCCEED;
1017 : }
1018 :
1019 : static str
1020 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)
1021 : {
1022 0 : BAT *lres = NULL, *rres = NULL, *l = NULL, *r = NULL, *ls = NULL, *rs = NULL;
1023 0 : BATiter l_iter, r_iter;
1024 0 : str msg = MAL_SUCCEED;
1025 0 : struct canditer l_ci, r_ci;
1026 0 : GEOSGeom l_geom, r_geom;
1027 0 : GEOSGeom *l_geoms = NULL, *r_geoms = NULL;
1028 0 : BUN est;
1029 :
1030 : //get the input BATs
1031 0 : l = BATdescriptor(*l_id);
1032 0 : r = BATdescriptor(*r_id);
1033 0 : if (l == NULL || r == NULL) {
1034 0 : if (l)
1035 0 : BBPunfix(l->batCacheid);
1036 0 : if (r)
1037 0 : BBPunfix(r->batCacheid);
1038 0 : throw(MAL, name, SQLSTATE(HY002) RUNTIME_OBJECT_MISSING);
1039 : }
1040 : //get the candidate lists
1041 0 : if (ls_id && !is_bat_nil(*ls_id) && !(ls = BATdescriptor(*ls_id)) && rs_id && !is_bat_nil(*rs_id) && !(rs = BATdescriptor(*rs_id))) {
1042 0 : msg = createException(MAL, name, SQLSTATE(HY002) RUNTIME_OBJECT_MISSING);
1043 0 : goto free;
1044 : }
1045 0 : canditer_init(&l_ci, l, ls);
1046 0 : canditer_init(&r_ci, r, rs);
1047 : //create new BATs for the output
1048 0 : est = is_lng_nil(*estimate) || *estimate == 0 ? l_ci.ncand : (BUN) *estimate;
1049 0 : if ((lres = COLnew(0, ATOMindex("oid"), est, TRANSIENT)) == NULL) {
1050 0 : msg = createException(MAL, name, SQLSTATE(HY013) MAL_MALLOC_FAIL);
1051 0 : goto free;
1052 : }
1053 0 : if ((rres = COLnew(0, ATOMindex("oid"), est, TRANSIENT)) == NULL) {
1054 0 : msg = createException(MAL, name, SQLSTATE(HY013) MAL_MALLOC_FAIL);
1055 0 : goto free;
1056 : }
1057 :
1058 : //Allocate arrays for reutilizing GEOS type conversion
1059 0 : if ((l_geoms = GDKmalloc(l_ci.ncand * sizeof(GEOSGeometry *))) == NULL || (r_geoms = GDKmalloc(r_ci.ncand * sizeof(GEOSGeometry *))) == NULL) {
1060 0 : msg = createException(MAL, name, SQLSTATE(HY013) MAL_MALLOC_FAIL);
1061 0 : goto free;
1062 : }
1063 :
1064 0 : l_iter = bat_iterator(l);
1065 0 : r_iter = bat_iterator(r);
1066 :
1067 : //Convert wkb to GEOS only once
1068 0 : for (BUN i = 0; i < l_ci.ncand; i++) {
1069 0 : oid l_oid = canditer_next(&l_ci);
1070 0 : l_geoms[i] = wkb2geos((const wkb*) BUNtvar(l_iter, l_oid - l->hseqbase));
1071 : }
1072 0 : for (BUN j = 0; j < r_ci.ncand; j++) {
1073 0 : oid r_oid = canditer_next(&r_ci);
1074 0 : r_geoms[j] = wkb2geos((const wkb*)BUNtvar(r_iter, r_oid - r->hseqbase));
1075 : }
1076 :
1077 0 : canditer_reset(&l_ci);
1078 0 : for (BUN i = 0; i < l_ci.ncand; i++) {
1079 0 : oid l_oid = canditer_next(&l_ci);
1080 0 : l_geom = l_geoms[i];
1081 0 : if (!nil_matches && l_geom == NULL)
1082 0 : continue;
1083 0 : canditer_reset(&r_ci);
1084 0 : for (BUN j = 0; j < r_ci.ncand; j++) {
1085 0 : oid r_oid = canditer_next(&r_ci);
1086 0 : r_geom = r_geoms[j];
1087 : //Null handling
1088 0 : if (r_geom == NULL) {
1089 0 : if (nil_matches && l_geom == NULL) {
1090 0 : if (BUNappend(lres, &l_oid, false) != GDK_SUCCEED || BUNappend(rres, &r_oid, false) != GDK_SUCCEED) {
1091 0 : msg = createException(MAL, name, SQLSTATE(HY013) MAL_MALLOC_FAIL);
1092 0 : bat_iterator_end(&l_iter);
1093 0 : bat_iterator_end(&r_iter);
1094 0 : goto free;
1095 : }
1096 : }
1097 : else
1098 0 : continue;
1099 : }
1100 : //TODO Do we need to do this check for every element?
1101 0 : if (GEOSGetSRID_r(geoshandle, l_geom) != GEOSGetSRID_r(geoshandle, r_geom)) {
1102 0 : msg = createException(MAL, name, SQLSTATE(38000) "Geometries of different SRID");
1103 0 : bat_iterator_end(&l_iter);
1104 0 : bat_iterator_end(&r_iter);
1105 0 : goto free;
1106 : }
1107 : //Apply the (Geom, Geom) -> bit function
1108 0 : bit cond = (*func)(l_geom, r_geom, double_flag);
1109 0 : if (cond != anti) {
1110 0 : if (BUNappend(lres, &l_oid, false) != GDK_SUCCEED || BUNappend(rres, &r_oid, false) != GDK_SUCCEED) {
1111 0 : msg = createException(MAL, name, SQLSTATE(HY013) MAL_MALLOC_FAIL);
1112 0 : bat_iterator_end(&l_iter);
1113 0 : bat_iterator_end(&r_iter);
1114 0 : goto free;
1115 : }
1116 : }
1117 : }
1118 : }
1119 : if (l_geoms) {
1120 0 : for (BUN i = 0; i < l_ci.ncand; i++) {
1121 0 : GEOSGeom_destroy_r(geoshandle, l_geoms[i]);
1122 : }
1123 0 : GDKfree(l_geoms);
1124 : }
1125 0 : if (r_geoms) {
1126 0 : for (BUN i = 0; i < r_ci.ncand; i++) {
1127 0 : GEOSGeom_destroy_r(geoshandle, r_geoms[i]);
1128 : }
1129 0 : GDKfree(r_geoms);
1130 : }
1131 0 : bat_iterator_end(&l_iter);
1132 0 : bat_iterator_end(&r_iter);
1133 0 : BBPunfix(l->batCacheid);
1134 0 : BBPunfix(r->batCacheid);
1135 0 : if (ls)
1136 0 : BBPunfix(ls->batCacheid);
1137 0 : if (rs)
1138 0 : BBPunfix(rs->batCacheid);
1139 0 : *lres_id = lres->batCacheid;
1140 0 : BBPkeepref(lres);
1141 0 : *rres_id = rres->batCacheid;
1142 0 : BBPkeepref(rres);
1143 0 : return MAL_SUCCEED;
1144 0 : free:
1145 0 : if (l_geoms) {
1146 0 : for (BUN i = 0; i < l_ci.ncand; i++) {
1147 0 : GEOSGeom_destroy_r(geoshandle, l_geoms[i]);
1148 : }
1149 0 : GDKfree(l_geoms);
1150 : }
1151 0 : if (r_geoms) {
1152 0 : for (BUN i = 0; i < r_ci.ncand; i++) {
1153 0 : GEOSGeom_destroy_r(geoshandle, r_geoms[i]);
1154 : }
1155 0 : GDKfree(r_geoms);
1156 : }
1157 0 : BBPunfix(l->batCacheid);
1158 0 : BBPunfix(r->batCacheid);
1159 0 : if (ls)
1160 0 : BBPunfix(ls->batCacheid);
1161 0 : if (rs)
1162 0 : BBPunfix(rs->batCacheid);
1163 0 : if (lres)
1164 0 : BBPreclaim(lres);
1165 0 : if (rres)
1166 0 : BBPreclaim(rres);
1167 : return msg;
1168 : }
1169 :
1170 : str
1171 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) {
1172 0 : double distance_within = 0;
1173 0 : BAT *d = NULL;
1174 : //Get the distance BAT and get the double value
1175 0 : if ((d = BATdescriptor(*d_id)) == NULL) {
1176 0 : throw(MAL, "geom.wkbDWithinGeographicJoin", SQLSTATE(HY002) RUNTIME_OBJECT_MISSING);
1177 : }
1178 0 : if (BATcount(d) == 1) {
1179 0 : distance_within = *((double*) Tloc(d, 0));
1180 : }
1181 0 : BBPunfix(d->batCacheid);
1182 0 : return filterJoinGeomGeomDoubleToBit(lres_id,rres_id,l_id,r_id,distance_within,ls_id,rs_id,*nil_matches,estimate,*anti,geosDistanceWithin,"geom.wkbDWithinGeographicJoin");
1183 : }
1184 :
1185 : str
1186 0 : wkbDWithinGeographicSelect(bat* outid, const bat *bid , const bat *sid, wkb * const *wkb_const, const dbl *distance_within, const bit *anti) {
1187 0 : return filterSelectGeomGeomDoubleToBit(outid,bid,sid,*wkb_const,*distance_within,*anti,geosDistanceWithin,"geom.wkbDWithinGeographicSelect");
1188 : }
1189 :
1190 : str
1191 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) {
1192 0 : return filterJoinGeomGeomDoubleToBit(lres_id,rres_id,l_id,r_id,0,ls_id,rs_id,*nil_matches,estimate,*anti,geosDistanceWithin,"geom.wkbIntersectsGeographicJoin");
1193 : }
1194 :
1195 : str
1196 0 : wkbIntersectsGeographicSelect(bat* outid, const bat *bid , const bat *sid, wkb * const *wkb_const, const bit *anti) {
1197 0 : return filterSelectGeomGeomDoubleToBit(outid,bid,sid,*wkb_const,0,*anti,geosDistanceWithin,"geom.wkbIntersectsGeographicSelect");
1198 : }
1199 :
1200 : static inline CartPoint3D
1201 0 : crossProduct (const CartPoint3D *p1, const CartPoint3D *p2)
1202 : {
1203 0 : CartPoint3D p3;
1204 0 : p3.x = p1->y * p2->z - p1->z * p2->y;
1205 0 : p3.y = p1->z * p2->x - p1->x * p2->z;
1206 0 : p3.z = p1->x * p2->y - p1->y * p2->x;
1207 0 : return p3;
1208 : }
1209 :
1210 : static inline double
1211 0 : dotProduct (const CartPoint3D *p1, const CartPoint3D *p2)
1212 : {
1213 0 : return (p1->x * p2->x) + (p1->y * p2->y) + (p1->z * p2->z);
1214 : }
1215 :
1216 : //
1217 : static inline int
1218 0 : angleRotation (const CartPoint2D p1, const CartPoint2D p2, const CartPoint2D p3)
1219 : {
1220 : // vector P = P1P2
1221 : // vector Q = P1P3
1222 : // side = || Q x P ||
1223 : // Q and P are 2D vectors, so z component is 0
1224 0 : double side = ((p3.x - p1.x) * (p2.y - p1.y) - (p2.x - p1.x) * (p3.y - p1.y));
1225 : //
1226 0 : if (side < 0)
1227 : return -1;
1228 0 : else if (side > 0)
1229 : return 1;
1230 : else
1231 0 : return 0;
1232 : }
1233 :
1234 : static void
1235 0 : normalize2D (CartPoint2D *p) {
1236 0 : double c = sqrt(p->x * p->x + p->y * p->y);
1237 0 : if (c == 0) {
1238 0 : p->x = p->y = 0;
1239 0 : return;
1240 : }
1241 0 : p->x = p->x / c;
1242 0 : p->y = p->y / c;
1243 : }
1244 :
1245 : static void
1246 0 : normalize3D (CartPoint3D *p) {
1247 0 : double c = sqrt(p->x * p->x + p->y * p->y + p->z * p->z);
1248 0 : if (c == 0) {
1249 0 : p->x = p->y = p->z = 0;
1250 0 : return;
1251 : }
1252 0 : p->x = p->x / c;
1253 0 : p->y = p->y / c;
1254 0 : p->z = p->z / c;
1255 : }
1256 :
1257 : static inline bool
1258 0 : FP_EQUALS (double x, double y)
1259 : {
1260 0 : return fabs(x-y) < 1e-12;
1261 : }
1262 :
1263 : str
1264 0 : geodeticEdgeBoundingBox(const CartPoint3D* p1, const CartPoint3D* p2, BoundingBox* mbox)
1265 : {
1266 0 : CartPoint3D p3, pn, ep3d;
1267 0 : CartPoint3D e[6];
1268 0 : CartPoint2D s1, s2, ep, o;
1269 0 : int rotation_to_origin, rotation_to_ep;
1270 :
1271 : // check coinciding points
1272 0 : if (FP_EQUALS(p1->x,p2->x) && FP_EQUALS(p1->y,p2->y) && FP_EQUALS(p1->z,p2->z))
1273 : return MAL_SUCCEED;
1274 : // check antipodal points
1275 0 : if (FP_EQUALS(p1->x,-p2->x) && FP_EQUALS(p1->y,-p2->y) && FP_EQUALS(p1->z,-p2->z))
1276 0 : throw(MAL, "geom.geodeticEdgeBoundingBox", SQLSTATE(38000) "Antipodal edge");
1277 :
1278 : // create the great circle plane coord system (p1, p3)
1279 : // pn = p1 x p2
1280 : // p3 = pn x p1
1281 : // TODO handle the narrow and wide angle cases
1282 0 : pn = crossProduct(p1,p2);
1283 0 : normalize3D(&pn);
1284 0 : p3 = crossProduct(&pn,p1);
1285 0 : normalize3D(&p3);
1286 :
1287 : // represent p1, p2 with (s1, s2) 2-D space
1288 : // s1.x = 1, s1.y = 0
1289 : // s2.x = p2 * p1, s2.y = p2 * p3
1290 0 : s1.x = 1;
1291 0 : s1.y = 0;
1292 0 : s2.x = dotProduct(p2, p1);
1293 0 : s2.y = dotProduct(p2, &p3);
1294 : // 2-D space origin
1295 : // O.x = 0, O.y = 0
1296 0 : o.x = 0;
1297 0 : o.y = 0;
1298 :
1299 : // create 3D endpoints E.x, E.-x, ...
1300 : // E.x = (1, 0, 0), E.-x = (-1, 0, 0) ...
1301 0 : memset(e, 0, sizeof(CartPoint2D) * 6);
1302 0 : e[0].x = e[1].y = e[2].z = 1;
1303 0 : e[3].x = e[4].y = e[5].z = -1;
1304 :
1305 : // find the rotation between s1->s2 and s1->O
1306 : // rot = norm( vec(s1,s2) x vec(s1,0))
1307 0 : rotation_to_origin = angleRotation(s1,s2,o);
1308 :
1309 : // for every endpoint E
1310 0 : for (int i = 0; i < 6; i++) {
1311 : // project the endpoint in the 2-D space
1312 0 : ep.x = dotProduct(&e[i],p1);
1313 0 : ep.y = dotProduct(&e[i],&p3);
1314 : // re-normalize it e.g. EP (for endpoint_projection)
1315 0 : normalize2D(&ep);
1316 : // ep_rot = norm( vec(s1,s2) x vec(s1,EP_end) )
1317 0 : rotation_to_ep = angleRotation(s1,s2,ep);
1318 0 : if (rotation_to_origin != rotation_to_ep) {
1319 : // convert the 2-D EP into 3-D space
1320 0 : ep3d.x = ep.x * p1->x + ep.y * p3.x;
1321 0 : ep3d.y = ep.x * p1->y + ep.y * p3.y;
1322 0 : ep3d.z = ep.x * p1->z + ep.y * p3.z;
1323 : // expand the mbox in order to include 3-D representation of EP
1324 0 : boundingBoxAddPoint(mbox,ep3d);
1325 : }
1326 : }
1327 : return MAL_SUCCEED;
1328 : }
|
