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