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 "geom_atoms.h"
14 :
15 : /***********************************************/
16 : /************* wkb type functions **************/
17 : /***********************************************/
18 :
19 : /* Creates the string representation (WKT) of a WKB */
20 : /* return length of resulting string. */
21 : ssize_t
22 982 : wkbTOSTR(char **geomWKT, size_t *len, const void *GEOMWKB, bool external)
23 : {
24 982 : const wkb *geomWKB = GEOMWKB;
25 982 : char *wkt = NULL;
26 982 : size_t dstStrLen = 5; /* "nil" */
27 :
28 : /* from WKB to GEOSGeometry */
29 982 : GEOSGeom geosGeometry = wkb2geos(geomWKB);
30 :
31 982 : if (geosGeometry) {
32 982 : size_t l;
33 982 : GEOSWKTWriter *WKT_wr = GEOSWKTWriter_create_r(geoshandle);
34 : //set the number of dimensions in the writer so that it can
35 : //read correctly the geometry coordinates
36 982 : GEOSWKTWriter_setOutputDimension_r(geoshandle, WKT_wr, GEOSGeom_getCoordinateDimension_r(geoshandle, geosGeometry));
37 982 : GEOSWKTWriter_setTrim_r(geoshandle, WKT_wr, 1);
38 982 : wkt = GEOSWKTWriter_write_r(geoshandle, WKT_wr, geosGeometry);
39 982 : if (wkt == NULL) {
40 0 : GDKerror("GEOSWKTWriter_write failed\n");
41 0 : return -1;
42 : }
43 982 : GEOSWKTWriter_destroy_r(geoshandle, WKT_wr);
44 982 : GEOSGeom_destroy_r(geoshandle, geosGeometry);
45 :
46 982 : l = strlen(wkt);
47 982 : dstStrLen = l;
48 982 : if (external)
49 831 : dstStrLen += 2; /* add quotes */
50 982 : if (*len < dstStrLen + 1 || *geomWKT == NULL) {
51 100 : *len = dstStrLen + 1;
52 100 : GDKfree(*geomWKT);
53 100 : if ((*geomWKT = GDKmalloc(*len)) == NULL) {
54 0 : GEOSFree_r(geoshandle, wkt);
55 0 : return -1;
56 : }
57 : }
58 982 : if (external)
59 831 : snprintf(*geomWKT, *len, "\"%s\"", wkt);
60 : else
61 151 : strcpy(*geomWKT, wkt);
62 982 : GEOSFree_r(geoshandle, wkt);
63 :
64 982 : return (ssize_t) dstStrLen;
65 : }
66 :
67 : /* geosGeometry == NULL */
68 0 : if (*len < 4 || *geomWKT == NULL) {
69 0 : GDKfree(*geomWKT);
70 0 : if ((*geomWKT = GDKmalloc(*len = 4)) == NULL)
71 : return -1;
72 : }
73 0 : if (external) {
74 0 : strcpy(*geomWKT, "nil");
75 0 : return 3;
76 : }
77 0 : strcpy(*geomWKT, str_nil);
78 0 : return 1;
79 : }
80 :
81 : ssize_t
82 2 : wkbFROMSTR(const char *geomWKT, size_t *len, void **GEOMWKB, bool external)
83 : {
84 2 : wkb **geomWKB = (wkb **) GEOMWKB;
85 2 : size_t parsedBytes;
86 2 : str err;
87 :
88 2 : if (external && strncmp(geomWKT, "nil", 3) == 0) {
89 0 : *geomWKB = wkbNULLcopy();
90 0 : if (*geomWKB == NULL)
91 : return -1;
92 : return 3;
93 : }
94 2 : err = wkbFROMSTR_withSRID(geomWKT, len, geomWKB, 0, &parsedBytes);
95 2 : if (err != MAL_SUCCEED) {
96 1 : GDKerror("%s", getExceptionMessageAndState(err));
97 1 : freeException(err);
98 1 : return -1;
99 : }
100 1 : return (ssize_t) parsedBytes;
101 : }
102 :
103 : BUN
104 9001 : wkbHASH(const void *W)
105 : {
106 9001 : const wkb *w = W;
107 9001 : int i;
108 9001 : BUN h = 0;
109 :
110 158111 : for (i = 0; i < (w->len - 1); i += 2) {
111 149110 : BUN a = ((unsigned char *) w->data)[i];
112 149110 : BUN b = ((unsigned char *) w->data)[i + 1];
113 : #if '\377' < 0 /* char is signed? */
114 : /* maybe sign extend */
115 149110 : if (a & 0x80)
116 27436 : a |= ~(BUN)0x7f;
117 149110 : if (b & 0x80)
118 27188 : b |= ~(BUN)0x7f;
119 : #endif
120 149110 : h = (h << 3) ^ (h >> 11) ^ (h >> 17) ^ (b << 8) ^ a;
121 : }
122 9001 : return h;
123 : }
124 :
125 : /* returns a pointer to a null wkb */
126 : const void *
127 340 : wkbNULL(void)
128 : {
129 340 : return &wkb_nil;
130 : }
131 :
132 : int
133 8494 : wkbCOMP(const void *L, const void *R)
134 : {
135 8494 : const wkb *l = L, *r = R;
136 :
137 8494 : if (l->srid != r->srid)
138 : return -1;
139 :
140 8082 : int len = l->len;
141 :
142 8082 : if (len != r->len)
143 5145 : return len - r->len;
144 :
145 2937 : if (len == ~(int) 0)
146 : return (0);
147 :
148 2601 : return memcmp(l->data, r->data, len);
149 : }
150 :
151 : /* read wkb from log */
152 : void *
153 469 : wkbREAD(void *A, size_t *dstlen, stream *s, size_t cnt)
154 : {
155 469 : wkb *a = A;
156 469 : int len;
157 469 : int srid;
158 :
159 469 : (void) cnt;
160 469 : assert(cnt == 1);
161 469 : if (mnstr_readInt(s, &len) != 1)
162 : return NULL;
163 469 : if (mnstr_readInt(s, &srid) != 1)
164 : return NULL;
165 469 : size_t wkblen = (size_t) wkb_size(len);
166 469 : if (a == NULL || *dstlen < wkblen) {
167 0 : if ((a = GDKrealloc(a, wkblen)) == NULL)
168 : return NULL;
169 0 : *dstlen = wkblen;
170 : }
171 469 : a->len = len;
172 469 : a->srid = srid;
173 469 : if (len > 0 && mnstr_read(s, (char *) a->data, len, 1) != 1) {
174 0 : GDKfree(a);
175 0 : return NULL;
176 : }
177 : return a;
178 : }
179 :
180 : /* write wkb to log */
181 : gdk_return
182 473 : wkbWRITE(const void *A, stream *s, size_t cnt)
183 : {
184 473 : const wkb *a = A;
185 473 : int len = a->len;
186 473 : int srid = a->srid;
187 :
188 473 : (void) cnt;
189 473 : assert(cnt == 1);
190 473 : if (!mnstr_writeInt(s, len)) /* 64bit: check for overflow */
191 : return GDK_FAIL;
192 473 : if (!mnstr_writeInt(s, srid)) /* 64bit: check for overflow */
193 : return GDK_FAIL;
194 943 : if (len > 0 && /* 64bit: check for overflow */
195 470 : mnstr_write(s, (char *) a->data, len, 1) < 0)
196 : return GDK_FAIL;
197 : return GDK_SUCCEED;
198 : }
199 :
200 : var_t
201 1011 : wkbPUT(BAT *b, var_t *bun, const void *VAL)
202 : {
203 1011 : const wkb *val = VAL;
204 1011 : char *base;
205 :
206 1011 : *bun = HEAP_malloc(b, wkb_size(val->len));
207 1009 : base = b->tvheap->base;
208 1009 : if (*bun != (var_t) -1) {
209 1010 : memcpy(&base[*bun], val, wkb_size(val->len));
210 1010 : b->tvheap->dirty = true;
211 : }
212 1009 : return *bun;
213 : }
214 :
215 : void
216 0 : wkbDEL(Heap *h, var_t *index)
217 : {
218 0 : HEAP_free(h, *index);
219 0 : }
220 :
221 : size_t
222 1370 : wkbLENGTH(const void *P)
223 : {
224 1370 : const wkb *p = P;
225 1370 : var_t len = wkb_size(p->len);
226 1370 : assert(len <= GDK_int_max);
227 1370 : return (size_t) len;
228 : }
229 :
230 : gdk_return
231 1074 : wkbHEAP(Heap *heap, size_t capacity)
232 : {
233 1074 : return HEAP_initialize(heap, capacity, 0, (int) sizeof(var_t));
234 : }
235 :
236 : /* Non-atom WKB functions */
237 : wkb *
238 3 : wkbNULLcopy(void)
239 : {
240 3 : wkb *n = GDKmalloc(sizeof(wkb_nil));
241 3 : if (n)
242 3 : *n = wkb_nil;
243 3 : return n;
244 : }
245 :
246 : wkb *
247 1 : wkbCopy(const wkb* src)
248 : {
249 1 : wkb *n = GDKmalloc(wkb_size(src->len));
250 1 : if (n) {
251 1 : n->len = src->len;
252 1 : n->srid = src->srid;
253 1 : memcpy(n->data, src->data, src->len);
254 : }
255 1 : return n;
256 : }
257 :
258 : /* returns the size of variable-sized atom wkb */
259 : var_t
260 6233 : wkb_size(size_t len)
261 : {
262 6233 : if (len == ~(size_t) 0)
263 180 : len = 0;
264 6233 : assert(offsetof(wkb, data) + len <= VAR_MAX);
265 6233 : return (var_t) (offsetof(wkb, data) + len);
266 : }
267 :
268 : /* Creates WKB representation (including srid) from WKT representation */
269 : /* return number of parsed characters. */
270 : str
271 1027 : wkbFROMSTR_withSRID(const char *geomWKT, size_t *len, wkb **geomWKB, int srid, size_t *nread)
272 : {
273 1027 : GEOSGeom geosGeometry = NULL; /* The geometry object that is parsed from the src string. */
274 1027 : GEOSWKTReader *WKT_reader;
275 1027 : static const char polyhedralSurface[] = "POLYHEDRALSURFACE";
276 1027 : static const char multiPolygon[] = "MULTIPOLYGON";
277 1027 : char *geomWKT_new = NULL;
278 1027 : size_t parsedCharacters = 0;
279 :
280 1027 : *nread = 0;
281 :
282 : /* we always allocate new memory */
283 1027 : GDKfree(*geomWKB);
284 1026 : *len = 0;
285 1026 : *geomWKB = NULL;
286 :
287 1026 : if (strNil(geomWKT)) {
288 0 : *geomWKB = wkbNULLcopy();
289 0 : if (*geomWKB == NULL)
290 0 : throw(MAL, "wkb.FromText", SQLSTATE(HY013) MAL_MALLOC_FAIL);
291 0 : *len = sizeof(wkb_nil);
292 0 : return MAL_SUCCEED;
293 : }
294 : //check whether the representation is binary (hex)
295 1026 : if (geomWKT[0] == '0') {
296 289 : str ret = wkbFromBinary(geomWKB, &geomWKT);
297 :
298 289 : if (ret != MAL_SUCCEED)
299 : return ret;
300 289 : *nread = strlen(geomWKT);
301 289 : *len = (size_t) wkb_size((*geomWKB)->len);
302 289 : return MAL_SUCCEED;
303 : }
304 : //check whether the geometry type is polyhedral surface
305 : //geos cannot handle this type of geometry but since it is
306 : //a special type of multipolygon I just change the type before
307 : //continuing. Of course this means that isValid for example does
308 : //not work correctly.
309 737 : if (strncasecmp(geomWKT, polyhedralSurface, strlen(polyhedralSurface)) == 0) {
310 2 : size_t sizeOfInfo = strlen(geomWKT) - strlen(polyhedralSurface) + strlen(multiPolygon) + 1;
311 2 : geomWKT_new = GDKmalloc(sizeOfInfo);
312 2 : if (geomWKT_new == NULL)
313 0 : throw(MAL, "wkb.FromText", SQLSTATE(HY013) MAL_MALLOC_FAIL);
314 2 : snprintf(geomWKT_new, sizeOfInfo, "%s%s", multiPolygon, geomWKT + strlen(polyhedralSurface));
315 2 : geomWKT = geomWKT_new;
316 : }
317 : ////////////////////////// UP TO HERE ///////////////////////////
318 :
319 737 : WKT_reader = GEOSWKTReader_create_r(geoshandle);
320 744 : if (WKT_reader == NULL) {
321 0 : if (geomWKT_new)
322 0 : GDKfree(geomWKT_new);
323 0 : throw(MAL, "wkb.FromText", SQLSTATE(38000) "Geos operation GEOSWKTReader_create failed");
324 : }
325 744 : geosGeometry = GEOSWKTReader_read_r(geoshandle, WKT_reader, geomWKT);
326 746 : GEOSWKTReader_destroy_r(geoshandle, WKT_reader);
327 :
328 740 : if (geosGeometry == NULL) {
329 55 : if (geomWKT_new)
330 0 : GDKfree(geomWKT_new);
331 55 : throw(MAL, "wkb.FromText", SQLSTATE(38000) "Geos operation GEOSWKTReader_read failed");
332 : }
333 :
334 685 : if (GEOSGeomTypeId_r(geoshandle, geosGeometry) == -1) {
335 0 : if (geomWKT_new)
336 0 : GDKfree(geomWKT_new);
337 0 : GEOSGeom_destroy_r(geoshandle, geosGeometry);
338 0 : throw(MAL, "wkb.FromText", SQLSTATE(38000) "Geos operation GEOSGeomTypeId failed");
339 : }
340 :
341 687 : GEOSSetSRID_r(geoshandle, geosGeometry, srid);
342 : /* the srid was lost with the transformation of the GEOSGeom to wkb
343 : * so we decided to store it in the wkb */
344 :
345 : /* we have a GEOSGeometry with number of coordinates and SRID and we
346 : * want to get the wkb out of it */
347 688 : *geomWKB = geos2wkb(geosGeometry);
348 684 : GEOSGeom_destroy_r(geoshandle, geosGeometry);
349 689 : if (*geomWKB == NULL) {
350 0 : if (geomWKT_new)
351 0 : GDKfree(geomWKT_new);
352 0 : throw(MAL, "wkb.FromText", SQLSTATE(38000) "Geos operation geos2wkb failed");
353 : }
354 :
355 689 : *len = (size_t) wkb_size((*geomWKB)->len);
356 684 : parsedCharacters = strlen(geomWKT);
357 684 : assert(parsedCharacters <= GDK_int_max);
358 :
359 684 : GDKfree(geomWKT_new);
360 :
361 681 : *nread = parsedCharacters;
362 681 : return MAL_SUCCEED;
363 : }
364 :
365 : /***********************************************/
366 : /************* mbr type functions **************/
367 : /***********************************************/
368 :
369 : #define MBR_WKTLEN 256
370 :
371 : /* TOSTR: print atom in a string. */
372 : /* return length of resulting string. */
373 : ssize_t
374 91 : mbrTOSTR(char **dst, size_t *len, const void *ATOM, bool external)
375 : {
376 91 : const mbr *atom = ATOM;
377 91 : char tempWkt[MBR_WKTLEN];
378 91 : size_t dstStrLen;
379 :
380 91 : if (!is_mbr_nil(atom)) {
381 91 : dstStrLen = (size_t) snprintf(tempWkt, MBR_WKTLEN,
382 : "BOX (%f %f, %f %f)",
383 91 : atom->xmin, atom->ymin,
384 91 : atom->xmax, atom->ymax);
385 : } else {
386 0 : tempWkt[0] = 0; /* not used */
387 0 : dstStrLen = 0;
388 : }
389 :
390 91 : if (*len < dstStrLen + 4 || *dst == NULL) {
391 0 : GDKfree(*dst);
392 0 : if ((*dst = GDKmalloc(*len = dstStrLen + 4)) == NULL)
393 : return -1;
394 : }
395 :
396 91 : if (dstStrLen > 4) {
397 91 : if (external) {
398 75 : snprintf(*dst, *len, "\"%s\"", tempWkt);
399 75 : dstStrLen += 2;
400 : } else {
401 16 : strcpy(*dst, tempWkt);
402 : }
403 0 : } else if (external) {
404 0 : strcpy(*dst, "nil");
405 0 : dstStrLen = 3;
406 : } else {
407 0 : strcpy(*dst, str_nil);
408 0 : dstStrLen = 1;
409 : }
410 91 : return (ssize_t) dstStrLen;
411 : }
412 :
413 : /* FROMSTR: parse string to mbr. */
414 : /* return number of parsed characters. */
415 : ssize_t
416 2 : mbrFROMSTR(const char *src, size_t *len, void **ATOM, bool external)
417 : {
418 2 : mbr **atom = (mbr **) ATOM;
419 2 : size_t nchars = 0; /* The number of characters parsed; the return value. */
420 2 : GEOSGeom geosMbr = NULL; /* The geometry object that is parsed from the src string. */
421 2 : double xmin = 0, ymin = 0, xmax = 0, ymax = 0;
422 2 : const char *c;
423 :
424 2 : if (*len < sizeof(mbr) || *atom == NULL) {
425 2 : GDKfree(*atom);
426 2 : if ((*atom = GDKmalloc(*len = sizeof(mbr))) == NULL)
427 : return -1;
428 : }
429 2 : if (external && strncmp(src, "nil", 3) == 0) {
430 0 : **atom = mbrNIL;
431 0 : return 3;
432 : }
433 2 : if (strNil(src)) {
434 0 : **atom = mbrNIL;
435 0 : return 1;
436 : }
437 :
438 2 : if ((strstr(src, "mbr") == src || strstr(src, "MBR") == src
439 0 : || strstr(src, "box") == src || strstr(src, "BOX") == src)
440 2 : && (c = strstr(src, "(")) != NULL) {
441 : /* Parse the mbr */
442 2 : if ((c - src) != 3 && (c - src) != 4) {
443 0 : GDKerror("ParseException: Expected a string like 'MBR(0 0,1 1)' or 'MBR (0 0,1 1)'\n");
444 0 : return -1;
445 : }
446 :
447 2 : if (sscanf(c, "(%lf %lf,%lf %lf)", &xmin, &ymin, &xmax, &ymax) != 4) {
448 0 : GDKerror("ParseException: Not enough coordinates.\n");
449 0 : return -1;
450 : }
451 0 : } else if ((geosMbr = GEOSGeomFromWKT(src)) == NULL) {
452 0 : GDKerror("GEOSGeomFromWKT failed\n");
453 0 : return -1;
454 : }
455 :
456 2 : if (geosMbr == NULL) {
457 2 : assert(GDK_flt_min <= xmin && xmin <= GDK_flt_max);
458 2 : assert(GDK_flt_min <= xmax && xmax <= GDK_flt_max);
459 2 : assert(GDK_flt_min <= ymin && ymin <= GDK_flt_max);
460 2 : assert(GDK_flt_min <= ymax && ymax <= GDK_flt_max);
461 2 : (*atom)->xmin = (float) xmin;
462 2 : (*atom)->ymin = (float) ymin;
463 2 : (*atom)->xmax = (float) xmax;
464 2 : (*atom)->ymax = (float) ymax;
465 2 : nchars = strlen(src);
466 : }
467 2 : if (geosMbr)
468 0 : GEOSGeom_destroy_r(geoshandle, geosMbr);
469 2 : assert(nchars <= GDK_int_max);
470 2 : return (ssize_t) nchars;
471 : }
472 :
473 : /* HASH: compute a hash value. */
474 : /* returns a positive integer hash value */
475 : BUN
476 10 : mbrHASH(const void *ATOM)
477 : {
478 10 : const mbr *atom = ATOM;
479 10 : return ATOMhash(TYPE_flt, &atom->xmin) ^ ATOMhash(TYPE_flt, &atom->ymin) ^
480 10 : ATOMhash(TYPE_flt, &atom->xmax) ^ ATOMhash(TYPE_flt, &atom->ymax);
481 : }
482 :
483 : const void *
484 340 : mbrNULL(void)
485 : {
486 340 : return &mbrNIL;
487 : }
488 :
489 : /* COMP: compare two mbrs. */
490 : /* returns int <0 if l<r, 0 if l==r, >0 else */
491 : int
492 1716 : mbrCOMP(const void *L, const void *R)
493 : {
494 : /* simple lexicographical ordering on (x,y) */
495 1716 : const mbr *l = L, *r = R;
496 1716 : int res;
497 1716 : if (is_mbr_nil(l))
498 156 : return -!is_mbr_nil(r);
499 1562 : if (is_mbr_nil(r))
500 : return 1;
501 958 : if (l->xmin == r->xmin)
502 512 : res = (l->ymin < r->ymin) ? -1 : (l->ymin != r->ymin);
503 : else
504 446 : res = (l->xmin < r->xmin) ? -1 : 1;
505 357 : if (res == 0) {
506 346 : if (l->xmax == r->xmax)
507 252 : res = (l->ymax < r->ymax) ? -1 : (l->ymax != r->ymax);
508 : else
509 94 : res = (l->xmax < r->xmax) ? -1 : 1;
510 : }
511 : return res;
512 : }
513 :
514 : /* read mbr from log */
515 : void *
516 10 : mbrREAD(void *A, size_t *dstlen, stream *s, size_t cnt)
517 : {
518 10 : mbr *a = A;
519 10 : mbr *c;
520 10 : size_t i;
521 10 : int v[4];
522 10 : flt vals[4];
523 :
524 10 : if (a == NULL || *dstlen < cnt * sizeof(mbr)) {
525 0 : if ((a = GDKrealloc(a, cnt * sizeof(mbr))) == NULL)
526 : return NULL;
527 0 : *dstlen = cnt * sizeof(mbr);
528 : }
529 24 : for (i = 0, c = a; i < cnt; i++, c++) {
530 14 : if (!mnstr_readIntArray(s, v, 4)) {
531 0 : if (a != A)
532 0 : GDKfree(a);
533 0 : return NULL;
534 : }
535 14 : memcpy(vals, v, 4 * sizeof(int));
536 14 : c->xmin = vals[0];
537 14 : c->ymin = vals[1];
538 14 : c->xmax = vals[2];
539 14 : c->ymax = vals[3];
540 : }
541 : return a;
542 : }
543 :
544 : /* write mbr to log */
545 : gdk_return
546 14 : mbrWRITE(const void *C, stream *s, size_t cnt)
547 : {
548 14 : const mbr *c = C;
549 14 : size_t i;
550 14 : flt vals[4];
551 14 : int v[4];
552 :
553 28 : for (i = 0; i < cnt; i++, c++) {
554 14 : vals[0] = c->xmin;
555 14 : vals[1] = c->ymin;
556 14 : vals[2] = c->xmax;
557 14 : vals[3] = c->ymax;
558 14 : memcpy(v, vals, 4 * sizeof(int));
559 14 : if (!mnstr_writeIntArray(s, v, 4))
560 : return GDK_FAIL;
561 : }
562 : return GDK_SUCCEED;
563 : }
564 :
565 : /* Non-atom mbr functions */
566 : /* Check if fixed-sized atom mbr is null */
567 : bool
568 5087 : is_mbr_nil(const mbr *m)
569 : {
570 5087 : return (m == NULL || is_flt_nil(m->xmin) || is_flt_nil(m->ymin) || is_flt_nil(m->xmax) || is_flt_nil(m->ymax));
571 : }
572 :
573 : /* MBR FUNCTIONS */
574 : /* MBR */
575 :
576 : /* Creates the mbr for the given geom_geometry. */
577 : str
578 675 : wkbMBR(mbr **geomMBR, wkb **geomWKB)
579 : {
580 675 : GEOSGeom geosGeometry;
581 675 : str ret = MAL_SUCCEED;
582 675 : bit empty;
583 :
584 : //check if the geometry is nil
585 675 : if (is_wkb_nil(*geomWKB)) {
586 2 : if ((*geomMBR = GDKmalloc(sizeof(mbr))) == NULL)
587 0 : throw(MAL, "geom.MBR", SQLSTATE(HY013) MAL_MALLOC_FAIL);
588 2 : **geomMBR = mbrNIL;
589 2 : return MAL_SUCCEED;
590 : }
591 : //check if the geometry is empty
592 670 : if ((ret = wkbIsEmpty(&empty, geomWKB)) != MAL_SUCCEED) {
593 : return ret;
594 : }
595 675 : if (empty) {
596 5 : if ((*geomMBR = GDKmalloc(sizeof(mbr))) == NULL)
597 0 : throw(MAL, "geom.MBR", SQLSTATE(HY013) MAL_MALLOC_FAIL);
598 5 : **geomMBR = mbrNIL;
599 5 : return MAL_SUCCEED;
600 : }
601 :
602 670 : geosGeometry = wkb2geos(*geomWKB);
603 675 : if (geosGeometry == NULL) {
604 0 : *geomMBR = NULL;
605 0 : throw(MAL, "geom.MBR", SQLSTATE(38000) "Geos problem converting GEOS to WKB");
606 : }
607 :
608 675 : *geomMBR = mbrFromGeos(geosGeometry);
609 :
610 675 : GEOSGeom_destroy_r(geoshandle, geosGeometry);
611 :
612 677 : if (*geomMBR == NULL || is_mbr_nil(*geomMBR)) {
613 0 : GDKfree(*geomMBR);
614 0 : *geomMBR = NULL;
615 0 : throw(MAL, "wkb.mbr", SQLSTATE(38000) "Geos failed to create mbr");
616 : }
617 :
618 : return MAL_SUCCEED;
619 : }
620 :
621 : str
622 23 : wkbBox2D(mbr **box, wkb **point1, wkb **point2)
623 : {
624 23 : GEOSGeom point1_geom, point2_geom;
625 23 : double xmin = 0.0, ymin = 0.0, xmax = 0.0, ymax = 0.0;
626 23 : str err = MAL_SUCCEED;
627 :
628 : //check null input
629 23 : if (is_wkb_nil(*point1) || is_wkb_nil(*point2)) {
630 1 : if ((*box = GDKmalloc(sizeof(mbr))) == NULL)
631 0 : throw(MAL, "geom.MakeBox2D", SQLSTATE(HY013) MAL_MALLOC_FAIL);
632 1 : **box = mbrNIL;
633 1 : return MAL_SUCCEED;
634 : }
635 : //check input not point geometries
636 22 : point1_geom = wkb2geos(*point1);
637 21 : point2_geom = wkb2geos(*point2);
638 22 : if (point1_geom == NULL || point2_geom == NULL) {
639 0 : if (point1_geom)
640 0 : GEOSGeom_destroy_r(geoshandle, point1_geom);
641 0 : if (point2_geom)
642 0 : GEOSGeom_destroy_r(geoshandle, point2_geom);
643 0 : throw(MAL, "geom.MakeBox2D", SQLSTATE(HY013) MAL_MALLOC_FAIL);
644 : }
645 43 : if (GEOSGeomTypeId_r(geoshandle, point1_geom) + 1 != wkbPoint_mdb ||
646 21 : GEOSGeomTypeId_r(geoshandle, point2_geom) + 1 != wkbPoint_mdb) {
647 1 : err = createException(MAL, "geom.MakeBox2D", SQLSTATE(38000) "Geometries should be points");
648 40 : } else if (GEOSGeomGetX_r(geoshandle, point1_geom, &xmin) == -1 ||
649 42 : GEOSGeomGetY_r(geoshandle, point1_geom, &ymin) == -1 ||
650 42 : GEOSGeomGetX_r(geoshandle, point2_geom, &xmax) == -1 ||
651 21 : GEOSGeomGetY_r(geoshandle, point2_geom, &ymax) == -1) {
652 :
653 0 : err = createException(MAL, "geom.MakeBox2D", SQLSTATE(38000) "Geos error in reading the points' coordinates");
654 : } else {
655 : //Assign the coordinates. Ensure that they are in correct order
656 20 : *box = GDKmalloc(sizeof(mbr));
657 21 : if (*box == NULL) {
658 0 : err = createException(MAL, "geom.MakeBox2D", SQLSTATE(HY013) MAL_MALLOC_FAIL);
659 : } else {
660 21 : (*box)->xmin = (float) (xmin < xmax ? xmin : xmax);
661 21 : (*box)->ymin = (float) (ymin < ymax ? ymin : ymax);
662 21 : (*box)->xmax = (float) (xmax > xmin ? xmax : xmin);
663 31 : (*box)->ymax = (float) (ymax > ymin ? ymax : ymin);
664 : }
665 : }
666 22 : GEOSGeom_destroy_r(geoshandle, point1_geom);
667 22 : GEOSGeom_destroy_r(geoshandle, point2_geom);
668 :
669 20 : return err;
670 : }
671 :
672 : static str
673 158 : mbrrelation_wkb(bit *out, wkb **geom1WKB, wkb **geom2WKB, str (*func)(bit *, mbr **, mbr **))
674 : {
675 158 : mbr *geom1MBR = NULL, *geom2MBR = NULL;
676 158 : str ret = MAL_SUCCEED;
677 :
678 158 : if (is_wkb_nil(*geom1WKB) || is_wkb_nil(*geom2WKB)) {
679 29 : *out = bit_nil;
680 29 : return MAL_SUCCEED;
681 : }
682 :
683 133 : ret = wkbMBR(&geom1MBR, geom1WKB);
684 132 : if (ret != MAL_SUCCEED) {
685 : return ret;
686 : }
687 :
688 132 : ret = wkbMBR(&geom2MBR, geom2WKB);
689 133 : if (ret != MAL_SUCCEED) {
690 0 : GDKfree(geom1MBR);
691 0 : return ret;
692 : }
693 :
694 133 : ret = (*func) (out, &geom1MBR, &geom2MBR);
695 :
696 132 : GDKfree(geom1MBR);
697 132 : GDKfree(geom2MBR);
698 :
699 132 : return ret;
700 : }
701 :
702 : /*returns true if the two mbrs overlap */
703 : str
704 14 : mbrOverlaps(bit *out, mbr **b1, mbr **b2)
705 : {
706 14 : if (is_mbr_nil(*b1) || is_mbr_nil(*b2))
707 0 : *out = bit_nil;
708 : else //they cannot overlap if b2 is left, right, above or below b1
709 28 : *out = !((*b2)->ymax < (*b1)->ymin ||
710 12 : (*b2)->ymin > (*b1)->ymax ||
711 12 : (*b2)->xmax < (*b1)->xmin ||
712 12 : (*b2)->xmin > (*b1)->xmax);
713 14 : return MAL_SUCCEED;
714 : }
715 :
716 : /*returns true if the mbrs of the two geometries overlap */
717 : str
718 6 : mbrOverlaps_wkb(bit *out, wkb **geom1WKB, wkb **geom2WKB)
719 : {
720 6 : return mbrrelation_wkb(out, geom1WKB, geom2WKB, mbrOverlaps);
721 : }
722 :
723 : /* returns true if b1 is above b2 */
724 : str
725 5 : mbrAbove(bit *out, mbr **b1, mbr **b2)
726 : {
727 5 : if (is_mbr_nil(*b1) || is_mbr_nil(*b2))
728 0 : *out = bit_nil;
729 : else
730 5 : *out = ((*b1)->ymin > (*b2)->ymax);
731 5 : return MAL_SUCCEED;
732 : }
733 :
734 : /*returns true if the mbrs of geom1 is above the mbr of geom2 */
735 : str
736 6 : mbrAbove_wkb(bit *out, wkb **geom1WKB, wkb **geom2WKB)
737 : {
738 6 : return mbrrelation_wkb(out, geom1WKB, geom2WKB, mbrAbove);
739 : }
740 :
741 : /* returns true if b1 is below b2 */
742 : str
743 5 : mbrBelow(bit *out, mbr **b1, mbr **b2)
744 : {
745 5 : if (is_mbr_nil(*b1) || is_mbr_nil(*b2))
746 0 : *out = bit_nil;
747 : else
748 5 : *out = ((*b1)->ymax < (*b2)->ymin);
749 5 : return MAL_SUCCEED;
750 : }
751 :
752 : /*returns true if the mbrs of geom1 is below the mbr of geom2 */
753 : str
754 6 : mbrBelow_wkb(bit *out, wkb **geom1WKB, wkb **geom2WKB)
755 : {
756 6 : return mbrrelation_wkb(out, geom1WKB, geom2WKB, mbrBelow);
757 : }
758 :
759 : /* returns true if box1 is left of box2 */
760 : str
761 65 : mbrLeft(bit *out, mbr **b1, mbr **b2)
762 : {
763 65 : if (is_mbr_nil(*b1) || is_mbr_nil(*b2))
764 0 : *out = bit_nil;
765 : else
766 65 : *out = ((*b1)->xmax < (*b2)->xmin);
767 65 : return MAL_SUCCEED;
768 : }
769 :
770 : /*returns true if the mbrs of geom1 is on the left of the mbr of geom2 */
771 : str
772 82 : mbrLeft_wkb(bit *out, wkb **geom1WKB, wkb **geom2WKB)
773 : {
774 82 : return mbrrelation_wkb(out, geom1WKB, geom2WKB, mbrLeft);
775 : }
776 :
777 : /* returns true if box1 is right of box2 */
778 : str
779 10 : mbrRight(bit *out, mbr **b1, mbr **b2)
780 : {
781 10 : if (is_mbr_nil(*b1) || is_mbr_nil(*b2))
782 0 : *out = bit_nil;
783 : else
784 11 : *out = ((*b1)->xmin > (*b2)->xmax);
785 11 : return MAL_SUCCEED;
786 : }
787 :
788 : /*returns true if the mbrs of geom1 is on the right of the mbr of geom2 */
789 : str
790 10 : mbrRight_wkb(bit *out, wkb **geom1WKB, wkb **geom2WKB)
791 : {
792 10 : return mbrrelation_wkb(out, geom1WKB, geom2WKB, mbrRight);
793 : }
794 :
795 : /* returns true if box1 overlaps or is above box2 when only the Y coordinate is considered*/
796 : str
797 0 : mbrOverlapOrAbove(bit *out, mbr **b1, mbr **b2)
798 : {
799 0 : if (is_mbr_nil(*b1) || is_mbr_nil(*b2))
800 0 : *out = bit_nil;
801 : else
802 0 : *out = ((*b1)->ymin >= (*b2)->ymin);
803 0 : return MAL_SUCCEED;
804 : }
805 :
806 : /*returns true if the mbrs of geom1 overlaps or is above the mbr of geom2 */
807 : str
808 0 : mbrOverlapOrAbove_wkb(bit *out, wkb **geom1WKB, wkb **geom2WKB)
809 : {
810 0 : return mbrrelation_wkb(out, geom1WKB, geom2WKB, mbrOverlapOrAbove);
811 : }
812 :
813 : /* returns true if box1 overlaps or is below box2 when only the Y coordinate is considered*/
814 : str
815 0 : mbrOverlapOrBelow(bit *out, mbr **b1, mbr **b2)
816 : {
817 0 : if (is_mbr_nil(*b1) || is_mbr_nil(*b2))
818 0 : *out = bit_nil;
819 : else
820 0 : *out = ((*b1)->ymax <= (*b2)->ymax);
821 0 : return MAL_SUCCEED;
822 : }
823 :
824 : /*returns true if the mbrs of geom1 overlaps or is below the mbr of geom2 */
825 : str
826 0 : mbrOverlapOrBelow_wkb(bit *out, wkb **geom1WKB, wkb **geom2WKB)
827 : {
828 0 : return mbrrelation_wkb(out, geom1WKB, geom2WKB, mbrOverlapOrBelow);
829 : }
830 :
831 : /* returns true if box1 overlaps or is left of box2 when only the X coordinate is considered*/
832 : str
833 10 : mbrOverlapOrLeft(bit *out, mbr **b1, mbr **b2)
834 : {
835 10 : if (is_mbr_nil(*b1) || is_mbr_nil(*b2))
836 0 : *out = bit_nil;
837 : else
838 10 : *out = ((*b1)->xmax <= (*b2)->xmax);
839 10 : return MAL_SUCCEED;
840 : }
841 :
842 : /*returns true if the mbrs of geom1 overlaps or is on the left of the mbr of geom2 */
843 : str
844 6 : mbrOverlapOrLeft_wkb(bit *out, wkb **geom1WKB, wkb **geom2WKB)
845 : {
846 6 : return mbrrelation_wkb(out, geom1WKB, geom2WKB, mbrOverlapOrLeft);
847 : }
848 :
849 : /* returns true if box1 overlaps or is right of box2 when only the X coordinate is considered*/
850 : str
851 10 : mbrOverlapOrRight(bit *out, mbr **b1, mbr **b2)
852 : {
853 10 : if (is_mbr_nil(*b1) || is_mbr_nil(*b2))
854 0 : *out = bit_nil;
855 : else
856 10 : *out = ((*b1)->xmin >= (*b2)->xmin);
857 10 : return MAL_SUCCEED;
858 : }
859 :
860 : /*returns true if the mbrs of geom1 overlaps or is on the right of the mbr of geom2 */
861 : str
862 6 : mbrOverlapOrRight_wkb(bit *out, wkb **geom1WKB, wkb **geom2WKB)
863 : {
864 6 : return mbrrelation_wkb(out, geom1WKB, geom2WKB, mbrOverlapOrRight);
865 : }
866 :
867 : /* returns true if b1 is contained in b2 */
868 : str
869 46 : mbrContained(bit *out, mbr **b1, mbr **b2)
870 : {
871 46 : if (is_mbr_nil(*b1) || is_mbr_nil(*b2))
872 0 : *out = bit_nil;
873 : else
874 75 : *out = (((*b1)->xmin >= (*b2)->xmin) && ((*b1)->xmax <= (*b2)->xmax) && ((*b1)->ymin >= (*b2)->ymin) && ((*b1)->ymax <= (*b2)->ymax));
875 46 : return MAL_SUCCEED;
876 : }
877 :
878 : /*returns true if the mbrs of geom1 is contained in the mbr of geom2 */
879 : str
880 18 : mbrContained_wkb(bit *out, wkb **geom1WKB, wkb **geom2WKB)
881 : {
882 18 : return mbrrelation_wkb(out, geom1WKB, geom2WKB, mbrContained);
883 : }
884 :
885 : /*returns true if b1 contains b2 */
886 : str
887 26 : mbrContains(bit *out, mbr **b1, mbr **b2)
888 : {
889 26 : return mbrContained(out, b2, b1);
890 : }
891 :
892 : /*returns true if the mbrs of geom1 contains the mbr of geom2 */
893 : str
894 6 : mbrContains_wkb(bit *out, wkb **geom1WKB, wkb **geom2WKB)
895 : {
896 6 : return mbrrelation_wkb(out, geom1WKB, geom2WKB, mbrContains);
897 : }
898 :
899 : /* returns true if the boxes are the same */
900 : str
901 31 : mbrEqual(bit *out, mbr **b1, mbr **b2)
902 : {
903 31 : if (is_mbr_nil(*b1) && is_mbr_nil(*b2))
904 0 : *out = 1;
905 31 : else if (is_mbr_nil(*b1) || is_mbr_nil(*b2))
906 0 : *out = 0;
907 : else
908 55 : *out = (((*b1)->xmin == (*b2)->xmin) && ((*b1)->xmax == (*b2)->xmax) && ((*b1)->ymin == (*b2)->ymin) && ((*b1)->ymax == (*b2)->ymax));
909 31 : return MAL_SUCCEED;
910 : }
911 :
912 : /*returns true if the mbrs of geom1 and the mbr of geom2 are the same */
913 : str
914 12 : mbrEqual_wkb(bit *out, wkb **geom1WKB, wkb **geom2WKB)
915 : {
916 12 : return mbrrelation_wkb(out, geom1WKB, geom2WKB, mbrEqual);
917 : }
918 :
919 : str
920 20 : mbrDiagonal(dbl *out, mbr **b)
921 : {
922 20 : double side_a = .0, side_b = .0;
923 :
924 20 : if (is_mbr_nil(*b)) {
925 0 : *out = dbl_nil;
926 0 : return MAL_SUCCEED;
927 : }
928 :
929 20 : side_a = (*b)->xmax - (*b)->xmin;
930 20 : side_b = (*b)->ymax - (*b)->ymin;
931 :
932 20 : *out = sqrt(pow(side_a, 2.0) + pow(side_b, 2.0));
933 :
934 20 : return MAL_SUCCEED;
935 : }
936 :
937 : /* returns the Euclidean distance of the centroids of the boxes */
938 : str
939 132 : mbrDistance(dbl *out, mbr **b1, mbr **b2)
940 : {
941 132 : double b1_Cx = 0.0, b1_Cy = 0.0, b2_Cx = 0.0, b2_Cy = 0.0;
942 :
943 132 : if (is_mbr_nil(*b1) || is_mbr_nil(*b2)) {
944 0 : *out = dbl_nil;
945 0 : return MAL_SUCCEED;
946 : }
947 : //compute the centroids of the two polygons
948 132 : b1_Cx = ((*b1)->xmin + (*b1)->xmax) / 2.0;
949 132 : b1_Cy = ((*b1)->ymin + (*b1)->ymax) / 2.0;
950 132 : b2_Cx = ((*b2)->xmin + (*b2)->xmax) / 2.0;
951 132 : b2_Cy = ((*b2)->ymin + (*b2)->ymax) / 2.0;
952 :
953 : //compute the euclidean distance
954 132 : *out = sqrt(pow(b2_Cx - b1_Cx, 2.0) + pow(b2_Cy - b1_Cy, 2.0));
955 :
956 132 : return MAL_SUCCEED;
957 : }
958 :
959 : /*returns the Euclidean distance of the centroids of the mbrs of the two geometries */
960 : str
961 166 : mbrDistance_wkb(dbl *out, wkb **geom1WKB, wkb **geom2WKB)
962 : {
963 166 : mbr *geom1MBR = NULL, *geom2MBR = NULL;
964 166 : str ret = MAL_SUCCEED;
965 :
966 166 : if (is_wkb_nil(*geom1WKB) || is_wkb_nil(*geom2WKB)) {
967 45 : *out = dbl_nil;
968 45 : return MAL_SUCCEED;
969 : }
970 :
971 124 : ret = wkbMBR(&geom1MBR, geom1WKB);
972 124 : if (ret != MAL_SUCCEED) {
973 : return ret;
974 : }
975 :
976 124 : ret = wkbMBR(&geom2MBR, geom2WKB);
977 122 : if (ret != MAL_SUCCEED) {
978 0 : GDKfree(geom1MBR);
979 0 : return ret;
980 : }
981 :
982 122 : ret = mbrDistance(out, &geom1MBR, &geom2MBR);
983 :
984 123 : GDKfree(geom1MBR);
985 125 : GDKfree(geom2MBR);
986 :
987 125 : return ret;
988 : }
989 :
990 : /* get Xmin, Ymin, Xmax, Ymax coordinates of mbr */
991 : str
992 205 : wkbCoordinateFromMBR(dbl *coordinateValue, mbr **geomMBR, int *coordinateIdx)
993 : {
994 : //check if the MBR is null
995 205 : if (is_mbr_nil(*geomMBR) || is_int_nil(*coordinateIdx)) {
996 8 : *coordinateValue = dbl_nil;
997 8 : return MAL_SUCCEED;
998 : }
999 :
1000 208 : switch (*coordinateIdx) {
1001 52 : case 1:
1002 52 : *coordinateValue = (*geomMBR)->xmin;
1003 52 : break;
1004 52 : case 2:
1005 52 : *coordinateValue = (*geomMBR)->ymin;
1006 52 : break;
1007 52 : case 3:
1008 52 : *coordinateValue = (*geomMBR)->xmax;
1009 52 : break;
1010 52 : case 4:
1011 52 : *coordinateValue = (*geomMBR)->ymax;
1012 52 : break;
1013 0 : default:
1014 0 : throw(MAL, "geom.coordinateFromMBR", SQLSTATE(38000) "Geos unrecognized coordinateIdx: %d\n", *coordinateIdx);
1015 : }
1016 :
1017 : return MAL_SUCCEED;
1018 : }
1019 :
1020 : str
1021 0 : wkbCoordinateFromWKB(dbl *coordinateValue, wkb **geomWKB, int *coordinateIdx)
1022 : {
1023 0 : mbr *geomMBR;
1024 0 : str ret = MAL_SUCCEED;
1025 0 : bit empty;
1026 :
1027 0 : if (is_wkb_nil(*geomWKB) || is_int_nil(*coordinateIdx)) {
1028 0 : *coordinateValue = dbl_nil;
1029 0 : return MAL_SUCCEED;
1030 : }
1031 :
1032 : //check if the geometry is empty
1033 0 : if ((ret = wkbIsEmpty(&empty, geomWKB)) != MAL_SUCCEED) {
1034 : return ret;
1035 : }
1036 :
1037 0 : if (empty) {
1038 0 : *coordinateValue = dbl_nil;
1039 0 : return MAL_SUCCEED;
1040 : }
1041 :
1042 0 : if ((ret = wkbMBR(&geomMBR, geomWKB)) != MAL_SUCCEED)
1043 : return ret;
1044 :
1045 0 : ret = wkbCoordinateFromMBR(coordinateValue, &geomMBR, coordinateIdx);
1046 :
1047 0 : GDKfree(geomMBR);
1048 :
1049 0 : return ret;
1050 : }
1051 :
1052 : str
1053 2 : mbrFromString(mbr **w, const char **src)
1054 : {
1055 2 : size_t len = *w ? sizeof(mbr) : 0;
1056 2 : char *errbuf;
1057 2 : str ex;
1058 :
1059 2 : if (mbrFROMSTR(*src, &len, (void **) w, false) >= 0)
1060 : return MAL_SUCCEED;
1061 0 : GDKfree(*w);
1062 0 : *w = NULL;
1063 0 : errbuf = GDKerrbuf;
1064 0 : if (errbuf) {
1065 0 : if (strncmp(errbuf, "!ERROR: ", 8) == 0)
1066 0 : errbuf += 8;
1067 : } else {
1068 : errbuf = "cannot parse string";
1069 : }
1070 :
1071 0 : ex = createException(MAL, "mbr.FromString", SQLSTATE(38000) "Geos %s", errbuf);
1072 :
1073 0 : GDKclrerr();
1074 :
1075 0 : return ex;
1076 : }
1077 :
1078 : /* COMMAND mbr
1079 : * Creates the mbr for the given geom_geometry.
1080 : */
1081 :
1082 : str
1083 0 : ordinatesMBR(mbr **res, flt *minX, flt *minY, flt *maxX, flt *maxY)
1084 : {
1085 0 : if ((*res = GDKmalloc(sizeof(mbr))) == NULL)
1086 0 : throw(MAL, "geom.mbr", SQLSTATE(HY013) MAL_MALLOC_FAIL);
1087 0 : if (is_flt_nil(*minX) || is_flt_nil(*minY) || is_flt_nil(*maxX) || is_flt_nil(*maxY))
1088 0 : **res = mbrNIL;
1089 : else {
1090 0 : (*res)->xmin = *minX;
1091 0 : (*res)->ymin = *minY;
1092 0 : (*res)->xmax = *maxX;
1093 0 : (*res)->ymax = *maxY;
1094 : }
1095 : return MAL_SUCCEED;
1096 : }
1097 :
1098 : str
1099 0 : mbrIntersects(bit *out, mbr** mbr1, mbr** mbr2) {
1100 0 : if (((*mbr1)->ymax < (*mbr2)->ymin) || ((*mbr1)->ymin > (*mbr2)->ymax))
1101 0 : (*out) = false;
1102 0 : else if (((*mbr1)->xmax < (*mbr2)->xmin) || ((*mbr1)->xmin > (*mbr2)->xmax))
1103 0 : (*out) = false;
1104 : else
1105 0 : (*out) = true;
1106 0 : return MAL_SUCCEED;
1107 : }
|
