PostGIS  2.3.8dev-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 #if PARANOIA_LEVEL > 0
395  assert(o1 < 4);
396  assert(o2 < 4);
397 #endif
398 
399  dp1 = ((double*)&p)+(unsigned)o1;
400  dp2 = ((double*)&p)+(unsigned)o2;
401  for (i=0 ; i < pa->npoints ; i++)
402  {
403  getPoint4d_p(pa, i, &p);
404  d = *dp2;
405  *dp2 = *dp1;
406  *dp1 = d;
407  ptarray_set_point4d(pa, i, &p);
408  }
409 }
410 
411 
419 POINTARRAY *
420 ptarray_segmentize2d(const POINTARRAY *ipa, double dist)
421 {
422  double segdist;
423  POINT4D p1, p2;
424  POINT4D pbuf;
425  POINTARRAY *opa;
426  int ipoff=0; /* input point offset */
427  int hasz = FLAGS_GET_Z(ipa->flags);
428  int hasm = FLAGS_GET_M(ipa->flags);
429 
430  pbuf.x = pbuf.y = pbuf.z = pbuf.m = 0;
431 
432  /* Initial storage */
433  opa = ptarray_construct_empty(hasz, hasm, ipa->npoints);
434 
435  /* Add first point */
436  getPoint4d_p(ipa, ipoff, &p1);
437  ptarray_append_point(opa, &p1, LW_FALSE);
438 
439  ipoff++;
440 
441  while (ipoff<ipa->npoints)
442  {
443  /*
444  * We use these pointers to avoid
445  * "strict-aliasing rules break" warning raised
446  * by gcc (3.3 and up).
447  *
448  * It looks that casting a variable address (also
449  * referred to as "type-punned pointer")
450  * breaks those "strict" rules.
451  *
452  */
453  POINT4D *p1ptr=&p1, *p2ptr=&p2;
454 
455  getPoint4d_p(ipa, ipoff, &p2);
456 
457  segdist = distance2d_pt_pt((POINT2D *)p1ptr, (POINT2D *)p2ptr);
458 
459  if (segdist > dist) /* add an intermediate point */
460  {
461  pbuf.x = p1.x + (p2.x-p1.x)/segdist * dist;
462  pbuf.y = p1.y + (p2.y-p1.y)/segdist * dist;
463  if( hasz )
464  pbuf.z = p1.z + (p2.z-p1.z)/segdist * dist;
465  if( hasm )
466  pbuf.m = p1.m + (p2.m-p1.m)/segdist * dist;
467  ptarray_append_point(opa, &pbuf, LW_FALSE);
468  p1 = pbuf;
469  }
470  else /* copy second point */
471  {
472  ptarray_append_point(opa, &p2, (ipa->npoints==2)?LW_TRUE:LW_FALSE);
473  p1 = p2;
474  ipoff++;
475  }
476 
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  if ( where == -1 ) where = pa->npoints;
543 
544  if ( where )
545  {
546  memcpy(getPoint_internal(ret, 0), getPoint_internal(pa, 0), ptsize*where);
547  }
548 
549  memcpy(getPoint_internal(ret, where), (uint8_t *)&pbuf, ptsize);
550 
551  if ( where+1 != ret->npoints )
552  {
553  memcpy(getPoint_internal(ret, where+1),
554  getPoint_internal(pa, where),
555  ptsize*(pa->npoints-where));
556  }
557 
558  return ret;
559 }
560 
561 POINTARRAY *
562 ptarray_removePoint(POINTARRAY *pa, uint32_t which)
563 {
564  POINTARRAY *ret;
565  size_t ptsize = ptarray_point_size(pa);
566 
567  LWDEBUGF(3, "pa %x which %d", pa, which);
568 
569 #if PARANOIA_LEVEL > 0
570  if ( which > pa->npoints-1 )
571  {
572  lwerror("ptarray_removePoint: offset (%d) out of range (%d..%d)",
573  which, 0, pa->npoints-1);
574  return NULL;
575  }
576 
577  if ( pa->npoints < 3 )
578  {
579  lwerror("ptarray_removePointe: can't remove a point from a 2-vertex POINTARRAY");
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  size_t size;
638 
639  LWDEBUG(3, "ptarray_clone_deep called.");
640 
641  out->flags = in->flags;
642  out->npoints = in->npoints;
643  out->maxpoints = in->npoints;
644 
645  FLAGS_SET_READONLY(out->flags, 0);
646 
647  size = in->npoints * ptarray_point_size(in);
648  out->serialized_pointlist = lwalloc(size);
649  memcpy(out->serialized_pointlist, in->serialized_pointlist, size);
650 
651  return out;
652 }
653 
657 POINTARRAY *
659 {
660  POINTARRAY *out = lwalloc(sizeof(POINTARRAY));
661 
662  LWDEBUG(3, "ptarray_clone_deep called.");
663 
664  out->flags = in->flags;
665  out->npoints = in->npoints;
666  out->maxpoints = in->maxpoints;
667 
668  FLAGS_SET_READONLY(out->flags, 1);
669 
671 
672  return out;
673 }
674 
679 int
681 {
682  if (!in)
683  {
684  lwerror("ptarray_is_closed: called with null point array");
685  return 0;
686  }
687  if (in->npoints <= 1 ) return in->npoints; /* single-point are closed, empty not closed */
688 
689  return 0 == memcmp(getPoint_internal(in, 0), getPoint_internal(in, in->npoints-1), ptarray_point_size(in));
690 }
691 
692 
693 int
695 {
696  if (!in)
697  {
698  lwerror("ptarray_is_closed_2d: called with null point array");
699  return 0;
700  }
701  if (in->npoints <= 1 ) return in->npoints; /* single-point are closed, empty not closed */
702 
703  return 0 == memcmp(getPoint_internal(in, 0), getPoint_internal(in, in->npoints-1), sizeof(POINT2D) );
704 }
705 
706 int
708 {
709  if (!in)
710  {
711  lwerror("ptarray_is_closed_3d: called with null point array");
712  return 0;
713  }
714  if (in->npoints <= 1 ) return in->npoints; /* single-point are closed, empty not closed */
715 
716  return 0 == memcmp(getPoint_internal(in, 0), getPoint_internal(in, in->npoints-1), sizeof(POINT3D) );
717 }
718 
719 int
721 {
722  if ( FLAGS_GET_Z(in->flags) )
723  return ptarray_is_closed_3d(in);
724  else
725  return ptarray_is_closed_2d(in);
726 }
727 
732 int
734 {
735  return ptarray_contains_point_partial(pa, pt, LW_TRUE, NULL);
736 }
737 
738 int
739 ptarray_contains_point_partial(const POINTARRAY *pa, const POINT2D *pt, int check_closed, int *winding_number)
740 {
741  int wn = 0;
742  int i;
743  double side;
744  const POINT2D *seg1;
745  const POINT2D *seg2;
746  double ymin, ymax;
747 
748  seg1 = getPoint2d_cp(pa, 0);
749  seg2 = getPoint2d_cp(pa, pa->npoints-1);
750  if ( check_closed && ! p2d_same(seg1, seg2) )
751  lwerror("ptarray_contains_point called on unclosed ring");
752 
753  for ( i=1; i < pa->npoints; i++ )
754  {
755  seg2 = getPoint2d_cp(pa, i);
756 
757  /* Zero length segments are ignored. */
758  if ( seg1->x == seg2->x && seg1->y == seg2->y )
759  {
760  seg1 = seg2;
761  continue;
762  }
763 
764  ymin = FP_MIN(seg1->y, seg2->y);
765  ymax = FP_MAX(seg1->y, seg2->y);
766 
767  /* Only test segments in our vertical range */
768  if ( pt->y > ymax || pt->y < ymin )
769  {
770  seg1 = seg2;
771  continue;
772  }
773 
774  side = lw_segment_side(seg1, seg2, pt);
775 
776  /*
777  * A point on the boundary of a ring is not contained.
778  * WAS: if (fabs(side) < 1e-12), see #852
779  */
780  if ( (side == 0) && lw_pt_in_seg(pt, seg1, seg2) )
781  {
782  return LW_BOUNDARY;
783  }
784 
785  /*
786  * If the point is to the left of the line, and it's rising,
787  * then the line is to the right of the point and
788  * circling counter-clockwise, so incremement.
789  */
790  if ( (side < 0) && (seg1->y <= pt->y) && (pt->y < seg2->y) )
791  {
792  wn++;
793  }
794 
795  /*
796  * If the point is to the right of the line, and it's falling,
797  * then the line is to the right of the point and circling
798  * clockwise, so decrement.
799  */
800  else if ( (side > 0) && (seg2->y <= pt->y) && (pt->y < seg1->y) )
801  {
802  wn--;
803  }
804 
805  seg1 = seg2;
806  }
807 
808  /* Sent out the winding number for calls that are building on this as a primitive */
809  if ( winding_number )
810  *winding_number = wn;
811 
812  /* Outside */
813  if (wn == 0)
814  {
815  return LW_OUTSIDE;
816  }
817 
818  /* Inside */
819  return LW_INSIDE;
820 }
821 
831 int
833 {
834  return ptarrayarc_contains_point_partial(pa, pt, LW_TRUE /* Check closed*/, NULL);
835 }
836 
837 int
838 ptarrayarc_contains_point_partial(const POINTARRAY *pa, const POINT2D *pt, int check_closed, int *winding_number)
839 {
840  int wn = 0;
841  int i, side;
842  const POINT2D *seg1;
843  const POINT2D *seg2;
844  const POINT2D *seg3;
845  GBOX gbox;
846 
847  /* Check for not an arc ring (always have odd # of points) */
848  if ( (pa->npoints % 2) == 0 )
849  {
850  lwerror("ptarrayarc_contains_point called with even number of points");
851  return LW_OUTSIDE;
852  }
853 
854  /* Check for not an arc ring (always have >= 3 points) */
855  if ( pa->npoints < 3 )
856  {
857  lwerror("ptarrayarc_contains_point called too-short pointarray");
858  return LW_OUTSIDE;
859  }
860 
861  /* Check for unclosed case */
862  seg1 = getPoint2d_cp(pa, 0);
863  seg3 = getPoint2d_cp(pa, pa->npoints-1);
864  if ( check_closed && ! p2d_same(seg1, seg3) )
865  {
866  lwerror("ptarrayarc_contains_point called on unclosed ring");
867  return LW_OUTSIDE;
868  }
869  /* OK, it's closed. Is it just one circle? */
870  else if ( p2d_same(seg1, seg3) && pa->npoints == 3 )
871  {
872  double radius, d;
873  POINT2D c;
874  seg2 = getPoint2d_cp(pa, 1);
875 
876  /* Wait, it's just a point, so it can't contain anything */
877  if ( lw_arc_is_pt(seg1, seg2, seg3) )
878  return LW_OUTSIDE;
879 
880  /* See if the point is within the circle radius */
881  radius = lw_arc_center(seg1, seg2, seg3, &c);
882  d = distance2d_pt_pt(pt, &c);
883  if ( FP_EQUALS(d, radius) )
884  return LW_BOUNDARY; /* Boundary of circle */
885  else if ( d < radius )
886  return LW_INSIDE; /* Inside circle */
887  else
888  return LW_OUTSIDE; /* Outside circle */
889  }
890  else if ( p2d_same(seg1, pt) || p2d_same(seg3, pt) )
891  {
892  return LW_BOUNDARY; /* Boundary case */
893  }
894 
895  /* Start on the ring */
896  seg1 = getPoint2d_cp(pa, 0);
897  for ( i=1; i < pa->npoints; i += 2 )
898  {
899  seg2 = getPoint2d_cp(pa, i);
900  seg3 = getPoint2d_cp(pa, i+1);
901 
902  /* Catch an easy boundary case */
903  if( p2d_same(seg3, pt) )
904  return LW_BOUNDARY;
905 
906  /* Skip arcs that have no size */
907  if ( lw_arc_is_pt(seg1, seg2, seg3) )
908  {
909  seg1 = seg3;
910  continue;
911  }
912 
913  /* Only test segments in our vertical range */
914  lw_arc_calculate_gbox_cartesian_2d(seg1, seg2, seg3, &gbox);
915  if ( pt->y > gbox.ymax || pt->y < gbox.ymin )
916  {
917  seg1 = seg3;
918  continue;
919  }
920 
921  /* Outside of horizontal range, and not between end points we also skip */
922  if ( (pt->x > gbox.xmax || pt->x < gbox.xmin) &&
923  (pt->y > FP_MAX(seg1->y, seg3->y) || pt->y < FP_MIN(seg1->y, seg3->y)) )
924  {
925  seg1 = seg3;
926  continue;
927  }
928 
929  side = lw_arc_side(seg1, seg2, seg3, pt);
930 
931  /* On the boundary */
932  if ( (side == 0) && lw_pt_in_arc(pt, seg1, seg2, seg3) )
933  {
934  return LW_BOUNDARY;
935  }
936 
937  /* Going "up"! Point to left of arc. */
938  if ( side < 0 && (seg1->y <= pt->y) && (pt->y < seg3->y) )
939  {
940  wn++;
941  }
942 
943  /* Going "down"! */
944  if ( side > 0 && (seg2->y <= pt->y) && (pt->y < seg1->y) )
945  {
946  wn--;
947  }
948 
949  /* Inside the arc! */
950  if ( pt->x <= gbox.xmax && pt->x >= gbox.xmin )
951  {
952  POINT2D C;
953  double radius = lw_arc_center(seg1, seg2, seg3, &C);
954  double d = distance2d_pt_pt(pt, &C);
955 
956  /* On the boundary! */
957  if ( d == radius )
958  return LW_BOUNDARY;
959 
960  /* Within the arc! */
961  if ( d < radius )
962  {
963  /* Left side, increment winding number */
964  if ( side < 0 )
965  wn++;
966  /* Right side, decrement winding number */
967  if ( side > 0 )
968  wn--;
969  }
970  }
971 
972  seg1 = seg3;
973  }
974 
975  /* Sent out the winding number for calls that are building on this as a primitive */
976  if ( winding_number )
977  *winding_number = wn;
978 
979  /* Outside */
980  if (wn == 0)
981  {
982  return LW_OUTSIDE;
983  }
984 
985  /* Inside */
986  return LW_INSIDE;
987 }
988 
994 double
996 {
997  const POINT2D *P1;
998  const POINT2D *P2;
999  const POINT2D *P3;
1000  double sum = 0.0;
1001  double x0, x, y1, y2;
1002  int i;
1003 
1004  if (! pa || pa->npoints < 3 )
1005  return 0.0;
1006 
1007  P1 = getPoint2d_cp(pa, 0);
1008  P2 = getPoint2d_cp(pa, 1);
1009  x0 = P1->x;
1010  for ( i = 2; i < pa->npoints; i++ )
1011  {
1012  P3 = getPoint2d_cp(pa, i);
1013  x = P2->x - x0;
1014  y1 = P3->y;
1015  y2 = P1->y;
1016  sum += x * (y2-y1);
1017 
1018  /* Move forwards! */
1019  P1 = P2;
1020  P2 = P3;
1021  }
1022  return sum / 2.0;
1023 }
1024 
1025 int
1027 {
1028  double area = 0;
1029  area = ptarray_signed_area(pa);
1030  if ( area > 0 ) return LW_FALSE;
1031  else return LW_TRUE;
1032 }
1033 
1034 POINTARRAY*
1035 ptarray_force_dims(const POINTARRAY *pa, int hasz, int hasm)
1036 {
1037  /* TODO handle zero-length point arrays */
1038  int i;
1039  int in_hasz = FLAGS_GET_Z(pa->flags);
1040  int in_hasm = FLAGS_GET_M(pa->flags);
1041  POINT4D pt;
1042  POINTARRAY *pa_out = ptarray_construct_empty(hasz, hasm, pa->npoints);
1043 
1044  for( i = 0; i < pa->npoints; i++ )
1045  {
1046  getPoint4d_p(pa, i, &pt);
1047  if( hasz && ! in_hasz )
1048  pt.z = 0.0;
1049  if( hasm && ! in_hasm )
1050  pt.m = 0.0;
1051  ptarray_append_point(pa_out, &pt, LW_TRUE);
1052  }
1053 
1054  return pa_out;
1055 }
1056 
1057 POINTARRAY *
1058 ptarray_substring(POINTARRAY *ipa, double from, double to, double tolerance)
1059 {
1060  POINTARRAY *dpa;
1061  POINT4D pt;
1062  POINT4D p1, p2;
1063  POINT4D *p1ptr=&p1; /* don't break strict-aliasing rule */
1064  POINT4D *p2ptr=&p2;
1065  int nsegs, i;
1066  double length, slength, tlength;
1067  int state = 0; /* 0=before, 1=inside */
1068 
1069  /*
1070  * Create a dynamic pointarray with an initial capacity
1071  * equal to full copy of input points
1072  */
1074 
1075  /* Compute total line length */
1076  length = ptarray_length_2d(ipa);
1077 
1078 
1079  LWDEBUGF(3, "Total length: %g", length);
1080 
1081 
1082  /* Get 'from' and 'to' lengths */
1083  from = length*from;
1084  to = length*to;
1085 
1086 
1087  LWDEBUGF(3, "From/To: %g/%g", from, to);
1088 
1089 
1090  tlength = 0;
1091  getPoint4d_p(ipa, 0, &p1);
1092  nsegs = ipa->npoints - 1;
1093  for ( i = 0; i < nsegs; i++ )
1094  {
1095  double dseg;
1096 
1097  getPoint4d_p(ipa, i+1, &p2);
1098 
1099 
1100  LWDEBUGF(3 ,"Segment %d: (%g,%g,%g,%g)-(%g,%g,%g,%g)",
1101  i, p1.x, p1.y, p1.z, p1.m, p2.x, p2.y, p2.z, p2.m);
1102 
1103 
1104  /* Find the length of this segment */
1105  slength = distance2d_pt_pt((POINT2D *)p1ptr, (POINT2D *)p2ptr);
1106 
1107  /*
1108  * We are before requested start.
1109  */
1110  if ( state == 0 ) /* before */
1111  {
1112 
1113  LWDEBUG(3, " Before start");
1114 
1115  if ( fabs ( from - ( tlength + slength ) ) <= tolerance )
1116  {
1117 
1118  LWDEBUG(3, " Second point is our start");
1119 
1120  /*
1121  * Second point is our start
1122  */
1123  ptarray_append_point(dpa, &p2, LW_FALSE);
1124  state=1; /* we're inside now */
1125  goto END;
1126  }
1127 
1128  else if ( fabs(from - tlength) <= tolerance )
1129  {
1130 
1131  LWDEBUG(3, " First point is our start");
1132 
1133  /*
1134  * First point is our start
1135  */
1136  ptarray_append_point(dpa, &p1, LW_FALSE);
1137 
1138  /*
1139  * We're inside now, but will check
1140  * 'to' point as well
1141  */
1142  state=1;
1143  }
1144 
1145  /*
1146  * Didn't reach the 'from' point,
1147  * nothing to do
1148  */
1149  else if ( from > tlength + slength ) goto END;
1150 
1151  else /* tlength < from < tlength+slength */
1152  {
1153 
1154  LWDEBUG(3, " Seg contains first point");
1155 
1156  /*
1157  * Our start is between first and
1158  * second point
1159  */
1160  dseg = (from - tlength) / slength;
1161 
1162  interpolate_point4d(&p1, &p2, &pt, dseg);
1163 
1164  ptarray_append_point(dpa, &pt, LW_FALSE);
1165 
1166  /*
1167  * We're inside now, but will check
1168  * 'to' point as well
1169  */
1170  state=1;
1171  }
1172  }
1173 
1174  if ( state == 1 ) /* inside */
1175  {
1176 
1177  LWDEBUG(3, " Inside");
1178 
1179  /*
1180  * 'to' point is our second point.
1181  */
1182  if ( fabs(to - ( tlength + slength ) ) <= tolerance )
1183  {
1184 
1185  LWDEBUG(3, " Second point is our end");
1186 
1187  ptarray_append_point(dpa, &p2, LW_FALSE);
1188  break; /* substring complete */
1189  }
1190 
1191  /*
1192  * 'to' point is our first point.
1193  * (should only happen if 'to' is 0)
1194  */
1195  else if ( fabs(to - tlength) <= tolerance )
1196  {
1197 
1198  LWDEBUG(3, " First point is our end");
1199 
1200  ptarray_append_point(dpa, &p1, LW_FALSE);
1201 
1202  break; /* substring complete */
1203  }
1204 
1205  /*
1206  * Didn't reach the 'end' point,
1207  * just copy second point
1208  */
1209  else if ( to > tlength + slength )
1210  {
1211  ptarray_append_point(dpa, &p2, LW_FALSE);
1212  goto END;
1213  }
1214 
1215  /*
1216  * 'to' point falls on this segment
1217  * Interpolate and break.
1218  */
1219  else if ( to < tlength + slength )
1220  {
1221 
1222  LWDEBUG(3, " Seg contains our end");
1223 
1224  dseg = (to - tlength) / slength;
1225  interpolate_point4d(&p1, &p2, &pt, dseg);
1226 
1227  ptarray_append_point(dpa, &pt, LW_FALSE);
1228 
1229  break;
1230  }
1231 
1232  else
1233  {
1234  LWDEBUG(3, "Unhandled case");
1235  }
1236  }
1237 
1238 
1239 END:
1240 
1241  tlength += slength;
1242  memcpy(&p1, &p2, sizeof(POINT4D));
1243  }
1244 
1245  LWDEBUGF(3, "Out of loop, ptarray has %d points", dpa->npoints);
1246 
1247  return dpa;
1248 }
1249 
1250 /*
1251  * Write into the *ret argument coordinates of the closes point on
1252  * the given segment to the reference input point.
1253  */
1254 void
1255 closest_point_on_segment(const POINT4D *p, const POINT4D *A, const POINT4D *B, POINT4D *ret)
1256 {
1257  double r;
1258 
1259  if ( FP_EQUALS(A->x, B->x) && FP_EQUALS(A->y, B->y) )
1260  {
1261  *ret = *A;
1262  return;
1263  }
1264 
1265  /*
1266  * We use comp.graphics.algorithms Frequently Asked Questions method
1267  *
1268  * (1) AC dot AB
1269  * r = ----------
1270  * ||AB||^2
1271  * r has the following meaning:
1272  * r=0 P = A
1273  * r=1 P = B
1274  * r<0 P is on the backward extension of AB
1275  * r>1 P is on the forward extension of AB
1276  * 0<r<1 P is interior to AB
1277  *
1278  */
1279  r = ( (p->x-A->x) * (B->x-A->x) + (p->y-A->y) * (B->y-A->y) )/( (B->x-A->x)*(B->x-A->x) +(B->y-A->y)*(B->y-A->y) );
1280 
1281  if (r<0)
1282  {
1283  *ret = *A;
1284  return;
1285  }
1286  if (r>1)
1287  {
1288  *ret = *B;
1289  return;
1290  }
1291 
1292  ret->x = A->x + ( (B->x - A->x) * r );
1293  ret->y = A->y + ( (B->y - A->y) * r );
1294  ret->z = A->z + ( (B->z - A->z) * r );
1295  ret->m = A->m + ( (B->m - A->m) * r );
1296 }
1297 
1298 /*
1299  * Given a point, returns the location of closest point on pointarray
1300  * and, optionally, it's actual distance from the point array.
1301  */
1302 double
1303 ptarray_locate_point(const POINTARRAY *pa, const POINT4D *p4d, double *mindistout, POINT4D *proj4d)
1304 {
1305  double mindist=-1;
1306  double tlen, plen;
1307  int t, seg=-1;
1308  POINT4D start4d, end4d, projtmp;
1309  POINT2D proj, p;
1310  const POINT2D *start = NULL, *end = NULL;
1311 
1312  /* Initialize our 2D copy of the input parameter */
1313  p.x = p4d->x;
1314  p.y = p4d->y;
1315 
1316  if ( ! proj4d ) proj4d = &projtmp;
1317 
1318  start = getPoint2d_cp(pa, 0);
1319 
1320  /* If the pointarray has only one point, the nearest point is */
1321  /* just that point */
1322  if ( pa->npoints == 1 )
1323  {
1324  getPoint4d_p(pa, 0, proj4d);
1325  if ( mindistout )
1326  *mindistout = distance2d_pt_pt(&p, start);
1327  return 0.0;
1328  }
1329 
1330  /* Loop through pointarray looking for nearest segment */
1331  for (t=1; t<pa->npoints; t++)
1332  {
1333  double dist;
1334  end = getPoint2d_cp(pa, t);
1335  dist = distance2d_pt_seg(&p, start, end);
1336 
1337  if (t==1 || dist < mindist )
1338  {
1339  mindist = dist;
1340  seg=t-1;
1341  }
1342 
1343  if ( mindist == 0 )
1344  {
1345  LWDEBUG(3, "Breaking on mindist=0");
1346  break;
1347  }
1348 
1349  start = end;
1350  }
1351 
1352  if ( mindistout ) *mindistout = mindist;
1353 
1354  LWDEBUGF(3, "Closest segment: %d", seg);
1355  LWDEBUGF(3, "mindist: %g", mindist);
1356 
1357  /*
1358  * We need to project the
1359  * point on the closest segment.
1360  */
1361  getPoint4d_p(pa, seg, &start4d);
1362  getPoint4d_p(pa, seg+1, &end4d);
1363  closest_point_on_segment(p4d, &start4d, &end4d, proj4d);
1364 
1365  /* Copy 4D values into 2D holder */
1366  proj.x = proj4d->x;
1367  proj.y = proj4d->y;
1368 
1369  LWDEBUGF(3, "Closest segment:%d, npoints:%d", seg, pa->npoints);
1370 
1371  /* For robustness, force 1 when closest point == endpoint */
1372  if ( (seg >= (pa->npoints-2)) && p2d_same(&proj, end) )
1373  {
1374  return 1.0;
1375  }
1376 
1377  LWDEBUGF(3, "Closest point on segment: %g,%g", proj.x, proj.y);
1378 
1379  tlen = ptarray_length_2d(pa);
1380 
1381  LWDEBUGF(3, "tlen %g", tlen);
1382 
1383  /* Location of any point on a zero-length line is 0 */
1384  /* See http://trac.osgeo.org/postgis/ticket/1772#comment:2 */
1385  if ( tlen == 0 ) return 0;
1386 
1387  plen=0;
1388  start = getPoint2d_cp(pa, 0);
1389  for (t=0; t<seg; t++, start=end)
1390  {
1391  end = getPoint2d_cp(pa, t+1);
1392  plen += distance2d_pt_pt(start, end);
1393 
1394  LWDEBUGF(4, "Segment %d made plen %g", t, plen);
1395  }
1396 
1397  plen+=distance2d_pt_pt(&proj, start);
1398 
1399  LWDEBUGF(3, "plen %g, tlen %g", plen, tlen);
1400 
1401  return plen/tlen;
1402 }
1403 
1413 void
1415 {
1416  int i;
1417  double x;
1418 
1419  for (i=0; i<pa->npoints; i++)
1420  {
1421  memcpy(&x, getPoint_internal(pa, i), sizeof(double));
1422  if ( x < 0 ) x+= 360;
1423  else if ( x > 180 ) x -= 360;
1424  memcpy(getPoint_internal(pa, i), &x, sizeof(double));
1425  }
1426 }
1427 
1428 
1429 /*
1430  * Returns a POINTARRAY with consecutive equal points
1431  * removed. Equality test on all dimensions of input.
1432  *
1433  * Always returns a newly allocated object.
1434  *
1435  */
1436 POINTARRAY *
1437 ptarray_remove_repeated_points_minpoints(const POINTARRAY *in, double tolerance, int minpoints)
1438 {
1439  POINTARRAY* out;
1440  size_t ptsize;
1441  size_t ipn, opn;
1442  const POINT2D *last_point, *this_point;
1443  double tolsq = tolerance * tolerance;
1444 
1445  if ( minpoints < 1 ) minpoints = 1;
1446 
1447  LWDEBUGF(3, "%s called", __func__);
1448 
1449  /* Single or zero point arrays can't have duplicates */
1450  if ( in->npoints < 3 ) return ptarray_clone_deep(in);
1451 
1452  ptsize = ptarray_point_size(in);
1453 
1454  LWDEBUGF(3, " ptsize: %d", ptsize);
1455 
1456  /* Allocate enough space for all points */
1457  out = ptarray_construct(FLAGS_GET_Z(in->flags),
1458  FLAGS_GET_M(in->flags), in->npoints);
1459 
1460  /* Now fill up the actual points (NOTE: could be optimized) */
1461 
1462  opn=1;
1463  /* Keep the first point */
1464  memcpy(getPoint_internal(out, 0), getPoint_internal(in, 0), ptsize);
1465  last_point = getPoint2d_cp(in, 0);
1466  LWDEBUGF(3, " first point copied, out points: %d", opn);
1467  for ( ipn = 1; ipn < in->npoints; ++ipn)
1468  {
1469  this_point = getPoint2d_cp(in, ipn);
1470  if ( ipn < in->npoints-minpoints+1 || opn >= minpoints ) /* need extra points to hit minponts */
1471  {
1472  if (
1473  (tolerance == 0 && memcmp(getPoint_internal(in, ipn-1), getPoint_internal(in, ipn), ptsize) == 0) || /* exact dupe */
1474  (tolerance > 0.0 && distance2d_sqr_pt_pt(last_point, this_point) <= tolsq) /* within the removal tolerance */
1475  ) continue;
1476  }
1477 
1478  /*
1479  * The point is different (see above) from the previous,
1480  * so we add it to output
1481  */
1482  memcpy(getPoint_internal(out, opn++), getPoint_internal(in, ipn), ptsize);
1483  last_point = this_point;
1484  LWDEBUGF(3, " Point %d differs from point %d. Out points: %d", ipn, ipn-1, opn);
1485  }
1486  /* Keep the last point */
1487  if ( memcmp(last_point, getPoint_internal(in, ipn-1), ptsize) != 0 )
1488  {
1489  memcpy(getPoint_internal(out, opn-1), getPoint_internal(in, ipn-1), ptsize);
1490  }
1491 
1492  LWDEBUGF(3, " in:%d out:%d", out->npoints, opn);
1493  out->npoints = opn;
1494 
1495  return out;
1496 }
1497 
1498 POINTARRAY *
1499 ptarray_remove_repeated_points(const POINTARRAY *in, double tolerance)
1500 {
1501  return ptarray_remove_repeated_points_minpoints(in, tolerance, 2);
1502 }
1503 
1504 static void
1505 ptarray_dp_findsplit(POINTARRAY *pts, int p1, int p2, int *split, double *dist)
1506 {
1507  int k;
1508  const POINT2D *pk, *pa, *pb;
1509  double tmp, d;
1510 
1511  LWDEBUG(4, "function called");
1512 
1513  *split = p1;
1514  d = -1;
1515 
1516  if (p1 + 1 < p2)
1517  {
1518 
1519  pa = getPoint2d_cp(pts, p1);
1520  pb = getPoint2d_cp(pts, p2);
1521 
1522  LWDEBUGF(4, "P%d(%f,%f) to P%d(%f,%f)",
1523  p1, pa->x, pa->y, p2, pb->x, pb->y);
1524 
1525  for (k=p1+1; k<p2; k++)
1526  {
1527  pk = getPoint2d_cp(pts, k);
1528 
1529  LWDEBUGF(4, "P%d(%f,%f)", k, pk->x, pk->y);
1530 
1531  /* distance computation */
1532  tmp = distance2d_sqr_pt_seg(pk, pa, pb);
1533 
1534  if (tmp > d)
1535  {
1536  d = tmp; /* record the maximum */
1537  *split = k;
1538 
1539  LWDEBUGF(4, "P%d is farthest (%g)", k, d);
1540  }
1541  }
1542  *dist = d;
1543 
1544  } /* length---should be redone if can == 0 */
1545  else
1546  {
1547  LWDEBUG(3, "segment too short, no split/no dist");
1548  *dist = -1;
1549  }
1550 
1551 }
1552 
1553 POINTARRAY *
1554 ptarray_simplify(POINTARRAY *inpts, double epsilon, unsigned int minpts)
1555 {
1556  int *stack; /* recursion stack */
1557  int sp=-1; /* recursion stack pointer */
1558  int p1, split;
1559  double dist;
1560  POINTARRAY *outpts;
1561  POINT4D pt;
1562 
1563  double eps_sqr = epsilon * epsilon;
1564 
1565  /* Allocate recursion stack */
1566  stack = lwalloc(sizeof(int)*inpts->npoints);
1567 
1568  p1 = 0;
1569  stack[++sp] = inpts->npoints-1;
1570 
1571  LWDEBUGF(2, "Input has %d pts and %d dims", inpts->npoints,
1572  FLAGS_NDIMS(inpts->flags));
1573 
1574  /* Allocate output POINTARRAY, and add first point. */
1575  outpts = ptarray_construct_empty(FLAGS_GET_Z(inpts->flags), FLAGS_GET_M(inpts->flags), inpts->npoints);
1576  getPoint4d_p(inpts, 0, &pt);
1577  ptarray_append_point(outpts, &pt, LW_FALSE);
1578 
1579  LWDEBUG(3, "Added P0 to simplified point array (size 1)");
1580 
1581  do
1582  {
1583 
1584  ptarray_dp_findsplit(inpts, p1, stack[sp], &split, &dist);
1585 
1586  LWDEBUGF(3, "Farthest point from P%d-P%d is P%d (dist. %g)", p1, stack[sp], split, dist);
1587 
1588  if (dist > eps_sqr || ( outpts->npoints+sp+1 < minpts && dist >= 0 ) )
1589  {
1590  LWDEBUGF(4, "Added P%d to stack (outpts:%d)", split, sp);
1591  stack[++sp] = split;
1592  }
1593  else
1594  {
1595  getPoint4d_p(inpts, stack[sp], &pt);
1596  LWDEBUGF(4, "npoints , minpoints %d %d", outpts->npoints, minpts);
1597  ptarray_append_point(outpts, &pt, LW_FALSE);
1598 
1599  LWDEBUGF(4, "Added P%d to simplified point array (size: %d)", stack[sp], outpts->npoints);
1600 
1601  p1 = stack[sp--];
1602  }
1603 
1604  LWDEBUGF(4, "stack pointer = %d", sp);
1605  }
1606  while (! (sp<0) );
1607 
1608  lwfree(stack);
1609  return outpts;
1610 }
1611 
1617 double
1619 {
1620  double dist = 0.0;
1621  int i;
1622  const POINT2D *a1;
1623  const POINT2D *a2;
1624  const POINT2D *a3;
1625 
1626  if ( pts->npoints % 2 != 1 )
1627  lwerror("arc point array with even number of points");
1628 
1629  a1 = getPoint2d_cp(pts, 0);
1630 
1631  for ( i=2; i < pts->npoints; i += 2 )
1632  {
1633  a2 = getPoint2d_cp(pts, i-1);
1634  a3 = getPoint2d_cp(pts, i);
1635  dist += lw_arc_length(a1, a2, a3);
1636  a1 = a3;
1637  }
1638  return dist;
1639 }
1640 
1644 double
1646 {
1647  double dist = 0.0;
1648  int i;
1649  const POINT2D *frm;
1650  const POINT2D *to;
1651 
1652  if ( pts->npoints < 2 ) return 0.0;
1653 
1654  frm = getPoint2d_cp(pts, 0);
1655 
1656  for ( i=1; i < pts->npoints; i++ )
1657  {
1658  to = getPoint2d_cp(pts, i);
1659 
1660  dist += sqrt( ((frm->x - to->x)*(frm->x - to->x)) +
1661  ((frm->y - to->y)*(frm->y - to->y)) );
1662 
1663  frm = to;
1664  }
1665  return dist;
1666 }
1667 
1672 double
1674 {
1675  double dist = 0.0;
1676  int i;
1677  POINT3DZ frm;
1678  POINT3DZ to;
1679 
1680  if ( pts->npoints < 2 ) return 0.0;
1681 
1682  /* compute 2d length if 3d is not available */
1683  if ( ! FLAGS_GET_Z(pts->flags) ) return ptarray_length_2d(pts);
1684 
1685  getPoint3dz_p(pts, 0, &frm);
1686  for ( i=1; i < pts->npoints; i++ )
1687  {
1688  getPoint3dz_p(pts, i, &to);
1689  dist += sqrt( ((frm.x - to.x)*(frm.x - to.x)) +
1690  ((frm.y - to.y)*(frm.y - to.y)) +
1691  ((frm.z - to.z)*(frm.z - to.z)) );
1692  frm = to;
1693  }
1694  return dist;
1695 }
1696 
1697 
1698 /*
1699  * Get a pointer to nth point of a POINTARRAY.
1700  *
1701  * Casting to returned pointer to POINT2D* should be safe,
1702  * as gserialized format always keeps the POINTARRAY pointer
1703  * aligned to double boundary.
1704  */
1705 uint8_t *
1707 {
1708  size_t size;
1709  uint8_t *ptr;
1710 
1711 #if PARANOIA_LEVEL > 0
1712  if ( pa == NULL )
1713  {
1714  lwerror("getPoint got NULL pointarray");
1715  return NULL;
1716  }
1717 
1718  LWDEBUGF(5, "(n=%d, pa.npoints=%d, pa.maxpoints=%d)",n,pa->npoints,pa->maxpoints);
1719 
1720  if ( ( n < 0 ) ||
1721  ( n > pa->npoints ) ||
1722  ( n >= pa->maxpoints ) )
1723  {
1724  lwerror("getPoint_internal called outside of ptarray range (n=%d, pa.npoints=%d, pa.maxpoints=%d)",n,pa->npoints,pa->maxpoints);
1725  return NULL; /*error */
1726  }
1727 #endif
1728 
1729  size = ptarray_point_size(pa);
1730 
1731  ptr = pa->serialized_pointlist + size * n;
1732  if ( FLAGS_NDIMS(pa->flags) == 2)
1733  {
1734  LWDEBUGF(5, "point = %g %g", *((double*)(ptr)), *((double*)(ptr+8)));
1735  }
1736  else if ( FLAGS_NDIMS(pa->flags) == 3)
1737  {
1738  LWDEBUGF(5, "point = %g %g %g", *((double*)(ptr)), *((double*)(ptr+8)), *((double*)(ptr+16)));
1739  }
1740  else if ( FLAGS_NDIMS(pa->flags) == 4)
1741  {
1742  LWDEBUGF(5, "point = %g %g %g %g", *((double*)(ptr)), *((double*)(ptr+8)), *((double*)(ptr+16)), *((double*)(ptr+24)));
1743  }
1744 
1745  return ptr;
1746 }
1747 
1748 
1752 void
1754 {
1755  int i;
1756  double x,y,z;
1757  POINT4D p4d;
1758 
1759  LWDEBUG(2, "lwgeom_affine_ptarray start");
1760 
1761  if ( FLAGS_GET_Z(pa->flags) )
1762  {
1763  LWDEBUG(3, " has z");
1764 
1765  for (i=0; i<pa->npoints; i++)
1766  {
1767  getPoint4d_p(pa, i, &p4d);
1768  x = p4d.x;
1769  y = p4d.y;
1770  z = p4d.z;
1771  p4d.x = a->afac * x + a->bfac * y + a->cfac * z + a->xoff;
1772  p4d.y = a->dfac * x + a->efac * y + a->ffac * z + a->yoff;
1773  p4d.z = a->gfac * x + a->hfac * y + a->ifac * z + a->zoff;
1774  ptarray_set_point4d(pa, i, &p4d);
1775 
1776  LWDEBUGF(3, " POINT %g %g %g => %g %g %g", x, y, z, p4d.x, p4d.y, p4d.z);
1777  }
1778  }
1779  else
1780  {
1781  LWDEBUG(3, " doesn't have z");
1782 
1783  for (i=0; i<pa->npoints; i++)
1784  {
1785  getPoint4d_p(pa, i, &p4d);
1786  x = p4d.x;
1787  y = p4d.y;
1788  p4d.x = a->afac * x + a->bfac * y + a->xoff;
1789  p4d.y = a->dfac * x + a->efac * y + a->yoff;
1790  ptarray_set_point4d(pa, i, &p4d);
1791 
1792  LWDEBUGF(3, " POINT %g %g => %g %g", x, y, p4d.x, p4d.y);
1793  }
1794  }
1795 
1796  LWDEBUG(3, "lwgeom_affine_ptarray end");
1797 
1798 }
1799 
1803 void
1805 {
1806  int i;
1807  POINT4D p4d;
1808 
1809  LWDEBUG(3, "ptarray_scale start");
1810 
1811  for (i=0; i<pa->npoints; ++i)
1812  {
1813  getPoint4d_p(pa, i, &p4d);
1814  p4d.x *= fact->x;
1815  p4d.y *= fact->y;
1816  p4d.z *= fact->z;
1817  p4d.m *= fact->m;
1818  ptarray_set_point4d(pa, i, &p4d);
1819  }
1820 
1821  LWDEBUG(3, "ptarray_scale end");
1822 
1823 }
1824 
1825 int
1827 {
1828  return getPoint4d_p(pa, 0, pt);
1829 }
1830 
1831 
1832 
1833 
1834 /*
1835  * Stick an array of points to the given gridspec.
1836  * Return "gridded" points in *outpts and their number in *outptsn.
1837  *
1838  * Two consecutive points falling on the same grid cell are collapsed
1839  * into one single point.
1840  *
1841  */
1842 POINTARRAY *
1843 ptarray_grid(const POINTARRAY *pa, const gridspec *grid)
1844 {
1845  POINT4D pt;
1846  int ipn; /* input point numbers */
1847  POINTARRAY *dpa;
1848 
1849  LWDEBUGF(2, "ptarray_grid called on %p", pa);
1850 
1852 
1853  for (ipn=0; ipn<pa->npoints; ++ipn)
1854  {
1855 
1856  getPoint4d_p(pa, ipn, &pt);
1857 
1858  if ( grid->xsize )
1859  pt.x = rint((pt.x - grid->ipx)/grid->xsize) *
1860  grid->xsize + grid->ipx;
1861 
1862  if ( grid->ysize )
1863  pt.y = rint((pt.y - grid->ipy)/grid->ysize) *
1864  grid->ysize + grid->ipy;
1865 
1866  if ( FLAGS_GET_Z(pa->flags) && grid->zsize )
1867  pt.z = rint((pt.z - grid->ipz)/grid->zsize) *
1868  grid->zsize + grid->ipz;
1869 
1870  if ( FLAGS_GET_M(pa->flags) && grid->msize )
1871  pt.m = rint((pt.m - grid->ipm)/grid->msize) *
1872  grid->msize + grid->ipm;
1873 
1874  ptarray_append_point(dpa, &pt, LW_FALSE);
1875 
1876  }
1877 
1878  return dpa;
1879 }
1880 
1881 int
1882 ptarray_npoints_in_rect(const POINTARRAY *pa, const GBOX *gbox)
1883 {
1884  const POINT2D *pt;
1885  int n = 0;
1886  int i;
1887  for ( i = 0; i < pa->npoints; i++ )
1888  {
1889  pt = getPoint2d_cp(pa, i);
1890  if ( gbox_contains_point2d(gbox, pt) )
1891  n++;
1892  }
1893  return n;
1894 }
1895 
1896 
void ptarray_set_point4d(POINTARRAY *pa, int n, const POINT4D *p4d)
Definition: lwgeom_api.c:549
double x
Definition: liblwgeom.h:351
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:1437
double z
Definition: liblwgeom.h:333
int ptarrayarc_contains_point_partial(const POINTARRAY *pa, const POINT2D *pt, int check_closed, int *winding_number)
Definition: ptarray.c:838
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:1058
double y
Definition: liblwgeom.h:333
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:351
uint8_t * serialized_pointlist
Definition: liblwgeom.h:365
double x
Definition: liblwgeom.h:333
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:603
void lwfree(void *mem)
Definition: lwutil.c:242
#define FLAGS_GET_READONLY(flags)
Definition: liblwgeom.h:143
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:1673
int npoints
Definition: liblwgeom.h:370
POINTARRAY * ptarray_clone_deep(const POINTARRAY *in)
Deep clone a pointarray (also clones serialized pointlist)
Definition: ptarray.c:634
POINTARRAY * ptarray_simplify(POINTARRAY *inpts, double epsilon, unsigned int minpts)
Definition: ptarray.c:1554
Datum area(PG_FUNCTION_ARGS)
double ptarray_locate_point(const POINTARRAY *pa, const POINT4D *p4d, double *mindistout, POINT4D *proj4d)
Definition: ptarray.c:1303
double xmax
Definition: liblwgeom.h:292
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:79
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:152
#define LWDEBUG(level, msg)
Definition: lwgeom_log.h:83
void ptarray_scale(POINTARRAY *pa, const POINT4D *fact)
Scale a pointarray.
Definition: ptarray.c:1804
int ptarray_is_closed_2d(const POINTARRAY *in)
Definition: ptarray.c:694
#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:1026
POINTARRAY * ptarray_removePoint(POINTARRAY *pa, uint32_t which)
Remove a point from a pointarray.
Definition: ptarray.c:562
double distance2d_pt_pt(const POINT2D *p1, const POINT2D *p2)
The old function nessecary for ptarray_segmentize2d in ptarray.c.
Definition: measures.c:2316
double zoff
Definition: liblwgeom.h:269
#define FP_MIN(A, B)
int ptarray_npoints_in_rect(const POINTARRAY *pa, const GBOX *gbox)
Definition: ptarray.c:1882
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:2341
void ptarray_longitude_shift(POINTARRAY *pa)
Longitude shift for a pointarray.
Definition: ptarray.c:1414
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:819
double ifac
Definition: liblwgeom.h:269
POINTARRAY * ptarray_clone(const POINTARRAY *in)
Clone a POINTARRAY object.
Definition: ptarray.c:658
double ffac
Definition: liblwgeom.h:269
double xoff
Definition: liblwgeom.h:269
double afac
Definition: liblwgeom.h:269
int ptarray_startpoint(const POINTARRAY *pa, POINT4D *pt)
Definition: ptarray.c:1826
#define LW_FAILURE
Definition: liblwgeom.h:78
double x
Definition: liblwgeom.h:327
int p2d_same(const POINT2D *p1, const POINT2D *p2)
Definition: lwalgorithm.c:49
double ymin
Definition: liblwgeom.h:293
POINTARRAY * ptarray_grid(const POINTARRAY *pa, const gridspec *grid)
Definition: ptarray.c:1843
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:332
int ptarray_contains_point_partial(const POINTARRAY *pa, const POINT2D *pt, int check_closed, int *winding_number)
Definition: ptarray.c:739
double xmin
Definition: liblwgeom.h:291
#define LW_FALSE
Definition: liblwgeom.h:76
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:485
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
uint8_t flags
Definition: liblwgeom.h:368
#define LW_TRUE
Return types for functions with status returns.
Definition: liblwgeom.h:75
void closest_point_on_segment(const POINT4D *p, const POINT4D *A, const POINT4D *B, POINT4D *ret)
Definition: ptarray.c:1255
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:720
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:269
double ymax
Definition: liblwgeom.h:294
enum LWORD_T LWORD
Ordinate names.
double y
Definition: liblwgeom.h:327
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:461
#define FLAGS_GET_Z(flags)
Macros for manipulating the &#39;flags&#39; byte.
Definition: liblwgeom.h:139
int ptarray_is_closed_3d(const POINTARRAY *in)
Definition: ptarray.c:707
double ptarray_length_2d(const POINTARRAY *pts)
Find the 2d length of the given POINTARRAY (even if it&#39;s 3d)
Definition: ptarray.c:1645
double dfac
Definition: liblwgeom.h:269
double z
Definition: liblwgeom.h:351
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:995
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:1618
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:235
double efac
Definition: liblwgeom.h:269
#define FLAGS_GET_M(flags)
Definition: liblwgeom.h:140
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:269
double distance2d_sqr_pt_pt(const POINT2D *p1, const POINT2D *p2)
Definition: measures.c:2326
int maxpoints
Definition: liblwgeom.h:371
#define LW_OUTSIDE
void * lwalloc(size_t size)
Definition: lwutil.c:227
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:420
double y
Definition: liblwgeom.h:351
int ptarrayarc_contains_point(const POINTARRAY *pa, const POINT2D *pt)
For POINTARRAYs representing CIRCULARSTRINGS.
Definition: ptarray.c:832
double distance2d_sqr_pt_seg(const POINT2D *p, const POINT2D *A, const POINT2D *B)
Definition: measures.c:2390
char ptarray_same(const POINTARRAY *pa1, const POINTARRAY *pa2)
Definition: ptarray.c:484
double gfac
Definition: liblwgeom.h:269
#define LWDEBUGF(level, msg,...)
Definition: lwgeom_log.h:88
#define FLAGS_NDIMS(flags)
Definition: liblwgeom.h:151
void ptarray_affine(POINTARRAY *pa, const AFFINE *a)
Affine transform a pointarray.
Definition: ptarray.c:1753
POINTARRAY * ptarray_force_dims(const POINTARRAY *pa, int hasz, int hasm)
Definition: ptarray.c:1035
uint8_t * getPoint_internal(const POINTARRAY *pa, int n)
Definition: ptarray.c:1706
int ptarray_is_closed(const POINTARRAY *in)
Check for ring closure using whatever dimensionality is declared on the pointarray.
Definition: ptarray.c:680
void lwerror(const char *fmt,...)
Write a notice out to the error handler.
Definition: lwutil.c:102
int gbox_contains_point2d(const GBOX *g, const POINT2D *p)
Definition: g_box.c:362
double hfac
Definition: liblwgeom.h:269
double bfac
Definition: liblwgeom.h:269
int getPoint4d_p(const POINTARRAY *pa, int n, POINT4D *point)
Definition: lwgeom_api.c:244
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:149
POINTARRAY * ptarray_remove_repeated_points(const POINTARRAY *in, double tolerance)
Definition: ptarray.c:1499
static void ptarray_dp_findsplit(POINTARRAY *pts, int p1, int p2, int *split, double *dist)
Definition: ptarray.c:1505
#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:733
Snap to grid.