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