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