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