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