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 (C) 2012-2021 Sandro Santilli <strk@kbt.io>
22 * Copyright (C) 2001-2006 Refractions Research Inc.
23 *
24 **********************************************************************/
25
26
27 #include <stdio.h>
28 #include <string.h>
29 #include <stdlib.h>
30
31 #include "../postgis_config.h"
32 /*#define POSTGIS_DEBUG_LEVEL 4*/
33 #include "liblwgeom_internal.h"
34 #include "lwgeom_log.h"
35
36 int
ptarray_has_z(const POINTARRAY * pa)37 ptarray_has_z(const POINTARRAY *pa)
38 {
39 if ( ! pa ) return LW_FALSE;
40 return FLAGS_GET_Z(pa->flags);
41 }
42
43 int
ptarray_has_m(const POINTARRAY * pa)44 ptarray_has_m(const POINTARRAY *pa)
45 {
46 if ( ! pa ) return LW_FALSE;
47 return FLAGS_GET_M(pa->flags);
48 }
49
50 POINTARRAY*
ptarray_construct(char hasz,char hasm,uint32_t npoints)51 ptarray_construct(char hasz, char hasm, uint32_t npoints)
52 {
53 POINTARRAY *pa = ptarray_construct_empty(hasz, hasm, npoints);
54 pa->npoints = npoints;
55 return pa;
56 }
57
58 POINTARRAY*
ptarray_construct_empty(char hasz,char hasm,uint32_t maxpoints)59 ptarray_construct_empty(char hasz, char hasm, uint32_t maxpoints)
60 {
61 POINTARRAY *pa = lwalloc(sizeof(POINTARRAY));
62 pa->serialized_pointlist = NULL;
63
64 /* Set our dimensionality info on the bitmap */
65 pa->flags = lwflags(hasz, hasm, 0);
66
67 /* We will be allocating a bit of room */
68 pa->npoints = 0;
69 pa->maxpoints = maxpoints;
70
71 /* Allocate the coordinate array */
72 if ( maxpoints > 0 )
73 pa->serialized_pointlist = lwalloc(maxpoints * ptarray_point_size(pa));
74 else
75 pa->serialized_pointlist = NULL;
76
77 return pa;
78 }
79
80 /*
81 * Add a point into a pointarray. Only adds as many dimensions as the
82 * pointarray supports.
83 */
84 int
ptarray_insert_point(POINTARRAY * pa,const POINT4D * p,uint32_t where)85 ptarray_insert_point(POINTARRAY *pa, const POINT4D *p, uint32_t where)
86 {
87 if (!pa || !p)
88 return LW_FAILURE;
89 size_t point_size = ptarray_point_size(pa);
90 LWDEBUGF(5,"pa = %p; p = %p; where = %d", pa, p, where);
91 LWDEBUGF(5,"pa->npoints = %d; pa->maxpoints = %d", pa->npoints, pa->maxpoints);
92
93 if ( FLAGS_GET_READONLY(pa->flags) )
94 {
95 lwerror("ptarray_insert_point: called on read-only point array");
96 return LW_FAILURE;
97 }
98
99 /* Error on invalid offset value */
100 if ( where > pa->npoints )
101 {
102 lwerror("ptarray_insert_point: offset out of range (%d)", where);
103 return LW_FAILURE;
104 }
105
106 /* If we have no storage, let's allocate some */
107 if( pa->maxpoints == 0 || ! pa->serialized_pointlist )
108 {
109 pa->maxpoints = 32;
110 pa->npoints = 0;
111 pa->serialized_pointlist = lwalloc(ptarray_point_size(pa) * pa->maxpoints);
112 }
113
114 /* Error out if we have a bad situation */
115 if ( pa->npoints > pa->maxpoints )
116 {
117 lwerror("npoints (%d) is greater than maxpoints (%d)", pa->npoints, pa->maxpoints);
118 return LW_FAILURE;
119 }
120
121 /* Check if we have enough storage, add more if necessary */
122 if( pa->npoints == pa->maxpoints )
123 {
124 pa->maxpoints *= 2;
125 pa->serialized_pointlist = lwrealloc(pa->serialized_pointlist, ptarray_point_size(pa) * pa->maxpoints);
126 }
127
128 /* Make space to insert the new point */
129 if( where < pa->npoints )
130 {
131 size_t copy_size = point_size * (pa->npoints - where);
132 memmove(getPoint_internal(pa, where+1), getPoint_internal(pa, where), copy_size);
133 LWDEBUGF(5,"copying %d bytes to start vertex %d from start vertex %d", copy_size, where+1, where);
134 }
135
136 /* We have one more point */
137 ++pa->npoints;
138
139 /* Copy the new point into the gap */
140 ptarray_set_point4d(pa, where, p);
141 LWDEBUGF(5,"copying new point to start vertex %d", point_size, where);
142
143 return LW_SUCCESS;
144 }
145
146 int
ptarray_append_point(POINTARRAY * pa,const POINT4D * pt,int repeated_points)147 ptarray_append_point(POINTARRAY *pa, const POINT4D *pt, int repeated_points)
148 {
149 /* Check for pathology */
150 if( ! pa || ! pt )
151 {
152 lwerror("ptarray_append_point: null input");
153 return LW_FAILURE;
154 }
155
156 /* Check for duplicate end point */
157 if ( repeated_points == LW_FALSE && pa->npoints > 0 )
158 {
159 POINT4D tmp;
160 getPoint4d_p(pa, pa->npoints-1, &tmp);
161 LWDEBUGF(4,"checking for duplicate end point (pt = POINT(%g %g) pa->npoints-q = POINT(%g %g))",pt->x,pt->y,tmp.x,tmp.y);
162
163 /* Return LW_SUCCESS and do nothing else if previous point in list is equal to this one */
164 if ( (pt->x == tmp.x) && (pt->y == tmp.y) &&
165 (FLAGS_GET_Z(pa->flags) ? pt->z == tmp.z : 1) &&
166 (FLAGS_GET_M(pa->flags) ? pt->m == tmp.m : 1) )
167 {
168 return LW_SUCCESS;
169 }
170 }
171
172 /* Append is just a special case of insert */
173 return ptarray_insert_point(pa, pt, pa->npoints);
174 }
175
176 int
ptarray_append_ptarray(POINTARRAY * pa1,POINTARRAY * pa2,double gap_tolerance)177 ptarray_append_ptarray(POINTARRAY *pa1, POINTARRAY *pa2, double gap_tolerance)
178 {
179 unsigned int poff = 0;
180 unsigned int npoints;
181 unsigned int ncap;
182 unsigned int ptsize;
183
184 /* Check for pathology */
185 if( ! pa1 || ! pa2 )
186 {
187 lwerror("ptarray_append_ptarray: null input");
188 return LW_FAILURE;
189 }
190
191 npoints = pa2->npoints;
192
193 if ( ! npoints ) return LW_SUCCESS; /* nothing more to do */
194
195 if( FLAGS_GET_READONLY(pa1->flags) )
196 {
197 lwerror("ptarray_append_ptarray: target pointarray is read-only");
198 return LW_FAILURE;
199 }
200
201 if( FLAGS_GET_ZM(pa1->flags) != FLAGS_GET_ZM(pa2->flags) )
202 {
203 lwerror("ptarray_append_ptarray: appending mixed dimensionality is not allowed");
204 return LW_FAILURE;
205 }
206
207 ptsize = ptarray_point_size(pa1);
208
209 /* Check for duplicate end point */
210 if ( pa1->npoints )
211 {
212 POINT2D tmp1, tmp2;
213 getPoint2d_p(pa1, pa1->npoints-1, &tmp1);
214 getPoint2d_p(pa2, 0, &tmp2);
215
216 /* If the end point and start point are the same, then don't copy start point */
217 if (p2d_same(&tmp1, &tmp2)) {
218 poff = 1;
219 --npoints;
220 }
221 else if ( gap_tolerance == 0 || ( gap_tolerance > 0 &&
222 distance2d_pt_pt(&tmp1, &tmp2) > gap_tolerance ) )
223 {
224 lwerror("Second line start point too far from first line end point");
225 return LW_FAILURE;
226 }
227 }
228
229 /* Check if we need extra space */
230 ncap = pa1->npoints + npoints;
231 if ( pa1->maxpoints < ncap )
232 {
233 pa1->maxpoints = ncap > pa1->maxpoints*2 ?
234 ncap : pa1->maxpoints*2;
235 pa1->serialized_pointlist = lwrealloc(pa1->serialized_pointlist, ptsize * pa1->maxpoints);
236 }
237
238 memcpy(getPoint_internal(pa1, pa1->npoints),
239 getPoint_internal(pa2, poff), ptsize * npoints);
240
241 pa1->npoints = ncap;
242
243 return LW_SUCCESS;
244 }
245
246 /*
247 * Add a point into a pointarray. Only adds as many dimensions as the
248 * pointarray supports.
249 */
250 int
ptarray_remove_point(POINTARRAY * pa,uint32_t where)251 ptarray_remove_point(POINTARRAY *pa, uint32_t where)
252 {
253 /* Check for pathology */
254 if( ! pa )
255 {
256 lwerror("ptarray_remove_point: null input");
257 return LW_FAILURE;
258 }
259
260 /* Error on invalid offset value */
261 if ( where >= pa->npoints )
262 {
263 lwerror("ptarray_remove_point: offset out of range (%d)", where);
264 return LW_FAILURE;
265 }
266
267 /* If the point is any but the last, we need to copy the data back one point */
268 if (where < pa->npoints - 1)
269 memmove(getPoint_internal(pa, where),
270 getPoint_internal(pa, where + 1),
271 ptarray_point_size(pa) * (pa->npoints - where - 1));
272
273 /* We have one less point */
274 pa->npoints--;
275
276 return LW_SUCCESS;
277 }
278
279 /**
280 * Build a new #POINTARRAY, but on top of someone else's ordinate array.
281 * Flag as read-only, so that ptarray_free() does not free the serialized_ptlist
282 */
ptarray_construct_reference_data(char hasz,char hasm,uint32_t npoints,uint8_t * ptlist)283 POINTARRAY* ptarray_construct_reference_data(char hasz, char hasm, uint32_t npoints, uint8_t *ptlist)
284 {
285 POINTARRAY *pa = lwalloc(sizeof(POINTARRAY));
286 LWDEBUGF(5, "hasz = %d, hasm = %d, npoints = %d, ptlist = %p", hasz, hasm, npoints, ptlist);
287 pa->flags = lwflags(hasz, hasm, 0);
288 FLAGS_SET_READONLY(pa->flags, 1); /* We don't own this memory, so we can't alter or free it. */
289 pa->npoints = npoints;
290 pa->maxpoints = npoints;
291 pa->serialized_pointlist = ptlist;
292 return pa;
293 }
294
295
296 POINTARRAY*
ptarray_construct_copy_data(char hasz,char hasm,uint32_t npoints,const uint8_t * ptlist)297 ptarray_construct_copy_data(char hasz, char hasm, uint32_t npoints, const uint8_t *ptlist)
298 {
299 POINTARRAY *pa = lwalloc(sizeof(POINTARRAY));
300
301 pa->flags = lwflags(hasz, hasm, 0);
302 pa->npoints = npoints;
303 pa->maxpoints = npoints;
304
305 if ( npoints > 0 )
306 {
307 pa->serialized_pointlist = lwalloc(ptarray_point_size(pa) * npoints);
308 memcpy(pa->serialized_pointlist, ptlist, ptarray_point_size(pa) * npoints);
309 }
310 else
311 {
312 pa->serialized_pointlist = NULL;
313 }
314
315 return pa;
316 }
317
318 void
ptarray_free(POINTARRAY * pa)319 ptarray_free(POINTARRAY *pa)
320 {
321 if (pa)
322 {
323 if (pa->serialized_pointlist && (!FLAGS_GET_READONLY(pa->flags)))
324 lwfree(pa->serialized_pointlist);
325 lwfree(pa);
326 }
327 }
328
329
330 void
ptarray_reverse_in_place(POINTARRAY * pa)331 ptarray_reverse_in_place(POINTARRAY *pa)
332 {
333 if (!pa->npoints)
334 return;
335 uint32_t i;
336 uint32_t last = pa->npoints - 1;
337 uint32_t mid = pa->npoints / 2;
338
339 double *d = (double*)(pa->serialized_pointlist);
340 int j;
341 int ndims = FLAGS_NDIMS(pa->flags);
342 for (i = 0; i < mid; i++)
343 {
344 for (j = 0; j < ndims; j++)
345 {
346 double buf;
347 buf = d[i*ndims+j];
348 d[i*ndims+j] = d[(last-i)*ndims+j];
349 d[(last-i)*ndims+j] = buf;
350 }
351 }
352 return;
353 }
354
355
356 /**
357 * Reverse X and Y axis on a given POINTARRAY
358 */
359 POINTARRAY*
ptarray_flip_coordinates(POINTARRAY * pa)360 ptarray_flip_coordinates(POINTARRAY *pa)
361 {
362 uint32_t i;
363 double d;
364 POINT4D p;
365
366 for (i=0 ; i < pa->npoints ; i++)
367 {
368 getPoint4d_p(pa, i, &p);
369 d = p.y;
370 p.y = p.x;
371 p.x = d;
372 ptarray_set_point4d(pa, i, &p);
373 }
374
375 return pa;
376 }
377
378 void
ptarray_swap_ordinates(POINTARRAY * pa,LWORD o1,LWORD o2)379 ptarray_swap_ordinates(POINTARRAY *pa, LWORD o1, LWORD o2)
380 {
381 uint32_t i;
382 double d, *dp1, *dp2;
383 POINT4D p;
384
385 dp1 = ((double*)&p)+(unsigned)o1;
386 dp2 = ((double*)&p)+(unsigned)o2;
387 for (i=0 ; i < pa->npoints ; i++)
388 {
389 getPoint4d_p(pa, i, &p);
390 d = *dp2;
391 *dp2 = *dp1;
392 *dp1 = d;
393 ptarray_set_point4d(pa, i, &p);
394 }
395 }
396
397 /**
398 * @brief Returns a modified #POINTARRAY so that no segment is
399 * longer than the given distance (computed using 2d).
400 *
401 * Every input point is kept.
402 * Z and M values for added points (if needed) are set proportionally.
403 */
404 POINTARRAY *
ptarray_segmentize2d(const POINTARRAY * ipa,double dist)405 ptarray_segmentize2d(const POINTARRAY *ipa, double dist)
406 {
407 double segdist;
408 POINT4D p1, p2;
409 POINT4D pbuf;
410 POINTARRAY *opa;
411 uint32_t i, j, nseg;
412 int hasz = FLAGS_GET_Z(ipa->flags);
413 int hasm = FLAGS_GET_M(ipa->flags);
414
415 pbuf.x = pbuf.y = pbuf.z = pbuf.m = 0;
416
417 /* Initial storage */
418 opa = ptarray_construct_empty(hasz, hasm, ipa->npoints);
419
420 /* Add first point */
421 getPoint4d_p(ipa, 0, &p1);
422 ptarray_append_point(opa, &p1, LW_FALSE);
423
424 /* Loop on all other input points */
425 for (i = 1; i < ipa->npoints; i++)
426 {
427 /*
428 * We use these pointers to avoid
429 * "strict-aliasing rules break" warning raised
430 * by gcc (3.3 and up).
431 *
432 * It looks that casting a variable address (also
433 * referred to as "type-punned pointer")
434 * breaks those "strict" rules.
435 */
436 POINT4D *p1ptr=&p1, *p2ptr=&p2;
437 double segments;
438
439 getPoint4d_p(ipa, i, &p2);
440
441 segdist = distance2d_pt_pt((POINT2D *)p1ptr, (POINT2D *)p2ptr);
442 /* Split input segment into shorter even chunks */
443 segments = ceil(segdist / dist);
444
445 /* Uses INT32_MAX instead of UINT32_MAX to be safe that it fits */
446 if (segments >= INT32_MAX)
447 {
448 lwnotice("%s:%d - %s: Too many segments required (%e)",
449 __FILE__, __LINE__,__func__, segments);
450 ptarray_free(opa);
451 return NULL;
452 }
453 nseg = segments;
454
455 for (j = 1; j < nseg; j++)
456 {
457 pbuf.x = p1.x + (p2.x - p1.x) * j / nseg;
458 pbuf.y = p1.y + (p2.y - p1.y) * j / nseg;
459 if (hasz)
460 pbuf.z = p1.z + (p2.z - p1.z) * j / nseg;
461 if (hasm)
462 pbuf.m = p1.m + (p2.m - p1.m) * j / nseg;
463 ptarray_append_point(opa, &pbuf, LW_FALSE);
464 LW_ON_INTERRUPT(ptarray_free(opa); return NULL);
465 }
466
467 ptarray_append_point(opa, &p2, (ipa->npoints == 2) ? LW_TRUE : LW_FALSE);
468 p1 = p2;
469 LW_ON_INTERRUPT(ptarray_free(opa); return NULL);
470 }
471
472 return opa;
473 }
474
475 char
ptarray_same(const POINTARRAY * pa1,const POINTARRAY * pa2)476 ptarray_same(const POINTARRAY *pa1, const POINTARRAY *pa2)
477 {
478 uint32_t i;
479 size_t ptsize;
480
481 if ( FLAGS_GET_ZM(pa1->flags) != FLAGS_GET_ZM(pa2->flags) ) return LW_FALSE;
482 LWDEBUG(5,"dimensions are the same");
483
484 if ( pa1->npoints != pa2->npoints ) return LW_FALSE;
485 LWDEBUG(5,"npoints are the same");
486
487 ptsize = ptarray_point_size(pa1);
488 LWDEBUGF(5, "ptsize = %d", ptsize);
489
490 for (i=0; i<pa1->npoints; i++)
491 {
492 if ( memcmp(getPoint_internal(pa1, i), getPoint_internal(pa2, i), ptsize) )
493 return LW_FALSE;
494 LWDEBUGF(5,"point #%d is the same",i);
495 }
496
497 return LW_TRUE;
498 }
499
500 POINTARRAY *
ptarray_addPoint(const POINTARRAY * pa,uint8_t * p,size_t pdims,uint32_t where)501 ptarray_addPoint(const POINTARRAY *pa, uint8_t *p, size_t pdims, uint32_t where)
502 {
503 POINTARRAY *ret;
504 POINT4D pbuf;
505 size_t ptsize = ptarray_point_size(pa);
506
507 LWDEBUGF(3, "pa %x p %x size %d where %d",
508 pa, p, pdims, where);
509
510 if ( pdims < 2 || pdims > 4 )
511 {
512 lwerror("ptarray_addPoint: point dimension out of range (%d)",
513 pdims);
514 return NULL;
515 }
516
517 if ( where > pa->npoints )
518 {
519 lwerror("ptarray_addPoint: offset out of range (%d)",
520 where);
521 return NULL;
522 }
523
524 LWDEBUG(3, "called with a %dD point");
525
526 pbuf.x = pbuf.y = pbuf.z = pbuf.m = 0.0;
527 memcpy((uint8_t *)&pbuf, p, pdims*sizeof(double));
528
529 LWDEBUG(3, "initialized point buffer");
530
531 ret = ptarray_construct(FLAGS_GET_Z(pa->flags),
532 FLAGS_GET_M(pa->flags), pa->npoints+1);
533
534
535 if ( where )
536 {
537 memcpy(getPoint_internal(ret, 0), getPoint_internal(pa, 0), ptsize*where);
538 }
539
540 memcpy(getPoint_internal(ret, where), (uint8_t *)&pbuf, ptsize);
541
542 if ( where+1 != ret->npoints )
543 {
544 memcpy(getPoint_internal(ret, where+1),
545 getPoint_internal(pa, where),
546 ptsize*(pa->npoints-where));
547 }
548
549 return ret;
550 }
551
552 POINTARRAY *
ptarray_removePoint(POINTARRAY * pa,uint32_t which)553 ptarray_removePoint(POINTARRAY *pa, uint32_t which)
554 {
555 POINTARRAY *ret;
556 size_t ptsize = ptarray_point_size(pa);
557
558 LWDEBUGF(3, "pa %x which %d", pa, which);
559
560 #if PARANOIA_LEVEL > 0
561 if ( which > pa->npoints-1 )
562 {
563 lwerror("%s [%d] offset (%d) out of range (%d..%d)", __FILE__, __LINE__,
564 which, 0, pa->npoints-1);
565 return NULL;
566 }
567
568 if ( pa->npoints < 3 )
569 {
570 lwerror("%s [%d] can't remove a point from a 2-vertex POINTARRAY", __FILE__, __LINE__);
571 return NULL;
572 }
573 #endif
574
575 ret = ptarray_construct(FLAGS_GET_Z(pa->flags),
576 FLAGS_GET_M(pa->flags), pa->npoints-1);
577
578 /* copy initial part */
579 if ( which )
580 {
581 memcpy(getPoint_internal(ret, 0), getPoint_internal(pa, 0), ptsize*which);
582 }
583
584 /* copy final part */
585 if ( which < pa->npoints-1 )
586 {
587 memcpy(getPoint_internal(ret, which), getPoint_internal(pa, which+1),
588 ptsize*(pa->npoints-which-1));
589 }
590
591 return ret;
592 }
593
594 POINTARRAY *
ptarray_merge(POINTARRAY * pa1,POINTARRAY * pa2)595 ptarray_merge(POINTARRAY *pa1, POINTARRAY *pa2)
596 {
597 POINTARRAY *pa;
598 size_t ptsize = ptarray_point_size(pa1);
599
600 if (FLAGS_GET_ZM(pa1->flags) != FLAGS_GET_ZM(pa2->flags))
601 lwerror("ptarray_cat: Mixed dimension");
602
603 pa = ptarray_construct( FLAGS_GET_Z(pa1->flags),
604 FLAGS_GET_M(pa1->flags),
605 pa1->npoints + pa2->npoints);
606
607 memcpy( getPoint_internal(pa, 0),
608 getPoint_internal(pa1, 0),
609 ptsize*(pa1->npoints));
610
611 memcpy( getPoint_internal(pa, pa1->npoints),
612 getPoint_internal(pa2, 0),
613 ptsize*(pa2->npoints));
614
615 ptarray_free(pa1);
616 ptarray_free(pa2);
617
618 return pa;
619 }
620
621
622 /**
623 * @brief Deep clone a pointarray (also clones serialized pointlist)
624 */
625 POINTARRAY *
ptarray_clone_deep(const POINTARRAY * in)626 ptarray_clone_deep(const POINTARRAY *in)
627 {
628 POINTARRAY *out = lwalloc(sizeof(POINTARRAY));
629
630 LWDEBUG(3, "ptarray_clone_deep called.");
631
632 out->flags = in->flags;
633 out->npoints = in->npoints;
634 out->maxpoints = in->npoints;
635
636 FLAGS_SET_READONLY(out->flags, 0);
637
638 if (!in->npoints)
639 {
640 // Avoid calling lwalloc of 0 bytes
641 out->serialized_pointlist = NULL;
642 }
643 else
644 {
645 size_t size = in->npoints * ptarray_point_size(in);
646 out->serialized_pointlist = lwalloc(size);
647 memcpy(out->serialized_pointlist, in->serialized_pointlist, size);
648 }
649
650 return out;
651 }
652
653 /**
654 * @brief Clone a POINTARRAY object. Serialized pointlist is not copied.
655 */
656 POINTARRAY *
ptarray_clone(const POINTARRAY * in)657 ptarray_clone(const POINTARRAY *in)
658 {
659 POINTARRAY *out = lwalloc(sizeof(POINTARRAY));
660
661 LWDEBUG(3, "ptarray_clone called.");
662
663 out->flags = in->flags;
664 out->npoints = in->npoints;
665 out->maxpoints = in->maxpoints;
666
667 FLAGS_SET_READONLY(out->flags, 1);
668
669 out->serialized_pointlist = in->serialized_pointlist;
670
671 return out;
672 }
673
674 /**
675 * Check for ring closure using whatever dimensionality is declared on the
676 * pointarray.
677 */
678 int
ptarray_is_closed(const POINTARRAY * in)679 ptarray_is_closed(const POINTARRAY *in)
680 {
681 if (!in)
682 {
683 lwerror("ptarray_is_closed: called with null point array");
684 return 0;
685 }
686 if (in->npoints <= 1 ) return in->npoints; /* single-point are closed, empty not closed */
687
688 return 0 == memcmp(getPoint_internal(in, 0), getPoint_internal(in, in->npoints-1), ptarray_point_size(in));
689 }
690
691
692 int
ptarray_is_closed_2d(const POINTARRAY * in)693 ptarray_is_closed_2d(const POINTARRAY *in)
694 {
695 if (!in)
696 {
697 lwerror("ptarray_is_closed_2d: called with null point array");
698 return 0;
699 }
700 if (in->npoints <= 1 ) return in->npoints; /* single-point are closed, empty not closed */
701
702 return 0 == memcmp(getPoint_internal(in, 0), getPoint_internal(in, in->npoints-1), sizeof(POINT2D) );
703 }
704
705 int
ptarray_is_closed_3d(const POINTARRAY * in)706 ptarray_is_closed_3d(const POINTARRAY *in)
707 {
708 if (!in)
709 {
710 lwerror("ptarray_is_closed_3d: called with null point array");
711 return 0;
712 }
713 if (in->npoints <= 1 ) return in->npoints; /* single-point are closed, empty not closed */
714
715 return 0 == memcmp(getPoint_internal(in, 0), getPoint_internal(in, in->npoints-1), sizeof(POINT3D) );
716 }
717
718 int
ptarray_is_closed_z(const POINTARRAY * in)719 ptarray_is_closed_z(const POINTARRAY *in)
720 {
721 if ( FLAGS_GET_Z(in->flags) )
722 return ptarray_is_closed_3d(in);
723 else
724 return ptarray_is_closed_2d(in);
725 }
726
727 /**
728 * Return 1 if the point is inside the POINTARRAY, -1 if it is outside,
729 * and 0 if it is on the boundary.
730 */
731 int
ptarray_contains_point(const POINTARRAY * pa,const POINT2D * pt)732 ptarray_contains_point(const POINTARRAY *pa, const POINT2D *pt)
733 {
734 return ptarray_contains_point_partial(pa, pt, LW_TRUE, NULL);
735 }
736
737 int
ptarray_contains_point_partial(const POINTARRAY * pa,const POINT2D * pt,int check_closed,int * winding_number)738 ptarray_contains_point_partial(const POINTARRAY *pa, const POINT2D *pt, int check_closed, int *winding_number)
739 {
740 int wn = 0;
741 uint32_t i;
742 double side;
743 const POINT2D *seg1;
744 const POINT2D *seg2;
745 double ymin, ymax;
746
747 seg1 = getPoint2d_cp(pa, 0);
748 seg2 = getPoint2d_cp(pa, pa->npoints-1);
749 if ( check_closed && ! p2d_same(seg1, seg2) )
750 lwerror("ptarray_contains_point called on unclosed ring");
751
752 for ( i=1; i < pa->npoints; i++ )
753 {
754 seg2 = getPoint2d_cp(pa, i);
755
756 /* Zero length segments are ignored. */
757 if ( seg1->x == seg2->x && seg1->y == seg2->y )
758 {
759 seg1 = seg2;
760 continue;
761 }
762
763 ymin = FP_MIN(seg1->y, seg2->y);
764 ymax = FP_MAX(seg1->y, seg2->y);
765
766 /* Only test segments in our vertical range */
767 if ( pt->y > ymax || pt->y < ymin )
768 {
769 seg1 = seg2;
770 continue;
771 }
772
773 side = lw_segment_side(seg1, seg2, pt);
774
775 /*
776 * A point on the boundary of a ring is not contained.
777 * WAS: if (fabs(side) < 1e-12), see #852
778 */
779 if ( (side == 0) && lw_pt_in_seg(pt, seg1, seg2) )
780 {
781 return LW_BOUNDARY;
782 }
783
784 /*
785 * If the point is to the left of the line, and it's rising,
786 * then the line is to the right of the point and
787 * circling counter-clockwise, so increment.
788 */
789 if ( (side < 0) && (seg1->y <= pt->y) && (pt->y < seg2->y) )
790 {
791 wn++;
792 }
793
794 /*
795 * If the point is to the right of the line, and it's falling,
796 * then the line is to the right of the point and circling
797 * clockwise, so decrement.
798 */
799 else if ( (side > 0) && (seg2->y <= pt->y) && (pt->y < seg1->y) )
800 {
801 wn--;
802 }
803
804 seg1 = seg2;
805 }
806
807 /* Sent out the winding number for calls that are building on this as a primitive */
808 if ( winding_number )
809 *winding_number = wn;
810
811 /* Outside */
812 if (wn == 0)
813 {
814 return LW_OUTSIDE;
815 }
816
817 /* Inside */
818 return LW_INSIDE;
819 }
820
821 /**
822 * For POINTARRAYs representing CIRCULARSTRINGS. That is, linked triples
823 * with each triple being control points of a circular arc. Such
824 * POINTARRAYs have an odd number of vertices.
825 *
826 * Return 1 if the point is inside the POINTARRAY, -1 if it is outside,
827 * and 0 if it is on the boundary.
828 */
829
830 int
ptarrayarc_contains_point(const POINTARRAY * pa,const POINT2D * pt)831 ptarrayarc_contains_point(const POINTARRAY *pa, const POINT2D *pt)
832 {
833 return ptarrayarc_contains_point_partial(pa, pt, LW_TRUE /* Check closed*/, NULL);
834 }
835
836 int
ptarrayarc_contains_point_partial(const POINTARRAY * pa,const POINT2D * pt,int check_closed,int * winding_number)837 ptarrayarc_contains_point_partial(const POINTARRAY *pa, const POINT2D *pt, int check_closed, int *winding_number)
838 {
839 int wn = 0;
840 uint32_t i;
841 int side;
842 const POINT2D *seg1;
843 const POINT2D *seg2;
844 const POINT2D *seg3;
845 GBOX gbox;
846
847 /* Check for not an arc ring (always have odd # of points) */
848 if ( (pa->npoints % 2) == 0 )
849 {
850 lwerror("ptarrayarc_contains_point called with even number of points");
851 return LW_OUTSIDE;
852 }
853
854 /* Check for not an arc ring (always have >= 3 points) */
855 if ( pa->npoints < 3 )
856 {
857 lwerror("ptarrayarc_contains_point called too-short pointarray");
858 return LW_OUTSIDE;
859 }
860
861 /* Check for unclosed case */
862 seg1 = getPoint2d_cp(pa, 0);
863 seg3 = getPoint2d_cp(pa, pa->npoints-1);
864 if ( check_closed && ! p2d_same(seg1, seg3) )
865 {
866 lwerror("ptarrayarc_contains_point called on unclosed ring");
867 return LW_OUTSIDE;
868 }
869 /* OK, it's closed. Is it just one circle? */
870 else if ( p2d_same(seg1, seg3) && pa->npoints == 3 )
871 {
872 double radius, d;
873 POINT2D c;
874 seg2 = getPoint2d_cp(pa, 1);
875
876 /* Wait, it's just a point, so it can't contain anything */
877 if ( lw_arc_is_pt(seg1, seg2, seg3) )
878 return LW_OUTSIDE;
879
880 /* See if the point is within the circle radius */
881 radius = lw_arc_center(seg1, seg2, seg3, &c);
882 d = distance2d_pt_pt(pt, &c);
883 if ( FP_EQUALS(d, radius) )
884 return LW_BOUNDARY; /* Boundary of circle */
885 else if ( d < radius )
886 return LW_INSIDE; /* Inside circle */
887 else
888 return LW_OUTSIDE; /* Outside circle */
889 }
890 else if ( p2d_same(seg1, pt) || p2d_same(seg3, pt) )
891 {
892 return LW_BOUNDARY; /* Boundary case */
893 }
894
895 /* Start on the ring */
896 seg1 = getPoint2d_cp(pa, 0);
897 for ( i=1; i < pa->npoints; i += 2 )
898 {
899 seg2 = getPoint2d_cp(pa, i);
900 seg3 = getPoint2d_cp(pa, i+1);
901
902 /* Catch an easy boundary case */
903 if( p2d_same(seg3, pt) )
904 return LW_BOUNDARY;
905
906 /* Skip arcs that have no size */
907 if ( lw_arc_is_pt(seg1, seg2, seg3) )
908 {
909 seg1 = seg3;
910 continue;
911 }
912
913 /* Only test segments in our vertical range */
914 lw_arc_calculate_gbox_cartesian_2d(seg1, seg2, seg3, &gbox);
915 if ( pt->y > gbox.ymax || pt->y < gbox.ymin )
916 {
917 seg1 = seg3;
918 continue;
919 }
920
921 /* Outside of horizontal range, and not between end points we also skip */
922 if ( (pt->x > gbox.xmax || pt->x < gbox.xmin) &&
923 (pt->y > FP_MAX(seg1->y, seg3->y) || pt->y < FP_MIN(seg1->y, seg3->y)) )
924 {
925 seg1 = seg3;
926 continue;
927 }
928
929 side = lw_arc_side(seg1, seg2, seg3, pt);
930
931 /* On the boundary */
932 if ( (side == 0) && lw_pt_in_arc(pt, seg1, seg2, seg3) )
933 {
934 return LW_BOUNDARY;
935 }
936
937 /* Going "up"! Point to left of arc. */
938 if ( side < 0 && (seg1->y <= pt->y) && (pt->y < seg3->y) )
939 {
940 wn++;
941 }
942
943 /* Going "down"! */
944 if ( side > 0 && (seg3->y <= pt->y) && (pt->y < seg1->y) )
945 {
946 wn--;
947 }
948
949 /* Inside the arc! */
950 if ( pt->x <= gbox.xmax && pt->x >= gbox.xmin )
951 {
952 POINT2D C;
953 double radius = lw_arc_center(seg1, seg2, seg3, &C);
954 double d = distance2d_pt_pt(pt, &C);
955
956 /* On the boundary! */
957 if ( d == radius )
958 return LW_BOUNDARY;
959
960 /* Within the arc! */
961 if ( d < radius )
962 {
963 /* Left side, increment winding number */
964 if ( side < 0 )
965 wn++;
966 /* Right side, decrement winding number */
967 if ( side > 0 )
968 wn--;
969 }
970 }
971
972 seg1 = seg3;
973 }
974
975 /* Sent out the winding number for calls that are building on this as a primitive */
976 if ( winding_number )
977 *winding_number = wn;
978
979 /* Outside */
980 if (wn == 0)
981 {
982 return LW_OUTSIDE;
983 }
984
985 /* Inside */
986 return LW_INSIDE;
987 }
988
989 /**
990 * Returns the area in cartesian units. Area is negative if ring is oriented CCW,
991 * positive if it is oriented CW and zero if the ring is degenerate or flat.
992 * http://en.wikipedia.org/wiki/Shoelace_formula
993 */
994 double
ptarray_signed_area(const POINTARRAY * pa)995 ptarray_signed_area(const POINTARRAY *pa)
996 {
997 const POINT2D *P1;
998 const POINT2D *P2;
999 const POINT2D *P3;
1000 double sum = 0.0;
1001 double x0, x, y1, y2;
1002 uint32_t i;
1003
1004 if (! pa || pa->npoints < 3 )
1005 return 0.0;
1006
1007 P1 = getPoint2d_cp(pa, 0);
1008 P2 = getPoint2d_cp(pa, 1);
1009 x0 = P1->x;
1010 for ( i = 2; i < pa->npoints; i++ )
1011 {
1012 P3 = getPoint2d_cp(pa, i);
1013 x = P2->x - x0;
1014 y1 = P3->y;
1015 y2 = P1->y;
1016 sum += x * (y2-y1);
1017
1018 /* Move forwards! */
1019 P1 = P2;
1020 P2 = P3;
1021 }
1022 return sum / 2.0;
1023 }
1024
1025 int
ptarray_isccw(const POINTARRAY * pa)1026 ptarray_isccw(const POINTARRAY *pa)
1027 {
1028 double area = 0;
1029 area = ptarray_signed_area(pa);
1030 if ( area > 0 ) return LW_FALSE;
1031 else return LW_TRUE;
1032 }
1033
1034 POINTARRAY*
ptarray_force_dims(const POINTARRAY * pa,int hasz,int hasm,double zval,double mval)1035 ptarray_force_dims(const POINTARRAY *pa, int hasz, int hasm, double zval, double mval)
1036 {
1037 /* TODO handle zero-length point arrays */
1038 uint32_t i;
1039 int in_hasz = FLAGS_GET_Z(pa->flags);
1040 int in_hasm = FLAGS_GET_M(pa->flags);
1041 POINT4D pt;
1042 POINTARRAY *pa_out = ptarray_construct_empty(hasz, hasm, pa->npoints);
1043
1044 for( i = 0; i < pa->npoints; i++ )
1045 {
1046 getPoint4d_p(pa, i, &pt);
1047 if( hasz && ! in_hasz )
1048 pt.z = zval;
1049 if( hasm && ! in_hasm )
1050 pt.m = mval;
1051 ptarray_append_point(pa_out, &pt, LW_TRUE);
1052 }
1053
1054 return pa_out;
1055 }
1056
1057 POINTARRAY *
ptarray_substring(POINTARRAY * ipa,double from,double to,double tolerance)1058 ptarray_substring(POINTARRAY *ipa, double from, double to, double tolerance)
1059 {
1060 POINTARRAY *dpa;
1061 POINT4D pt;
1062 POINT4D p1, p2;
1063 POINT4D *p1ptr=&p1; /* don't break strict-aliasing rule */
1064 POINT4D *p2ptr=&p2;
1065 int nsegs, i;
1066 double length, slength, tlength;
1067 int state = 0; /* 0=before, 1=inside */
1068
1069 /*
1070 * Create a dynamic pointarray with an initial capacity
1071 * equal to full copy of input points
1072 */
1073 dpa = ptarray_construct_empty(FLAGS_GET_Z(ipa->flags), FLAGS_GET_M(ipa->flags), ipa->npoints);
1074
1075 /* Compute total line length */
1076 length = ptarray_length_2d(ipa);
1077
1078
1079 LWDEBUGF(3, "Total length: %g", length);
1080
1081
1082 /* Get 'from' and 'to' lengths */
1083 from = length*from;
1084 to = length*to;
1085
1086
1087 LWDEBUGF(3, "From/To: %g/%g", from, to);
1088
1089
1090 tlength = 0;
1091 getPoint4d_p(ipa, 0, &p1);
1092 nsegs = ipa->npoints - 1;
1093 for ( i = 0; i < nsegs; i++ )
1094 {
1095 double dseg;
1096
1097 getPoint4d_p(ipa, i+1, &p2);
1098
1099
1100 LWDEBUGF(3 ,"Segment %d: (%g,%g,%g,%g)-(%g,%g,%g,%g)",
1101 i, p1.x, p1.y, p1.z, p1.m, p2.x, p2.y, p2.z, p2.m);
1102
1103
1104 /* Find the length of this segment */
1105 slength = distance2d_pt_pt((POINT2D *)p1ptr, (POINT2D *)p2ptr);
1106
1107 /*
1108 * We are before requested start.
1109 */
1110 if ( state == 0 ) /* before */
1111 {
1112
1113 LWDEBUG(3, " Before start");
1114
1115 if ( fabs ( from - ( tlength + slength ) ) <= tolerance )
1116 {
1117
1118 LWDEBUG(3, " Second point is our start");
1119
1120 /*
1121 * Second point is our start
1122 */
1123 ptarray_append_point(dpa, &p2, LW_FALSE);
1124 state=1; /* we're inside now */
1125 goto END;
1126 }
1127
1128 else if ( fabs(from - tlength) <= tolerance )
1129 {
1130
1131 LWDEBUG(3, " First point is our start");
1132
1133 /*
1134 * First point is our start
1135 */
1136 ptarray_append_point(dpa, &p1, LW_FALSE);
1137
1138 /*
1139 * We're inside now, but will check
1140 * 'to' point as well
1141 */
1142 state=1;
1143 }
1144
1145 /*
1146 * Didn't reach the 'from' point,
1147 * nothing to do
1148 */
1149 else if ( from > tlength + slength ) goto END;
1150
1151 else /* tlength < from < tlength+slength */
1152 {
1153
1154 LWDEBUG(3, " Seg contains first point");
1155
1156 /*
1157 * Our start is between first and
1158 * second point
1159 */
1160 dseg = (from - tlength) / slength;
1161
1162 interpolate_point4d(&p1, &p2, &pt, dseg);
1163
1164 ptarray_append_point(dpa, &pt, LW_FALSE);
1165
1166 /*
1167 * We're inside now, but will check
1168 * 'to' point as well
1169 */
1170 state=1;
1171 }
1172 }
1173
1174 if ( state == 1 ) /* inside */
1175 {
1176
1177 LWDEBUG(3, " Inside");
1178
1179 /*
1180 * 'to' point is our second point.
1181 */
1182 if ( fabs(to - ( tlength + slength ) ) <= tolerance )
1183 {
1184
1185 LWDEBUG(3, " Second point is our end");
1186
1187 ptarray_append_point(dpa, &p2, LW_FALSE);
1188 break; /* substring complete */
1189 }
1190
1191 /*
1192 * 'to' point is our first point.
1193 * (should only happen if 'to' is 0)
1194 */
1195 else if ( fabs(to - tlength) <= tolerance )
1196 {
1197
1198 LWDEBUG(3, " First point is our end");
1199
1200 ptarray_append_point(dpa, &p1, LW_FALSE);
1201
1202 break; /* substring complete */
1203 }
1204
1205 /*
1206 * Didn't reach the 'end' point,
1207 * just copy second point
1208 */
1209 else if ( to > tlength + slength )
1210 {
1211 ptarray_append_point(dpa, &p2, LW_FALSE);
1212 goto END;
1213 }
1214
1215 /*
1216 * 'to' point falls on this segment
1217 * Interpolate and break.
1218 */
1219 else if ( to < tlength + slength )
1220 {
1221
1222 LWDEBUG(3, " Seg contains our end");
1223
1224 dseg = (to - tlength) / slength;
1225 interpolate_point4d(&p1, &p2, &pt, dseg);
1226
1227 ptarray_append_point(dpa, &pt, LW_FALSE);
1228
1229 break;
1230 }
1231
1232 else
1233 {
1234 LWDEBUG(3, "Unhandled case");
1235 }
1236 }
1237
1238
1239 END:
1240
1241 tlength += slength;
1242 memcpy(&p1, &p2, sizeof(POINT4D));
1243 }
1244
1245 LWDEBUGF(3, "Out of loop, ptarray has %d points", dpa->npoints);
1246
1247 return dpa;
1248 }
1249
1250 /*
1251 * Write into the *ret argument coordinates of the closes point on
1252 * the given segment to the reference input point.
1253 */
1254 void
closest_point_on_segment(const POINT4D * p,const POINT4D * A,const POINT4D * B,POINT4D * ret)1255 closest_point_on_segment(const POINT4D *p, const POINT4D *A, const POINT4D *B, POINT4D *ret)
1256 {
1257 double r;
1258
1259 if ( FP_EQUALS(A->x, B->x) && FP_EQUALS(A->y, B->y) )
1260 {
1261 *ret = *A;
1262 return;
1263 }
1264
1265 /*
1266 * We use comp.graphics.algorithms Frequently Asked Questions method
1267 *
1268 * (1) AC dot AB
1269 * r = ----------
1270 * ||AB||^2
1271 * r has the following meaning:
1272 * r=0 P = A
1273 * r=1 P = B
1274 * r<0 P is on the backward extension of AB
1275 * r>1 P is on the forward extension of AB
1276 * 0<r<1 P is interior to AB
1277 *
1278 */
1279 r = ( (p->x-A->x) * (B->x-A->x) + (p->y-A->y) * (B->y-A->y) )/( (B->x-A->x)*(B->x-A->x) +(B->y-A->y)*(B->y-A->y) );
1280
1281 if (r<=0)
1282 {
1283 *ret = *A;
1284 return;
1285 }
1286 if (r>=1)
1287 {
1288 *ret = *B;
1289 return;
1290 }
1291
1292 ret->x = A->x + ( (B->x - A->x) * r );
1293 ret->y = A->y + ( (B->y - A->y) * r );
1294 ret->z = A->z + ( (B->z - A->z) * r );
1295 ret->m = A->m + ( (B->m - A->m) * r );
1296 }
1297
1298 int
ptarray_closest_segment_2d(const POINTARRAY * pa,const POINT2D * qp,double * dist)1299 ptarray_closest_segment_2d(const POINTARRAY *pa, const POINT2D *qp, double *dist)
1300 {
1301 const POINT2D *start = getPoint2d_cp(pa, 0), *end = NULL;
1302 uint32_t t, seg=0;
1303 double mindist=DBL_MAX;
1304
1305 /* Loop through pointarray looking for nearest segment */
1306 for (t=1; t<pa->npoints; t++)
1307 {
1308 double dist_sqr;
1309 end = getPoint2d_cp(pa, t);
1310 dist_sqr = distance2d_sqr_pt_seg(qp, start, end);
1311
1312 if (dist_sqr < mindist)
1313 {
1314 mindist = dist_sqr;
1315 seg=t-1;
1316 if ( mindist == 0 )
1317 {
1318 LWDEBUG(3, "Breaking on mindist=0");
1319 break;
1320 }
1321 }
1322
1323 start = end;
1324 }
1325
1326 if ( dist ) *dist = sqrt(mindist);
1327 return seg;
1328 }
1329
1330
1331 int
ptarray_closest_vertex_2d(const POINTARRAY * pa,const POINT2D * qp,double * dist)1332 ptarray_closest_vertex_2d(const POINTARRAY *pa, const POINT2D *qp, double *dist)
1333 {
1334 uint32_t t, pn=0;
1335 const POINT2D *p;
1336 double mindist = DBL_MAX;
1337
1338 /* Loop through pointarray looking for nearest segment */
1339 for (t=0; t<pa->npoints; t++)
1340 {
1341 double dist_sqr;
1342 p = getPoint2d_cp(pa, t);
1343 dist_sqr = distance2d_sqr_pt_pt(p, qp);
1344
1345 if (dist_sqr < mindist)
1346 {
1347 mindist = dist_sqr;
1348 pn = t;
1349 if ( mindist == 0 )
1350 {
1351 LWDEBUG(3, "Breaking on mindist=0");
1352 break;
1353 }
1354 }
1355 }
1356 if ( dist ) *dist = sqrt(mindist);
1357 return pn;
1358 }
1359
1360 /*
1361 * Given a point, returns the location of closest point on pointarray
1362 * and, optionally, it's actual distance from the point array.
1363 */
1364 double
ptarray_locate_point(const POINTARRAY * pa,const POINT4D * p4d,double * mindistout,POINT4D * proj4d)1365 ptarray_locate_point(const POINTARRAY *pa, const POINT4D *p4d, double *mindistout, POINT4D *proj4d)
1366 {
1367 double mindist=DBL_MAX;
1368 double tlen, plen;
1369 uint32_t t, seg=0;
1370 POINT4D start4d, end4d, projtmp;
1371 POINT2D proj, p;
1372 const POINT2D *start = NULL, *end = NULL;
1373
1374 /* Initialize our 2D copy of the input parameter */
1375 p.x = p4d->x;
1376 p.y = p4d->y;
1377
1378 if ( ! proj4d ) proj4d = &projtmp;
1379
1380 /* Check for special cases (length 0 and 1) */
1381 if ( pa->npoints <= 1 )
1382 {
1383 if ( pa->npoints == 1 )
1384 {
1385 getPoint4d_p(pa, 0, proj4d);
1386 if ( mindistout )
1387 *mindistout = distance2d_pt_pt(&p, getPoint2d_cp(pa, 0));
1388 }
1389 return 0.0;
1390 }
1391
1392 start = getPoint2d_cp(pa, 0);
1393 /* Loop through pointarray looking for nearest segment */
1394 for (t=1; t<pa->npoints; t++)
1395 {
1396 double dist_sqr;
1397 end = getPoint2d_cp(pa, t);
1398 dist_sqr = distance2d_sqr_pt_seg(&p, start, end);
1399
1400 if (dist_sqr < mindist)
1401 {
1402 mindist = dist_sqr;
1403 seg=t-1;
1404 if ( mindist == 0 )
1405 {
1406 LWDEBUG(3, "Breaking on mindist=0");
1407 break;
1408 }
1409 }
1410
1411 start = end;
1412 }
1413 mindist = sqrt(mindist);
1414
1415 if ( mindistout ) *mindistout = mindist;
1416
1417 LWDEBUGF(3, "Closest segment: %d", seg);
1418 LWDEBUGF(3, "mindist: %g", mindist);
1419
1420 /*
1421 * We need to project the
1422 * point on the closest segment.
1423 */
1424 getPoint4d_p(pa, seg, &start4d);
1425 getPoint4d_p(pa, seg+1, &end4d);
1426 closest_point_on_segment(p4d, &start4d, &end4d, proj4d);
1427
1428 /* Copy 4D values into 2D holder */
1429 proj.x = proj4d->x;
1430 proj.y = proj4d->y;
1431
1432 LWDEBUGF(3, "Closest segment:%d, npoints:%d", seg, pa->npoints);
1433
1434 /* For robustness, force 1 when closest point == endpoint */
1435 if ( (seg >= (pa->npoints-2)) && p2d_same(&proj, end) )
1436 {
1437 return 1.0;
1438 }
1439
1440 LWDEBUGF(3, "Closest point on segment: %g,%g", proj.x, proj.y);
1441
1442 tlen = ptarray_length_2d(pa);
1443
1444 LWDEBUGF(3, "tlen %g", tlen);
1445
1446 /* Location of any point on a zero-length line is 0 */
1447 /* See http://trac.osgeo.org/postgis/ticket/1772#comment:2 */
1448 if ( tlen == 0 ) return 0;
1449
1450 plen=0;
1451 start = getPoint2d_cp(pa, 0);
1452 for (t=0; t<seg; t++, start=end)
1453 {
1454 end = getPoint2d_cp(pa, t+1);
1455 plen += distance2d_pt_pt(start, end);
1456
1457 LWDEBUGF(4, "Segment %d made plen %g", t, plen);
1458 }
1459
1460 plen+=distance2d_pt_pt(&proj, start);
1461
1462 LWDEBUGF(3, "plen %g, tlen %g", plen, tlen);
1463
1464 return plen/tlen;
1465 }
1466
1467 /**
1468 * @brief Longitude shift for a pointarray.
1469 * Y remains the same
1470 * X is converted:
1471 * from -180..180 to 0..360
1472 * from 0..360 to -180..180
1473 * X < 0 becomes X + 360
1474 * X > 180 becomes X - 360
1475 */
1476 void
ptarray_longitude_shift(POINTARRAY * pa)1477 ptarray_longitude_shift(POINTARRAY *pa)
1478 {
1479 uint32_t i;
1480 double x;
1481
1482 for (i=0; i<pa->npoints; i++)
1483 {
1484 memcpy(&x, getPoint_internal(pa, i), sizeof(double));
1485 if ( x < 0 ) x+= 360;
1486 else if ( x > 180 ) x -= 360;
1487 memcpy(getPoint_internal(pa, i), &x, sizeof(double));
1488 }
1489 }
1490
1491
1492 /*
1493 * Returns a POINTARRAY with consecutive equal points
1494 * removed. Equality test on all dimensions of input.
1495 *
1496 * Always returns a newly allocated object.
1497 */
1498 static POINTARRAY *
ptarray_remove_repeated_points_minpoints(const POINTARRAY * in,double tolerance,int minpoints)1499 ptarray_remove_repeated_points_minpoints(const POINTARRAY *in, double tolerance, int minpoints)
1500 {
1501 POINTARRAY *out = ptarray_clone_deep(in);
1502 ptarray_remove_repeated_points_in_place(out, tolerance, minpoints);
1503 return out;
1504 }
1505
1506 POINTARRAY *
ptarray_remove_repeated_points(const POINTARRAY * in,double tolerance)1507 ptarray_remove_repeated_points(const POINTARRAY *in, double tolerance)
1508 {
1509 return ptarray_remove_repeated_points_minpoints(in, tolerance, 2);
1510 }
1511
1512
1513 void
ptarray_remove_repeated_points_in_place(POINTARRAY * pa,double tolerance,uint32_t min_points)1514 ptarray_remove_repeated_points_in_place(POINTARRAY *pa, double tolerance, uint32_t min_points)
1515 {
1516 uint32_t i;
1517 double tolsq = tolerance * tolerance;
1518 const POINT2D *last = NULL;
1519 const POINT2D *pt;
1520 uint32_t n_points = pa->npoints;
1521 uint32_t n_points_out = 1;
1522 size_t pt_size = ptarray_point_size(pa);
1523
1524 double dsq = FLT_MAX;
1525
1526 /* No-op on short inputs */
1527 if ( n_points <= min_points ) return;
1528
1529 last = getPoint2d_cp(pa, 0);
1530 void *p_to = ((char *)last) + pt_size;
1531 for (i = 1; i < n_points; i++)
1532 {
1533 int last_point = (i == n_points - 1);
1534
1535 /* Look straight into the abyss */
1536 pt = getPoint2d_cp(pa, i);
1537
1538 /* Don't drop points if we are running short of points */
1539 if (n_points + n_points_out > min_points + i)
1540 {
1541 if (tolerance > 0.0)
1542 {
1543 /* Only drop points that are within our tolerance */
1544 dsq = distance2d_sqr_pt_pt(last, pt);
1545 /* Allow any point but the last one to be dropped */
1546 if (!last_point && dsq <= tolsq)
1547 {
1548 continue;
1549 }
1550 }
1551 else
1552 {
1553 /* At tolerance zero, only skip exact dupes */
1554 if (memcmp((char*)pt, (char*)last, pt_size) == 0)
1555 continue;
1556 }
1557
1558 /* Got to last point, and it's not very different from */
1559 /* the point that preceded it. We want to keep the last */
1560 /* point, not the second-to-last one, so we pull our write */
1561 /* index back one value */
1562 if (last_point && n_points_out > 1 && tolerance > 0.0 && dsq <= tolsq)
1563 {
1564 n_points_out--;
1565 p_to = (char*)p_to - pt_size;
1566 }
1567 }
1568
1569 /* Compact all remaining values to front of array */
1570 memcpy(p_to, pt, pt_size);
1571 n_points_out++;
1572 p_to = (char*)p_to + pt_size;
1573 last = pt;
1574 }
1575 /* Adjust array length */
1576 pa->npoints = n_points_out;
1577 return;
1578 }
1579
1580 /* Out of the points in pa [itfist .. itlast], finds the one that's farthest away from
1581 * the segment determined by pts[itfist] and pts[itlast].
1582 * Returns itfirst if no point was found futher away than max_distance_sqr
1583 */
1584 static uint32_t
ptarray_dp_findsplit_in_place(const POINTARRAY * pts,uint32_t it_first,uint32_t it_last,double max_distance_sqr)1585 ptarray_dp_findsplit_in_place(const POINTARRAY *pts, uint32_t it_first, uint32_t it_last, double max_distance_sqr)
1586 {
1587 uint32_t split = it_first;
1588 if ((it_first - it_last) < 2)
1589 return it_first;
1590
1591 const POINT2D *A = getPoint2d_cp(pts, it_first);
1592 const POINT2D *B = getPoint2d_cp(pts, it_last);
1593
1594 if (distance2d_sqr_pt_pt(A, B) < DBL_EPSILON)
1595 {
1596 /* If p1 == p2, we can just calculate the distance from each point to A */
1597 for (uint32_t itk = it_first + 1; itk < it_last; itk++)
1598 {
1599 const POINT2D *pk = getPoint2d_cp(pts, itk);
1600 double distance_sqr = distance2d_sqr_pt_pt(pk, A);
1601 if (distance_sqr > max_distance_sqr)
1602 {
1603 split = itk;
1604 max_distance_sqr = distance_sqr;
1605 }
1606 }
1607 return split;
1608 }
1609
1610 /* This is based on distance2d_sqr_pt_seg, but heavily inlined here to avoid recalculations */
1611 double ba_x = (B->x - A->x);
1612 double ba_y = (B->y - A->y);
1613 double ab_length_sqr = (ba_x * ba_x + ba_y * ba_y);
1614 /* To avoid the division by ab_length_sqr in the 3rd path, we normalize here
1615 * and multiply in the first two paths [(dot_ac_ab < 0) and (> ab_length_sqr)] */
1616 max_distance_sqr *= ab_length_sqr;
1617 for (uint32_t itk = it_first + 1; itk < it_last; itk++)
1618 {
1619 const POINT2D *C = getPoint2d_cp(pts, itk);
1620 double distance_sqr;
1621 double ca_x = (C->x - A->x);
1622 double ca_y = (C->y - A->y);
1623 double dot_ac_ab = (ca_x * ba_x + ca_y * ba_y);
1624
1625 if (dot_ac_ab <= 0.0)
1626 {
1627 distance_sqr = distance2d_sqr_pt_pt(C, A) * ab_length_sqr;
1628 }
1629 else if (dot_ac_ab >= ab_length_sqr)
1630 {
1631 distance_sqr = distance2d_sqr_pt_pt(C, B) * ab_length_sqr;
1632 }
1633 else
1634 {
1635 double s_numerator = ca_x * ba_y - ca_y * ba_x;
1636 distance_sqr = s_numerator * s_numerator; /* Missing division by ab_length_sqr on purpose */
1637 }
1638
1639 if (distance_sqr > max_distance_sqr)
1640 {
1641 split = itk;
1642 max_distance_sqr = distance_sqr;
1643 }
1644 }
1645 return split;
1646 }
1647
1648 /* O(N) simplification for tolearnce = 0 */
1649 static void
ptarray_simplify_in_place_tolerance0(POINTARRAY * pa)1650 ptarray_simplify_in_place_tolerance0(POINTARRAY *pa)
1651 {
1652 uint32_t kept_it = 0;
1653 uint32_t last_it = pa->npoints - 1;
1654 const POINT2D *kept_pt = getPoint2d_cp(pa, 0);
1655 const size_t pt_size = ptarray_point_size(pa);
1656
1657 for (uint32_t i = 1; i < last_it; i++)
1658 {
1659 const POINT2D *curr_pt = getPoint2d_cp(pa, i);
1660 const POINT2D *next_pt = getPoint2d_cp(pa, i + 1);
1661
1662 double ba_x = next_pt->x - kept_pt->x;
1663 double ba_y = next_pt->y - kept_pt->y;
1664 double ab_length_sqr = ba_x * ba_x + ba_y * ba_y;
1665
1666 double ca_x = curr_pt->x - kept_pt->x;
1667 double ca_y = curr_pt->y - kept_pt->y;
1668 double dot_ac_ab = ca_x * ba_x + ca_y * ba_y;
1669 double s_numerator = ca_x * ba_y - ca_y * ba_x;
1670
1671 if (dot_ac_ab < 0.0 || dot_ac_ab > ab_length_sqr || s_numerator != 0)
1672 {
1673 kept_it++;
1674 kept_pt = curr_pt;
1675 if (kept_it != i)
1676 memcpy(pa->serialized_pointlist + pt_size * kept_it,
1677 pa->serialized_pointlist + pt_size * i,
1678 pt_size);
1679 }
1680 }
1681
1682 /* Append last point */
1683 kept_it++;
1684 if (kept_it != last_it)
1685 memcpy(pa->serialized_pointlist + pt_size * kept_it,
1686 pa->serialized_pointlist + pt_size * last_it,
1687 pt_size);
1688 pa->npoints = kept_it + 1;
1689 }
1690
1691 void
ptarray_simplify_in_place(POINTARRAY * pa,double tolerance,uint32_t minpts)1692 ptarray_simplify_in_place(POINTARRAY *pa, double tolerance, uint32_t minpts)
1693 {
1694 /* Do not try to simplify really short things */
1695 if (pa->npoints < 3 || pa->npoints <= minpts)
1696 return;
1697
1698 if (tolerance == 0 && minpts <= 2)
1699 {
1700 ptarray_simplify_in_place_tolerance0(pa);
1701 return;
1702 }
1703
1704 /* We use this array to keep track of the points we are keeping, so
1705 * we store just TRUE / FALSE in their position */
1706 uint8_t *kept_points = lwalloc(sizeof(uint8_t) * pa->npoints);
1707 memset(kept_points, LW_FALSE, sizeof(uint8_t) * pa->npoints);
1708 kept_points[0] = LW_TRUE;
1709 kept_points[pa->npoints - 1] = LW_TRUE;
1710 uint32_t keptn = 2;
1711
1712 /* We use this array as a stack to store the iterators that we are going to need
1713 * in the following steps.
1714 * This is ~10% faster than iterating over @kept_points looking for them
1715 */
1716 uint32_t *iterator_stack = lwalloc(sizeof(uint32_t) * pa->npoints);
1717 iterator_stack[0] = 0;
1718 uint32_t iterator_stack_size = 1;
1719
1720 uint32_t it_first = 0;
1721 uint32_t it_last = pa->npoints - 1;
1722
1723 const double tolerance_sqr = tolerance * tolerance;
1724 /* For the first @minpts points we ignore the tolerance */
1725 double it_tol = keptn >= minpts ? tolerance_sqr : -1.0;
1726
1727 while (iterator_stack_size)
1728 {
1729 uint32_t split = ptarray_dp_findsplit_in_place(pa, it_first, it_last, it_tol);
1730 if (split == it_first)
1731 {
1732 it_first = it_last;
1733 it_last = iterator_stack[--iterator_stack_size];
1734 }
1735 else
1736 {
1737 kept_points[split] = LW_TRUE;
1738 keptn++;
1739
1740 iterator_stack[iterator_stack_size++] = it_last;
1741 it_last = split;
1742 it_tol = keptn >= minpts ? tolerance_sqr : -1.0;
1743 }
1744 }
1745
1746 const size_t pt_size = ptarray_point_size(pa);
1747 /* The first point is already in place, so we don't need to copy it */
1748 size_t kept_it = 1;
1749 if (keptn == 2)
1750 {
1751 /* If there are 2 points remaining, it has to be first and last as
1752 * we added those at the start */
1753 memcpy(pa->serialized_pointlist + pt_size * kept_it,
1754 pa->serialized_pointlist + pt_size * (pa->npoints - 1),
1755 pt_size);
1756 }
1757 else if (pa->npoints != keptn) /* We don't need to move any points if we are keeping them all */
1758 {
1759 for (uint32_t i = 1; i < pa->npoints; i++)
1760 {
1761 if (kept_points[i])
1762 {
1763 memcpy(pa->serialized_pointlist + pt_size * kept_it,
1764 pa->serialized_pointlist + pt_size * i,
1765 pt_size);
1766 kept_it++;
1767 }
1768 }
1769 }
1770 pa->npoints = keptn;
1771
1772 lwfree(kept_points);
1773 lwfree(iterator_stack);
1774 }
1775
1776 /************************************************************************/
1777
1778 /**
1779 * Find the 2d length of the given #POINTARRAY, using circular
1780 * arc interpolation between each coordinate triple.
1781 * Length(A1, A2, A3, A4, A5) = Length(A1, A2, A3)+Length(A3, A4, A5)
1782 */
1783 double
ptarray_arc_length_2d(const POINTARRAY * pts)1784 ptarray_arc_length_2d(const POINTARRAY *pts)
1785 {
1786 double dist = 0.0;
1787 uint32_t i;
1788 const POINT2D *a1;
1789 const POINT2D *a2;
1790 const POINT2D *a3;
1791
1792 if ( pts->npoints % 2 != 1 )
1793 lwerror("arc point array with even number of points");
1794
1795 a1 = getPoint2d_cp(pts, 0);
1796
1797 for ( i=2; i < pts->npoints; i += 2 )
1798 {
1799 a2 = getPoint2d_cp(pts, i-1);
1800 a3 = getPoint2d_cp(pts, i);
1801 dist += lw_arc_length(a1, a2, a3);
1802 a1 = a3;
1803 }
1804 return dist;
1805 }
1806
1807 /**
1808 * Find the 2d length of the given #POINTARRAY (even if it's 3d)
1809 */
1810 double
ptarray_length_2d(const POINTARRAY * pts)1811 ptarray_length_2d(const POINTARRAY *pts)
1812 {
1813 double dist = 0.0;
1814 uint32_t i;
1815 const POINT2D *frm;
1816 const POINT2D *to;
1817
1818 if ( pts->npoints < 2 ) return 0.0;
1819
1820 frm = getPoint2d_cp(pts, 0);
1821
1822 for ( i=1; i < pts->npoints; i++ )
1823 {
1824 to = getPoint2d_cp(pts, i);
1825
1826 dist += sqrt( ((frm->x - to->x)*(frm->x - to->x)) +
1827 ((frm->y - to->y)*(frm->y - to->y)) );
1828
1829 frm = to;
1830 }
1831 return dist;
1832 }
1833
1834 /**
1835 * Find the 3d/2d length of the given #POINTARRAY
1836 * (depending on its dimensionality)
1837 */
1838 double
ptarray_length(const POINTARRAY * pts)1839 ptarray_length(const POINTARRAY *pts)
1840 {
1841 double dist = 0.0;
1842 uint32_t i;
1843 POINT3DZ frm;
1844 POINT3DZ to;
1845
1846 if ( pts->npoints < 2 ) return 0.0;
1847
1848 /* compute 2d length if 3d is not available */
1849 if ( ! FLAGS_GET_Z(pts->flags) ) return ptarray_length_2d(pts);
1850
1851 getPoint3dz_p(pts, 0, &frm);
1852 for ( i=1; i < pts->npoints; i++ )
1853 {
1854 getPoint3dz_p(pts, i, &to);
1855 dist += sqrt( ((frm.x - to.x)*(frm.x - to.x)) +
1856 ((frm.y - to.y)*(frm.y - to.y)) +
1857 ((frm.z - to.z)*(frm.z - to.z)) );
1858 frm = to;
1859 }
1860 return dist;
1861 }
1862
1863
1864
1865 /**
1866 * Affine transform a pointarray.
1867 */
1868 void
ptarray_affine(POINTARRAY * pa,const AFFINE * a)1869 ptarray_affine(POINTARRAY *pa, const AFFINE *a)
1870 {
1871 if (FLAGS_GET_Z(pa->flags))
1872 {
1873 for (uint32_t i = 0; i < pa->npoints; i++)
1874 {
1875 POINT4D *p4d = (POINT4D *)(getPoint_internal(pa, i));
1876 double x = p4d->x;
1877 double y = p4d->y;
1878 double z = p4d->z;
1879 p4d->x = a->afac * x + a->bfac * y + a->cfac * z + a->xoff;
1880 p4d->y = a->dfac * x + a->efac * y + a->ffac * z + a->yoff;
1881 p4d->z = a->gfac * x + a->hfac * y + a->ifac * z + a->zoff;
1882 }
1883 }
1884 else
1885 {
1886 for (uint32_t i = 0; i < pa->npoints; i++)
1887 {
1888 POINT2D *pt = (POINT2D *)(getPoint_internal(pa, i));
1889 double x = pt->x;
1890 double y = pt->y;
1891 pt->x = a->afac * x + a->bfac * y + a->xoff;
1892 pt->y = a->dfac * x + a->efac * y + a->yoff;
1893 }
1894 }
1895 }
1896
1897 /**
1898 * WARNING, make sure you send in only 16-member double arrays
1899 * or obviously things will go pear-shaped fast.
1900 */
1901 #if 0
1902 static int gluInvertMatrix(const double *m, double *invOut)
1903 {
1904 double inv[16], det;
1905 int i;
1906
1907 inv[0] = m[5] * m[10] * m[15] -
1908 m[5] * m[11] * m[14] -
1909 m[9] * m[6] * m[15] +
1910 m[9] * m[7] * m[14] +
1911 m[13] * m[6] * m[11] -
1912 m[13] * m[7] * m[10];
1913
1914 inv[4] = -m[4] * m[10] * m[15] +
1915 m[4] * m[11] * m[14] +
1916 m[8] * m[6] * m[15] -
1917 m[8] * m[7] * m[14] -
1918 m[12] * m[6] * m[11] +
1919 m[12] * m[7] * m[10];
1920
1921 inv[8] = m[4] * m[9] * m[15] -
1922 m[4] * m[11] * m[13] -
1923 m[8] * m[5] * m[15] +
1924 m[8] * m[7] * m[13] +
1925 m[12] * m[5] * m[11] -
1926 m[12] * m[7] * m[9];
1927
1928 inv[12] = -m[4] * m[9] * m[14] +
1929 m[4] * m[10] * m[13] +
1930 m[8] * m[5] * m[14] -
1931 m[8] * m[6] * m[13] -
1932 m[12] * m[5] * m[10] +
1933 m[12] * m[6] * m[9];
1934
1935 inv[1] = -m[1] * m[10] * m[15] +
1936 m[1] * m[11] * m[14] +
1937 m[9] * m[2] * m[15] -
1938 m[9] * m[3] * m[14] -
1939 m[13] * m[2] * m[11] +
1940 m[13] * m[3] * m[10];
1941
1942 inv[5] = m[0] * m[10] * m[15] -
1943 m[0] * m[11] * m[14] -
1944 m[8] * m[2] * m[15] +
1945 m[8] * m[3] * m[14] +
1946 m[12] * m[2] * m[11] -
1947 m[12] * m[3] * m[10];
1948
1949 inv[9] = -m[0] * m[9] * m[15] +
1950 m[0] * m[11] * m[13] +
1951 m[8] * m[1] * m[15] -
1952 m[8] * m[3] * m[13] -
1953 m[12] * m[1] * m[11] +
1954 m[12] * m[3] * m[9];
1955
1956 inv[13] = m[0] * m[9] * m[14] -
1957 m[0] * m[10] * m[13] -
1958 m[8] * m[1] * m[14] +
1959 m[8] * m[2] * m[13] +
1960 m[12] * m[1] * m[10] -
1961 m[12] * m[2] * m[9];
1962
1963 inv[2] = m[1] * m[6] * m[15] -
1964 m[1] * m[7] * m[14] -
1965 m[5] * m[2] * m[15] +
1966 m[5] * m[3] * m[14] +
1967 m[13] * m[2] * m[7] -
1968 m[13] * m[3] * m[6];
1969
1970 inv[6] = -m[0] * m[6] * m[15] +
1971 m[0] * m[7] * m[14] +
1972 m[4] * m[2] * m[15] -
1973 m[4] * m[3] * m[14] -
1974 m[12] * m[2] * m[7] +
1975 m[12] * m[3] * m[6];
1976
1977 inv[10] = m[0] * m[5] * m[15] -
1978 m[0] * m[7] * m[13] -
1979 m[4] * m[1] * m[15] +
1980 m[4] * m[3] * m[13] +
1981 m[12] * m[1] * m[7] -
1982 m[12] * m[3] * m[5];
1983
1984 inv[14] = -m[0] * m[5] * m[14] +
1985 m[0] * m[6] * m[13] +
1986 m[4] * m[1] * m[14] -
1987 m[4] * m[2] * m[13] -
1988 m[12] * m[1] * m[6] +
1989 m[12] * m[2] * m[5];
1990
1991 inv[3] = -m[1] * m[6] * m[11] +
1992 m[1] * m[7] * m[10] +
1993 m[5] * m[2] * m[11] -
1994 m[5] * m[3] * m[10] -
1995 m[9] * m[2] * m[7] +
1996 m[9] * m[3] * m[6];
1997
1998 inv[7] = m[0] * m[6] * m[11] -
1999 m[0] * m[7] * m[10] -
2000 m[4] * m[2] * m[11] +
2001 m[4] * m[3] * m[10] +
2002 m[8] * m[2] * m[7] -
2003 m[8] * m[3] * m[6];
2004
2005 inv[11] = -m[0] * m[5] * m[11] +
2006 m[0] * m[7] * m[9] +
2007 m[4] * m[1] * m[11] -
2008 m[4] * m[3] * m[9] -
2009 m[8] * m[1] * m[7] +
2010 m[8] * m[3] * m[5];
2011
2012 inv[15] = m[0] * m[5] * m[10] -
2013 m[0] * m[6] * m[9] -
2014 m[4] * m[1] * m[10] +
2015 m[4] * m[2] * m[9] +
2016 m[8] * m[1] * m[6] -
2017 m[8] * m[2] * m[5];
2018
2019 det = m[0] * inv[0] + m[1] * inv[4] + m[2] * inv[8] + m[3] * inv[12];
2020
2021 if (det == 0)
2022 return LW_FALSE;
2023
2024 det = 1.0 / det;
2025
2026 for (i = 0; i < 16; i++)
2027 invOut[i] = inv[i] * det;
2028
2029 return LW_TRUE;
2030 }
2031 #endif
2032
2033 /**
2034 * Scale a pointarray.
2035 */
2036 void
ptarray_scale(POINTARRAY * pa,const POINT4D * fact)2037 ptarray_scale(POINTARRAY *pa, const POINT4D *fact)
2038 {
2039 uint32_t i;
2040 POINT4D p4d;
2041 LWDEBUG(3, "ptarray_scale start");
2042 for (i=0; i<pa->npoints; i++)
2043 {
2044 getPoint4d_p(pa, i, &p4d);
2045 p4d.x *= fact->x;
2046 p4d.y *= fact->y;
2047 p4d.z *= fact->z;
2048 p4d.m *= fact->m;
2049 ptarray_set_point4d(pa, i, &p4d);
2050 }
2051 LWDEBUG(3, "ptarray_scale end");
2052 }
2053
2054 int
ptarray_startpoint(const POINTARRAY * pa,POINT4D * pt)2055 ptarray_startpoint(const POINTARRAY *pa, POINT4D *pt)
2056 {
2057 return getPoint4d_p(pa, 0, pt);
2058 }
2059
2060
2061 /*
2062 * Stick an array of points to the given gridspec.
2063 * Return "gridded" points in *outpts and their number in *outptsn.
2064 *
2065 * Two consecutive points falling on the same grid cell are collapsed
2066 * into one single point.
2067 *
2068 */
2069 void
ptarray_grid_in_place(POINTARRAY * pa,const gridspec * grid)2070 ptarray_grid_in_place(POINTARRAY *pa, const gridspec *grid)
2071 {
2072 uint32_t j = 0;
2073 POINT4D *p, *p_out = NULL;
2074 double x, y, z = 0, m = 0;
2075 uint32_t ndims = FLAGS_NDIMS(pa->flags);
2076 uint32_t has_z = FLAGS_GET_Z(pa->flags);
2077 uint32_t has_m = FLAGS_GET_M(pa->flags);
2078
2079 for (uint32_t i = 0; i < pa->npoints; i++)
2080 {
2081 /* Look straight into the abyss */
2082 p = (POINT4D *)(getPoint_internal(pa, i));
2083 x = p->x;
2084 y = p->y;
2085 if (ndims > 2)
2086 z = p->z;
2087 if (ndims > 3)
2088 m = p->m;
2089
2090 if (grid->xsize > 0)
2091 x = rint((x - grid->ipx) / grid->xsize) * grid->xsize + grid->ipx;
2092
2093 if (grid->ysize > 0)
2094 y = rint((y - grid->ipy) / grid->ysize) * grid->ysize + grid->ipy;
2095
2096 /* Read and round this point */
2097 /* Z is always in third position */
2098 if (has_z && grid->zsize > 0)
2099 z = rint((z - grid->ipz) / grid->zsize) * grid->zsize + grid->ipz;
2100
2101 /* M might be in 3rd or 4th position */
2102 if (has_m && grid->msize > 0)
2103 {
2104 /* In POINT ZM, M is in 4th position, in POINT M, M is in 3rd position which is Z in POINT4D */
2105 if (has_z)
2106 m = rint((m - grid->ipm) / grid->msize) * grid->msize + grid->ipm;
2107 else
2108 z = rint((z - grid->ipm) / grid->msize) * grid->msize + grid->ipm;
2109 }
2110
2111 /* Skip duplicates */
2112 if (p_out && p_out->x == x && p_out->y == y && (ndims > 2 ? p_out->z == z : 1) &&
2113 (ndims > 3 ? p_out->m == m : 1))
2114 continue;
2115
2116 /* Write rounded values into the next available point */
2117 p_out = (POINT4D *)(getPoint_internal(pa, j++));
2118 p_out->x = x;
2119 p_out->y = y;
2120 if (ndims > 2)
2121 p_out->z = z;
2122 if (ndims > 3)
2123 p_out->m = m;
2124 }
2125
2126 /* Update output ptarray length */
2127 pa->npoints = j;
2128 return;
2129 }
2130
2131 int
ptarray_npoints_in_rect(const POINTARRAY * pa,const GBOX * gbox)2132 ptarray_npoints_in_rect(const POINTARRAY *pa, const GBOX *gbox)
2133 {
2134 const POINT2D *pt;
2135 int n = 0;
2136 uint32_t i;
2137 for ( i = 0; i < pa->npoints; i++ )
2138 {
2139 pt = getPoint2d_cp(pa, i);
2140 if ( gbox_contains_point2d(gbox, pt) )
2141 n++;
2142 }
2143 return n;
2144 }
2145
2146
2147 /*
2148 * Reorder the vertices of a closed pointarray so that the
2149 * given point is the first/last one.
2150 *
2151 * Error out if pointarray is not closed or it does not
2152 * contain the given point.
2153 */
2154 int
ptarray_scroll_in_place(POINTARRAY * pa,const POINT4D * pt)2155 ptarray_scroll_in_place(POINTARRAY *pa, const POINT4D *pt)
2156 {
2157 POINTARRAY *tmp;
2158 int found;
2159 uint32_t it;
2160 int ptsize;
2161
2162 if ( ! ptarray_is_closed_2d(pa) )
2163 {
2164 lwerror("ptarray_scroll_in_place: input POINTARRAY is not closed");
2165 return LW_FAILURE;
2166 }
2167
2168 ptsize = ptarray_point_size(pa);
2169
2170 /* Find the point in the array */
2171 found = 0;
2172 for ( it = 0; it < pa->npoints; ++it )
2173 {
2174 if ( ! memcmp(getPoint_internal(pa, it), pt, ptsize) )
2175 {
2176 found = 1;
2177 break;
2178 }
2179 }
2180
2181 if ( ! found )
2182 {
2183 lwerror("ptarray_scroll_in_place: input POINTARRAY does not contain the given point");
2184 return LW_FAILURE;
2185 }
2186
2187 if ( 0 == it )
2188 {
2189 /* Point is already the start/end point, just clone the input */
2190 return LW_SUCCESS;
2191 }
2192
2193 /* TODO: reduce allocations */
2194 tmp = ptarray_construct(FLAGS_GET_Z(pa->flags), FLAGS_GET_M(pa->flags), pa->npoints);
2195
2196 bzero(getPoint_internal(tmp, 0), ptsize * pa->npoints);
2197 /* Copy the block from found point to last point into the output array */
2198 memcpy(
2199 getPoint_internal(tmp, 0),
2200 getPoint_internal(pa, it),
2201 ptsize * ( pa->npoints - it )
2202 );
2203
2204 /* Copy the block from second point to the found point into the last portion of the
2205 * return */
2206 memcpy(
2207 getPoint_internal(tmp, pa->npoints - it),
2208 getPoint_internal(pa, 1),
2209 ptsize * ( it )
2210 );
2211
2212 /* Copy the resulting pointarray back to source one */
2213 memcpy(
2214 getPoint_internal(pa, 0),
2215 getPoint_internal(tmp, 0),
2216 ptsize * ( pa->npoints )
2217 );
2218
2219 ptarray_free(tmp);
2220
2221 return LW_SUCCESS;
2222 }
2223