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 2001-2006 Refractions Research Inc.
22 * Copyright 2017 Darafei Praliaskouski <me@komzpa.net>
23 *
24 **********************************************************************/
25
26
27
28 #include "liblwgeom_internal.h"
29 #include "lwgeom_log.h"
30
31 #include <stdio.h>
32 #include <errno.h>
33 #include <assert.h>
34 #include "../postgis_svn_revision.h"
35
36 const char *
lwgeom_version()37 lwgeom_version()
38 {
39 static char *ptr = NULL;
40 static char buf[256];
41 if ( ! ptr )
42 {
43 ptr = buf;
44 snprintf(ptr, 256, LIBLWGEOM_VERSION" r%d", POSTGIS_SVN_REVISION);
45 }
46
47 return ptr;
48 }
49
50
51 inline float
next_float_down(double d)52 next_float_down(double d)
53 {
54 float result;
55 if (d > (double)FLT_MAX)
56 return FLT_MAX;
57 if (d <= (double)-FLT_MAX)
58 return -FLT_MAX;
59 result = d;
60
61 if ( ((double)result) <=d )
62 return result;
63
64 return nextafterf(result, -1*FLT_MAX);
65
66 }
67
68 /*
69 * Returns the float that's very close to the input, but >=.
70 * handles the funny differences in float4 and float8 reps.
71 */
72 inline float
next_float_up(double d)73 next_float_up(double d)
74 {
75 float result;
76 if (d >= (double)FLT_MAX)
77 return FLT_MAX;
78 if (d < (double)-FLT_MAX)
79 return -FLT_MAX;
80 result = d;
81
82 if ( ((double)result) >=d )
83 return result;
84
85 return nextafterf(result, FLT_MAX);
86 }
87
88
89
90
91 /************************************************************************
92 * POINTARRAY support functions
93 *
94 * TODO: should be moved to ptarray.c probably
95 *
96 ************************************************************************/
97
98 /*
99 * Copies a point from the point array into the parameter point
100 * will set point's z=NO_Z_VALUE if pa is 2d
101 * will set point's m=NO_M_VALUE if pa is 3d or 2d
102 *
103 * NOTE: point is a real POINT3D *not* a pointer
104 */
105 POINT4D
getPoint4d(const POINTARRAY * pa,uint32_t n)106 getPoint4d(const POINTARRAY *pa, uint32_t n)
107 {
108 POINT4D result;
109 getPoint4d_p(pa, n, &result);
110 return result;
111 }
112
113 /*
114 * Copies a point from the point array into the parameter point
115 * will set point's z=NO_Z_VALUE if pa is 2d
116 * will set point's m=NO_M_VALUE if pa is 3d or 2d
117 *
118 * NOTE: this will modify the point4d pointed to by 'point'.
119 *
120 * @return 0 on error, 1 on success
121 */
122 int
getPoint4d_p(const POINTARRAY * pa,uint32_t n,POINT4D * op)123 getPoint4d_p(const POINTARRAY *pa, uint32_t n, POINT4D *op)
124 {
125 uint8_t *ptr;
126 int zmflag;
127
128 if ( ! pa )
129 {
130 lwerror("%s [%d] NULL POINTARRAY input", __FILE__, __LINE__);
131 return 0;
132 }
133
134 if ( n>=pa->npoints )
135 {
136 lwnotice("%s [%d] called with n=%d and npoints=%d", __FILE__, __LINE__, n, pa->npoints);
137 return 0;
138 }
139
140 LWDEBUG(4, "getPoint4d_p called.");
141
142 /* Get a pointer to nth point offset and zmflag */
143 ptr=getPoint_internal(pa, n);
144 zmflag=FLAGS_GET_ZM(pa->flags);
145
146 LWDEBUGF(4, "ptr %p, zmflag %d", ptr, zmflag);
147
148 switch (zmflag)
149 {
150 case 0: /* 2d */
151 memcpy(op, ptr, sizeof(POINT2D));
152 op->m=NO_M_VALUE;
153 op->z=NO_Z_VALUE;
154 break;
155
156 case 3: /* ZM */
157 memcpy(op, ptr, sizeof(POINT4D));
158 break;
159
160 case 2: /* Z */
161 memcpy(op, ptr, sizeof(POINT3DZ));
162 op->m=NO_M_VALUE;
163 break;
164
165 case 1: /* M */
166 memcpy(op, ptr, sizeof(POINT3DM));
167 op->m=op->z; /* we use Z as temporary storage */
168 op->z=NO_Z_VALUE;
169 break;
170
171 default:
172 lwerror("Unknown ZM flag ??");
173 return 0;
174 }
175 return 1;
176
177 }
178
179
180
181 /*
182 * Copy a point from the point array into the parameter point
183 * will set point's z=NO_Z_VALUE if pa is 2d
184 * NOTE: point is a real POINT3DZ *not* a pointer
185 */
186 POINT3DZ
getPoint3dz(const POINTARRAY * pa,uint32_t n)187 getPoint3dz(const POINTARRAY *pa, uint32_t n)
188 {
189 POINT3DZ result;
190 getPoint3dz_p(pa, n, &result);
191 return result;
192 }
193
194 /*
195 * Copy a point from the point array into the parameter point
196 * will set point's z=NO_Z_VALUE if pa is 2d
197 *
198 * NOTE: point is a real POINT3DZ *not* a pointer
199 */
200 POINT3DM
getPoint3dm(const POINTARRAY * pa,uint32_t n)201 getPoint3dm(const POINTARRAY *pa, uint32_t n)
202 {
203 POINT3DM result;
204 getPoint3dm_p(pa, n, &result);
205 return result;
206 }
207
208 /*
209 * Copy a point from the point array into the parameter point
210 * will set point's z=NO_Z_VALUE if pa is 2d
211 *
212 * NOTE: this will modify the point3dz pointed to by 'point'.
213 */
214 int
getPoint3dz_p(const POINTARRAY * pa,uint32_t n,POINT3DZ * op)215 getPoint3dz_p(const POINTARRAY *pa, uint32_t n, POINT3DZ *op)
216 {
217 uint8_t *ptr;
218
219 if ( ! pa )
220 {
221 lwerror("%s [%d] NULL POINTARRAY input", __FILE__, __LINE__);
222 return 0;
223 }
224
225 if ( n>=pa->npoints )
226 {
227 lwnotice("%s [%d] called with n=%d and npoints=%d", __FILE__, __LINE__, n, pa->npoints);
228 return 0;
229 }
230
231 LWDEBUGF(2, "getPoint3dz_p called on array of %d-dimensions / %u pts",
232 FLAGS_NDIMS(pa->flags), pa->npoints);
233
234 /* Get a pointer to nth point offset */
235 ptr=getPoint_internal(pa, n);
236
237 /*
238 * if input POINTARRAY has the Z, it is always
239 * at third position so make a single copy
240 */
241 if ( FLAGS_GET_Z(pa->flags) )
242 {
243 memcpy(op, ptr, sizeof(POINT3DZ));
244 }
245
246 /*
247 * Otherwise copy the 2d part and initialize
248 * Z to NO_Z_VALUE
249 */
250 else
251 {
252 memcpy(op, ptr, sizeof(POINT2D));
253 op->z=NO_Z_VALUE;
254 }
255
256 return 1;
257
258 }
259
260 /*
261 * Copy a point from the point array into the parameter point
262 * will set point's m=NO_Z_VALUE if pa has no M
263 *
264 * NOTE: this will modify the point3dm pointed to by 'point'.
265 */
266 int
getPoint3dm_p(const POINTARRAY * pa,uint32_t n,POINT3DM * op)267 getPoint3dm_p(const POINTARRAY *pa, uint32_t n, POINT3DM *op)
268 {
269 uint8_t *ptr;
270 int zmflag;
271
272 if ( ! pa )
273 {
274 lwerror("%s [%d] NULL POINTARRAY input", __FILE__, __LINE__);
275 return 0;
276 }
277
278 if ( n>=pa->npoints )
279 {
280 lwnotice("%s [%d] called with n=%d and npoints=%d", __FILE__, __LINE__, n, pa->npoints);
281 return 0;
282 }
283
284 LWDEBUGF(2, "getPoint3dm_p(%d) called on array of %d-dimensions / %u pts",
285 n, FLAGS_NDIMS(pa->flags), pa->npoints);
286
287
288 /* Get a pointer to nth point offset and zmflag */
289 ptr=getPoint_internal(pa, n);
290 zmflag=FLAGS_GET_ZM(pa->flags);
291
292 /*
293 * if input POINTARRAY has the M and NO Z,
294 * we can issue a single memcpy
295 */
296 if ( zmflag == 1 )
297 {
298 memcpy(op, ptr, sizeof(POINT3DM));
299 return 1;
300 }
301
302 /*
303 * Otherwise copy the 2d part and
304 * initialize M to NO_M_VALUE
305 */
306 memcpy(op, ptr, sizeof(POINT2D));
307
308 /*
309 * Then, if input has Z skip it and
310 * copy next double, otherwise initialize
311 * M to NO_M_VALUE
312 */
313 if ( zmflag == 3 )
314 {
315 ptr+=sizeof(POINT3DZ);
316 memcpy(&(op->m), ptr, sizeof(double));
317 }
318 else
319 {
320 op->m=NO_M_VALUE;
321 }
322
323 return 1;
324 }
325
326
327 /*
328 * Copy a point from the point array into the parameter point
329 * z value (if present) is not returned.
330 *
331 * NOTE: point is a real POINT2D *not* a pointer
332 */
333 POINT2D
getPoint2d(const POINTARRAY * pa,uint32_t n)334 getPoint2d(const POINTARRAY *pa, uint32_t n)
335 {
336 const POINT2D *result;
337 result = getPoint2d_cp(pa, n);
338 return *result;
339 }
340
341 /*
342 * Copy a point from the point array into the parameter point
343 * z value (if present) is not returned.
344 *
345 * NOTE: this will modify the point2d pointed to by 'point'.
346 */
347 int
getPoint2d_p(const POINTARRAY * pa,uint32_t n,POINT2D * point)348 getPoint2d_p(const POINTARRAY *pa, uint32_t n, POINT2D *point)
349 {
350 if ( ! pa )
351 {
352 lwerror("%s [%d] NULL POINTARRAY input", __FILE__, __LINE__);
353 return 0;
354 }
355
356 if ( n>=pa->npoints )
357 {
358 lwnotice("%s [%d] called with n=%d and npoints=%d", __FILE__, __LINE__, n, pa->npoints);
359 return 0;
360 }
361
362 /* this does x,y */
363 memcpy(point, getPoint_internal(pa, n), sizeof(POINT2D));
364 return 1;
365 }
366
367 /**
368 * Returns a pointer into the POINTARRAY serialized_ptlist,
369 * suitable for reading from. This is very high performance
370 * and declared const because you aren't allowed to muck with the
371 * values, only read them.
372 */
373 const POINT2D*
getPoint2d_cp(const POINTARRAY * pa,uint32_t n)374 getPoint2d_cp(const POINTARRAY *pa, uint32_t n)
375 {
376 if ( ! pa ) return 0;
377
378 if ( n>=pa->npoints )
379 {
380 lwerror("getPoint2d_cp: point offset out of range");
381 return 0; /*error */
382 }
383
384 return (const POINT2D*)getPoint_internal(pa, n);
385 }
386
387 const POINT3DZ*
getPoint3dz_cp(const POINTARRAY * pa,uint32_t n)388 getPoint3dz_cp(const POINTARRAY *pa, uint32_t n)
389 {
390 if ( ! pa ) return 0;
391
392 if ( ! FLAGS_GET_Z(pa->flags) )
393 {
394 lwerror("getPoint3dz_cp: no Z coordinates in point array");
395 return 0; /*error */
396 }
397
398 if ( n>=pa->npoints )
399 {
400 lwerror("getPoint3dz_cp: point offset out of range");
401 return 0; /*error */
402 }
403
404 return (const POINT3DZ*)getPoint_internal(pa, n);
405 }
406
407 const POINT4D*
getPoint4d_cp(const POINTARRAY * pa,uint32_t n)408 getPoint4d_cp(const POINTARRAY* pa, uint32_t n)
409 {
410 if (!pa) return 0;
411
412 if (!(FLAGS_GET_Z(pa->flags) && FLAGS_GET_M(pa->flags)))
413 {
414 lwerror("getPoint4d_cp: no Z and M coordinates in point array");
415 return 0; /*error */
416 }
417
418 if (n >= pa->npoints)
419 {
420 lwerror("getPoint4d_cp: point offset out of range");
421 return 0; /*error */
422 }
423
424 return (const POINT4D*)getPoint_internal(pa, n);
425 }
426
427 /*
428 * set point N to the given value
429 * NOTE that the pointarray can be of any
430 * dimension, the appropriate ordinate values
431 * will be extracted from it
432 *
433 */
434 void
ptarray_set_point4d(POINTARRAY * pa,uint32_t n,const POINT4D * p4d)435 ptarray_set_point4d(POINTARRAY *pa, uint32_t n, const POINT4D *p4d)
436 {
437 uint8_t *ptr;
438 assert(n < pa->npoints);
439 ptr=getPoint_internal(pa, n);
440 switch ( FLAGS_GET_ZM(pa->flags) )
441 {
442 case 3:
443 memcpy(ptr, p4d, sizeof(POINT4D));
444 break;
445 case 2:
446 memcpy(ptr, p4d, sizeof(POINT3DZ));
447 break;
448 case 1:
449 memcpy(ptr, p4d, sizeof(POINT2D));
450 ptr+=sizeof(POINT2D);
451 memcpy(ptr, &(p4d->m), sizeof(double));
452 break;
453 case 0:
454 memcpy(ptr, p4d, sizeof(POINT2D));
455 break;
456 }
457 }
458
459 void
ptarray_copy_point(POINTARRAY * pa,uint32_t from,uint32_t to)460 ptarray_copy_point(POINTARRAY *pa, uint32_t from, uint32_t to)
461 {
462 int ndims = FLAGS_NDIMS(pa->flags);
463 switch (ndims)
464 {
465 case 2:
466 {
467 POINT2D *p_from = (POINT2D*)(getPoint_internal(pa, from));
468 POINT2D *p_to = (POINT2D*)(getPoint_internal(pa, to));
469 *p_to = *p_from;
470 return;
471 }
472 case 3:
473 {
474 POINT3D *p_from = (POINT3D*)(getPoint_internal(pa, from));
475 POINT3D *p_to = (POINT3D*)(getPoint_internal(pa, to));
476 *p_to = *p_from;
477 return;
478 }
479 case 4:
480 {
481 POINT4D *p_from = (POINT4D*)(getPoint_internal(pa, from));
482 POINT4D *p_to = (POINT4D*)(getPoint_internal(pa, to));
483 *p_to = *p_from;
484 return;
485 }
486 default:
487 {
488 lwerror("%s: unsupported number of dimensions - %d", __func__, ndims);
489 return;
490 }
491 }
492 return;
493 }
494
495
496 /************************************************
497 * debugging routines
498 ************************************************/
499
printBOX3D(BOX3D * box)500 void printBOX3D(BOX3D *box)
501 {
502 lwnotice("BOX3D: %g %g, %g %g", box->xmin, box->ymin,
503 box->xmax, box->ymax);
504 }
505
printPA(POINTARRAY * pa)506 void printPA(POINTARRAY *pa)
507 {
508 uint32_t t;
509 POINT4D pt;
510 char *mflag;
511
512
513 if ( FLAGS_GET_M(pa->flags) ) mflag = "M";
514 else mflag = "";
515
516 lwnotice(" POINTARRAY%s{", mflag);
517 lwnotice(" ndims=%i, ptsize=%i",
518 FLAGS_NDIMS(pa->flags), ptarray_point_size(pa));
519 lwnotice(" npoints = %i", pa->npoints);
520
521 for (t =0; t<pa->npoints; t++)
522 {
523 getPoint4d_p(pa, t, &pt);
524 if (FLAGS_NDIMS(pa->flags) == 2)
525 {
526 lwnotice(" %i : %lf,%lf",t,pt.x,pt.y);
527 }
528 if (FLAGS_NDIMS(pa->flags) == 3)
529 {
530 lwnotice(" %i : %lf,%lf,%lf",t,pt.x,pt.y,pt.z);
531 }
532 if (FLAGS_NDIMS(pa->flags) == 4)
533 {
534 lwnotice(" %i : %lf,%lf,%lf,%lf",t,pt.x,pt.y,pt.z,pt.m);
535 }
536 }
537
538 lwnotice(" }");
539 }
540
541
542 /**
543 * Given a string with at least 2 chars in it, convert them to
544 * a byte value. No error checking done!
545 */
546 uint8_t
parse_hex(char * str)547 parse_hex(char *str)
548 {
549 /* do this a little brute force to make it faster */
550
551 uint8_t result_high = 0;
552 uint8_t result_low = 0;
553
554 switch (str[0])
555 {
556 case '0' :
557 result_high = 0;
558 break;
559 case '1' :
560 result_high = 1;
561 break;
562 case '2' :
563 result_high = 2;
564 break;
565 case '3' :
566 result_high = 3;
567 break;
568 case '4' :
569 result_high = 4;
570 break;
571 case '5' :
572 result_high = 5;
573 break;
574 case '6' :
575 result_high = 6;
576 break;
577 case '7' :
578 result_high = 7;
579 break;
580 case '8' :
581 result_high = 8;
582 break;
583 case '9' :
584 result_high = 9;
585 break;
586 case 'A' :
587 case 'a' :
588 result_high = 10;
589 break;
590 case 'B' :
591 case 'b' :
592 result_high = 11;
593 break;
594 case 'C' :
595 case 'c' :
596 result_high = 12;
597 break;
598 case 'D' :
599 case 'd' :
600 result_high = 13;
601 break;
602 case 'E' :
603 case 'e' :
604 result_high = 14;
605 break;
606 case 'F' :
607 case 'f' :
608 result_high = 15;
609 break;
610 }
611 switch (str[1])
612 {
613 case '0' :
614 result_low = 0;
615 break;
616 case '1' :
617 result_low = 1;
618 break;
619 case '2' :
620 result_low = 2;
621 break;
622 case '3' :
623 result_low = 3;
624 break;
625 case '4' :
626 result_low = 4;
627 break;
628 case '5' :
629 result_low = 5;
630 break;
631 case '6' :
632 result_low = 6;
633 break;
634 case '7' :
635 result_low = 7;
636 break;
637 case '8' :
638 result_low = 8;
639 break;
640 case '9' :
641 result_low = 9;
642 break;
643 case 'A' :
644 case 'a' :
645 result_low = 10;
646 break;
647 case 'B' :
648 case 'b' :
649 result_low = 11;
650 break;
651 case 'C' :
652 case 'c' :
653 result_low = 12;
654 break;
655 case 'D' :
656 case 'd' :
657 result_low = 13;
658 break;
659 case 'E' :
660 case 'e' :
661 result_low = 14;
662 break;
663 case 'F' :
664 case 'f' :
665 result_low = 15;
666 break;
667 }
668 return (uint8_t) ((result_high<<4) + result_low);
669 }
670
671
672 /**
673 * Given one byte, populate result with two byte representing
674 * the hex number.
675 *
676 * Ie. deparse_hex( 255, mystr)
677 * -> mystr[0] = 'F' and mystr[1] = 'F'
678 *
679 * No error checking done
680 */
681 void
deparse_hex(uint8_t str,char * result)682 deparse_hex(uint8_t str, char *result)
683 {
684 int input_high;
685 int input_low;
686 static char outchr[]=
687 {
688 "0123456789ABCDEF"
689 };
690
691 input_high = (str>>4);
692 input_low = (str & 0x0F);
693
694 result[0] = outchr[input_high];
695 result[1] = outchr[input_low];
696
697 }
698
699
700 /**
701 * Find interpolation point I
702 * between point A and point B
703 * so that the len(AI) == len(AB)*F
704 * and I falls on AB segment.
705 *
706 * Example:
707 *
708 * F=0.5 : A----I----B
709 * F=1 : A---------B==I
710 * F=0 : A==I---------B
711 * F=.2 : A-I-------B
712 */
713 void
interpolate_point4d(const POINT4D * A,const POINT4D * B,POINT4D * I,double F)714 interpolate_point4d(const POINT4D *A, const POINT4D *B, POINT4D *I, double F)
715 {
716 #if PARANOIA_LEVEL > 0
717 if (F < 0 || F > 1) lwerror("interpolate_point4d: invalid F (%g)", F);
718 #endif
719 I->x=A->x+((B->x-A->x)*F);
720 I->y=A->y+((B->y-A->y)*F);
721 I->z=A->z+((B->z-A->z)*F);
722 I->m=A->m+((B->m-A->m)*F);
723 }
724
725
726 int _lwgeom_interrupt_requested = 0;
727 void
lwgeom_request_interrupt()728 lwgeom_request_interrupt() {
729 _lwgeom_interrupt_requested = 1;
730 }
731 void
lwgeom_cancel_interrupt()732 lwgeom_cancel_interrupt() {
733 _lwgeom_interrupt_requested = 0;
734 }
735
736 lwinterrupt_callback *_lwgeom_interrupt_callback = 0;
737 lwinterrupt_callback *
lwgeom_register_interrupt_callback(lwinterrupt_callback * cb)738 lwgeom_register_interrupt_callback(lwinterrupt_callback *cb) {
739 lwinterrupt_callback *old = _lwgeom_interrupt_callback;
740 _lwgeom_interrupt_callback = cb;
741 return old;
742 }
743