PostGIS  2.4.9dev-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 
30 #include "../postgis_config.h"
31 /*#define POSTGIS_DEBUG_LEVEL 4*/
32 #include "liblwgeom_internal.h"
33 #include "lwgeom_log.h"
34 
35 int
37 {
38  if ( ! pa ) return LW_FALSE;
39  return FLAGS_GET_Z(pa->flags);
40 }
41 
42 int
44 {
45  if ( ! pa ) return LW_FALSE;
46  return FLAGS_GET_M(pa->flags);
47 }
48 
49 /*
50  * Size of point represeneted in the POINTARRAY
51  * 16 for 2d, 24 for 3d, 32 for 4d
52  */
53 int inline
55 {
56  LWDEBUGF(5, "ptarray_point_size: FLAGS_NDIMS(pa->flags)=%x",FLAGS_NDIMS(pa->flags));
57 
58  return sizeof(double)*FLAGS_NDIMS(pa->flags);
59 }
60 
62 ptarray_construct(char hasz, char hasm, uint32_t npoints)
63 {
64  POINTARRAY *pa = ptarray_construct_empty(hasz, hasm, npoints);
65  pa->npoints = npoints;
66  return pa;
67 }
68 
70 ptarray_construct_empty(char hasz, char hasm, uint32_t maxpoints)
71 {
72  POINTARRAY *pa = lwalloc(sizeof(POINTARRAY));
73  pa->serialized_pointlist = NULL;
74 
75  /* Set our dimsionality info on the bitmap */
76  pa->flags = gflags(hasz, hasm, 0);
77 
78  /* We will be allocating a bit of room */
79  pa->npoints = 0;
80  pa->maxpoints = maxpoints;
81 
82  /* Allocate the coordinate array */
83  if ( maxpoints > 0 )
84  pa->serialized_pointlist = lwalloc(maxpoints * ptarray_point_size(pa));
85  else
86  pa->serialized_pointlist = NULL;
87 
88  return pa;
89 }
90 
91 /*
92 * Add a point into a pointarray. Only adds as many dimensions as the
93 * pointarray supports.
94 */
95 int
96 ptarray_insert_point(POINTARRAY *pa, const POINT4D *p, int where)
97 {
98  size_t point_size = ptarray_point_size(pa);
99  LWDEBUGF(5,"pa = %p; p = %p; where = %d", pa, p, where);
100  LWDEBUGF(5,"pa->npoints = %d; pa->maxpoints = %d", pa->npoints, pa->maxpoints);
101 
102  if ( FLAGS_GET_READONLY(pa->flags) )
103  {
104  lwerror("ptarray_insert_point: called on read-only point array");
105  return LW_FAILURE;
106  }
107 
108  /* Error on invalid offset value */
109  if ( where > pa->npoints || where < 0)
110  {
111  lwerror("ptarray_insert_point: offset out of range (%d)", where);
112  return LW_FAILURE;
113  }
114 
115  /* If we have no storage, let's allocate some */
116  if( pa->maxpoints == 0 || ! pa->serialized_pointlist )
117  {
118  pa->maxpoints = 32;
119  pa->npoints = 0;
121  }
122 
123  /* Error out if we have a bad situation */
124  if ( pa->npoints > pa->maxpoints )
125  {
126  lwerror("npoints (%d) is greated than maxpoints (%d)", pa->npoints, pa->maxpoints);
127  return LW_FAILURE;
128  }
129 
130  /* Check if we have enough storage, add more if necessary */
131  if( pa->npoints == pa->maxpoints )
132  {
133  pa->maxpoints *= 2;
135  }
136 
137  /* Make space to insert the new point */
138  if( where < pa->npoints )
139  {
140  size_t copy_size = point_size * (pa->npoints - where);
141  memmove(getPoint_internal(pa, where+1), getPoint_internal(pa, where), copy_size);
142  LWDEBUGF(5,"copying %d bytes to start vertex %d from start vertex %d", copy_size, where+1, where);
143  }
144 
145  /* We have one more point */
146  ++pa->npoints;
147 
148  /* Copy the new point into the gap */
149  ptarray_set_point4d(pa, where, p);
150  LWDEBUGF(5,"copying new point to start vertex %d", point_size, where);
151 
152  return LW_SUCCESS;
153 }
154 
155 int
156 ptarray_append_point(POINTARRAY *pa, const POINT4D *pt, int repeated_points)
157 {
158 
159  /* Check for pathology */
160  if( ! pa || ! pt )
161  {
162  lwerror("ptarray_append_point: null input");
163  return LW_FAILURE;
164  }
165 
166  /* Check for duplicate end point */
167  if ( repeated_points == LW_FALSE && pa->npoints > 0 )
168  {
169  POINT4D tmp;
170  getPoint4d_p(pa, pa->npoints-1, &tmp);
171  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);
172 
173  /* Return LW_SUCCESS and do nothing else if previous point in list is equal to this one */
174  if ( (pt->x == tmp.x) && (pt->y == tmp.y) &&
175  (FLAGS_GET_Z(pa->flags) ? pt->z == tmp.z : 1) &&
176  (FLAGS_GET_M(pa->flags) ? pt->m == tmp.m : 1) )
177  {
178  return LW_SUCCESS;
179  }
180  }
181 
182  /* Append is just a special case of insert */
183  return ptarray_insert_point(pa, pt, pa->npoints);
184 }
185 
186 int
187 ptarray_append_ptarray(POINTARRAY *pa1, POINTARRAY *pa2, double gap_tolerance)
188 {
189  unsigned int poff = 0;
190  unsigned int npoints;
191  unsigned int ncap;
192  unsigned int ptsize;
193 
194  /* Check for pathology */
195  if( ! pa1 || ! pa2 )
196  {
197  lwerror("ptarray_append_ptarray: null input");
198  return LW_FAILURE;
199  }
200 
201  npoints = pa2->npoints;
202 
203  if ( ! npoints ) return LW_SUCCESS; /* nothing more to do */
204 
205  if( FLAGS_GET_READONLY(pa1->flags) )
206  {
207  lwerror("ptarray_append_ptarray: target pointarray is read-only");
208  return LW_FAILURE;
209  }
210 
211  if( FLAGS_GET_ZM(pa1->flags) != FLAGS_GET_ZM(pa2->flags) )
212  {
213  lwerror("ptarray_append_ptarray: appending mixed dimensionality is not allowed");
214  return LW_FAILURE;
215  }
216 
217  ptsize = ptarray_point_size(pa1);
218 
219  /* Check for duplicate end point */
220  if ( pa1->npoints )
221  {
222  POINT2D tmp1, tmp2;
223  getPoint2d_p(pa1, pa1->npoints-1, &tmp1);
224  getPoint2d_p(pa2, 0, &tmp2);
225 
226  /* If the end point and start point are the same, then don't copy start point */
227  if (p2d_same(&tmp1, &tmp2)) {
228  poff = 1;
229  --npoints;
230  }
231  else if ( gap_tolerance == 0 || ( gap_tolerance > 0 &&
232  distance2d_pt_pt(&tmp1, &tmp2) > gap_tolerance ) )
233  {
234  lwerror("Second line start point too far from first line end point");
235  return LW_FAILURE;
236  }
237  }
238 
239  /* Check if we need extra space */
240  ncap = pa1->npoints + npoints;
241  if ( pa1->maxpoints < ncap )
242  {
243  pa1->maxpoints = ncap > pa1->maxpoints*2 ?
244  ncap : pa1->maxpoints*2;
245  pa1->serialized_pointlist = lwrealloc(pa1->serialized_pointlist, ptsize * pa1->maxpoints);
246  }
247 
248  memcpy(getPoint_internal(pa1, pa1->npoints),
249  getPoint_internal(pa2, poff), ptsize * npoints);
250 
251  pa1->npoints = ncap;
252 
253  return LW_SUCCESS;
254 }
255 
256 /*
257 * Add a point into a pointarray. Only adds as many dimensions as the
258 * pointarray supports.
259 */
260 int
262 {
263  size_t ptsize = ptarray_point_size(pa);
264 
265  /* Check for pathology */
266  if( ! pa )
267  {
268  lwerror("ptarray_remove_point: null input");
269  return LW_FAILURE;
270  }
271 
272  /* Error on invalid offset value */
273  if ( where >= pa->npoints || where < 0)
274  {
275  lwerror("ptarray_remove_point: offset out of range (%d)", where);
276  return LW_FAILURE;
277  }
278 
279  /* If the point is any but the last, we need to copy the data back one point */
280  if( where < pa->npoints - 1 )
281  {
282  memmove(getPoint_internal(pa, where), getPoint_internal(pa, where+1), ptsize * (pa->npoints - where - 1));
283  }
284 
285  /* We have one less point */
286  pa->npoints--;
287 
288  return LW_SUCCESS;
289 }
290 
295 POINTARRAY* ptarray_construct_reference_data(char hasz, char hasm, uint32_t npoints, uint8_t *ptlist)
296 {
297  POINTARRAY *pa = lwalloc(sizeof(POINTARRAY));
298  LWDEBUGF(5, "hasz = %d, hasm = %d, npoints = %d, ptlist = %p", hasz, hasm, npoints, ptlist);
299  pa->flags = gflags(hasz, hasm, 0);
300  FLAGS_SET_READONLY(pa->flags, 1); /* We don't own this memory, so we can't alter or free it. */
301  pa->npoints = npoints;
302  pa->maxpoints = npoints;
303  pa->serialized_pointlist = ptlist;
304  return pa;
305 }
306 
307 
308 POINTARRAY*
309 ptarray_construct_copy_data(char hasz, char hasm, uint32_t npoints, const uint8_t *ptlist)
310 {
311  POINTARRAY *pa = lwalloc(sizeof(POINTARRAY));
312 
313  pa->flags = gflags(hasz, hasm, 0);
314  pa->npoints = npoints;
315  pa->maxpoints = npoints;
316 
317  if ( npoints > 0 )
318  {
319  pa->serialized_pointlist = lwalloc(ptarray_point_size(pa) * npoints);
320  memcpy(pa->serialized_pointlist, ptlist, ptarray_point_size(pa) * npoints);
321  }
322  else
323  {
324  pa->serialized_pointlist = NULL;
325  }
326 
327  return pa;
328 }
329 
331 {
332  if(pa)
333  {
334  if(pa->serialized_pointlist && ( ! FLAGS_GET_READONLY(pa->flags) ) )
336  lwfree(pa);
337  LWDEBUG(5,"Freeing a PointArray");
338  }
339 }
340 
341 
342 void
344 {
345  /* TODO change this to double array operations once point array is double aligned */
346  POINT4D pbuf;
347  uint32_t i;
348  int ptsize = ptarray_point_size(pa);
349  int last = pa->npoints-1;
350  int mid = pa->npoints/2;
351 
352  for (i=0; i<mid; i++)
353  {
354  uint8_t *from, *to;
355  from = getPoint_internal(pa, i);
356  to = getPoint_internal(pa, (last-i));
357  memcpy((uint8_t *)&pbuf, to, ptsize);
358  memcpy(to, from, ptsize);
359  memcpy(from, (uint8_t *)&pbuf, ptsize);
360  }
361 
362 }
363 
364 
368 POINTARRAY*
370 {
371  int i;
372  double d;
373  POINT4D p;
374 
375  for (i=0 ; i < pa->npoints ; i++)
376  {
377  getPoint4d_p(pa, i, &p);
378  d = p.y;
379  p.y = p.x;
380  p.x = d;
381  ptarray_set_point4d(pa, i, &p);
382  }
383 
384  return pa;
385 }
386 
387 void
389 {
390  int i;
391  double d, *dp1, *dp2;
392  POINT4D p;
393 
394  dp1 = ((double*)&p)+(unsigned)o1;
395  dp2 = ((double*)&p)+(unsigned)o2;
396  for (i=0 ; i < pa->npoints ; i++)
397  {
398  getPoint4d_p(pa, i, &p);
399  d = *dp2;
400  *dp2 = *dp1;
401  *dp1 = d;
402  ptarray_set_point4d(pa, i, &p);
403  }
404 }
405 
406 
414 POINTARRAY *
415 ptarray_segmentize2d(const POINTARRAY *ipa, double dist)
416 {
417  double segdist;
418  POINT4D p1, p2;
419  POINT4D pbuf;
420  POINTARRAY *opa;
421  int ipoff=0; /* input point offset */
422  int hasz = FLAGS_GET_Z(ipa->flags);
423  int hasm = FLAGS_GET_M(ipa->flags);
424 
425  pbuf.x = pbuf.y = pbuf.z = pbuf.m = 0;
426 
427  /* Initial storage */
428  opa = ptarray_construct_empty(hasz, hasm, ipa->npoints);
429 
430  /* Add first point */
431  getPoint4d_p(ipa, ipoff, &p1);
432  ptarray_append_point(opa, &p1, LW_FALSE);
433 
434  ipoff++;
435 
436  while (ipoff<ipa->npoints)
437  {
438  /*
439  * We use these pointers to avoid
440  * "strict-aliasing rules break" warning raised
441  * by gcc (3.3 and up).
442  *
443  * It looks that casting a variable address (also
444  * referred to as "type-punned pointer")
445  * breaks those "strict" rules.
446  *
447  */
448  POINT4D *p1ptr=&p1, *p2ptr=&p2;
449 
450  getPoint4d_p(ipa, ipoff, &p2);
451 
452  segdist = distance2d_pt_pt((POINT2D *)p1ptr, (POINT2D *)p2ptr);
453 
454  if (segdist > dist) /* add an intermediate point */
455  {
456  pbuf.x = p1.x + (p2.x-p1.x)/segdist * dist;
457  pbuf.y = p1.y + (p2.y-p1.y)/segdist * dist;
458  if( hasz )
459  pbuf.z = p1.z + (p2.z-p1.z)/segdist * dist;
460  if( hasm )
461  pbuf.m = p1.m + (p2.m-p1.m)/segdist * dist;
462  ptarray_append_point(opa, &pbuf, LW_FALSE);
463  p1 = pbuf;
464  }
465  else /* copy second point */
466  {
467  ptarray_append_point(opa, &p2, (ipa->npoints==2)?LW_TRUE:LW_FALSE);
468  p1 = p2;
469  ipoff++;
470  }
471 
472  LW_ON_INTERRUPT(ptarray_free(opa); return NULL);
473  }
474 
475  return opa;
476 }
477 
478 char
479 ptarray_same(const POINTARRAY *pa1, const POINTARRAY *pa2)
480 {
481  uint32_t i;
482  size_t ptsize;
483 
484  if ( FLAGS_GET_ZM(pa1->flags) != FLAGS_GET_ZM(pa2->flags) ) return LW_FALSE;
485  LWDEBUG(5,"dimensions are the same");
486 
487  if ( pa1->npoints != pa2->npoints ) return LW_FALSE;
488  LWDEBUG(5,"npoints are the same");
489 
490  ptsize = ptarray_point_size(pa1);
491  LWDEBUGF(5, "ptsize = %d", ptsize);
492 
493  for (i=0; i<pa1->npoints; i++)
494  {
495  if ( memcmp(getPoint_internal(pa1, i), getPoint_internal(pa2, i), ptsize) )
496  return LW_FALSE;
497  LWDEBUGF(5,"point #%d is the same",i);
498  }
499 
500  return LW_TRUE;
501 }
502 
503 POINTARRAY *
504 ptarray_addPoint(const POINTARRAY *pa, uint8_t *p, size_t pdims, uint32_t where)
505 {
506  POINTARRAY *ret;
507  POINT4D pbuf;
508  size_t ptsize = ptarray_point_size(pa);
509 
510  LWDEBUGF(3, "pa %x p %x size %d where %d",
511  pa, p, pdims, where);
512 
513  if ( pdims < 2 || pdims > 4 )
514  {
515  lwerror("ptarray_addPoint: point dimension out of range (%d)",
516  pdims);
517  return NULL;
518  }
519 
520  if ( where > pa->npoints )
521  {
522  lwerror("ptarray_addPoint: offset out of range (%d)",
523  where);
524  return NULL;
525  }
526 
527  LWDEBUG(3, "called with a %dD point");
528 
529  pbuf.x = pbuf.y = pbuf.z = pbuf.m = 0.0;
530  memcpy((uint8_t *)&pbuf, p, pdims*sizeof(double));
531 
532  LWDEBUG(3, "initialized point buffer");
533 
535  FLAGS_GET_M(pa->flags), pa->npoints+1);
536 
537  if ( where == -1 ) where = pa->npoints;
538 
539  if ( where )
540  {
541  memcpy(getPoint_internal(ret, 0), getPoint_internal(pa, 0), ptsize*where);
542  }
543 
544  memcpy(getPoint_internal(ret, where), (uint8_t *)&pbuf, ptsize);
545 
546  if ( where+1 != ret->npoints )
547  {
548  memcpy(getPoint_internal(ret, where+1),
549  getPoint_internal(pa, where),
550  ptsize*(pa->npoints-where));
551  }
552 
553  return ret;
554 }
555 
556 POINTARRAY *
558 {
559  POINTARRAY *ret;
560  size_t ptsize = ptarray_point_size(pa);
561 
562  LWDEBUGF(3, "pa %x which %d", pa, which);
563 
564 #if PARANOIA_LEVEL > 0
565  if ( which > pa->npoints-1 )
566  {
567  lwerror("%s [%d] offset (%d) out of range (%d..%d)", __FILE__, __LINE__,
568  which, 0, pa->npoints-1);
569  return NULL;
570  }
571 
572  if ( pa->npoints < 3 )
573  {
574  lwerror("%s [%d] can't remove a point from a 2-vertex POINTARRAY", __FILE__, __LINE__);
575  return NULL;
576  }
577 #endif
578 
580  FLAGS_GET_M(pa->flags), pa->npoints-1);
581 
582  /* copy initial part */
583  if ( which )
584  {
585  memcpy(getPoint_internal(ret, 0), getPoint_internal(pa, 0), ptsize*which);
586  }
587 
588  /* copy final part */
589  if ( which < pa->npoints-1 )
590  {
591  memcpy(getPoint_internal(ret, which), getPoint_internal(pa, which+1),
592  ptsize*(pa->npoints-which-1));
593  }
594 
595  return ret;
596 }
597 
598 POINTARRAY *
600 {
601  POINTARRAY *pa;
602  size_t ptsize = ptarray_point_size(pa1);
603 
604  if (FLAGS_GET_ZM(pa1->flags) != FLAGS_GET_ZM(pa2->flags))
605  lwerror("ptarray_cat: Mixed dimension");
606 
607  pa = ptarray_construct( FLAGS_GET_Z(pa1->flags),
608  FLAGS_GET_M(pa1->flags),
609  pa1->npoints + pa2->npoints);
610 
611  memcpy( getPoint_internal(pa, 0),
612  getPoint_internal(pa1, 0),
613  ptsize*(pa1->npoints));
614 
615  memcpy( getPoint_internal(pa, pa1->npoints),
616  getPoint_internal(pa2, 0),
617  ptsize*(pa2->npoints));
618 
619  ptarray_free(pa1);
620  ptarray_free(pa2);
621 
622  return pa;
623 }
624 
625 
629 POINTARRAY *
631 {
632  POINTARRAY *out = lwalloc(sizeof(POINTARRAY));
633 
634  LWDEBUG(3, "ptarray_clone_deep called.");
635 
636  out->flags = in->flags;
637  out->npoints = in->npoints;
638  out->maxpoints = in->npoints;
639 
640  FLAGS_SET_READONLY(out->flags, 0);
641 
642  if (!in->npoints)
643  {
644  // Avoid calling lwalloc of 0 bytes
645  out->serialized_pointlist = NULL;
646  }
647  else
648  {
649  size_t size = in->npoints * ptarray_point_size(in);
650  out->serialized_pointlist = lwalloc(size);
651  memcpy(out->serialized_pointlist, in->serialized_pointlist, size);
652  }
653 
654  return out;
655 }
656 
660 POINTARRAY *
662 {
663  POINTARRAY *out = lwalloc(sizeof(POINTARRAY));
664 
665  LWDEBUG(3, "ptarray_clone_deep called.");
666 
667  out->flags = in->flags;
668  out->npoints = in->npoints;
669  out->maxpoints = in->maxpoints;
670 
671  FLAGS_SET_READONLY(out->flags, 1);
672 
674 
675  return out;
676 }
677 
682 int
684 {
685  if (!in)
686  {
687  lwerror("ptarray_is_closed: called with null point array");
688  return 0;
689  }
690  if (in->npoints <= 1 ) return in->npoints; /* single-point are closed, empty not closed */
691 
692  return 0 == memcmp(getPoint_internal(in, 0), getPoint_internal(in, in->npoints-1), ptarray_point_size(in));
693 }
694 
695 
696 int
698 {
699  if (!in)
700  {
701  lwerror("ptarray_is_closed_2d: called with null point array");
702  return 0;
703  }
704  if (in->npoints <= 1 ) return in->npoints; /* single-point are closed, empty not closed */
705 
706  return 0 == memcmp(getPoint_internal(in, 0), getPoint_internal(in, in->npoints-1), sizeof(POINT2D) );
707 }
708 
709 int
711 {
712  if (!in)
713  {
714  lwerror("ptarray_is_closed_3d: called with null point array");
715  return 0;
716  }
717  if (in->npoints <= 1 ) return in->npoints; /* single-point are closed, empty not closed */
718 
719  return 0 == memcmp(getPoint_internal(in, 0), getPoint_internal(in, in->npoints-1), sizeof(POINT3D) );
720 }
721 
722 int
724 {
725  if ( FLAGS_GET_Z(in->flags) )
726  return ptarray_is_closed_3d(in);
727  else
728  return ptarray_is_closed_2d(in);
729 }
730 
735 int
737 {
738  return ptarray_contains_point_partial(pa, pt, LW_TRUE, NULL);
739 }
740 
741 int
742 ptarray_contains_point_partial(const POINTARRAY *pa, const POINT2D *pt, int check_closed, int *winding_number)
743 {
744  int wn = 0;
745  int i;
746  double side;
747  const POINT2D *seg1;
748  const POINT2D *seg2;
749  double ymin, ymax;
750 
751  seg1 = getPoint2d_cp(pa, 0);
752  seg2 = getPoint2d_cp(pa, pa->npoints-1);
753  if ( check_closed && ! p2d_same(seg1, seg2) )
754  lwerror("ptarray_contains_point called on unclosed ring");
755 
756  for ( i=1; i < pa->npoints; i++ )
757  {
758  seg2 = getPoint2d_cp(pa, i);
759 
760  /* Zero length segments are ignored. */
761  if ( seg1->x == seg2->x && seg1->y == seg2->y )
762  {
763  seg1 = seg2;
764  continue;
765  }
766 
767  ymin = FP_MIN(seg1->y, seg2->y);
768  ymax = FP_MAX(seg1->y, seg2->y);
769 
770  /* Only test segments in our vertical range */
771  if ( pt->y > ymax || pt->y < ymin )
772  {
773  seg1 = seg2;
774  continue;
775  }
776 
777  side = lw_segment_side(seg1, seg2, pt);
778 
779  /*
780  * A point on the boundary of a ring is not contained.
781  * WAS: if (fabs(side) < 1e-12), see #852
782  */
783  if ( (side == 0) && lw_pt_in_seg(pt, seg1, seg2) )
784  {
785  return LW_BOUNDARY;
786  }
787 
788  /*
789  * If the point is to the left of the line, and it's rising,
790  * then the line is to the right of the point and
791  * circling counter-clockwise, so incremement.
792  */
793  if ( (side < 0) && (seg1->y <= pt->y) && (pt->y < seg2->y) )
794  {
795  wn++;
796  }
797 
798  /*
799  * If the point is to the right of the line, and it's falling,
800  * then the line is to the right of the point and circling
801  * clockwise, so decrement.
802  */
803  else if ( (side > 0) && (seg2->y <= pt->y) && (pt->y < seg1->y) )
804  {
805  wn--;
806  }
807 
808  seg1 = seg2;
809  }
810 
811  /* Sent out the winding number for calls that are building on this as a primitive */
812  if ( winding_number )
813  *winding_number = wn;
814 
815  /* Outside */
816  if (wn == 0)
817  {
818  return LW_OUTSIDE;
819  }
820 
821  /* Inside */
822  return LW_INSIDE;
823 }
824 
834 int
836 {
837  return ptarrayarc_contains_point_partial(pa, pt, LW_TRUE /* Check closed*/, NULL);
838 }
839 
840 int
841 ptarrayarc_contains_point_partial(const POINTARRAY *pa, const POINT2D *pt, int check_closed, int *winding_number)
842 {
843  int wn = 0;
844  int i, side;
845  const POINT2D *seg1;
846  const POINT2D *seg2;
847  const POINT2D *seg3;
848  GBOX gbox;
849 
850  /* Check for not an arc ring (always have odd # of points) */
851  if ( (pa->npoints % 2) == 0 )
852  {
853  lwerror("ptarrayarc_contains_point called with even number of points");
854  return LW_OUTSIDE;
855  }
856 
857  /* Check for not an arc ring (always have >= 3 points) */
858  if ( pa->npoints < 3 )
859  {
860  lwerror("ptarrayarc_contains_point called too-short pointarray");
861  return LW_OUTSIDE;
862  }
863 
864  /* Check for unclosed case */
865  seg1 = getPoint2d_cp(pa, 0);
866  seg3 = getPoint2d_cp(pa, pa->npoints-1);
867  if ( check_closed && ! p2d_same(seg1, seg3) )
868  {
869  lwerror("ptarrayarc_contains_point called on unclosed ring");
870  return LW_OUTSIDE;
871  }
872  /* OK, it's closed. Is it just one circle? */
873  else if ( p2d_same(seg1, seg3) && pa->npoints == 3 )
874  {
875  double radius, d;
876  POINT2D c;
877  seg2 = getPoint2d_cp(pa, 1);
878 
879  /* Wait, it's just a point, so it can't contain anything */
880  if ( lw_arc_is_pt(seg1, seg2, seg3) )
881  return LW_OUTSIDE;
882 
883  /* See if the point is within the circle radius */
884  radius = lw_arc_center(seg1, seg2, seg3, &c);
885  d = distance2d_pt_pt(pt, &c);
886  if ( FP_EQUALS(d, radius) )
887  return LW_BOUNDARY; /* Boundary of circle */
888  else if ( d < radius )
889  return LW_INSIDE; /* Inside circle */
890  else
891  return LW_OUTSIDE; /* Outside circle */
892  }
893  else if ( p2d_same(seg1, pt) || p2d_same(seg3, pt) )
894  {
895  return LW_BOUNDARY; /* Boundary case */
896  }
897 
898  /* Start on the ring */
899  seg1 = getPoint2d_cp(pa, 0);
900  for ( i=1; i < pa->npoints; i += 2 )
901  {
902  seg2 = getPoint2d_cp(pa, i);
903  seg3 = getPoint2d_cp(pa, i+1);
904 
905  /* Catch an easy boundary case */
906  if( p2d_same(seg3, pt) )
907  return LW_BOUNDARY;
908 
909  /* Skip arcs that have no size */
910  if ( lw_arc_is_pt(seg1, seg2, seg3) )
911  {
912  seg1 = seg3;
913  continue;
914  }
915 
916  /* Only test segments in our vertical range */
917  lw_arc_calculate_gbox_cartesian_2d(seg1, seg2, seg3, &gbox);
918  if ( pt->y > gbox.ymax || pt->y < gbox.ymin )
919  {
920  seg1 = seg3;
921  continue;
922  }
923 
924  /* Outside of horizontal range, and not between end points we also skip */
925  if ( (pt->x > gbox.xmax || pt->x < gbox.xmin) &&
926  (pt->y > FP_MAX(seg1->y, seg3->y) || pt->y < FP_MIN(seg1->y, seg3->y)) )
927  {
928  seg1 = seg3;
929  continue;
930  }
931 
932  side = lw_arc_side(seg1, seg2, seg3, pt);
933 
934  /* On the boundary */
935  if ( (side == 0) && lw_pt_in_arc(pt, seg1, seg2, seg3) )
936  {
937  return LW_BOUNDARY;
938  }
939 
940  /* Going "up"! Point to left of arc. */
941  if ( side < 0 && (seg1->y <= pt->y) && (pt->y < seg3->y) )
942  {
943  wn++;
944  }
945 
946  /* Going "down"! */
947  if ( side > 0 && (seg2->y <= pt->y) && (pt->y < seg1->y) )
948  {
949  wn--;
950  }
951 
952  /* Inside the arc! */
953  if ( pt->x <= gbox.xmax && pt->x >= gbox.xmin )
954  {
955  POINT2D C;
956  double radius = lw_arc_center(seg1, seg2, seg3, &C);
957  double d = distance2d_pt_pt(pt, &C);
958 
959  /* On the boundary! */
960  if ( d == radius )
961  return LW_BOUNDARY;
962 
963  /* Within the arc! */
964  if ( d < radius )
965  {
966  /* Left side, increment winding number */
967  if ( side < 0 )
968  wn++;
969  /* Right side, decrement winding number */
970  if ( side > 0 )
971  wn--;
972  }
973  }
974 
975  seg1 = seg3;
976  }
977 
978  /* Sent out the winding number for calls that are building on this as a primitive */
979  if ( winding_number )
980  *winding_number = wn;
981 
982  /* Outside */
983  if (wn == 0)
984  {
985  return LW_OUTSIDE;
986  }
987 
988  /* Inside */
989  return LW_INSIDE;
990 }
991 
997 double
999 {
1000  const POINT2D *P1;
1001  const POINT2D *P2;
1002  const POINT2D *P3;
1003  double sum = 0.0;
1004  double x0, x, y1, y2;
1005  int i;
1006 
1007  if (! pa || pa->npoints < 3 )
1008  return 0.0;
1009 
1010  P1 = getPoint2d_cp(pa, 0);
1011  P2 = getPoint2d_cp(pa, 1);
1012  x0 = P1->x;
1013  for ( i = 2; i < pa->npoints; i++ )
1014  {
1015  P3 = getPoint2d_cp(pa, i);
1016  x = P2->x - x0;
1017  y1 = P3->y;
1018  y2 = P1->y;
1019  sum += x * (y2-y1);
1020 
1021  /* Move forwards! */
1022  P1 = P2;
1023  P2 = P3;
1024  }
1025  return sum / 2.0;
1026 }
1027 
1028 int
1030 {
1031  double area = 0;
1032  area = ptarray_signed_area(pa);
1033  if ( area > 0 ) return LW_FALSE;
1034  else return LW_TRUE;
1035 }
1036 
1037 POINTARRAY*
1038 ptarray_force_dims(const POINTARRAY *pa, int hasz, int hasm)
1039 {
1040  /* TODO handle zero-length point arrays */
1041  int i;
1042  int in_hasz = FLAGS_GET_Z(pa->flags);
1043  int in_hasm = FLAGS_GET_M(pa->flags);
1044  POINT4D pt;
1045  POINTARRAY *pa_out = ptarray_construct_empty(hasz, hasm, pa->npoints);
1046 
1047  for( i = 0; i < pa->npoints; i++ )
1048  {
1049  getPoint4d_p(pa, i, &pt);
1050  if( hasz && ! in_hasz )
1051  pt.z = 0.0;
1052  if( hasm && ! in_hasm )
1053  pt.m = 0.0;
1054  ptarray_append_point(pa_out, &pt, LW_TRUE);
1055  }
1056 
1057  return pa_out;
1058 }
1059 
1060 POINTARRAY *
1061 ptarray_substring(POINTARRAY *ipa, double from, double to, double tolerance)
1062 {
1063  POINTARRAY *dpa;
1064  POINT4D pt;
1065  POINT4D p1, p2;
1066  POINT4D *p1ptr=&p1; /* don't break strict-aliasing rule */
1067  POINT4D *p2ptr=&p2;
1068  int nsegs, i;
1069  double length, slength, tlength;
1070  int state = 0; /* 0=before, 1=inside */
1071 
1072  /*
1073  * Create a dynamic pointarray with an initial capacity
1074  * equal to full copy of input points
1075  */
1077 
1078  /* Compute total line length */
1079  length = ptarray_length_2d(ipa);
1080 
1081 
1082  LWDEBUGF(3, "Total length: %g", length);
1083 
1084 
1085  /* Get 'from' and 'to' lengths */
1086  from = length*from;
1087  to = length*to;
1088 
1089 
1090  LWDEBUGF(3, "From/To: %g/%g", from, to);
1091 
1092 
1093  tlength = 0;
1094  getPoint4d_p(ipa, 0, &p1);
1095  nsegs = ipa->npoints - 1;
1096  for ( i = 0; i < nsegs; i++ )
1097  {
1098  double dseg;
1099 
1100  getPoint4d_p(ipa, i+1, &p2);
1101 
1102 
1103  LWDEBUGF(3 ,"Segment %d: (%g,%g,%g,%g)-(%g,%g,%g,%g)",
1104  i, p1.x, p1.y, p1.z, p1.m, p2.x, p2.y, p2.z, p2.m);
1105 
1106 
1107  /* Find the length of this segment */
1108  slength = distance2d_pt_pt((POINT2D *)p1ptr, (POINT2D *)p2ptr);
1109 
1110  /*
1111  * We are before requested start.
1112  */
1113  if ( state == 0 ) /* before */
1114  {
1115 
1116  LWDEBUG(3, " Before start");
1117 
1118  if ( fabs ( from - ( tlength + slength ) ) <= tolerance )
1119  {
1120 
1121  LWDEBUG(3, " Second point is our start");
1122 
1123  /*
1124  * Second point is our start
1125  */
1126  ptarray_append_point(dpa, &p2, LW_FALSE);
1127  state=1; /* we're inside now */
1128  goto END;
1129  }
1130 
1131  else if ( fabs(from - tlength) <= tolerance )
1132  {
1133 
1134  LWDEBUG(3, " First point is our start");
1135 
1136  /*
1137  * First point is our start
1138  */
1139  ptarray_append_point(dpa, &p1, LW_FALSE);
1140 
1141  /*
1142  * We're inside now, but will check
1143  * 'to' point as well
1144  */
1145  state=1;
1146  }
1147 
1148  /*
1149  * Didn't reach the 'from' point,
1150  * nothing to do
1151  */
1152  else if ( from > tlength + slength ) goto END;
1153 
1154  else /* tlength < from < tlength+slength */
1155  {
1156 
1157  LWDEBUG(3, " Seg contains first point");
1158 
1159  /*
1160  * Our start is between first and
1161  * second point
1162  */
1163  dseg = (from - tlength) / slength;
1164 
1165  interpolate_point4d(&p1, &p2, &pt, dseg);
1166 
1167  ptarray_append_point(dpa, &pt, LW_FALSE);
1168 
1169  /*
1170  * We're inside now, but will check
1171  * 'to' point as well
1172  */
1173  state=1;
1174  }
1175  }
1176 
1177  if ( state == 1 ) /* inside */
1178  {
1179 
1180  LWDEBUG(3, " Inside");
1181 
1182  /*
1183  * 'to' point is our second point.
1184  */
1185  if ( fabs(to - ( tlength + slength ) ) <= tolerance )
1186  {
1187 
1188  LWDEBUG(3, " Second point is our end");
1189 
1190  ptarray_append_point(dpa, &p2, LW_FALSE);
1191  break; /* substring complete */
1192  }
1193 
1194  /*
1195  * 'to' point is our first point.
1196  * (should only happen if 'to' is 0)
1197  */
1198  else if ( fabs(to - tlength) <= tolerance )
1199  {
1200 
1201  LWDEBUG(3, " First point is our end");
1202 
1203  ptarray_append_point(dpa, &p1, LW_FALSE);
1204 
1205  break; /* substring complete */
1206  }
1207 
1208  /*
1209  * Didn't reach the 'end' point,
1210  * just copy second point
1211  */
1212  else if ( to > tlength + slength )
1213  {
1214  ptarray_append_point(dpa, &p2, LW_FALSE);
1215  goto END;
1216  }
1217 
1218  /*
1219  * 'to' point falls on this segment
1220  * Interpolate and break.
1221  */
1222  else if ( to < tlength + slength )
1223  {
1224 
1225  LWDEBUG(3, " Seg contains our end");
1226 
1227  dseg = (to - tlength) / slength;
1228  interpolate_point4d(&p1, &p2, &pt, dseg);
1229 
1230  ptarray_append_point(dpa, &pt, LW_FALSE);
1231 
1232  break;
1233  }
1234 
1235  else
1236  {
1237  LWDEBUG(3, "Unhandled case");
1238  }
1239  }
1240 
1241 
1242 END:
1243 
1244  tlength += slength;
1245  memcpy(&p1, &p2, sizeof(POINT4D));
1246  }
1247 
1248  LWDEBUGF(3, "Out of loop, ptarray has %d points", dpa->npoints);
1249 
1250  return dpa;
1251 }
1252 
1253 /*
1254  * Write into the *ret argument coordinates of the closes point on
1255  * the given segment to the reference input point.
1256  */
1257 void
1258 closest_point_on_segment(const POINT4D *p, const POINT4D *A, const POINT4D *B, POINT4D *ret)
1259 {
1260  double r;
1261 
1262  if ( FP_EQUALS(A->x, B->x) && FP_EQUALS(A->y, B->y) )
1263  {
1264  *ret = *A;
1265  return;
1266  }
1267 
1268  /*
1269  * We use comp.graphics.algorithms Frequently Asked Questions method
1270  *
1271  * (1) AC dot AB
1272  * r = ----------
1273  * ||AB||^2
1274  * r has the following meaning:
1275  * r=0 P = A
1276  * r=1 P = B
1277  * r<0 P is on the backward extension of AB
1278  * r>1 P is on the forward extension of AB
1279  * 0<r<1 P is interior to AB
1280  *
1281  */
1282  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) );
1283 
1284  if (r<0)
1285  {
1286  *ret = *A;
1287  return;
1288  }
1289  if (r>1)
1290  {
1291  *ret = *B;
1292  return;
1293  }
1294 
1295  ret->x = A->x + ( (B->x - A->x) * r );
1296  ret->y = A->y + ( (B->y - A->y) * r );
1297  ret->z = A->z + ( (B->z - A->z) * r );
1298  ret->m = A->m + ( (B->m - A->m) * r );
1299 }
1300 
1301 /*
1302  * Given a point, returns the location of closest point on pointarray
1303  * and, optionally, it's actual distance from the point array.
1304  */
1305 double
1306 ptarray_locate_point(const POINTARRAY *pa, const POINT4D *p4d, double *mindistout, POINT4D *proj4d)
1307 {
1308  double mindist=-1;
1309  double tlen, plen;
1310  int t, seg=-1;
1311  POINT4D start4d, end4d, projtmp;
1312  POINT2D proj, p;
1313  const POINT2D *start = NULL, *end = NULL;
1314 
1315  /* Initialize our 2D copy of the input parameter */
1316  p.x = p4d->x;
1317  p.y = p4d->y;
1318 
1319  if ( ! proj4d ) proj4d = &projtmp;
1320 
1321  start = getPoint2d_cp(pa, 0);
1322 
1323  /* If the pointarray has only one point, the nearest point is */
1324  /* just that point */
1325  if ( pa->npoints == 1 )
1326  {
1327  getPoint4d_p(pa, 0, proj4d);
1328  if ( mindistout )
1329  *mindistout = distance2d_pt_pt(&p, start);
1330  return 0.0;
1331  }
1332 
1333  /* Loop through pointarray looking for nearest segment */
1334  for (t=1; t<pa->npoints; t++)
1335  {
1336  double dist;
1337  end = getPoint2d_cp(pa, t);
1338  dist = distance2d_pt_seg(&p, start, end);
1339 
1340  if (t==1 || dist < mindist )
1341  {
1342  mindist = dist;
1343  seg=t-1;
1344  }
1345 
1346  if ( mindist == 0 )
1347  {
1348  LWDEBUG(3, "Breaking on mindist=0");
1349  break;
1350  }
1351 
1352  start = end;
1353  }
1354 
1355  if ( mindistout ) *mindistout = mindist;
1356 
1357  LWDEBUGF(3, "Closest segment: %d", seg);
1358  LWDEBUGF(3, "mindist: %g", mindist);
1359 
1360  /*
1361  * We need to project the
1362  * point on the closest segment.
1363  */
1364  getPoint4d_p(pa, seg, &start4d);
1365  getPoint4d_p(pa, seg+1, &end4d);
1366  closest_point_on_segment(p4d, &start4d, &end4d, proj4d);
1367 
1368  /* Copy 4D values into 2D holder */
1369  proj.x = proj4d->x;
1370  proj.y = proj4d->y;
1371 
1372  LWDEBUGF(3, "Closest segment:%d, npoints:%d", seg, pa->npoints);
1373 
1374  /* For robustness, force 1 when closest point == endpoint */
1375  if ( (seg >= (pa->npoints-2)) && p2d_same(&proj, end) )
1376  {
1377  return 1.0;
1378  }
1379 
1380  LWDEBUGF(3, "Closest point on segment: %g,%g", proj.x, proj.y);
1381 
1382  tlen = ptarray_length_2d(pa);
1383 
1384  LWDEBUGF(3, "tlen %g", tlen);
1385 
1386  /* Location of any point on a zero-length line is 0 */
1387  /* See http://trac.osgeo.org/postgis/ticket/1772#comment:2 */
1388  if ( tlen == 0 ) return 0;
1389 
1390  plen=0;
1391  start = getPoint2d_cp(pa, 0);
1392  for (t=0; t<seg; t++, start=end)
1393  {
1394  end = getPoint2d_cp(pa, t+1);
1395  plen += distance2d_pt_pt(start, end);
1396 
1397  LWDEBUGF(4, "Segment %d made plen %g", t, plen);
1398  }
1399 
1400  plen+=distance2d_pt_pt(&proj, start);
1401 
1402  LWDEBUGF(3, "plen %g, tlen %g", plen, tlen);
1403 
1404  return plen/tlen;
1405 }
1406 
1416 void
1418 {
1419  int i;
1420  double x;
1421 
1422  for (i=0; i<pa->npoints; i++)
1423  {
1424  memcpy(&x, getPoint_internal(pa, i), sizeof(double));
1425  if ( x < 0 ) x+= 360;
1426  else if ( x > 180 ) x -= 360;
1427  memcpy(getPoint_internal(pa, i), &x, sizeof(double));
1428  }
1429 }
1430 
1431 static void
1433 {
1434  int ndims = FLAGS_NDIMS(pa->flags);
1435  switch (ndims)
1436  {
1437  case 2:
1438  {
1439  POINT2D *p_from = (POINT2D*)(getPoint_internal(pa, from));
1440  POINT2D *p_to = (POINT2D*)(getPoint_internal(pa, to));
1441  *p_to = *p_from;
1442  return;
1443  }
1444  case 3:
1445  {
1446  POINT3D *p_from = (POINT3D*)(getPoint_internal(pa, from));
1447  POINT3D *p_to = (POINT3D*)(getPoint_internal(pa, to));
1448  *p_to = *p_from;
1449  return;
1450  }
1451  case 4:
1452  {
1453  POINT4D *p_from = (POINT4D*)(getPoint_internal(pa, from));
1454  POINT4D *p_to = (POINT4D*)(getPoint_internal(pa, to));
1455  *p_to = *p_from;
1456  return;
1457  }
1458  default:
1459  {
1460  lwerror("%s: unsupported number of dimensions - %d", __func__, ndims);
1461  return;
1462  }
1463  }
1464  return;
1465 }
1466 
1467 static void
1469 {
1470  uint32_t i;
1471  double tolsq = tolerance * tolerance;
1472  const POINT2D *last = NULL;
1473  const POINT2D *pt;
1474  uint32_t n_points = pa->npoints;
1475  uint32_t n_points_out = 1;
1476  size_t pt_size = ptarray_point_size(pa);
1477 
1478  double dsq = FLT_MAX;
1479 
1480  /* No-op on short inputs */
1481  if ( n_points <= min_points ) return;
1482 
1483  last = getPoint2d_cp(pa, 0);
1484  for (i = 1; i < n_points; i++)
1485  {
1486  int last_point = (i == n_points-1);
1487 
1488  /* Look straight into the abyss */
1489  pt = getPoint2d_cp(pa, i);
1490 
1491  /* Don't drop points if we are running short of points */
1492  if (n_points + n_points_out > min_points + i)
1493  {
1494  if (tolerance > 0.0)
1495  {
1496  /* Only drop points that are within our tolerance */
1497  dsq = distance2d_sqr_pt_pt(last, pt);
1498  /* Allow any point but the last one to be dropped */
1499  if (!last_point && dsq <= tolsq)
1500  {
1501  continue;
1502  }
1503  }
1504  else
1505  {
1506  /* At tolerance zero, only skip exact dupes */
1507  if (memcmp((char*)pt, (char*)last, pt_size) == 0)
1508  continue;
1509  }
1510 
1511  /* Got to last point, and it's not very different from */
1512  /* the point that preceded it. We want to keep the last */
1513  /* point, not the second-to-last one, so we pull our write */
1514  /* index back one value */
1515  if (last_point && n_points_out > 1 && tolerance > 0.0 && dsq <= tolsq)
1516  {
1517  n_points_out--;
1518  }
1519  }
1520 
1521  /* Compact all remaining values to front of array */
1522  ptarray_copy_point(pa, i, n_points_out++);
1523  last = pt;
1524  }
1525  /* Adjust array length */
1526  pa->npoints = n_points_out;
1527  return;
1528 }
1529 
1530 /*
1531  * Returns a POINTARRAY with consecutive equal points
1532  * removed. Equality test on all dimensions of input.
1533  *
1534  * Always returns a newly allocated object.
1535  *
1536  */
1537 POINTARRAY *
1538 ptarray_remove_repeated_points_minpoints(const POINTARRAY *in, double tolerance, int minpoints)
1539 {
1540  POINTARRAY *out = ptarray_clone_deep(in);
1541  ptarray_remove_repeated_points_in_place(out, tolerance, minpoints);
1542  return out;
1543 }
1544 
1545 POINTARRAY *
1546 ptarray_remove_repeated_points(const POINTARRAY *in, double tolerance)
1547 {
1548  return ptarray_remove_repeated_points_minpoints(in, tolerance, 2);
1549 }
1550 
1551 static void
1552 ptarray_dp_findsplit(POINTARRAY *pts, int p1, int p2, int *split, double *dist)
1553 {
1554  int k;
1555  const POINT2D *pk, *pa, *pb;
1556  double tmp, d;
1557 
1558  LWDEBUG(4, "function called");
1559 
1560  *split = p1;
1561  d = -1;
1562 
1563  if (p1 + 1 < p2)
1564  {
1565 
1566  pa = getPoint2d_cp(pts, p1);
1567  pb = getPoint2d_cp(pts, p2);
1568 
1569  LWDEBUGF(4, "P%d(%f,%f) to P%d(%f,%f)",
1570  p1, pa->x, pa->y, p2, pb->x, pb->y);
1571 
1572  for (k=p1+1; k<p2; k++)
1573  {
1574  pk = getPoint2d_cp(pts, k);
1575 
1576  LWDEBUGF(4, "P%d(%f,%f)", k, pk->x, pk->y);
1577 
1578  /* distance computation */
1579  tmp = distance2d_sqr_pt_seg(pk, pa, pb);
1580 
1581  if (tmp > d)
1582  {
1583  d = tmp; /* record the maximum */
1584  *split = k;
1585 
1586  LWDEBUGF(4, "P%d is farthest (%g)", k, d);
1587  }
1588  }
1589  *dist = d;
1590 
1591  } /* length---should be redone if can == 0 */
1592  else
1593  {
1594  LWDEBUG(3, "segment too short, no split/no dist");
1595  *dist = -1;
1596  }
1597 
1598 }
1599 
1600 POINTARRAY *
1601 ptarray_simplify(POINTARRAY *inpts, double epsilon, unsigned int minpts)
1602 {
1603  int *stack; /* recursion stack */
1604  int sp=-1; /* recursion stack pointer */
1605  int p1, split;
1606  double dist;
1607  POINTARRAY *outpts;
1608  POINT4D pt;
1609 
1610  double eps_sqr = epsilon * epsilon;
1611 
1612  /* Allocate recursion stack */
1613  stack = lwalloc(sizeof(int)*inpts->npoints);
1614 
1615  p1 = 0;
1616  stack[++sp] = inpts->npoints-1;
1617 
1618  LWDEBUGF(2, "Input has %d pts and %d dims", inpts->npoints,
1619  FLAGS_NDIMS(inpts->flags));
1620 
1621  /* Allocate output POINTARRAY, and add first point. */
1622  outpts = ptarray_construct_empty(FLAGS_GET_Z(inpts->flags), FLAGS_GET_M(inpts->flags), inpts->npoints);
1623  getPoint4d_p(inpts, 0, &pt);
1624  ptarray_append_point(outpts, &pt, LW_FALSE);
1625 
1626  LWDEBUG(3, "Added P0 to simplified point array (size 1)");
1627 
1628  do
1629  {
1630 
1631  ptarray_dp_findsplit(inpts, p1, stack[sp], &split, &dist);
1632 
1633  LWDEBUGF(3, "Farthest point from P%d-P%d is P%d (dist. %g)", p1, stack[sp], split, dist);
1634 
1635  if (dist > eps_sqr || ( outpts->npoints+sp+1 < minpts && dist >= 0 ) )
1636  {
1637  LWDEBUGF(4, "Added P%d to stack (outpts:%d)", split, sp);
1638  stack[++sp] = split;
1639  }
1640  else
1641  {
1642  getPoint4d_p(inpts, stack[sp], &pt);
1643  LWDEBUGF(4, "npoints , minpoints %d %d", outpts->npoints, minpts);
1644  ptarray_append_point(outpts, &pt, LW_FALSE);
1645 
1646  LWDEBUGF(4, "Added P%d to simplified point array (size: %d)", stack[sp], outpts->npoints);
1647 
1648  p1 = stack[sp--];
1649  }
1650 
1651  LWDEBUGF(4, "stack pointer = %d", sp);
1652  }
1653  while (! (sp<0) );
1654 
1655  lwfree(stack);
1656  return outpts;
1657 }
1658 
1664 double
1666 {
1667  double dist = 0.0;
1668  int i;
1669  const POINT2D *a1;
1670  const POINT2D *a2;
1671  const POINT2D *a3;
1672 
1673  if ( pts->npoints % 2 != 1 )
1674  lwerror("arc point array with even number of points");
1675 
1676  a1 = getPoint2d_cp(pts, 0);
1677 
1678  for ( i=2; i < pts->npoints; i += 2 )
1679  {
1680  a2 = getPoint2d_cp(pts, i-1);
1681  a3 = getPoint2d_cp(pts, i);
1682  dist += lw_arc_length(a1, a2, a3);
1683  a1 = a3;
1684  }
1685  return dist;
1686 }
1687 
1691 double
1693 {
1694  double dist = 0.0;
1695  int i;
1696  const POINT2D *frm;
1697  const POINT2D *to;
1698 
1699  if ( pts->npoints < 2 ) return 0.0;
1700 
1701  frm = getPoint2d_cp(pts, 0);
1702 
1703  for ( i=1; i < pts->npoints; i++ )
1704  {
1705  to = getPoint2d_cp(pts, i);
1706 
1707  dist += sqrt( ((frm->x - to->x)*(frm->x - to->x)) +
1708  ((frm->y - to->y)*(frm->y - to->y)) );
1709 
1710  frm = to;
1711  }
1712  return dist;
1713 }
1714 
1719 double
1721 {
1722  double dist = 0.0;
1723  int i;
1724  POINT3DZ frm;
1725  POINT3DZ to;
1726 
1727  if ( pts->npoints < 2 ) return 0.0;
1728 
1729  /* compute 2d length if 3d is not available */
1730  if ( ! FLAGS_GET_Z(pts->flags) ) return ptarray_length_2d(pts);
1731 
1732  getPoint3dz_p(pts, 0, &frm);
1733  for ( i=1; i < pts->npoints; i++ )
1734  {
1735  getPoint3dz_p(pts, i, &to);
1736  dist += sqrt( ((frm.x - to.x)*(frm.x - to.x)) +
1737  ((frm.y - to.y)*(frm.y - to.y)) +
1738  ((frm.z - to.z)*(frm.z - to.z)) );
1739  frm = to;
1740  }
1741  return dist;
1742 }
1743 
1744 
1745 /*
1746  * Get a pointer to nth point of a POINTARRAY.
1747  *
1748  * Casting to returned pointer to POINT2D* should be safe,
1749  * as gserialized format always keeps the POINTARRAY pointer
1750  * aligned to double boundary.
1751  */
1752 uint8_t *
1754 {
1755  size_t size;
1756  uint8_t *ptr;
1757 
1758 #if PARANOIA_LEVEL > 0
1759  if ( pa == NULL )
1760  {
1761  lwerror("%s [%d] got NULL pointarray", __FILE__, __LINE__);
1762  return NULL;
1763  }
1764 
1765  LWDEBUGF(5, "(n=%d, pa.npoints=%d, pa.maxpoints=%d)",n,pa->npoints,pa->maxpoints);
1766 
1767  if ( ( n < 0 ) ||
1768  ( n > pa->npoints ) ||
1769  ( n >= pa->maxpoints ) )
1770  {
1771  lwerror("%s [%d] called outside of ptarray range (n=%d, pa.npoints=%d, pa.maxpoints=%d)", __FILE__, __LINE__, n, pa->npoints, pa->maxpoints);
1772  return NULL; /*error */
1773  }
1774 #endif
1775 
1776  size = ptarray_point_size(pa);
1777 
1778  ptr = pa->serialized_pointlist + size * n;
1779  if ( FLAGS_NDIMS(pa->flags) == 2)
1780  {
1781  LWDEBUGF(5, "point = %g %g", *((double*)(ptr)), *((double*)(ptr+8)));
1782  }
1783  else if ( FLAGS_NDIMS(pa->flags) == 3)
1784  {
1785  LWDEBUGF(5, "point = %g %g %g", *((double*)(ptr)), *((double*)(ptr+8)), *((double*)(ptr+16)));
1786  }
1787  else if ( FLAGS_NDIMS(pa->flags) == 4)
1788  {
1789  LWDEBUGF(5, "point = %g %g %g %g", *((double*)(ptr)), *((double*)(ptr+8)), *((double*)(ptr+16)), *((double*)(ptr+24)));
1790  }
1791 
1792  return ptr;
1793 }
1794 
1795 
1799 void
1801 {
1802  int i;
1803  double x,y,z;
1804  POINT4D p4d;
1805 
1806  LWDEBUG(2, "lwgeom_affine_ptarray start");
1807 
1808  if ( FLAGS_GET_Z(pa->flags) )
1809  {
1810  LWDEBUG(3, " has z");
1811 
1812  for (i=0; i<pa->npoints; i++)
1813  {
1814  getPoint4d_p(pa, i, &p4d);
1815  x = p4d.x;
1816  y = p4d.y;
1817  z = p4d.z;
1818  p4d.x = a->afac * x + a->bfac * y + a->cfac * z + a->xoff;
1819  p4d.y = a->dfac * x + a->efac * y + a->ffac * z + a->yoff;
1820  p4d.z = a->gfac * x + a->hfac * y + a->ifac * z + a->zoff;
1821  ptarray_set_point4d(pa, i, &p4d);
1822 
1823  LWDEBUGF(3, " POINT %g %g %g => %g %g %g", x, y, z, p4d.x, p4d.y, p4d.z);
1824  }
1825  }
1826  else
1827  {
1828  LWDEBUG(3, " doesn't have z");
1829 
1830  for (i=0; i<pa->npoints; i++)
1831  {
1832  getPoint4d_p(pa, i, &p4d);
1833  x = p4d.x;
1834  y = p4d.y;
1835  p4d.x = a->afac * x + a->bfac * y + a->xoff;
1836  p4d.y = a->dfac * x + a->efac * y + a->yoff;
1837  ptarray_set_point4d(pa, i, &p4d);
1838 
1839  LWDEBUGF(3, " POINT %g %g => %g %g", x, y, p4d.x, p4d.y);
1840  }
1841  }
1842 
1843  LWDEBUG(3, "lwgeom_affine_ptarray end");
1844 
1845 }
1846 
1850 void
1852 {
1853  int i;
1854  POINT4D p4d;
1855 
1856  LWDEBUG(3, "ptarray_scale start");
1857 
1858  for (i=0; i<pa->npoints; ++i)
1859  {
1860  getPoint4d_p(pa, i, &p4d);
1861  p4d.x *= fact->x;
1862  p4d.y *= fact->y;
1863  p4d.z *= fact->z;
1864  p4d.m *= fact->m;
1865  ptarray_set_point4d(pa, i, &p4d);
1866  }
1867 
1868  LWDEBUG(3, "ptarray_scale end");
1869 
1870 }
1871 
1872 int
1874 {
1875  return getPoint4d_p(pa, 0, pt);
1876 }
1877 
1878 
1879 
1880 
1881 /*
1882  * Stick an array of points to the given gridspec.
1883  * Return "gridded" points in *outpts and their number in *outptsn.
1884  *
1885  * Two consecutive points falling on the same grid cell are collapsed
1886  * into one single point.
1887  *
1888  */
1889 POINTARRAY *
1890 ptarray_grid(const POINTARRAY *pa, const gridspec *grid)
1891 {
1892  POINT4D pt;
1893  int ipn; /* input point numbers */
1894  POINTARRAY *dpa;
1895 
1896  LWDEBUGF(2, "ptarray_grid called on %p", pa);
1897 
1899 
1900  for (ipn=0; ipn<pa->npoints; ++ipn)
1901  {
1902 
1903  getPoint4d_p(pa, ipn, &pt);
1904 
1905  if ( grid->xsize )
1906  pt.x = rint((pt.x - grid->ipx)/grid->xsize) *
1907  grid->xsize + grid->ipx;
1908 
1909  if ( grid->ysize )
1910  pt.y = rint((pt.y - grid->ipy)/grid->ysize) *
1911  grid->ysize + grid->ipy;
1912 
1913  if ( FLAGS_GET_Z(pa->flags) && grid->zsize )
1914  pt.z = rint((pt.z - grid->ipz)/grid->zsize) *
1915  grid->zsize + grid->ipz;
1916 
1917  if ( FLAGS_GET_M(pa->flags) && grid->msize )
1918  pt.m = rint((pt.m - grid->ipm)/grid->msize) *
1919  grid->msize + grid->ipm;
1920 
1921  ptarray_append_point(dpa, &pt, LW_FALSE);
1922 
1923  }
1924 
1925  return dpa;
1926 }
1927 
1928 int
1929 ptarray_npoints_in_rect(const POINTARRAY *pa, const GBOX *gbox)
1930 {
1931  const POINT2D *pt;
1932  int n = 0;
1933  int i;
1934  for ( i = 0; i < pa->npoints; i++ )
1935  {
1936  pt = getPoint2d_cp(pa, i);
1937  if ( gbox_contains_point2d(gbox, pt) )
1938  n++;
1939  }
1940  return n;
1941 }
1942 
1943 
void ptarray_set_point4d(POINTARRAY *pa, int n, const POINT4D *p4d)
Definition: lwgeom_api.c:437
double x
Definition: liblwgeom.h:352
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:227
POINTARRAY * ptarray_remove_repeated_points_minpoints(const POINTARRAY *in, double tolerance, int minpoints)
Definition: ptarray.c:1538
double z
Definition: liblwgeom.h:334
int ptarrayarc_contains_point_partial(const POINTARRAY *pa, const POINT2D *pt, int check_closed, int *winding_number)
Definition: ptarray.c:841
POINTARRAY * ptarray_substring(POINTARRAY *ipa, double from, double to, double tolerance)
start location (distance from start / total distance) end location (distance from start / total dist...
Definition: ptarray.c:1061
double y
Definition: liblwgeom.h:334
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:105
double m
Definition: liblwgeom.h:352
uint8_t * serialized_pointlist
Definition: liblwgeom.h:366
double x
Definition: liblwgeom.h:334
char * r
Definition: cu_in_wkt.c:24
int ptarray_point_size(const POINTARRAY *pa)
Definition: ptarray.c:54
POINTARRAY * ptarray_merge(POINTARRAY *pa1, POINTARRAY *pa2)
Merge two given POINTARRAY and returns a pointer on the new aggregate one.
Definition: ptarray.c:599
void lwfree(void *mem)
Definition: lwutil.c:244
#define FLAGS_GET_READONLY(flags)
Definition: liblwgeom.h:144
int lw_arc_calculate_gbox_cartesian_2d(const POINT2D *A1, const POINT2D *A2, const POINT2D *A3, GBOX *gbox)
Definition: g_box.c:464
double ptarray_length(const POINTARRAY *pts)
Find the 3d/2d length of the given POINTARRAY (depending on its dimensionality)
Definition: ptarray.c:1720
int npoints
Definition: liblwgeom.h:371
POINTARRAY * ptarray_clone_deep(const POINTARRAY *in)
Deep clone a pointarray (also clones serialized pointlist)
Definition: ptarray.c:630
POINTARRAY * ptarray_simplify(POINTARRAY *inpts, double epsilon, unsigned int minpts)
Definition: ptarray.c:1601
Datum area(PG_FUNCTION_ARGS)
double ptarray_locate_point(const POINTARRAY *pa, const POINT4D *p4d, double *mindistout, POINT4D *proj4d)
Definition: ptarray.c:1306
double xmax
Definition: liblwgeom.h:293
POINTARRAY * ptarray_flip_coordinates(POINTARRAY *pa)
Reverse X and Y axis on a given POINTARRAY.
Definition: ptarray.c:369
void ptarray_reverse(POINTARRAY *pa)
Definition: ptarray.c:343
#define LW_SUCCESS
Definition: liblwgeom.h:80
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&#39;s ordinate array.
Definition: ptarray.c:295
#define FLAGS_GET_ZM(flags)
Definition: liblwgeom.h:153
#define LWDEBUG(level, msg)
Definition: lwgeom_log.h:83
void ptarray_scale(POINTARRAY *pa, const POINT4D *fact)
Scale a pointarray.
Definition: ptarray.c:1851
int ptarray_is_closed_2d(const POINTARRAY *in)
Definition: ptarray.c:697
#define LW_ON_INTERRUPT(x)
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, then a duplicate point will not be added.
Definition: ptarray.c:156
void ptarray_swap_ordinates(POINTARRAY *pa, LWORD o1, LWORD o2)
Swap ordinate values o1 and o2 on a given POINTARRAY.
Definition: ptarray.c:388
int ptarray_isccw(const POINTARRAY *pa)
Definition: ptarray.c:1029
POINTARRAY * ptarray_removePoint(POINTARRAY *pa, uint32_t which)
Remove a point from a pointarray.
Definition: ptarray.c:557
double distance2d_pt_pt(const POINT2D *p1, const POINT2D *p2)
The old function nessecary for ptarray_segmentize2d in ptarray.c.
Definition: measures.c:2317
double zoff
Definition: liblwgeom.h:270
#define FP_MIN(A, B)
int ptarray_npoints_in_rect(const POINTARRAY *pa, const GBOX *gbox)
Definition: ptarray.c:1929
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:309
double distance2d_pt_seg(const POINT2D *p, const POINT2D *A, const POINT2D *B)
The old function nessecary for ptarray_segmentize2d in ptarray.c.
Definition: measures.c:2342
void ptarray_longitude_shift(POINTARRAY *pa)
Longitude shift for a pointarray.
Definition: ptarray.c:1417
void interpolate_point4d(POINT4D *A, 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:680
double ifac
Definition: liblwgeom.h:270
POINTARRAY * ptarray_clone(const POINTARRAY *in)
Clone a POINTARRAY object.
Definition: ptarray.c:661
double ffac
Definition: liblwgeom.h:270
double xoff
Definition: liblwgeom.h:270
double afac
Definition: liblwgeom.h:270
int ptarray_startpoint(const POINTARRAY *pa, POINT4D *pt)
Definition: ptarray.c:1873
#define LW_FAILURE
Definition: liblwgeom.h:79
unsigned int uint32_t
Definition: uthash.h:78
double x
Definition: liblwgeom.h:328
int p2d_same(const POINT2D *p1, const POINT2D *p2)
Definition: lwalgorithm.c:49
double ymin
Definition: liblwgeom.h:294
POINTARRAY * ptarray_grid(const POINTARRAY *pa, const gridspec *grid)
Definition: ptarray.c:1890
double lw_arc_length(const POINT2D *A1, const POINT2D *A2, const POINT2D *A3)
Returns the length of a circular arc segment.
Definition: lwalgorithm.c:118
int getPoint3dz_p(const POINTARRAY *pa, int n, POINT3DZ *point)
Definition: lwgeom_api.c:214
int ptarray_contains_point_partial(const POINTARRAY *pa, const POINT2D *pt, int check_closed, int *winding_number)
Definition: ptarray.c:742
double xmin
Definition: liblwgeom.h:292
#define LW_FALSE
Definition: liblwgeom.h:77
const POINT2D * getPoint2d_cp(const POINTARRAY *pa, int n)
Returns a POINT2D pointer into the POINTARRAY serialized_ptlist, suitable for reading from...
Definition: lwgeom_api.c:373
POINTARRAY * ptarray_addPoint(const POINTARRAY *pa, uint8_t *p, size_t pdims, uint32_t where)
Add a point in a pointarray.
Definition: ptarray.c:504
uint8_t flags
Definition: liblwgeom.h:369
#define LW_TRUE
Return types for functions with status returns.
Definition: liblwgeom.h:76
void closest_point_on_segment(const POINT4D *p, const POINT4D *A, const POINT4D *B, POINT4D *ret)
Definition: ptarray.c:1258
int ptarray_insert_point(POINTARRAY *pa, const POINT4D *p, int where)
Insert a point into an existing POINTARRAY.
Definition: ptarray.c:96
int ptarray_is_closed_z(const POINTARRAY *in)
Definition: ptarray.c:723
static void ptarray_copy_point(POINTARRAY *pa, uint32_t from, uint32_t to)
Definition: ptarray.c:1432
int ptarray_has_z(const POINTARRAY *pa)
Definition: ptarray.c:36
#define LW_INSIDE
Constants for point-in-polygon return values.
int ptarray_has_m(const POINTARRAY *pa)
Definition: ptarray.c:43
double cfac
Definition: liblwgeom.h:270
double ymax
Definition: liblwgeom.h:295
enum LWORD_T LWORD
Ordinate names.
double y
Definition: liblwgeom.h:328
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:62
int getPoint2d_p(const POINTARRAY *pa, int n, POINT2D *point)
Definition: lwgeom_api.c:347
#define FLAGS_GET_Z(flags)
Macros for manipulating the &#39;flags&#39; byte.
Definition: liblwgeom.h:140
int ptarray_is_closed_3d(const POINTARRAY *in)
Definition: ptarray.c:710
double ptarray_length_2d(const POINTARRAY *pts)
Find the 2d length of the given POINTARRAY (even if it&#39;s 3d)
Definition: ptarray.c:1692
double dfac
Definition: liblwgeom.h:270
double z
Definition: liblwgeom.h:352
static void ptarray_remove_repeated_points_in_place(POINTARRAY *pa, double tolerance, uint32_t min_points)
Definition: ptarray.c:1468
uint8_t gflags(int hasz, int hasm, int geodetic)
Construct a new flags char.
Definition: g_util.c:145
double ptarray_signed_area(const POINTARRAY *pa)
Returns the area in cartesian units.
Definition: ptarray.c:998
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:1665
int ptarray_remove_point(POINTARRAY *pa, int where)
Remove a point from an existing POINTARRAY.
Definition: ptarray.c:261
POINTARRAY * ptarray_construct_empty(char hasz, char hasm, uint32_t maxpoints)
Create a new POINTARRAY with no points.
Definition: ptarray.c:70
int lw_arc_side(const POINT2D *A1, const POINT2D *A2, const POINT2D *A3, const POINT2D *Q)
Definition: lwalgorithm.c:178
#define LW_BOUNDARY
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:85
void * lwrealloc(void *mem, size_t size)
Definition: lwutil.c:237
double efac
Definition: liblwgeom.h:270
#define FLAGS_GET_M(flags)
Definition: liblwgeom.h:141
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:95
#define FP_EQUALS(A, B)
double yoff
Definition: liblwgeom.h:270
double distance2d_sqr_pt_pt(const POINT2D *p1, const POINT2D *p2)
Definition: measures.c:2327
int maxpoints
Definition: liblwgeom.h:372
#define LW_OUTSIDE
void * lwalloc(size_t size)
Definition: lwutil.c:229
int lw_segment_side(const POINT2D *p1, const POINT2D *p2, const POINT2D *q)
lw_segment_side()
Definition: lwalgorithm.c:64
void ptarray_free(POINTARRAY *pa)
Definition: ptarray.c:330
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:415
double y
Definition: liblwgeom.h:352
int ptarrayarc_contains_point(const POINTARRAY *pa, const POINT2D *pt)
For POINTARRAYs representing CIRCULARSTRINGS.
Definition: ptarray.c:835
double distance2d_sqr_pt_seg(const POINT2D *p, const POINT2D *A, const POINT2D *B)
Definition: measures.c:2391
char ptarray_same(const POINTARRAY *pa1, const POINTARRAY *pa2)
Definition: ptarray.c:479
double gfac
Definition: liblwgeom.h:270
#define LWDEBUGF(level, msg,...)
Definition: lwgeom_log.h:88
#define FLAGS_NDIMS(flags)
Definition: liblwgeom.h:152
unsigned char uint8_t
Definition: uthash.h:79
void ptarray_affine(POINTARRAY *pa, const AFFINE *a)
Affine transform a pointarray.
Definition: ptarray.c:1800
POINTARRAY * ptarray_force_dims(const POINTARRAY *pa, int hasz, int hasm)
Definition: ptarray.c:1038
uint8_t * getPoint_internal(const POINTARRAY *pa, int n)
Definition: ptarray.c:1753
int ptarray_is_closed(const POINTARRAY *in)
Check for ring closure using whatever dimensionality is declared on the pointarray.
Definition: ptarray.c:683
void lwerror(const char *fmt,...)
Write a notice out to the error handler.
Definition: lwutil.c:190
int gbox_contains_point2d(const GBOX *g, const POINT2D *p)
Definition: g_box.c:362
double hfac
Definition: liblwgeom.h:270
double bfac
Definition: liblwgeom.h:270
int getPoint4d_p(const POINTARRAY *pa, int n, POINT4D *point)
Definition: lwgeom_api.c:122
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:187
#define FLAGS_SET_READONLY(flags, value)
Definition: liblwgeom.h:150
POINTARRAY * ptarray_remove_repeated_points(const POINTARRAY *in, double tolerance)
Definition: ptarray.c:1546
static void ptarray_dp_findsplit(POINTARRAY *pts, int p1, int p2, int *split, double *dist)
Definition: ptarray.c:1552
#define FP_MAX(A, B)
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:736
Snap to grid.