1 /**********************************************************************
2 *
3 * PostGIS - Spatial Types for PostgreSQL
4 * http://postgis.net
5 *
6 * PostGIS is free software: you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation, either version 2 of the License, or
9 * (at your option) any later version.
10 *
11 * PostGIS is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with PostGIS. If not, see <http://www.gnu.org/licenses/>.
18 *
19 **********************************************************************
20 *
21 * Copyright 2009 Paul Ramsey <pramsey@cleverelephant.ca>
22 *
23 **********************************************************************/
24
25 #include "liblwgeom_internal.h"
26 #include "lwgeodetic.h"
27 #include "lwgeom_log.h"
28 #include <stdlib.h>
29 #include <math.h>
30
31
gbox_new(lwflags_t flags)32 GBOX* gbox_new(lwflags_t flags)
33 {
34 GBOX *g = (GBOX*)lwalloc(sizeof(GBOX));
35 gbox_init(g);
36 g->flags = flags;
37 return g;
38 }
39
gbox_init(GBOX * gbox)40 void gbox_init(GBOX *gbox)
41 {
42 memset(gbox, 0, sizeof(GBOX));
43 }
44
gbox_clone(const GBOX * gbox)45 GBOX* gbox_clone(const GBOX *gbox)
46 {
47 GBOX *g = lwalloc(sizeof(GBOX));
48 memcpy(g, gbox, sizeof(GBOX));
49 return g;
50 }
51
52 /* TODO to be removed */
box3d_from_gbox(const GBOX * gbox)53 BOX3D* box3d_from_gbox(const GBOX *gbox)
54 {
55 BOX3D *b;
56 assert(gbox);
57
58 b = lwalloc(sizeof(BOX3D));
59
60 b->xmin = gbox->xmin;
61 b->xmax = gbox->xmax;
62 b->ymin = gbox->ymin;
63 b->ymax = gbox->ymax;
64
65 if ( FLAGS_GET_Z(gbox->flags) )
66 {
67 b->zmin = gbox->zmin;
68 b->zmax = gbox->zmax;
69 }
70 else
71 {
72 b->zmin = b->zmax = 0.0;
73 }
74
75 b->srid = SRID_UNKNOWN;
76 return b;
77 }
78
79 /* TODO to be removed */
box3d_to_gbox(const BOX3D * b3d)80 GBOX* box3d_to_gbox(const BOX3D *b3d)
81 {
82 GBOX *b;
83 assert(b3d);
84
85 b = lwalloc(sizeof(GBOX));
86
87 b->xmin = b3d->xmin;
88 b->xmax = b3d->xmax;
89 b->ymin = b3d->ymin;
90 b->ymax = b3d->ymax;
91 b->zmin = b3d->zmin;
92 b->zmax = b3d->zmax;
93
94 return b;
95 }
96
gbox_expand(GBOX * g,double d)97 void gbox_expand(GBOX *g, double d)
98 {
99 g->xmin -= d;
100 g->xmax += d;
101 g->ymin -= d;
102 g->ymax += d;
103 if (FLAGS_GET_Z(g->flags) || FLAGS_GET_GEODETIC(g->flags))
104 {
105 g->zmin -= d;
106 g->zmax += d;
107 }
108 if (FLAGS_GET_M(g->flags))
109 {
110 g->mmin -= d;
111 g->mmax += d;
112 }
113 }
114
gbox_expand_xyzm(GBOX * g,double dx,double dy,double dz,double dm)115 void gbox_expand_xyzm(GBOX *g, double dx, double dy, double dz, double dm)
116 {
117 g->xmin -= dx;
118 g->xmax += dx;
119 g->ymin -= dy;
120 g->ymax += dy;
121
122 if (FLAGS_GET_Z(g->flags))
123 {
124 g->zmin -= dz;
125 g->zmax += dz;
126 }
127
128 if (FLAGS_GET_M(g->flags))
129 {
130 g->mmin -= dm;
131 g->mmax += dm;
132 }
133 }
134
gbox_union(const GBOX * g1,const GBOX * g2,GBOX * gout)135 int gbox_union(const GBOX *g1, const GBOX *g2, GBOX *gout)
136 {
137 if ( ( ! g1 ) && ( ! g2 ) )
138 return LW_FALSE;
139 else if (!g1)
140 {
141 memcpy(gout, g2, sizeof(GBOX));
142 return LW_TRUE;
143 }
144 else if (!g2)
145 {
146 memcpy(gout, g1, sizeof(GBOX));
147 return LW_TRUE;
148 }
149
150 gout->flags = g1->flags;
151
152 gout->xmin = FP_MIN(g1->xmin, g2->xmin);
153 gout->xmax = FP_MAX(g1->xmax, g2->xmax);
154
155 gout->ymin = FP_MIN(g1->ymin, g2->ymin);
156 gout->ymax = FP_MAX(g1->ymax, g2->ymax);
157
158 gout->zmin = FP_MIN(g1->zmin, g2->zmin);
159 gout->zmax = FP_MAX(g1->zmax, g2->zmax);
160
161 return LW_TRUE;
162 }
163
gbox_same(const GBOX * g1,const GBOX * g2)164 int gbox_same(const GBOX *g1, const GBOX *g2)
165 {
166 if (FLAGS_GET_ZM(g1->flags) != FLAGS_GET_ZM(g2->flags))
167 return LW_FALSE;
168
169 if (!gbox_same_2d(g1, g2)) return LW_FALSE;
170
171 if (FLAGS_GET_Z(g1->flags) && (g1->zmin != g2->zmin || g1->zmax != g2->zmax))
172 return LW_FALSE;
173 if (FLAGS_GET_M(g1->flags) && (g1->mmin != g2->mmin || g1->mmax != g2->mmax))
174 return LW_FALSE;
175
176 return LW_TRUE;
177 }
178
gbox_same_2d(const GBOX * g1,const GBOX * g2)179 int gbox_same_2d(const GBOX *g1, const GBOX *g2)
180 {
181 if (g1->xmin == g2->xmin && g1->ymin == g2->ymin &&
182 g1->xmax == g2->xmax && g1->ymax == g2->ymax)
183 return LW_TRUE;
184 return LW_FALSE;
185 }
186
gbox_same_2d_float(const GBOX * g1,const GBOX * g2)187 int gbox_same_2d_float(const GBOX *g1, const GBOX *g2)
188 {
189 if ((g1->xmax == g2->xmax || next_float_up(g1->xmax) == next_float_up(g2->xmax)) &&
190 (g1->ymax == g2->ymax || next_float_up(g1->ymax) == next_float_up(g2->ymax)) &&
191 (g1->xmin == g2->xmin || next_float_down(g1->xmin) == next_float_down(g1->xmin)) &&
192 (g1->ymin == g2->ymin || next_float_down(g2->ymin) == next_float_down(g2->ymin)))
193 return LW_TRUE;
194 return LW_FALSE;
195 }
196
gbox_is_valid(const GBOX * gbox)197 int gbox_is_valid(const GBOX *gbox)
198 {
199 /* X */
200 if ( ! isfinite(gbox->xmin) || isnan(gbox->xmin) ||
201 ! isfinite(gbox->xmax) || isnan(gbox->xmax) )
202 return LW_FALSE;
203
204 /* Y */
205 if ( ! isfinite(gbox->ymin) || isnan(gbox->ymin) ||
206 ! isfinite(gbox->ymax) || isnan(gbox->ymax) )
207 return LW_FALSE;
208
209 /* Z */
210 if ( FLAGS_GET_GEODETIC(gbox->flags) || FLAGS_GET_Z(gbox->flags) )
211 {
212 if ( ! isfinite(gbox->zmin) || isnan(gbox->zmin) ||
213 ! isfinite(gbox->zmax) || isnan(gbox->zmax) )
214 return LW_FALSE;
215 }
216
217 /* M */
218 if ( FLAGS_GET_M(gbox->flags) )
219 {
220 if ( ! isfinite(gbox->mmin) || isnan(gbox->mmin) ||
221 ! isfinite(gbox->mmax) || isnan(gbox->mmax) )
222 return LW_FALSE;
223 }
224
225 return LW_TRUE;
226 }
227
gbox_merge_point3d(const POINT3D * p,GBOX * gbox)228 int gbox_merge_point3d(const POINT3D *p, GBOX *gbox)
229 {
230 if ( gbox->xmin > p->x ) gbox->xmin = p->x;
231 if ( gbox->ymin > p->y ) gbox->ymin = p->y;
232 if ( gbox->zmin > p->z ) gbox->zmin = p->z;
233 if ( gbox->xmax < p->x ) gbox->xmax = p->x;
234 if ( gbox->ymax < p->y ) gbox->ymax = p->y;
235 if ( gbox->zmax < p->z ) gbox->zmax = p->z;
236 return LW_SUCCESS;
237 }
238
gbox_init_point3d(const POINT3D * p,GBOX * gbox)239 int gbox_init_point3d(const POINT3D *p, GBOX *gbox)
240 {
241 gbox->xmin = gbox->xmax = p->x;
242 gbox->ymin = gbox->ymax = p->y;
243 gbox->zmin = gbox->zmax = p->z;
244 return LW_SUCCESS;
245 }
246
gbox_contains_point3d(const GBOX * gbox,const POINT3D * pt)247 int gbox_contains_point3d(const GBOX *gbox, const POINT3D *pt)
248 {
249 if ( gbox->xmin > pt->x || gbox->ymin > pt->y || gbox->zmin > pt->z ||
250 gbox->xmax < pt->x || gbox->ymax < pt->y || gbox->zmax < pt->z )
251 {
252 return LW_FALSE;
253 }
254 return LW_TRUE;
255 }
256
gbox_merge(const GBOX * new_box,GBOX * merge_box)257 int gbox_merge(const GBOX *new_box, GBOX *merge_box)
258 {
259 assert(merge_box);
260
261 if ( FLAGS_GET_ZM(merge_box->flags) != FLAGS_GET_ZM(new_box->flags) )
262 return LW_FAILURE;
263
264 if ( new_box->xmin < merge_box->xmin) merge_box->xmin = new_box->xmin;
265 if ( new_box->ymin < merge_box->ymin) merge_box->ymin = new_box->ymin;
266 if ( new_box->xmax > merge_box->xmax) merge_box->xmax = new_box->xmax;
267 if ( new_box->ymax > merge_box->ymax) merge_box->ymax = new_box->ymax;
268
269 if ( FLAGS_GET_Z(merge_box->flags) || FLAGS_GET_GEODETIC(merge_box->flags) )
270 {
271 if ( new_box->zmin < merge_box->zmin) merge_box->zmin = new_box->zmin;
272 if ( new_box->zmax > merge_box->zmax) merge_box->zmax = new_box->zmax;
273 }
274 if ( FLAGS_GET_M(merge_box->flags) )
275 {
276 if ( new_box->mmin < merge_box->mmin) merge_box->mmin = new_box->mmin;
277 if ( new_box->mmax > merge_box->mmax) merge_box->mmax = new_box->mmax;
278 }
279
280 return LW_SUCCESS;
281 }
282
gbox_overlaps(const GBOX * g1,const GBOX * g2)283 int gbox_overlaps(const GBOX *g1, const GBOX *g2)
284 {
285
286 /* Make sure our boxes are consistent */
287 if ( FLAGS_GET_GEODETIC(g1->flags) != FLAGS_GET_GEODETIC(g2->flags) )
288 lwerror("gbox_overlaps: cannot compare geodetic and non-geodetic boxes");
289
290 /* Check X/Y first */
291 if ( g1->xmax < g2->xmin || g1->ymax < g2->ymin ||
292 g1->xmin > g2->xmax || g1->ymin > g2->ymax )
293 return LW_FALSE;
294
295 /* Deal with the geodetic case special: we only compare the geodetic boxes (x/y/z) */
296 /* Never the M dimension */
297 if ( FLAGS_GET_GEODETIC(g1->flags) && FLAGS_GET_GEODETIC(g2->flags) )
298 {
299 if ( g1->zmax < g2->zmin || g1->zmin > g2->zmax )
300 return LW_FALSE;
301 else
302 return LW_TRUE;
303 }
304
305 /* If both geodetic or both have Z, check Z */
306 if ( FLAGS_GET_Z(g1->flags) && FLAGS_GET_Z(g2->flags) )
307 {
308 if ( g1->zmax < g2->zmin || g1->zmin > g2->zmax )
309 return LW_FALSE;
310 }
311
312 /* If both have M, check M */
313 if ( FLAGS_GET_M(g1->flags) && FLAGS_GET_M(g2->flags) )
314 {
315 if ( g1->mmax < g2->mmin || g1->mmin > g2->mmax )
316 return LW_FALSE;
317 }
318
319 return LW_TRUE;
320 }
321
322 int
gbox_overlaps_2d(const GBOX * g1,const GBOX * g2)323 gbox_overlaps_2d(const GBOX *g1, const GBOX *g2)
324 {
325
326 /* Make sure our boxes are consistent */
327 if ( FLAGS_GET_GEODETIC(g1->flags) != FLAGS_GET_GEODETIC(g2->flags) )
328 lwerror("gbox_overlaps: cannot compare geodetic and non-geodetic boxes");
329
330 /* Check X/Y first */
331 if ( g1->xmax < g2->xmin || g1->ymax < g2->ymin ||
332 g1->xmin > g2->xmax || g1->ymin > g2->ymax )
333 return LW_FALSE;
334
335 return LW_TRUE;
336 }
337
338 int
gbox_contains_2d(const GBOX * g1,const GBOX * g2)339 gbox_contains_2d(const GBOX *g1, const GBOX *g2)
340 {
341 if ( ( g2->xmin < g1->xmin ) || ( g2->xmax > g1->xmax ) ||
342 ( g2->ymin < g1->ymin ) || ( g2->ymax > g1->ymax ) )
343 {
344 return LW_FALSE;
345 }
346 return LW_TRUE;
347 }
348
349 int
gbox_contains_point2d(const GBOX * g,const POINT2D * p)350 gbox_contains_point2d(const GBOX *g, const POINT2D *p)
351 {
352 if ( ( g->xmin <= p->x ) && ( g->xmax >= p->x ) &&
353 ( g->ymin <= p->y ) && ( g->ymax >= p->y ) )
354 {
355 return LW_TRUE;
356 }
357 return LW_FALSE;
358 }
359
360 /**
361 * Warning, this function is only good for x/y/z boxes, used
362 * in unit testing of geodetic box generation.
363 */
gbox_from_string(const char * str)364 GBOX* gbox_from_string(const char *str)
365 {
366 const char *ptr = str;
367 char *nextptr;
368 char *gbox_start = strstr(str, "GBOX((");
369 GBOX *gbox = gbox_new(lwflags(0,0,1));
370 if ( ! gbox_start ) return NULL; /* No header found */
371 ptr += 6;
372 gbox->xmin = strtod(ptr, &nextptr);
373 if ( ptr == nextptr ) return NULL; /* No double found */
374 ptr = nextptr + 1;
375 gbox->ymin = strtod(ptr, &nextptr);
376 if ( ptr == nextptr ) return NULL; /* No double found */
377 ptr = nextptr + 1;
378 gbox->zmin = strtod(ptr, &nextptr);
379 if ( ptr == nextptr ) return NULL; /* No double found */
380 ptr = nextptr + 3;
381 gbox->xmax = strtod(ptr, &nextptr);
382 if ( ptr == nextptr ) return NULL; /* No double found */
383 ptr = nextptr + 1;
384 gbox->ymax = strtod(ptr, &nextptr);
385 if ( ptr == nextptr ) return NULL; /* No double found */
386 ptr = nextptr + 1;
387 gbox->zmax = strtod(ptr, &nextptr);
388 if ( ptr == nextptr ) return NULL; /* No double found */
389 return gbox;
390 }
391
gbox_to_string(const GBOX * gbox)392 char* gbox_to_string(const GBOX *gbox)
393 {
394 const size_t sz = 138;
395 char *str = NULL;
396
397 if ( ! gbox )
398 return lwstrdup("NULL POINTER");
399
400 str = (char*)lwalloc(sz);
401
402 if ( FLAGS_GET_GEODETIC(gbox->flags) )
403 {
404 snprintf(str, sz, "GBOX((%.8g,%.8g,%.8g),(%.8g,%.8g,%.8g))", gbox->xmin, gbox->ymin, gbox->zmin, gbox->xmax, gbox->ymax, gbox->zmax);
405 return str;
406 }
407 if ( FLAGS_GET_Z(gbox->flags) && FLAGS_GET_M(gbox->flags) )
408 {
409 snprintf(str, sz, "GBOX((%.8g,%.8g,%.8g,%.8g),(%.8g,%.8g,%.8g,%.8g))", gbox->xmin, gbox->ymin, gbox->zmin, gbox->mmin, gbox->xmax, gbox->ymax, gbox->zmax, gbox->mmax);
410 return str;
411 }
412 if ( FLAGS_GET_Z(gbox->flags) )
413 {
414 snprintf(str, sz, "GBOX((%.8g,%.8g,%.8g),(%.8g,%.8g,%.8g))", gbox->xmin, gbox->ymin, gbox->zmin, gbox->xmax, gbox->ymax, gbox->zmax);
415 return str;
416 }
417 if ( FLAGS_GET_M(gbox->flags) )
418 {
419 snprintf(str, sz, "GBOX((%.8g,%.8g,%.8g),(%.8g,%.8g,%.8g))", gbox->xmin, gbox->ymin, gbox->mmin, gbox->xmax, gbox->ymax, gbox->mmax);
420 return str;
421 }
422 snprintf(str, sz, "GBOX((%.8g,%.8g),(%.8g,%.8g))", gbox->xmin, gbox->ymin, gbox->xmax, gbox->ymax);
423 return str;
424 }
425
gbox_copy(const GBOX * box)426 GBOX* gbox_copy(const GBOX *box)
427 {
428 GBOX *copy = (GBOX*)lwalloc(sizeof(GBOX));
429 memcpy(copy, box, sizeof(GBOX));
430 return copy;
431 }
432
gbox_duplicate(const GBOX * original,GBOX * duplicate)433 void gbox_duplicate(const GBOX *original, GBOX *duplicate)
434 {
435 assert(duplicate);
436 assert(original);
437 memcpy(duplicate, original, sizeof(GBOX));
438 }
439
gbox_serialized_size(lwflags_t flags)440 size_t gbox_serialized_size(lwflags_t flags)
441 {
442 if (FLAGS_GET_GEODETIC(flags))
443 return 6 * sizeof(float);
444 else
445 return 2 * FLAGS_NDIMS(flags) * sizeof(float);
446 }
447
448
449 /* ********************************************************************************
450 ** Compute cartesian bounding GBOX boxes from LWGEOM.
451 */
452
lw_arc_calculate_gbox_cartesian_2d(const POINT2D * A1,const POINT2D * A2,const POINT2D * A3,GBOX * gbox)453 int lw_arc_calculate_gbox_cartesian_2d(const POINT2D *A1, const POINT2D *A2, const POINT2D *A3, GBOX *gbox)
454 {
455 POINT2D xmin, ymin, xmax, ymax;
456 POINT2D C;
457 int A2_side;
458 double radius_A;
459
460 LWDEBUG(2, "lw_arc_calculate_gbox_cartesian_2d called.");
461
462 radius_A = lw_arc_center(A1, A2, A3, &C);
463
464 /* Negative radius signals straight line, p1/p2/p3 are collinear */
465 if (radius_A < 0.0)
466 {
467 gbox->xmin = FP_MIN(A1->x, A3->x);
468 gbox->ymin = FP_MIN(A1->y, A3->y);
469 gbox->xmax = FP_MAX(A1->x, A3->x);
470 gbox->ymax = FP_MAX(A1->y, A3->y);
471 return LW_SUCCESS;
472 }
473
474 /* Matched start/end points imply circle */
475 if ( A1->x == A3->x && A1->y == A3->y )
476 {
477 gbox->xmin = C.x - radius_A;
478 gbox->ymin = C.y - radius_A;
479 gbox->xmax = C.x + radius_A;
480 gbox->ymax = C.y + radius_A;
481 return LW_SUCCESS;
482 }
483
484 /* First approximation, bounds of start/end points */
485 gbox->xmin = FP_MIN(A1->x, A3->x);
486 gbox->ymin = FP_MIN(A1->y, A3->y);
487 gbox->xmax = FP_MAX(A1->x, A3->x);
488 gbox->ymax = FP_MAX(A1->y, A3->y);
489
490 /* Create points for the possible extrema */
491 xmin.x = C.x - radius_A;
492 xmin.y = C.y;
493 ymin.x = C.x;
494 ymin.y = C.y - radius_A;
495 xmax.x = C.x + radius_A;
496 xmax.y = C.y;
497 ymax.x = C.x;
498 ymax.y = C.y + radius_A;
499
500 /* Divide the circle into two parts, one on each side of a line
501 joining p1 and p3. The circle extrema on the same side of that line
502 as p2 is on, are also the extrema of the bbox. */
503
504 A2_side = lw_segment_side(A1, A3, A2);
505
506 if ( A2_side == lw_segment_side(A1, A3, &xmin) )
507 gbox->xmin = xmin.x;
508
509 if ( A2_side == lw_segment_side(A1, A3, &ymin) )
510 gbox->ymin = ymin.y;
511
512 if ( A2_side == lw_segment_side(A1, A3, &xmax) )
513 gbox->xmax = xmax.x;
514
515 if ( A2_side == lw_segment_side(A1, A3, &ymax) )
516 gbox->ymax = ymax.y;
517
518 return LW_SUCCESS;
519 }
520
521
lw_arc_calculate_gbox_cartesian(const POINT4D * p1,const POINT4D * p2,const POINT4D * p3,GBOX * gbox)522 static int lw_arc_calculate_gbox_cartesian(const POINT4D *p1, const POINT4D *p2, const POINT4D *p3, GBOX *gbox)
523 {
524 int rv;
525
526 LWDEBUG(2, "lw_arc_calculate_gbox_cartesian called.");
527
528 rv = lw_arc_calculate_gbox_cartesian_2d((POINT2D*)p1, (POINT2D*)p2, (POINT2D*)p3, gbox);
529 gbox->zmin = FP_MIN(p1->z, p3->z);
530 gbox->mmin = FP_MIN(p1->m, p3->m);
531 gbox->zmax = FP_MAX(p1->z, p3->z);
532 gbox->mmax = FP_MAX(p1->m, p3->m);
533 return rv;
534 }
535
536 static void
ptarray_calculate_gbox_cartesian_2d(const POINTARRAY * pa,GBOX * gbox)537 ptarray_calculate_gbox_cartesian_2d(const POINTARRAY *pa, GBOX *gbox)
538 {
539 const POINT2D *p = getPoint2d_cp(pa, 0);
540
541 gbox->xmax = gbox->xmin = p->x;
542 gbox->ymax = gbox->ymin = p->y;
543
544 for (uint32_t i = 1; i < pa->npoints; i++)
545 {
546 p = getPoint2d_cp(pa, i);
547 gbox->xmin = FP_MIN(gbox->xmin, p->x);
548 gbox->xmax = FP_MAX(gbox->xmax, p->x);
549 gbox->ymin = FP_MIN(gbox->ymin, p->y);
550 gbox->ymax = FP_MAX(gbox->ymax, p->y);
551 }
552 }
553
554 /* Works with X/Y/Z. Needs to be adjusted after if X/Y/M was required */
555 static void
ptarray_calculate_gbox_cartesian_3d(const POINTARRAY * pa,GBOX * gbox)556 ptarray_calculate_gbox_cartesian_3d(const POINTARRAY *pa, GBOX *gbox)
557 {
558 const POINT3D *p = getPoint3d_cp(pa, 0);
559
560 gbox->xmax = gbox->xmin = p->x;
561 gbox->ymax = gbox->ymin = p->y;
562 gbox->zmax = gbox->zmin = p->z;
563
564 for (uint32_t i = 1; i < pa->npoints; i++)
565 {
566 p = getPoint3d_cp(pa, i);
567 gbox->xmin = FP_MIN(gbox->xmin, p->x);
568 gbox->xmax = FP_MAX(gbox->xmax, p->x);
569 gbox->ymin = FP_MIN(gbox->ymin, p->y);
570 gbox->ymax = FP_MAX(gbox->ymax, p->y);
571 gbox->zmin = FP_MIN(gbox->zmin, p->z);
572 gbox->zmax = FP_MAX(gbox->zmax, p->z);
573 }
574 }
575
576 static void
ptarray_calculate_gbox_cartesian_4d(const POINTARRAY * pa,GBOX * gbox)577 ptarray_calculate_gbox_cartesian_4d(const POINTARRAY *pa, GBOX *gbox)
578 {
579 const POINT4D *p = getPoint4d_cp(pa, 0);
580
581 gbox->xmax = gbox->xmin = p->x;
582 gbox->ymax = gbox->ymin = p->y;
583 gbox->zmax = gbox->zmin = p->z;
584 gbox->mmax = gbox->mmin = p->m;
585
586 for (uint32_t i = 1; i < pa->npoints; i++)
587 {
588 p = getPoint4d_cp(pa, i);
589 gbox->xmin = FP_MIN(gbox->xmin, p->x);
590 gbox->xmax = FP_MAX(gbox->xmax, p->x);
591 gbox->ymin = FP_MIN(gbox->ymin, p->y);
592 gbox->ymax = FP_MAX(gbox->ymax, p->y);
593 gbox->zmin = FP_MIN(gbox->zmin, p->z);
594 gbox->zmax = FP_MAX(gbox->zmax, p->z);
595 gbox->mmin = FP_MIN(gbox->mmin, p->m);
596 gbox->mmax = FP_MAX(gbox->mmax, p->m);
597 }
598 }
599
600 int
ptarray_calculate_gbox_cartesian(const POINTARRAY * pa,GBOX * gbox)601 ptarray_calculate_gbox_cartesian(const POINTARRAY *pa, GBOX *gbox)
602 {
603 if (!pa || pa->npoints == 0)
604 return LW_FAILURE;
605 if (!gbox)
606 return LW_FAILURE;
607
608 int has_z = FLAGS_GET_Z(pa->flags);
609 int has_m = FLAGS_GET_M(pa->flags);
610 gbox->flags = lwflags(has_z, has_m, 0);
611 LWDEBUGF(4, "ptarray_calculate_gbox Z: %d M: %d", has_z, has_m);
612 int coordinates = 2 + has_z + has_m;
613
614 switch (coordinates)
615 {
616 case 2:
617 {
618 ptarray_calculate_gbox_cartesian_2d(pa, gbox);
619 break;
620 }
621 case 3:
622 {
623 if (has_z)
624 {
625 ptarray_calculate_gbox_cartesian_3d(pa, gbox);
626 }
627 else
628 {
629 double zmin = gbox->zmin;
630 double zmax = gbox->zmax;
631 ptarray_calculate_gbox_cartesian_3d(pa, gbox);
632 gbox->mmin = gbox->zmin;
633 gbox->mmax = gbox->zmax;
634 gbox->zmin = zmin;
635 gbox->zmax = zmax;
636 }
637 break;
638 }
639 default:
640 {
641 ptarray_calculate_gbox_cartesian_4d(pa, gbox);
642 break;
643 }
644 }
645 return LW_SUCCESS;
646 }
647
lwcircstring_calculate_gbox_cartesian(LWCIRCSTRING * curve,GBOX * gbox)648 static int lwcircstring_calculate_gbox_cartesian(LWCIRCSTRING *curve, GBOX *gbox)
649 {
650 GBOX tmp;
651 POINT4D p1, p2, p3;
652 uint32_t i;
653
654 if (!curve) return LW_FAILURE;
655 if (curve->points->npoints < 3) return LW_FAILURE;
656
657 tmp.flags =
658 lwflags(FLAGS_GET_Z(curve->flags), FLAGS_GET_M(curve->flags), 0);
659
660 /* Initialize */
661 gbox->xmin = gbox->ymin = gbox->zmin = gbox->mmin = FLT_MAX;
662 gbox->xmax = gbox->ymax = gbox->zmax = gbox->mmax = -1*FLT_MAX;
663
664 for ( i = 2; i < curve->points->npoints; i += 2 )
665 {
666 getPoint4d_p(curve->points, i-2, &p1);
667 getPoint4d_p(curve->points, i-1, &p2);
668 getPoint4d_p(curve->points, i, &p3);
669
670 if (lw_arc_calculate_gbox_cartesian(&p1, &p2, &p3, &tmp) == LW_FAILURE)
671 continue;
672
673 gbox_merge(&tmp, gbox);
674 }
675
676 return LW_SUCCESS;
677 }
678
lwpoint_calculate_gbox_cartesian(LWPOINT * point,GBOX * gbox)679 static int lwpoint_calculate_gbox_cartesian(LWPOINT *point, GBOX *gbox)
680 {
681 if ( ! point ) return LW_FAILURE;
682 return ptarray_calculate_gbox_cartesian( point->point, gbox );
683 }
684
lwline_calculate_gbox_cartesian(LWLINE * line,GBOX * gbox)685 static int lwline_calculate_gbox_cartesian(LWLINE *line, GBOX *gbox)
686 {
687 if ( ! line ) return LW_FAILURE;
688 return ptarray_calculate_gbox_cartesian( line->points, gbox );
689 }
690
lwtriangle_calculate_gbox_cartesian(LWTRIANGLE * triangle,GBOX * gbox)691 static int lwtriangle_calculate_gbox_cartesian(LWTRIANGLE *triangle, GBOX *gbox)
692 {
693 if ( ! triangle ) return LW_FAILURE;
694 return ptarray_calculate_gbox_cartesian( triangle->points, gbox );
695 }
696
lwpoly_calculate_gbox_cartesian(LWPOLY * poly,GBOX * gbox)697 static int lwpoly_calculate_gbox_cartesian(LWPOLY *poly, GBOX *gbox)
698 {
699 if ( ! poly ) return LW_FAILURE;
700 if ( poly->nrings == 0 ) return LW_FAILURE;
701 /* Just need to check outer ring */
702 return ptarray_calculate_gbox_cartesian( poly->rings[0], gbox );
703 }
704
lwcollection_calculate_gbox_cartesian(LWCOLLECTION * coll,GBOX * gbox)705 static int lwcollection_calculate_gbox_cartesian(LWCOLLECTION *coll, GBOX *gbox)
706 {
707 GBOX subbox;
708 uint32_t i;
709 int result = LW_FAILURE;
710 int first = LW_TRUE;
711 assert(coll);
712 if ( (coll->ngeoms == 0) || !gbox)
713 return LW_FAILURE;
714
715 subbox.flags = coll->flags;
716
717 for ( i = 0; i < coll->ngeoms; i++ )
718 {
719 if ( lwgeom_calculate_gbox_cartesian((LWGEOM*)(coll->geoms[i]), &subbox) == LW_SUCCESS )
720 {
721 /* Keep a copy of the sub-bounding box for later
722 if ( coll->geoms[i]->bbox )
723 lwfree(coll->geoms[i]->bbox);
724 coll->geoms[i]->bbox = gbox_copy(&subbox); */
725 if ( first )
726 {
727 gbox_duplicate(&subbox, gbox);
728 first = LW_FALSE;
729 }
730 else
731 {
732 gbox_merge(&subbox, gbox);
733 }
734 result = LW_SUCCESS;
735 }
736 }
737 return result;
738 }
739
lwgeom_calculate_gbox_cartesian(const LWGEOM * lwgeom,GBOX * gbox)740 int lwgeom_calculate_gbox_cartesian(const LWGEOM *lwgeom, GBOX *gbox)
741 {
742 if ( ! lwgeom ) return LW_FAILURE;
743 LWDEBUGF(4, "lwgeom_calculate_gbox got type (%d) - %s", lwgeom->type, lwtype_name(lwgeom->type));
744
745 switch (lwgeom->type)
746 {
747 case POINTTYPE:
748 return lwpoint_calculate_gbox_cartesian((LWPOINT *)lwgeom, gbox);
749 case LINETYPE:
750 return lwline_calculate_gbox_cartesian((LWLINE *)lwgeom, gbox);
751 case CIRCSTRINGTYPE:
752 return lwcircstring_calculate_gbox_cartesian((LWCIRCSTRING *)lwgeom, gbox);
753 case POLYGONTYPE:
754 return lwpoly_calculate_gbox_cartesian((LWPOLY *)lwgeom, gbox);
755 case TRIANGLETYPE:
756 return lwtriangle_calculate_gbox_cartesian((LWTRIANGLE *)lwgeom, gbox);
757 case COMPOUNDTYPE:
758 case CURVEPOLYTYPE:
759 case MULTIPOINTTYPE:
760 case MULTILINETYPE:
761 case MULTICURVETYPE:
762 case MULTIPOLYGONTYPE:
763 case MULTISURFACETYPE:
764 case POLYHEDRALSURFACETYPE:
765 case TINTYPE:
766 case COLLECTIONTYPE:
767 return lwcollection_calculate_gbox_cartesian((LWCOLLECTION *)lwgeom, gbox);
768 }
769 /* Never get here, please. */
770 lwerror("unsupported type (%d) - %s", lwgeom->type, lwtype_name(lwgeom->type));
771 return LW_FAILURE;
772 }
773
gbox_float_round(GBOX * gbox)774 void gbox_float_round(GBOX *gbox)
775 {
776 gbox->xmin = next_float_down(gbox->xmin);
777 gbox->xmax = next_float_up(gbox->xmax);
778
779 gbox->ymin = next_float_down(gbox->ymin);
780 gbox->ymax = next_float_up(gbox->ymax);
781
782 if ( FLAGS_GET_M(gbox->flags) )
783 {
784 gbox->mmin = next_float_down(gbox->mmin);
785 gbox->mmax = next_float_up(gbox->mmax);
786 }
787
788 if ( FLAGS_GET_Z(gbox->flags) )
789 {
790 gbox->zmin = next_float_down(gbox->zmin);
791 gbox->zmax = next_float_up(gbox->zmax);
792 }
793 }
794
795 inline static uint64_t
uint64_interleave_2(uint64_t x,uint64_t y)796 uint64_interleave_2(uint64_t x, uint64_t y)
797 {
798 x = (x | (x << 16)) & 0x0000FFFF0000FFFFULL;
799 x = (x | (x << 8)) & 0x00FF00FF00FF00FFULL;
800 x = (x | (x << 4)) & 0x0F0F0F0F0F0F0F0FULL;
801 x = (x | (x << 2)) & 0x3333333333333333ULL;
802 x = (x | (x << 1)) & 0x5555555555555555ULL;
803
804 y = (y | (y << 16)) & 0x0000FFFF0000FFFFULL;
805 y = (y | (y << 8)) & 0x00FF00FF00FF00FFULL;
806 y = (y | (y << 4)) & 0x0F0F0F0F0F0F0F0FULL;
807 y = (y | (y << 2)) & 0x3333333333333333ULL;
808 y = (y | (y << 1)) & 0x5555555555555555ULL;
809
810 return x | (y << 1);
811 }
812
813 /* Based on https://github.com/rawrunprotected/hilbert_curves Public Domain code */
814 inline static uint64_t
uint32_hilbert(uint32_t px,uint32_t py)815 uint32_hilbert(uint32_t px, uint32_t py)
816 {
817 uint64_t x = px;
818 uint64_t y = py;
819
820 uint64_t A, B, C, D;
821
822 // Initial prefix scan round, prime with x and y
823 {
824 uint64_t a = x ^ y;
825 uint64_t b = 0xFFFFFFFFULL ^ a;
826 uint64_t c = 0xFFFFFFFFULL ^ (x | y);
827 uint64_t d = x & (y ^ 0xFFFFFFFFULL);
828
829 A = a | (b >> 1);
830 B = (a >> 1) ^ a;
831 C = ((c >> 1) ^ (b & (d >> 1))) ^ c;
832 D = ((a & (c >> 1)) ^ (d >> 1)) ^ d;
833 }
834
835 {
836 uint64_t a = A;
837 uint64_t b = B;
838 uint64_t c = C;
839 uint64_t d = D;
840
841 A = ((a & (a >> 2)) ^ (b & (b >> 2)));
842 B = ((a & (b >> 2)) ^ (b & ((a ^ b) >> 2)));
843 C ^= ((a & (c >> 2)) ^ (b & (d >> 2)));
844 D ^= ((b & (c >> 2)) ^ ((a ^ b) & (d >> 2)));
845 }
846
847 {
848 uint64_t a = A;
849 uint64_t b = B;
850 uint64_t c = C;
851 uint64_t d = D;
852
853 A = ((a & (a >> 4)) ^ (b & (b >> 4)));
854 B = ((a & (b >> 4)) ^ (b & ((a ^ b) >> 4)));
855 C ^= ((a & (c >> 4)) ^ (b & (d >> 4)));
856 D ^= ((b & (c >> 4)) ^ ((a ^ b) & (d >> 4)));
857 }
858
859 {
860 uint64_t a = A;
861 uint64_t b = B;
862 uint64_t c = C;
863 uint64_t d = D;
864
865 A = ((a & (a >> 8)) ^ (b & (b >> 8)));
866 B = ((a & (b >> 8)) ^ (b & ((a ^ b) >> 8)));
867 C ^= ((a & (c >> 8)) ^ (b & (d >> 8)));
868 D ^= ((b & (c >> 8)) ^ ((a ^ b) & (d >> 8)));
869 }
870
871 {
872 uint64_t a = A;
873 uint64_t b = B;
874 uint64_t c = C;
875 uint64_t d = D;
876
877 C ^= ((a & (c >> 16)) ^ (b & (d >> 16)));
878 D ^= ((b & (c >> 16)) ^ ((a ^ b) & (d >> 16)));
879 }
880
881 // Undo transformation prefix scan
882 uint64_t a = C ^ (C >> 1);
883 uint64_t b = D ^ (D >> 1);
884
885 // Recover index bits
886 uint64_t i0 = x ^ y;
887 uint64_t i1 = b | (0xFFFFFFFFULL ^ (i0 | a));
888
889 return uint64_interleave_2(i0, i1);
890 }
891
892 uint64_t
gbox_get_sortable_hash(const GBOX * g,const int32_t srid)893 gbox_get_sortable_hash(const GBOX *g, const int32_t srid)
894 {
895 union floatuint {
896 uint32_t u;
897 float f;
898 };
899
900 union floatuint x, y;
901
902 /*
903 * Since in theory the bitwise representation of an IEEE
904 * float is sortable (exponents come before mantissa, etc)
905 * we just copy the bits directly into an int and then
906 * interleave those ints.
907 */
908 if (FLAGS_GET_GEODETIC(g->flags))
909 {
910 GEOGRAPHIC_POINT gpt;
911 POINT3D p;
912 p.x = (g->xmax + g->xmin) / 2.0;
913 p.y = (g->ymax + g->ymin) / 2.0;
914 p.z = (g->zmax + g->zmin) / 2.0;
915 normalize(&p);
916 cart2geog(&p, &gpt);
917 /* We know range for geography, so build the curve taking it into account */
918 x.f = 1.5 + gpt.lon / 512.0;
919 y.f = 1.5 + gpt.lat / 256.0;
920 }
921 else
922 {
923 x.f = (g->xmax + g->xmin) / 2;
924 y.f = (g->ymax + g->ymin) / 2;
925 /*
926 * Tweak for popular SRID values: push floating point values into 1..2 range,
927 * a region where exponent is constant and thus Hilbert curve
928 * doesn't have compression artifact when X or Y value is close to 0.
929 * If someone has out of bounds value it will still expose the arifact but not crash.
930 * TODO: reconsider when we will have machinery to properly get bounds by SRID.
931 */
932 if (srid == 3857 || srid == 3395)
933 {
934 x.f = 1.5 + x.f / 67108864.0;
935 y.f = 1.5 + y.f / 67108864.0;
936 }
937 else if (srid == 4326)
938 {
939 x.f = 1.5 + x.f / 512.0;
940 y.f = 1.5 + y.f / 256.0;
941 }
942 }
943
944 return uint32_hilbert(y.u, x.u);
945 }
946