PostGIS  3.2.2dev-r@@SVN_REVISION@@
lwgeom.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) 2001-2006 Refractions Research Inc.
22  * Copyright (C) 2017-2018 Daniel Baston <dbaston@gmail.com>
23  *
24  **********************************************************************/
25 
26 
27 #include <stdio.h>
28 #include <stdlib.h>
29 #include <stdarg.h>
30 
31 #include "liblwgeom_internal.h"
32 #include "lwgeom_log.h"
33 
34 #define out_stack_size 32
35 
37 void
39 {
40  LWCOLLECTION *coll;
41  uint32_t i;
42 
43  switch (lwgeom->type)
44  {
45  case POLYGONTYPE:
46  lwpoly_force_clockwise((LWPOLY *)lwgeom);
47  return;
48 
49  case TRIANGLETYPE:
51  return;
52 
53  /* Not handle POLYHEDRALSURFACE and TIN
54  as they are supposed to be well oriented */
55  case MULTIPOLYGONTYPE:
56  case COLLECTIONTYPE:
57  coll = (LWCOLLECTION *)lwgeom;
58  for (i=0; i<coll->ngeoms; i++)
59  lwgeom_force_clockwise(coll->geoms[i]);
60  return;
61  }
62 }
63 
65 int
67 {
68  switch (lwgeom->type)
69  {
70  case POLYGONTYPE:
71  return lwpoly_is_clockwise((LWPOLY *)lwgeom);
72 
73  case TRIANGLETYPE:
74  return lwtriangle_is_clockwise((LWTRIANGLE *)lwgeom);
75 
76  case MULTIPOLYGONTYPE:
77  case COLLECTIONTYPE:
78  {
79  uint32_t i;
80  LWCOLLECTION* coll = (LWCOLLECTION *)lwgeom;
81 
82  for (i=0; i < coll->ngeoms; i++)
83  if (!lwgeom_is_clockwise(coll->geoms[i]))
84  return LW_FALSE;
85  return LW_TRUE;
86  }
87  default:
88  return LW_TRUE;
89  return LW_FALSE;
90  }
91 }
92 
93 LWGEOM *
94 lwgeom_reverse(const LWGEOM *geom)
95 {
96  LWGEOM *geomout = lwgeom_clone_deep(geom);
97  lwgeom_reverse_in_place(geomout);
98  return geomout;
99 }
100 
102 void
104 {
105  uint32_t i;
106  LWCOLLECTION *col;
107  if (!geom)
108  return;
109 
110  switch (geom->type)
111  {
112  case MULTIPOINTTYPE:
113  case POINTTYPE:
114  {
115  return;
116  }
117  case TRIANGLETYPE:
118  case CIRCSTRINGTYPE:
119  case LINETYPE:
120  {
121  LWLINE *line = (LWLINE *)(geom);
123  return;
124  }
125  case POLYGONTYPE:
126  {
127  LWPOLY *poly = (LWPOLY *)(geom);
128  if (!poly->rings)
129  return;
130  uint32_t r;
131  for (r = 0; r < poly->nrings; r++)
133  return;
134  }
135  case MULTICURVETYPE:
136  case MULTILINETYPE:
137  case MULTIPOLYGONTYPE:
138  case MULTISURFACETYPE:
140  case TINTYPE:
141  case COLLECTIONTYPE:
142  case COMPOUNDTYPE:
143  case CURVEPOLYTYPE:
144  {
145  col = (LWCOLLECTION *)(geom);
146  if (!col->geoms)
147  return;
148  for (i=0; i<col->ngeoms; i++)
150  return;
151  }
152  default:
153  {
154  lwerror("%s: Unknown geometry type: %s", __func__, lwtype_name(geom->type));
155  return;
156  }
157 
158  }
159 }
160 
161 LWLINE *
162 lwgeom_as_lwline(const LWGEOM *lwgeom)
163 {
164  if ( lwgeom == NULL ) return NULL;
165  if ( lwgeom->type == LINETYPE )
166  return (LWLINE *)lwgeom;
167  else return NULL;
168 }
169 
170 LWCIRCSTRING *
172 {
173  if ( lwgeom == NULL ) return NULL;
174  if ( lwgeom->type == CIRCSTRINGTYPE )
175  return (LWCIRCSTRING *)lwgeom;
176  else return NULL;
177 }
178 
179 LWCOMPOUND *
181 {
182  if ( lwgeom == NULL ) return NULL;
183  if ( lwgeom->type == COMPOUNDTYPE )
184  return (LWCOMPOUND *)lwgeom;
185  else return NULL;
186 }
187 
188 LWCURVEPOLY *
190 {
191  if ( lwgeom == NULL ) return NULL;
192  if ( lwgeom->type == CURVEPOLYTYPE )
193  return (LWCURVEPOLY *)lwgeom;
194  else return NULL;
195 }
196 
197 LWPOLY *
198 lwgeom_as_lwpoly(const LWGEOM *lwgeom)
199 {
200  if ( lwgeom == NULL ) return NULL;
201  if ( lwgeom->type == POLYGONTYPE )
202  return (LWPOLY *)lwgeom;
203  else return NULL;
204 }
205 
206 LWTRIANGLE *
208 {
209  if ( lwgeom == NULL ) return NULL;
210  if ( lwgeom->type == TRIANGLETYPE )
211  return (LWTRIANGLE *)lwgeom;
212  else return NULL;
213 }
214 
215 LWCOLLECTION *
217 {
218  if ( lwgeom == NULL ) return NULL;
219  if ( lwgeom_is_collection(lwgeom) )
220  return (LWCOLLECTION*)lwgeom;
221  else return NULL;
222 }
223 
224 LWMPOINT *
226 {
227  if ( lwgeom == NULL ) return NULL;
228  if ( lwgeom->type == MULTIPOINTTYPE )
229  return (LWMPOINT *)lwgeom;
230  else return NULL;
231 }
232 
233 LWMLINE *
234 lwgeom_as_lwmline(const LWGEOM *lwgeom)
235 {
236  if ( lwgeom == NULL ) return NULL;
237  if ( lwgeom->type == MULTILINETYPE )
238  return (LWMLINE *)lwgeom;
239  else return NULL;
240 }
241 
242 LWMPOLY *
243 lwgeom_as_lwmpoly(const LWGEOM *lwgeom)
244 {
245  if ( lwgeom == NULL ) return NULL;
246  if ( lwgeom->type == MULTIPOLYGONTYPE )
247  return (LWMPOLY *)lwgeom;
248  else return NULL;
249 }
250 
251 LWPSURFACE *
253 {
254  if ( lwgeom->type == POLYHEDRALSURFACETYPE )
255  return (LWPSURFACE *)lwgeom;
256  else return NULL;
257 }
258 
259 LWTIN *
260 lwgeom_as_lwtin(const LWGEOM *lwgeom)
261 {
262  if ( lwgeom->type == TINTYPE )
263  return (LWTIN *)lwgeom;
264  else return NULL;
265 }
266 
268 {
269  return (LWGEOM *)obj;
270 }
271 
273 {
274  return (LWGEOM *)obj;
275 }
276 
278 {
279  if ( obj == NULL ) return NULL;
280  return (LWGEOM *)obj;
281 }
283 {
284  if ( obj == NULL ) return NULL;
285  return (LWGEOM *)obj;
286 }
288 {
289  if ( obj == NULL ) return NULL;
290  return (LWGEOM *)obj;
291 }
293 {
294  if ( obj == NULL ) return NULL;
295  return (LWGEOM *)obj;
296 }
298 {
299  if ( obj == NULL ) return NULL;
300  return (LWGEOM *)obj;
301 }
303 {
304  if ( obj == NULL ) return NULL;
305  return (LWGEOM *)obj;
306 }
308 {
309  if ( obj == NULL ) return NULL;
310  return (LWGEOM *)obj;
311 }
313 {
314  if ( obj == NULL ) return NULL;
315  return (LWGEOM *)obj;
316 }
318 {
319  if ( obj == NULL ) return NULL;
320  return (LWGEOM *)obj;
321 }
323 {
324  if ( obj == NULL ) return NULL;
325  return (LWGEOM *)obj;
326 }
328 {
329  if ( obj == NULL ) return NULL;
330  return (LWGEOM *)obj;
331 }
332 
333 
337 uint8_t MULTITYPE[NUMTYPES] =
338 {
339  0,
340  MULTIPOINTTYPE, /* 1 */
341  MULTILINETYPE, /* 2 */
342  MULTIPOLYGONTYPE, /* 3 */
343  0,0,0,0,
344  MULTICURVETYPE, /* 8 */
345  MULTICURVETYPE, /* 9 */
346  MULTISURFACETYPE, /* 10 */
347  POLYHEDRALSURFACETYPE, /* 11 */
348  0, 0,
349  TINTYPE, /* 14 */
350  0
351 };
352 
353 uint8_t lwtype_multitype(uint8_t type)
354 {
355  if (type > 15) return 0;
356  return MULTITYPE[type];
357 }
358 
362 LWGEOM *
363 lwgeom_as_multi(const LWGEOM *lwgeom)
364 {
365  LWGEOM **ogeoms;
366  LWGEOM *ogeom = NULL;
367  GBOX *box = NULL;
368  int type;
369 
370  type = lwgeom->type;
371 
372  if ( ! MULTITYPE[type] ) return lwgeom_clone(lwgeom);
373 
374  if( lwgeom_is_empty(lwgeom) )
375  {
377  MULTITYPE[type],
378  lwgeom->srid,
379  FLAGS_GET_Z(lwgeom->flags),
380  FLAGS_GET_M(lwgeom->flags)
381  );
382  }
383  else
384  {
385  ogeoms = lwalloc(sizeof(LWGEOM*));
386  ogeoms[0] = lwgeom_clone(lwgeom);
387 
388  /* Sub-geometries are not allowed to have bboxes or SRIDs, move the bbox to the collection */
389  box = ogeoms[0]->bbox;
390  ogeoms[0]->bbox = NULL;
391  ogeoms[0]->srid = SRID_UNKNOWN;
392 
393  ogeom = (LWGEOM *)lwcollection_construct(MULTITYPE[type], lwgeom->srid, box, 1, ogeoms);
394  }
395 
396  return ogeom;
397 }
398 
402 LWGEOM *
403 lwgeom_as_curve(const LWGEOM *lwgeom)
404 {
405  LWGEOM *ogeom;
406  int type = lwgeom->type;
407  /*
408  int hasz = FLAGS_GET_Z(lwgeom->flags);
409  int hasm = FLAGS_GET_M(lwgeom->flags);
410  int32_t srid = lwgeom->srid;
411  */
412 
413  switch(type)
414  {
415  case LINETYPE:
416  /* turn to COMPOUNDCURVE */
417  ogeom = (LWGEOM*)lwcompound_construct_from_lwline((LWLINE*)lwgeom);
418  break;
419  case POLYGONTYPE:
421  break;
422  case MULTILINETYPE:
423  /* turn to MULTICURVE */
424  ogeom = lwgeom_clone(lwgeom);
425  ogeom->type = MULTICURVETYPE;
426  break;
427  case MULTIPOLYGONTYPE:
428  /* turn to MULTISURFACE */
429  ogeom = lwgeom_clone(lwgeom);
430  ogeom->type = MULTISURFACETYPE;
431  break;
432  case COLLECTIONTYPE:
433  default:
434  ogeom = lwgeom_clone(lwgeom);
435  break;
436  }
437 
438  /* TODO: copy bbox from input geom ? */
439 
440  return ogeom;
441 }
442 
443 
450 void
452 {
453  if ( ! lwgeom )
454  lwerror("lwgeom_release: someone called on 0x0");
455 
456  LWDEBUGF(3, "releasing type %s", lwtype_name(lwgeom->type));
457 
458  /* Drop bounding box (always a copy) */
459  if ( lwgeom->bbox )
460  {
461  LWDEBUGF(3, "lwgeom_release: releasing bbox. %p", lwgeom->bbox);
462  lwfree(lwgeom->bbox);
463  }
464  lwfree(lwgeom);
465 
466 }
467 
468 
469 /* @brief Clone LWGEOM object. Serialized point lists are not copied.
470  *
471  * @see ptarray_clone
472  */
473 LWGEOM *
474 lwgeom_clone(const LWGEOM *lwgeom)
475 {
476  LWDEBUGF(2, "lwgeom_clone called with %p, %s",
477  lwgeom, lwtype_name(lwgeom->type));
478 
479  switch (lwgeom->type)
480  {
481  case POINTTYPE:
482  return (LWGEOM *)lwpoint_clone((LWPOINT *)lwgeom);
483  case LINETYPE:
484  return (LWGEOM *)lwline_clone((LWLINE *)lwgeom);
485  case CIRCSTRINGTYPE:
486  return (LWGEOM *)lwcircstring_clone((LWCIRCSTRING *)lwgeom);
487  case POLYGONTYPE:
488  return (LWGEOM *)lwpoly_clone((LWPOLY *)lwgeom);
489  case TRIANGLETYPE:
490  return (LWGEOM *)lwtriangle_clone((LWTRIANGLE *)lwgeom);
491  case COMPOUNDTYPE:
492  case CURVEPOLYTYPE:
493  case MULTICURVETYPE:
494  case MULTISURFACETYPE:
495  case MULTIPOINTTYPE:
496  case MULTILINETYPE:
497  case MULTIPOLYGONTYPE:
499  case TINTYPE:
500  case COLLECTIONTYPE:
501  return (LWGEOM *)lwcollection_clone((LWCOLLECTION *)lwgeom);
502  default:
503  lwerror("lwgeom_clone: Unknown geometry type: %s", lwtype_name(lwgeom->type));
504  return NULL;
505  }
506 }
507 
511 LWGEOM *
512 lwgeom_clone_deep(const LWGEOM *lwgeom)
513 {
514  LWDEBUGF(2, "lwgeom_clone called with %p, %s",
515  lwgeom, lwtype_name(lwgeom->type));
516 
517  switch (lwgeom->type)
518  {
519  case POINTTYPE:
520  case LINETYPE:
521  case CIRCSTRINGTYPE:
522  case TRIANGLETYPE:
523  return (LWGEOM *)lwline_clone_deep((LWLINE *)lwgeom);
524  case POLYGONTYPE:
525  return (LWGEOM *)lwpoly_clone_deep((LWPOLY *)lwgeom);
526  case COMPOUNDTYPE:
527  case CURVEPOLYTYPE:
528  case MULTICURVETYPE:
529  case MULTISURFACETYPE:
530  case MULTIPOINTTYPE:
531  case MULTILINETYPE:
532  case MULTIPOLYGONTYPE:
534  case TINTYPE:
535  case COLLECTIONTYPE:
536  return (LWGEOM *)lwcollection_clone_deep((LWCOLLECTION *)lwgeom);
537  default:
538  lwerror("lwgeom_clone_deep: Unknown geometry type: %s", lwtype_name(lwgeom->type));
539  return NULL;
540  }
541 }
542 
543 
547 char*
548 lwgeom_to_ewkt(const LWGEOM *lwgeom)
549 {
550  char* wkt = NULL;
551  size_t wkt_size = 0;
552 
553  wkt = lwgeom_to_wkt(lwgeom, WKT_EXTENDED, 12, &wkt_size);
554 
555  if ( ! wkt )
556  {
557  lwerror("Error writing geom %p to WKT", lwgeom);
558  }
559 
560  return wkt;
561 }
562 
573 char
574 lwgeom_same(const LWGEOM *lwgeom1, const LWGEOM *lwgeom2)
575 {
576  LWDEBUGF(2, "lwgeom_same(%s, %s) called",
577  lwtype_name(lwgeom1->type),
578  lwtype_name(lwgeom2->type));
579 
580  if ( lwgeom1->type != lwgeom2->type )
581  {
582  LWDEBUG(3, " type differ");
583 
584  return LW_FALSE;
585  }
586 
587  if ( FLAGS_GET_ZM(lwgeom1->flags) != FLAGS_GET_ZM(lwgeom2->flags) )
588  {
589  LWDEBUG(3, " ZM flags differ");
590 
591  return LW_FALSE;
592  }
593 
594  /* Check boxes if both already computed */
595  if ( lwgeom1->bbox && lwgeom2->bbox )
596  {
597  /*lwnotice("bbox1:%p, bbox2:%p", lwgeom1->bbox, lwgeom2->bbox);*/
598  if ( ! gbox_same(lwgeom1->bbox, lwgeom2->bbox) )
599  {
600  LWDEBUG(3, " bounding boxes differ");
601 
602  return LW_FALSE;
603  }
604  }
605 
606  /* geoms have same type, invoke type-specific function */
607  switch (lwgeom1->type)
608  {
609  case POINTTYPE:
610  return lwpoint_same((LWPOINT *)lwgeom1,
611  (LWPOINT *)lwgeom2);
612  case LINETYPE:
613  return lwline_same((LWLINE *)lwgeom1,
614  (LWLINE *)lwgeom2);
615  case POLYGONTYPE:
616  return lwpoly_same((LWPOLY *)lwgeom1,
617  (LWPOLY *)lwgeom2);
618  case TRIANGLETYPE:
619  return lwtriangle_same((LWTRIANGLE *)lwgeom1,
620  (LWTRIANGLE *)lwgeom2);
621  case CIRCSTRINGTYPE:
622  return lwcircstring_same((LWCIRCSTRING *)lwgeom1,
623  (LWCIRCSTRING *)lwgeom2);
624  case MULTIPOINTTYPE:
625  case MULTILINETYPE:
626  case MULTIPOLYGONTYPE:
627  case MULTICURVETYPE:
628  case MULTISURFACETYPE:
629  case COMPOUNDTYPE:
630  case CURVEPOLYTYPE:
632  case TINTYPE:
633  case COLLECTIONTYPE:
634  return lwcollection_same((LWCOLLECTION *)lwgeom1,
635  (LWCOLLECTION *)lwgeom2);
636  default:
637  lwerror("lwgeom_same: unsupported geometry type: %s",
638  lwtype_name(lwgeom1->type));
639  return LW_FALSE;
640  }
641 
642 }
643 
644 int
645 lwpoint_inside_circle(const LWPOINT *p, double cx, double cy, double rad)
646 {
647  const POINT2D *pt;
648  POINT2D center;
649 
650  if ( ! p || ! p->point )
651  return LW_FALSE;
652 
653  pt = getPoint2d_cp(p->point, 0);
654 
655  center.x = cx;
656  center.y = cy;
657 
658  if ( distance2d_pt_pt(pt, &center) < rad )
659  return LW_TRUE;
660 
661  return LW_FALSE;
662 }
663 
664 void
666 {
667  if ( lwgeom->bbox ) lwfree(lwgeom->bbox);
668  lwgeom->bbox = NULL;
669  FLAGS_SET_BBOX(lwgeom->flags, 0);
670 }
671 
677 void
679 {
680  /* an empty LWGEOM has no bbox */
681  if ( lwgeom_is_empty(lwgeom) ) return;
682 
683  if ( lwgeom->bbox ) return;
684  FLAGS_SET_BBOX(lwgeom->flags, 1);
685  lwgeom->bbox = gbox_new(lwgeom->flags);
686  lwgeom_calculate_gbox(lwgeom, lwgeom->bbox);
687 }
688 
689 void
691 {
692  lwgeom_drop_bbox(lwgeom);
693  lwgeom_add_bbox(lwgeom);
694 }
695 
696 void
698 {
699  if ( lwgeom_is_empty(lwgeom) ) return;
700 
701  FLAGS_SET_BBOX(lwgeom->flags, 1);
702 
703  if ( ! ( gbox || lwgeom->bbox ) )
704  {
705  lwgeom->bbox = gbox_new(lwgeom->flags);
706  lwgeom_calculate_gbox(lwgeom, lwgeom->bbox);
707  }
708  else if ( gbox && ! lwgeom->bbox )
709  {
710  lwgeom->bbox = gbox_clone(gbox);
711  }
712 
713  if ( lwgeom_is_collection(lwgeom) )
714  {
715  uint32_t i;
716  LWCOLLECTION *lwcol = (LWCOLLECTION*)lwgeom;
717 
718  for ( i = 0; i < lwcol->ngeoms; i++ )
719  {
720  lwgeom_add_bbox_deep(lwcol->geoms[i], lwgeom->bbox);
721  }
722  }
723 }
724 
725 const GBOX *
727 {
728  /* add it if not already there */
729  lwgeom_add_bbox((LWGEOM *)lwg);
730  return lwg->bbox;
731 }
732 
733 
738 int lwgeom_calculate_gbox(const LWGEOM *lwgeom, GBOX *gbox)
739 {
740  gbox->flags = lwgeom->flags;
741  if( FLAGS_GET_GEODETIC(lwgeom->flags) )
742  return lwgeom_calculate_gbox_geodetic(lwgeom, gbox);
743  else
744  return lwgeom_calculate_gbox_cartesian(lwgeom, gbox);
745 }
746 
747 void
749 {
750  lwgeom->srid = SRID_UNKNOWN; /* TODO: To be changed to SRID_UNKNOWN */
751 }
752 
753 LWGEOM *
754 lwgeom_segmentize2d(const LWGEOM *lwgeom, double dist)
755 {
756  switch (lwgeom->type)
757  {
758  case LINETYPE:
759  return (LWGEOM *)lwline_segmentize2d((LWLINE *)lwgeom,
760  dist);
761  case POLYGONTYPE:
762  return (LWGEOM *)lwpoly_segmentize2d((LWPOLY *)lwgeom,
763  dist);
764  case MULTILINETYPE:
765  case MULTIPOLYGONTYPE:
766  case COLLECTIONTYPE:
768  (LWCOLLECTION *)lwgeom, dist);
769 
770  default:
771  return lwgeom_clone(lwgeom);
772  }
773 }
774 
775 LWGEOM*
777 {
778  return lwgeom_force_dims(geom, 0, 0, 0, 0);
779 }
780 
781 LWGEOM*
782 lwgeom_force_3dz(const LWGEOM *geom, double zval)
783 {
784  return lwgeom_force_dims(geom, 1, 0, zval, 0);
785 }
786 
787 LWGEOM*
788 lwgeom_force_3dm(const LWGEOM *geom, double mval)
789 {
790  return lwgeom_force_dims(geom, 0, 1, 0, mval);
791 }
792 
793 LWGEOM*
794 lwgeom_force_4d(const LWGEOM *geom, double zval, double mval)
795 {
796  return lwgeom_force_dims(geom, 1, 1, zval, mval);
797 }
798 
799 LWGEOM*
800 lwgeom_force_dims(const LWGEOM *geom, int hasz, int hasm, double zval, double mval)
801 {
802  if (!geom)
803  return NULL;
804  switch(geom->type)
805  {
806  case POINTTYPE:
807  return lwpoint_as_lwgeom(lwpoint_force_dims((LWPOINT*)geom, hasz, hasm, zval, mval));
808  case CIRCSTRINGTYPE:
809  case LINETYPE:
810  case TRIANGLETYPE:
811  return lwline_as_lwgeom(lwline_force_dims((LWLINE*)geom, hasz, hasm, zval, mval));
812  case POLYGONTYPE:
813  return lwpoly_as_lwgeom(lwpoly_force_dims((LWPOLY*)geom, hasz, hasm, zval, mval));
814  case COMPOUNDTYPE:
815  case CURVEPOLYTYPE:
816  case MULTICURVETYPE:
817  case MULTISURFACETYPE:
818  case MULTIPOINTTYPE:
819  case MULTILINETYPE:
820  case MULTIPOLYGONTYPE:
822  case TINTYPE:
823  case COLLECTIONTYPE:
824  return lwcollection_as_lwgeom(lwcollection_force_dims((LWCOLLECTION*)geom, hasz, hasm, zval, mval));
825  default:
826  lwerror("lwgeom_force_2d: unsupported geom type: %s", lwtype_name(geom->type));
827  return NULL;
828  }
829 }
830 
831 LWGEOM*
832 lwgeom_force_sfs(LWGEOM *geom, int version)
833 {
834  LWCOLLECTION *col;
835  uint32_t i;
836  LWGEOM *g;
837 
838  /* SFS 1.2 version */
839  if (version == 120)
840  {
841  switch(geom->type)
842  {
843  /* SQL/MM types */
844  case CIRCSTRINGTYPE:
845  case COMPOUNDTYPE:
846  case CURVEPOLYTYPE:
847  case MULTICURVETYPE:
848  case MULTISURFACETYPE:
849  return lwgeom_stroke(geom, 32);
850 
851  case COLLECTIONTYPE:
852  col = (LWCOLLECTION*)geom;
853  for ( i = 0; i < col->ngeoms; i++ )
854  col->geoms[i] = lwgeom_force_sfs((LWGEOM*)col->geoms[i], version);
855 
856  return lwcollection_as_lwgeom((LWCOLLECTION*)geom);
857 
858  default:
859  return (LWGEOM *)geom;
860  }
861  }
862 
863 
864  /* SFS 1.1 version */
865  switch(geom->type)
866  {
867  /* SQL/MM types */
868  case CIRCSTRINGTYPE:
869  case COMPOUNDTYPE:
870  case CURVEPOLYTYPE:
871  case MULTICURVETYPE:
872  case MULTISURFACETYPE:
873  return lwgeom_stroke(geom, 32);
874 
875  /* SFS 1.2 types */
876  case TRIANGLETYPE:
877  g = lwpoly_as_lwgeom(lwpoly_from_lwlines((LWLINE*)geom, 0, NULL));
878  lwgeom_free(geom);
879  return g;
880 
881  case TINTYPE:
882  col = (LWCOLLECTION*) geom;
883  for ( i = 0; i < col->ngeoms; i++ )
884  {
885  g = lwpoly_as_lwgeom(lwpoly_from_lwlines((LWLINE*)col->geoms[i], 0, NULL));
886  lwgeom_free(col->geoms[i]);
887  col->geoms[i] = g;
888  }
889  col->type = COLLECTIONTYPE;
890  return lwmpoly_as_lwgeom((LWMPOLY*)geom);
891 
893  geom->type = COLLECTIONTYPE;
894  return (LWGEOM *)geom;
895 
896  /* Collection */
897  case COLLECTIONTYPE:
898  col = (LWCOLLECTION*)geom;
899  for ( i = 0; i < col->ngeoms; i++ )
900  col->geoms[i] = lwgeom_force_sfs((LWGEOM*)col->geoms[i], version);
901 
902  return lwcollection_as_lwgeom((LWCOLLECTION*)geom);
903 
904  default:
905  return (LWGEOM *)geom;
906  }
907 }
908 
909 int32_t
911 {
912  if ( ! geom ) return SRID_UNKNOWN;
913  return geom->srid;
914 }
915 
916 int
917 lwgeom_has_z(const LWGEOM *geom)
918 {
919  if ( ! geom ) return LW_FALSE;
920  return FLAGS_GET_Z(geom->flags);
921 }
922 
923 int
924 lwgeom_has_m(const LWGEOM *geom)
925 {
926  if ( ! geom ) return LW_FALSE;
927  return FLAGS_GET_M(geom->flags);
928 }
929 
930 int
932 {
933  if ( ! geom ) return LW_FALSE;
934  return FLAGS_GET_SOLID(geom->flags);
935 }
936 
937 int
938 lwgeom_ndims(const LWGEOM *geom)
939 {
940  if ( ! geom ) return 0;
941  return FLAGS_NDIMS(geom->flags);
942 }
943 
944 
945 
946 void
948 {
949  LWPOINT *pt;
950  LWLINE *ln;
951  LWPOLY *ply;
952  LWCOLLECTION *col;
953  uint32_t i;
954 
956  if ( geom->bbox )
958 
959  switch(geom->type)
960  {
961  case POINTTYPE:
962  pt = (LWPOINT*)geom;
963  if ( pt->point )
965  break;
966  case LINETYPE:
967  ln = (LWLINE*)geom;
968  if ( ln->points )
970  break;
971  case POLYGONTYPE:
972  ply = (LWPOLY*)geom;
973  for ( i = 0; i < ply->nrings; i++ )
974  FLAGS_SET_GEODETIC(ply->rings[i]->flags, value);
975  break;
976  case MULTIPOINTTYPE:
977  case MULTILINETYPE:
978  case MULTIPOLYGONTYPE:
979  case COLLECTIONTYPE:
980  col = (LWCOLLECTION*)geom;
981  for ( i = 0; i < col->ngeoms; i++ )
982  lwgeom_set_geodetic(col->geoms[i], value);
983  break;
984  default:
985  lwerror("lwgeom_set_geodetic: unsupported geom type: %s", lwtype_name(geom->type));
986  return;
987  }
988 }
989 
990 void
992 {
993  uint32_t i;
994  switch (lwgeom->type)
995  {
996  LWPOINT *point;
997  LWLINE *line;
998  LWPOLY *poly;
999  LWTRIANGLE *triangle;
1000  LWCOLLECTION *coll;
1001 
1002  case POINTTYPE:
1003  point = (LWPOINT *)lwgeom;
1005  return;
1006  case LINETYPE:
1007  line = (LWLINE *)lwgeom;
1009  return;
1010  case POLYGONTYPE:
1011  poly = (LWPOLY *)lwgeom;
1012  for (i=0; i<poly->nrings; i++)
1013  ptarray_longitude_shift(poly->rings[i]);
1014  return;
1015  case TRIANGLETYPE:
1016  triangle = (LWTRIANGLE *)lwgeom;
1017  ptarray_longitude_shift(triangle->points);
1018  return;
1019  case MULTIPOINTTYPE:
1020  case MULTILINETYPE:
1021  case MULTIPOLYGONTYPE:
1022  case POLYHEDRALSURFACETYPE:
1023  case TINTYPE:
1024  case COLLECTIONTYPE:
1025  coll = (LWCOLLECTION *)lwgeom;
1026  for (i=0; i<coll->ngeoms; i++)
1027  lwgeom_longitude_shift(coll->geoms[i]);
1028  return;
1029  default:
1030  lwerror("lwgeom_longitude_shift: unsupported geom type: %s",
1031  lwtype_name(lwgeom->type));
1032  }
1033 }
1034 
1035 int
1037 {
1038  int type = geom->type;
1039 
1040  if( lwgeom_is_empty(geom) )
1041  return LW_FALSE;
1042 
1043  /* Test linear types for closure */
1044  switch (type)
1045  {
1046  case LINETYPE:
1047  return lwline_is_closed((LWLINE*)geom);
1048  case POLYGONTYPE:
1049  return lwpoly_is_closed((LWPOLY*)geom);
1050  case CIRCSTRINGTYPE:
1051  return lwcircstring_is_closed((LWCIRCSTRING*)geom);
1052  case COMPOUNDTYPE:
1053  return lwcompound_is_closed((LWCOMPOUND*)geom);
1054  case TINTYPE:
1055  return lwtin_is_closed((LWTIN*)geom);
1056  case POLYHEDRALSURFACETYPE:
1057  return lwpsurface_is_closed((LWPSURFACE*)geom);
1058  }
1059 
1060  /* Recurse into collections and see if anything is not closed */
1061  if ( lwgeom_is_collection(geom) )
1062  {
1063  LWCOLLECTION *col = lwgeom_as_lwcollection(geom);
1064  uint32_t i;
1065  int closed;
1066  for ( i = 0; i < col->ngeoms; i++ )
1067  {
1068  closed = lwgeom_is_closed(col->geoms[i]);
1069  if ( ! closed )
1070  return LW_FALSE;
1071  }
1072  return LW_TRUE;
1073  }
1074 
1075  /* All non-linear non-collection types we will call closed */
1076  return LW_TRUE;
1077 }
1078 
1079 int
1081 {
1082  if( ! geom ) return LW_FALSE;
1083  return lwtype_is_collection(geom->type);
1084 }
1085 
1087 int
1089 {
1090  switch (type)
1091  {
1092  case MULTIPOINTTYPE:
1093  case MULTILINETYPE:
1094  case MULTIPOLYGONTYPE:
1095  case COLLECTIONTYPE:
1096  case CURVEPOLYTYPE:
1097  case COMPOUNDTYPE:
1098  case MULTICURVETYPE:
1099  case MULTISURFACETYPE:
1100  case POLYHEDRALSURFACETYPE:
1101  case TINTYPE:
1102  return LW_TRUE;
1103  break;
1104 
1105  default:
1106  return LW_FALSE;
1107  }
1108 }
1109 
1113 uint32_t
1115 {
1116  switch (type)
1117  {
1118  case POINTTYPE:
1119  return MULTIPOINTTYPE;
1120  case LINETYPE:
1121  return MULTILINETYPE;
1122  case POLYGONTYPE:
1123  return MULTIPOLYGONTYPE;
1124  case CIRCSTRINGTYPE:
1125  return MULTICURVETYPE;
1126  case COMPOUNDTYPE:
1127  return MULTICURVETYPE;
1128  case CURVEPOLYTYPE:
1129  return MULTISURFACETYPE;
1130  case TRIANGLETYPE:
1131  return TINTYPE;
1132  default:
1133  return COLLECTIONTYPE;
1134  }
1135 }
1136 
1137 
1138 void lwgeom_free(LWGEOM *lwgeom)
1139 {
1140 
1141  /* There's nothing here to free... */
1142  if( ! lwgeom ) return;
1143 
1144  LWDEBUGF(5,"freeing a %s",lwtype_name(lwgeom->type));
1145 
1146  switch (lwgeom->type)
1147  {
1148  case POINTTYPE:
1149  lwpoint_free((LWPOINT *)lwgeom);
1150  break;
1151  case LINETYPE:
1152  lwline_free((LWLINE *)lwgeom);
1153  break;
1154  case POLYGONTYPE:
1155  lwpoly_free((LWPOLY *)lwgeom);
1156  break;
1157  case CIRCSTRINGTYPE:
1158  lwcircstring_free((LWCIRCSTRING *)lwgeom);
1159  break;
1160  case TRIANGLETYPE:
1161  lwtriangle_free((LWTRIANGLE *)lwgeom);
1162  break;
1163  case MULTIPOINTTYPE:
1164  lwmpoint_free((LWMPOINT *)lwgeom);
1165  break;
1166  case MULTILINETYPE:
1167  lwmline_free((LWMLINE *)lwgeom);
1168  break;
1169  case MULTIPOLYGONTYPE:
1170  lwmpoly_free((LWMPOLY *)lwgeom);
1171  break;
1172  case POLYHEDRALSURFACETYPE:
1173  lwpsurface_free((LWPSURFACE *)lwgeom);
1174  break;
1175  case TINTYPE:
1176  lwtin_free((LWTIN *)lwgeom);
1177  break;
1178  case CURVEPOLYTYPE:
1179  case COMPOUNDTYPE:
1180  case MULTICURVETYPE:
1181  case MULTISURFACETYPE:
1182  case COLLECTIONTYPE:
1183  lwcollection_free((LWCOLLECTION *)lwgeom);
1184  break;
1185  default:
1186  lwerror("lwgeom_free called with unknown type (%d) %s", lwgeom->type, lwtype_name(lwgeom->type));
1187  }
1188  return;
1189 }
1190 
1191 int lwgeom_needs_bbox(const LWGEOM *geom)
1192 {
1193  assert(geom);
1194  if ( geom->type == POINTTYPE )
1195  {
1196  return LW_FALSE;
1197  }
1198  else if ( geom->type == LINETYPE )
1199  {
1200  if ( lwgeom_count_vertices(geom) <= 2 )
1201  return LW_FALSE;
1202  else
1203  return LW_TRUE;
1204  }
1205  else if ( geom->type == MULTIPOINTTYPE )
1206  {
1207  if ( ((LWCOLLECTION*)geom)->ngeoms == 1 )
1208  return LW_FALSE;
1209  else
1210  return LW_TRUE;
1211  }
1212  else if ( geom->type == MULTILINETYPE )
1213  {
1214  if ( ((LWCOLLECTION*)geom)->ngeoms == 1 && lwgeom_count_vertices(geom) <= 2 )
1215  return LW_FALSE;
1216  else
1217  return LW_TRUE;
1218  }
1219  else
1220  {
1221  return LW_TRUE;
1222  }
1223 }
1224 
1229 uint32_t lwgeom_count_vertices(const LWGEOM *geom)
1230 {
1231  int result = 0;
1232 
1233  /* Null? Zero. */
1234  if( ! geom ) return 0;
1235 
1236  LWDEBUGF(4, "lwgeom_count_vertices got type %s",
1237  lwtype_name(geom->type));
1238 
1239  /* Empty? Zero. */
1240  if( lwgeom_is_empty(geom) ) return 0;
1241 
1242  switch (geom->type)
1243  {
1244  case POINTTYPE:
1245  result = 1;
1246  break;
1247  case TRIANGLETYPE:
1248  case CIRCSTRINGTYPE:
1249  case LINETYPE:
1250  result = lwline_count_vertices((LWLINE *)geom);
1251  break;
1252  case POLYGONTYPE:
1253  result = lwpoly_count_vertices((LWPOLY *)geom);
1254  break;
1255  case COMPOUNDTYPE:
1256  case CURVEPOLYTYPE:
1257  case MULTICURVETYPE:
1258  case MULTISURFACETYPE:
1259  case MULTIPOINTTYPE:
1260  case MULTILINETYPE:
1261  case MULTIPOLYGONTYPE:
1262  case POLYHEDRALSURFACETYPE:
1263  case TINTYPE:
1264  case COLLECTIONTYPE:
1266  break;
1267  default:
1268  lwerror("%s: unsupported input geometry type: %s",
1269  __func__, lwtype_name(geom->type));
1270  break;
1271  }
1272  LWDEBUGF(3, "counted %d vertices", result);
1273  return result;
1274 }
1275 
1281 int lwgeom_dimension(const LWGEOM *geom)
1282 {
1283 
1284  /* Null? Zero. */
1285  if( ! geom ) return -1;
1286 
1287  LWDEBUGF(4, "lwgeom_dimension got type %s",
1288  lwtype_name(geom->type));
1289 
1290  /* Empty? Zero. */
1291  /* if( lwgeom_is_empty(geom) ) return 0; */
1292 
1293  switch (geom->type)
1294  {
1295  case POINTTYPE:
1296  case MULTIPOINTTYPE:
1297  return 0;
1298  case CIRCSTRINGTYPE:
1299  case LINETYPE:
1300  case COMPOUNDTYPE:
1301  case MULTICURVETYPE:
1302  case MULTILINETYPE:
1303  return 1;
1304  case TRIANGLETYPE:
1305  case POLYGONTYPE:
1306  case CURVEPOLYTYPE:
1307  case MULTISURFACETYPE:
1308  case MULTIPOLYGONTYPE:
1309  case TINTYPE:
1310  return 2;
1311  case POLYHEDRALSURFACETYPE:
1312  {
1313  /* A closed polyhedral surface contains a volume. */
1314  int closed = lwpsurface_is_closed((LWPSURFACE*)geom);
1315  return ( closed ? 3 : 2 );
1316  }
1317  case COLLECTIONTYPE:
1318  {
1319  int maxdim = 0;
1320  uint32_t i;
1321  LWCOLLECTION *col = (LWCOLLECTION*)geom;
1322  for( i = 0; i < col->ngeoms; i++ )
1323  {
1324  int dim = lwgeom_dimension(col->geoms[i]);
1325  maxdim = ( dim > maxdim ? dim : maxdim );
1326  }
1327  return maxdim;
1328  }
1329  default:
1330  lwerror("%s: unsupported input geometry type: %s",
1331  __func__, lwtype_name(geom->type));
1332  }
1333  return -1;
1334 }
1335 
1339 uint32_t lwgeom_count_rings(const LWGEOM *geom)
1340 {
1341  int result = 0;
1342 
1343  /* Null? Empty? Zero. */
1344  if( ! geom || lwgeom_is_empty(geom) )
1345  return 0;
1346 
1347  switch (geom->type)
1348  {
1349  case POINTTYPE:
1350  case CIRCSTRINGTYPE:
1351  case COMPOUNDTYPE:
1352  case MULTICURVETYPE:
1353  case MULTIPOINTTYPE:
1354  case MULTILINETYPE:
1355  case LINETYPE:
1356  result = 0;
1357  break;
1358  case TRIANGLETYPE:
1359  result = 1;
1360  break;
1361  case POLYGONTYPE:
1362  result = ((LWPOLY *)geom)->nrings;
1363  break;
1364  case CURVEPOLYTYPE:
1365  result = ((LWCURVEPOLY *)geom)->nrings;
1366  break;
1367  case MULTISURFACETYPE:
1368  case MULTIPOLYGONTYPE:
1369  case POLYHEDRALSURFACETYPE:
1370  case TINTYPE:
1371  case COLLECTIONTYPE:
1372  {
1373  LWCOLLECTION *col = (LWCOLLECTION*)geom;
1374  uint32_t i = 0;
1375  for( i = 0; i < col->ngeoms; i++ )
1376  result += lwgeom_count_rings(col->geoms[i]);
1377  break;
1378  }
1379  default:
1380  lwerror("lwgeom_count_rings: unsupported input geometry type: %s", lwtype_name(geom->type));
1381  break;
1382  }
1383  LWDEBUGF(3, "counted %d rings", result);
1384  return result;
1385 }
1386 
1387 int lwgeom_has_srid(const LWGEOM *geom)
1388 {
1389  if ( geom->srid != SRID_UNKNOWN )
1390  return LW_TRUE;
1391 
1392  return LW_FALSE;
1393 }
1394 
1395 
1397 {
1398  uint32_t i;
1399  int dimensionality = 0;
1400  for ( i = 0; i < col->ngeoms; i++ )
1401  {
1402  int d = lwgeom_dimensionality(col->geoms[i]);
1403  if ( d > dimensionality )
1404  dimensionality = d;
1405  }
1406  return dimensionality;
1407 }
1408 
1409 extern int lwgeom_dimensionality(const LWGEOM *geom)
1410 {
1411  int dim;
1412 
1413  LWDEBUGF(3, "lwgeom_dimensionality got type %s",
1414  lwtype_name(geom->type));
1415 
1416  switch (geom->type)
1417  {
1418  case POINTTYPE:
1419  case MULTIPOINTTYPE:
1420  return 0;
1421  break;
1422  case LINETYPE:
1423  case CIRCSTRINGTYPE:
1424  case MULTILINETYPE:
1425  case COMPOUNDTYPE:
1426  case MULTICURVETYPE:
1427  return 1;
1428  break;
1429  case POLYGONTYPE:
1430  case TRIANGLETYPE:
1431  case CURVEPOLYTYPE:
1432  case MULTIPOLYGONTYPE:
1433  case MULTISURFACETYPE:
1434  return 2;
1435  break;
1436 
1437  case POLYHEDRALSURFACETYPE:
1438  case TINTYPE:
1439  dim = lwgeom_is_closed(geom)?3:2;
1440  return dim;
1441  break;
1442 
1443  case COLLECTIONTYPE:
1444  return lwcollection_dimensionality((const LWCOLLECTION *)geom);
1445  break;
1446  default:
1447  lwerror("lwgeom_dimensionality: unsupported input geometry type: %s",
1448  lwtype_name(geom->type));
1449  break;
1450  }
1451  return 0;
1452 }
1453 
1454 extern LWGEOM* lwgeom_remove_repeated_points(const LWGEOM *in, double tolerance)
1455 {
1456  LWGEOM *out = lwgeom_clone_deep(in);
1458  return out;
1459 }
1460 
1462 {
1463  LWCOLLECTION *col;
1464  LWPOLY *poly;
1465  uint32_t i;
1466 
1467  if ( (!in) || lwgeom_is_empty(in) ) return;
1468 
1469  /* TODO: check for lwgeom NOT having the specified dimension ? */
1470 
1471  LWDEBUGF(4, "lwgeom_flip_coordinates, got type: %s",
1472  lwtype_name(in->type));
1473 
1474  switch (in->type)
1475  {
1476  case POINTTYPE:
1477  ptarray_swap_ordinates(lwgeom_as_lwpoint(in)->point, o1, o2);
1478  break;
1479 
1480  case LINETYPE:
1481  ptarray_swap_ordinates(lwgeom_as_lwline(in)->points, o1, o2);
1482  break;
1483 
1484  case CIRCSTRINGTYPE:
1485  ptarray_swap_ordinates(lwgeom_as_lwcircstring(in)->points, o1, o2);
1486  break;
1487 
1488  case POLYGONTYPE:
1489  poly = (LWPOLY *) in;
1490  for (i=0; i<poly->nrings; i++)
1491  {
1492  ptarray_swap_ordinates(poly->rings[i], o1, o2);
1493  }
1494  break;
1495 
1496  case TRIANGLETYPE:
1497  ptarray_swap_ordinates(lwgeom_as_lwtriangle(in)->points, o1, o2);
1498  break;
1499 
1500  case MULTIPOINTTYPE:
1501  case MULTILINETYPE:
1502  case MULTIPOLYGONTYPE:
1503  case COLLECTIONTYPE:
1504  case COMPOUNDTYPE:
1505  case CURVEPOLYTYPE:
1506  case MULTISURFACETYPE:
1507  case MULTICURVETYPE:
1508  case POLYHEDRALSURFACETYPE:
1509  case TINTYPE:
1510  col = (LWCOLLECTION *) in;
1511  for (i=0; i<col->ngeoms; i++)
1512  {
1513  lwgeom_swap_ordinates(col->geoms[i], o1, o2);
1514  }
1515  break;
1516 
1517  default:
1518  lwerror("lwgeom_swap_ordinates: unsupported geometry type: %s",
1519  lwtype_name(in->type));
1520  return;
1521  }
1522 
1523  /* only refresh bbox if X or Y changed */
1524  if ( in->bbox && (o1 < 2 || o2 < 2) )
1525  {
1526  lwgeom_refresh_bbox(in);
1527  }
1528 }
1529 
1530 void lwgeom_set_srid(LWGEOM *geom, int32_t srid)
1531 {
1532  uint32_t i;
1533 
1534  LWDEBUGF(4,"entered with srid=%d",srid);
1535 
1536  geom->srid = srid;
1537 
1538  if ( lwgeom_is_collection(geom) )
1539  {
1540  /* All the children are set to the same SRID value */
1541  LWCOLLECTION *col = lwgeom_as_lwcollection(geom);
1542  for ( i = 0; i < col->ngeoms; i++ )
1543  {
1544  lwgeom_set_srid(col->geoms[i], srid);
1545  }
1546  }
1547 }
1548 
1549 
1550 /**************************************************************/
1551 
1552 static int
1553 cmp_point_x(const void *pa, const void *pb)
1554 {
1555  LWPOINT *p1 = *((LWPOINT **)pa);
1556  LWPOINT *p2 = *((LWPOINT **)pb);
1557 
1558  const POINT2D *pt1 = getPoint2d_cp(p1->point, 0);
1559  const POINT2D *pt2 = getPoint2d_cp(p2->point, 0);
1560 
1561  if (pt1 && pt2)
1562  return (pt1->x > pt2->x) ? 1 : ((pt1->x < pt2->x) ? -1 : 0);
1563 
1564  if (pt1) return -1;
1565  if (pt2) return 1;
1566  return 0;
1567 }
1568 
1569 static int
1570 cmp_point_y(const void *pa, const void *pb)
1571 {
1572  LWPOINT *p1 = *((LWPOINT **)pa);
1573  LWPOINT *p2 = *((LWPOINT **)pb);
1574 
1575  const POINT2D *pt1 = getPoint2d_cp(p1->point, 0);
1576  const POINT2D *pt2 = getPoint2d_cp(p2->point, 0);
1577 
1578  if (pt1 && pt2)
1579  return (pt1->y > pt2->y) ? 1 : ((pt1->y < pt2->y) ? -1 : 0);
1580 
1581  if (pt1) return -1;
1582  if (pt2) return 1;
1583  return 0;
1584 }
1585 
1586 int
1588 {
1589  int geometry_modified = LW_FALSE;
1590  switch (geom->type)
1591  {
1592  /* No-op! Cannot remove points */
1593  case POINTTYPE:
1594  case TRIANGLETYPE:
1595  return geometry_modified;
1596  case LINETYPE: {
1597  LWLINE *g = (LWLINE *)(geom);
1598  POINTARRAY *pa = g->points;
1599  uint32_t npoints = pa->npoints;
1600  ptarray_remove_repeated_points_in_place(pa, tolerance, 2);
1601  geometry_modified = npoints != pa->npoints;
1602  /* Invalid input, discard the collapsed line */
1603  if (pa->npoints < 2)
1604  {
1605  pa->npoints = 0;
1606  geometry_modified = LW_TRUE;
1607  }
1608  break;
1609  }
1610  case POLYGONTYPE: {
1611  uint32_t j = 0;
1612  LWPOLY *g = (LWPOLY *)(geom);
1613  for (uint32_t i = 0; i < g->nrings; i++)
1614  {
1615  POINTARRAY *pa = g->rings[i];
1616  uint32_t npoints = pa->npoints;
1617  ptarray_remove_repeated_points_in_place(pa, tolerance, 4);
1618  geometry_modified |= npoints != pa->npoints;
1619  /* Drop collapsed rings */
1620  if (pa->npoints < 4)
1621  {
1622  geometry_modified = LW_TRUE;
1623  ptarray_free(pa);
1624  continue;
1625  }
1626  g->rings[j++] = pa;
1627  }
1628  /* Update ring count */
1629  g->nrings = j;
1630  break;
1631  }
1632  case MULTIPOINTTYPE: {
1633  double tolsq = tolerance * tolerance;
1634  LWMPOINT *mpt = (LWMPOINT *)geom;
1635 
1636  for (uint8_t dim = 0; dim < 2; dim++)
1637  {
1638  /* sort by y, then by x - this way the result is sorted by x */
1639  qsort(mpt->geoms, mpt->ngeoms, sizeof(LWPOINT *), dim ? cmp_point_x : cmp_point_y);
1640  for (uint32_t i = 0; i < mpt->ngeoms; i++)
1641  {
1642  if (!mpt->geoms[i])
1643  continue;
1644 
1645  const POINT2D *pti = getPoint2d_cp(mpt->geoms[i]->point, 0);
1646  if (!pti) continue;
1647 
1648  /* check upcoming points if they're within tolerance of current one */
1649  for (uint32_t j = i + 1; j < mpt->ngeoms; j++)
1650  {
1651  if (!mpt->geoms[j])
1652  continue;
1653 
1654  const POINT2D *ptj = getPoint2d_cp(mpt->geoms[j]->point, 0);
1655  if (!ptj) continue;
1656 
1657  /* check that the point is in the strip of tolerance around the point */
1658  if ((dim ? ptj->x - pti->x : ptj->y - pti->y) > tolerance)
1659  break;
1660 
1661  /* remove any upcoming point that is within tolerance circle */
1662  if (distance2d_sqr_pt_pt(pti, ptj) <= tolsq)
1663  {
1664  lwpoint_free(mpt->geoms[j]);
1665  mpt->geoms[j] = NULL;
1666  geometry_modified = LW_TRUE;
1667  }
1668  }
1669  }
1670 
1671  /* null out empties */
1672  for (uint32_t j = 0; j < mpt->ngeoms; j++)
1673  {
1674  if (mpt->geoms[j] && lwpoint_is_empty(mpt->geoms[j]))
1675  {
1676  lwpoint_free(mpt->geoms[j]);
1677  mpt->geoms[j] = NULL;
1678  geometry_modified = LW_TRUE;
1679  }
1680  }
1681 
1682  /* compactify array of points */
1683  uint32_t i = 0;
1684  for (uint32_t j = 0; j < mpt->ngeoms; j++)
1685  if (mpt->geoms[j])
1686  mpt->geoms[i++] = mpt->geoms[j];
1687  mpt->ngeoms = i;
1688  }
1689 
1690  break;
1691  }
1692 
1693  case CIRCSTRINGTYPE:
1694  /* Dunno how to handle these, will return untouched */
1695  return geometry_modified;
1696 
1697  /* Can process most multi* types as generic collection */
1698  case MULTILINETYPE:
1699  case MULTIPOLYGONTYPE:
1700  case TINTYPE:
1701  case COLLECTIONTYPE:
1702  /* Curve types we mostly ignore, but allow the linear */
1703  /* portions to be processed by recursing into them */
1704  case MULTICURVETYPE:
1705  case CURVEPOLYTYPE:
1706  case MULTISURFACETYPE:
1707  case COMPOUNDTYPE: {
1708  uint32_t i, j = 0;
1709  LWCOLLECTION *col = (LWCOLLECTION *)(geom);
1710  for (i = 0; i < col->ngeoms; i++)
1711  {
1712  LWGEOM *g = col->geoms[i];
1713  if (!g)
1714  continue;
1715  geometry_modified |= lwgeom_remove_repeated_points_in_place(g, tolerance);
1716  /* Drop zero'ed out geometries */
1717  if (lwgeom_is_empty(g))
1718  {
1719  lwgeom_free(g);
1720  continue;
1721  }
1722  col->geoms[j++] = g;
1723  }
1724  /* Update geometry count */
1725  col->ngeoms = j;
1726  break;
1727  }
1728  default: {
1729  lwerror("%s: unsupported geometry type: %s", __func__, lwtype_name(geom->type));
1730  break;
1731  }
1732  }
1733 
1734  if (geometry_modified)
1735  lwgeom_drop_bbox(geom);
1736  return geometry_modified;
1737 }
1738 
1739 
1740 /**************************************************************/
1741 
1742 int
1743 lwgeom_simplify_in_place(LWGEOM *geom, double epsilon, int preserve_collapsed)
1744 {
1745  int modified = LW_FALSE;
1746  switch (geom->type)
1747  {
1748  /* No-op! Cannot simplify points or triangles */
1749  case POINTTYPE:
1750  return modified;
1751  case TRIANGLETYPE:
1752  {
1753  if (preserve_collapsed)
1754  return modified;
1755  LWTRIANGLE *t = lwgeom_as_lwtriangle(geom);
1756  POINTARRAY *pa = t->points;
1757  ptarray_simplify_in_place(pa, epsilon, 0);
1758  if (pa->npoints < 3)
1759  {
1760  pa->npoints = 0;
1761  modified = LW_TRUE;
1762  }
1763  break;
1764  }
1765  case LINETYPE:
1766  {
1767  LWLINE *g = (LWLINE*)(geom);
1768  POINTARRAY *pa = g->points;
1769  uint32_t in_npoints = pa->npoints;
1770  ptarray_simplify_in_place(pa, epsilon, 2);
1771  modified = in_npoints != pa->npoints;
1772  /* Invalid output */
1773  if (pa->npoints == 1 && pa->maxpoints > 1)
1774  {
1775  /* Use first point as last point */
1776  if (preserve_collapsed)
1777  {
1778  pa->npoints = 2;
1779  ptarray_copy_point(pa, 0, 1);
1780  }
1781  /* Finish the collapse process */
1782  else
1783  {
1784  pa->npoints = 0;
1785  }
1786  }
1787  /* Duped output, force collapse */
1788  if (pa->npoints == 2 && !preserve_collapsed)
1789  {
1790  if (p2d_same(getPoint2d_cp(pa, 0), getPoint2d_cp(pa, 1)))
1791  pa->npoints = 0;
1792  }
1793  break;
1794  }
1795  case POLYGONTYPE:
1796  {
1797  uint32_t i, j = 0;
1798  LWPOLY *g = (LWPOLY*)(geom);
1799  for (i = 0; i < g->nrings; i++)
1800  {
1801  POINTARRAY *pa = g->rings[i];
1802  /* Only stop collapse on first ring */
1803  int minpoints = (preserve_collapsed && i == 0) ? 4 : 0;
1804  /* Skip zero'ed out rings */
1805  if(!pa)
1806  continue;
1807  uint32_t in_npoints = pa->npoints;
1808  ptarray_simplify_in_place(pa, epsilon, minpoints);
1809  modified |= in_npoints != pa->npoints;
1810  /* Drop collapsed rings */
1811  if(pa->npoints < 4)
1812  {
1813  if (i == 0)
1814  {
1815  /* If the outter ring is dropped, all can be dropped */
1816  for (i = 0; i < g->nrings; i++)
1817  {
1818  pa = g->rings[i];
1819  ptarray_free(pa);
1820  }
1821  break;
1822  }
1823  else
1824  {
1825  /* Drop this inner ring only */
1826  ptarray_free(pa);
1827  continue;
1828  }
1829  }
1830  g->rings[j++] = pa;
1831  }
1832  /* Update ring count */
1833  g->nrings = j;
1834  break;
1835  }
1836  /* Can process all multi* types as generic collection */
1837  case MULTIPOINTTYPE:
1838  case MULTILINETYPE:
1839  case MULTIPOLYGONTYPE:
1840  case TINTYPE:
1841  case COLLECTIONTYPE:
1842  {
1843  uint32_t i, j = 0;
1844  LWCOLLECTION *col = (LWCOLLECTION*)geom;
1845  for (i = 0; i < col->ngeoms; i++)
1846  {
1847  LWGEOM *g = col->geoms[i];
1848  if (!g) continue;
1849  modified |= lwgeom_simplify_in_place(g, epsilon, preserve_collapsed);
1850  /* Drop zero'ed out geometries */
1851  if(lwgeom_is_empty(g))
1852  {
1853  lwgeom_free(g);
1854  continue;
1855  }
1856  col->geoms[j++] = g;
1857  }
1858  /* Update geometry count */
1859  col->ngeoms = j;
1860  break;
1861  }
1862  default:
1863  {
1864  lwerror("%s: unsupported geometry type: %s", __func__, lwtype_name(geom->type));
1865  break;
1866  }
1867  }
1868 
1869  if (modified)
1870  {
1871  lwgeom_drop_bbox(geom);
1872  }
1873  return modified;
1874 }
1875 
1876 LWGEOM* lwgeom_simplify(const LWGEOM *igeom, double dist, int preserve_collapsed)
1877 {
1878  LWGEOM *lwgeom_out = lwgeom_clone_deep(igeom);
1879  lwgeom_simplify_in_place(lwgeom_out, dist, preserve_collapsed);
1880  if (lwgeom_is_empty(lwgeom_out))
1881  {
1882  lwgeom_free(lwgeom_out);
1883  return NULL;
1884  }
1885  return lwgeom_out;
1886 }
1887 
1888 /**************************************************************/
1889 
1890 
1891 double lwgeom_area(const LWGEOM *geom)
1892 {
1893  int type = geom->type;
1894 
1895  if ( type == POLYGONTYPE )
1896  return lwpoly_area((LWPOLY*)geom);
1897  else if ( type == CURVEPOLYTYPE )
1898  return lwcurvepoly_area((LWCURVEPOLY*)geom);
1899  else if (type == TRIANGLETYPE )
1900  return lwtriangle_area((LWTRIANGLE*)geom);
1901  else if ( lwgeom_is_collection(geom) )
1902  {
1903  double area = 0.0;
1904  uint32_t i;
1905  LWCOLLECTION *col = (LWCOLLECTION*)geom;
1906  for ( i = 0; i < col->ngeoms; i++ )
1907  area += lwgeom_area(col->geoms[i]);
1908  return area;
1909  }
1910  else
1911  return 0.0;
1912 }
1913 
1914 double lwgeom_perimeter(const LWGEOM *geom)
1915 {
1916  int type = geom->type;
1917  if ( type == POLYGONTYPE )
1918  return lwpoly_perimeter((LWPOLY*)geom);
1919  else if ( type == CURVEPOLYTYPE )
1920  return lwcurvepoly_perimeter((LWCURVEPOLY*)geom);
1921  else if ( type == TRIANGLETYPE )
1922  return lwtriangle_perimeter((LWTRIANGLE*)geom);
1923  else if ( lwgeom_is_collection(geom) )
1924  {
1925  double perimeter = 0.0;
1926  uint32_t i;
1927  LWCOLLECTION *col = (LWCOLLECTION*)geom;
1928  for ( i = 0; i < col->ngeoms; i++ )
1929  perimeter += lwgeom_perimeter(col->geoms[i]);
1930  return perimeter;
1931  }
1932  else
1933  return 0.0;
1934 }
1935 
1936 double lwgeom_perimeter_2d(const LWGEOM *geom)
1937 {
1938  int type = geom->type;
1939  if ( type == POLYGONTYPE )
1940  return lwpoly_perimeter_2d((LWPOLY*)geom);
1941  else if ( type == CURVEPOLYTYPE )
1942  return lwcurvepoly_perimeter_2d((LWCURVEPOLY*)geom);
1943  else if ( type == TRIANGLETYPE )
1944  return lwtriangle_perimeter_2d((LWTRIANGLE*)geom);
1945  else if ( lwgeom_is_collection(geom) )
1946  {
1947  double perimeter = 0.0;
1948  uint32_t i;
1949  LWCOLLECTION *col = (LWCOLLECTION*)geom;
1950  for ( i = 0; i < col->ngeoms; i++ )
1951  perimeter += lwgeom_perimeter_2d(col->geoms[i]);
1952  return perimeter;
1953  }
1954  else
1955  return 0.0;
1956 }
1957 
1958 double lwgeom_length(const LWGEOM *geom)
1959 {
1960  int type = geom->type;
1961  if ( type == LINETYPE )
1962  return lwline_length((LWLINE*)geom);
1963  else if ( type == CIRCSTRINGTYPE )
1964  return lwcircstring_length((LWCIRCSTRING*)geom);
1965  else if ( type == COMPOUNDTYPE )
1966  return lwcompound_length((LWCOMPOUND*)geom);
1967  else if ( lwgeom_is_collection(geom) )
1968  {
1969  double length = 0.0;
1970  uint32_t i;
1971  LWCOLLECTION *col = (LWCOLLECTION*)geom;
1972  for ( i = 0; i < col->ngeoms; i++ )
1973  length += lwgeom_length(col->geoms[i]);
1974  return length;
1975  }
1976  else
1977  return 0.0;
1978 }
1979 
1980 double lwgeom_length_2d(const LWGEOM *geom)
1981 {
1982  int type = geom->type;
1983  if ( type == LINETYPE )
1984  return lwline_length_2d((LWLINE*)geom);
1985  else if ( type == CIRCSTRINGTYPE )
1986  return lwcircstring_length_2d((LWCIRCSTRING*)geom);
1987  else if ( type == COMPOUNDTYPE )
1988  return lwcompound_length_2d((LWCOMPOUND*)geom);
1989  else if ( lwgeom_is_collection(geom) )
1990  {
1991  double length = 0.0;
1992  uint32_t i;
1993  LWCOLLECTION *col = (LWCOLLECTION*)geom;
1994  for ( i = 0; i < col->ngeoms; i++ )
1995  length += lwgeom_length_2d(col->geoms[i]);
1996  return length;
1997  }
1998  else
1999  return 0.0;
2000 }
2001 
2002 void
2003 lwgeom_affine(LWGEOM *geom, const AFFINE *affine)
2004 {
2005  int type = geom->type;
2006  uint32_t i;
2007 
2008  switch(type)
2009  {
2010  /* Take advantage of fact tht pt/ln/circ/tri have same memory structure */
2011  case POINTTYPE:
2012  case LINETYPE:
2013  case CIRCSTRINGTYPE:
2014  case TRIANGLETYPE:
2015  {
2016  LWLINE *l = (LWLINE*)geom;
2017  ptarray_affine(l->points, affine);
2018  break;
2019  }
2020  case POLYGONTYPE:
2021  {
2022  LWPOLY *p = (LWPOLY*)geom;
2023  for( i = 0; i < p->nrings; i++ )
2024  ptarray_affine(p->rings[i], affine);
2025  break;
2026  }
2027  case CURVEPOLYTYPE:
2028  {
2029  LWCURVEPOLY *c = (LWCURVEPOLY*)geom;
2030  for( i = 0; i < c->nrings; i++ )
2031  lwgeom_affine(c->rings[i], affine);
2032  break;
2033  }
2034  default:
2035  {
2036  if( lwgeom_is_collection(geom) )
2037  {
2038  LWCOLLECTION *c = (LWCOLLECTION*)geom;
2039  for( i = 0; i < c->ngeoms; i++ )
2040  {
2041  lwgeom_affine(c->geoms[i], affine);
2042  }
2043  }
2044  else
2045  {
2046  lwerror("lwgeom_affine: unable to handle type '%s'", lwtype_name(type));
2047  }
2048  }
2049  }
2050 
2051  /* Recompute bbox if needed */
2052  if (geom->bbox)
2053  lwgeom_refresh_bbox(geom);
2054 }
2055 
2056 void
2057 lwgeom_scale(LWGEOM *geom, const POINT4D *factor)
2058 {
2059  int type = geom->type;
2060  uint32_t i;
2061 
2062  switch(type)
2063  {
2064  /* Take advantage of fact tht pt/ln/circ/tri have same memory structure */
2065  case POINTTYPE:
2066  case LINETYPE:
2067  case CIRCSTRINGTYPE:
2068  case TRIANGLETYPE:
2069  {
2070  LWLINE *l = (LWLINE*)geom;
2071  ptarray_scale(l->points, factor);
2072  break;
2073  }
2074  case POLYGONTYPE:
2075  {
2076  LWPOLY *p = (LWPOLY*)geom;
2077  for( i = 0; i < p->nrings; i++ )
2078  ptarray_scale(p->rings[i], factor);
2079  break;
2080  }
2081  case CURVEPOLYTYPE:
2082  {
2083  LWCURVEPOLY *c = (LWCURVEPOLY*)geom;
2084  for( i = 0; i < c->nrings; i++ )
2085  lwgeom_scale(c->rings[i], factor);
2086  break;
2087  }
2088  default:
2089  {
2090  if( lwgeom_is_collection(geom) )
2091  {
2092  LWCOLLECTION *c = (LWCOLLECTION*)geom;
2093  for( i = 0; i < c->ngeoms; i++ )
2094  {
2095  lwgeom_scale(c->geoms[i], factor);
2096  }
2097  }
2098  else
2099  {
2100  lwerror("lwgeom_scale: unable to handle type '%s'", lwtype_name(type));
2101  }
2102  }
2103  }
2104 
2105  /* Recompute bbox if needed */
2106  if (geom->bbox)
2107  lwgeom_refresh_bbox(geom);
2108 }
2109 
2110 LWGEOM *
2111 lwgeom_construct_empty(uint8_t type, int32_t srid, char hasz, char hasm)
2112 {
2113  switch(type)
2114  {
2115  case POINTTYPE:
2116  return lwpoint_as_lwgeom(lwpoint_construct_empty(srid, hasz, hasm));
2117  case LINETYPE:
2118  return lwline_as_lwgeom(lwline_construct_empty(srid, hasz, hasm));
2119  case POLYGONTYPE:
2120  return lwpoly_as_lwgeom(lwpoly_construct_empty(srid, hasz, hasm));
2121  case CURVEPOLYTYPE:
2122  return lwcurvepoly_as_lwgeom(lwcurvepoly_construct_empty(srid, hasz, hasm));
2123  case CIRCSTRINGTYPE:
2124  return lwcircstring_as_lwgeom(lwcircstring_construct_empty(srid, hasz, hasm));
2125  case TRIANGLETYPE:
2126  return lwtriangle_as_lwgeom(lwtriangle_construct_empty(srid, hasz, hasm));
2127  case COMPOUNDTYPE:
2128  case MULTIPOINTTYPE:
2129  case MULTILINETYPE:
2130  case MULTIPOLYGONTYPE:
2131  case COLLECTIONTYPE:
2132  return lwcollection_as_lwgeom(lwcollection_construct_empty(type, srid, hasz, hasm));
2133  default:
2134  lwerror("lwgeom_construct_empty: unsupported geometry type: %s",
2135  lwtype_name(type));
2136  return NULL;
2137  }
2138 }
2139 
2140 int
2141 lwgeom_startpoint(const LWGEOM *lwgeom, POINT4D *pt)
2142 {
2143  if ( ! lwgeom || lwgeom_is_empty(lwgeom) )
2144  return LW_FAILURE;
2145 
2146  switch( lwgeom->type )
2147  {
2148  case POINTTYPE:
2149  return ptarray_startpoint(((LWPOINT*)lwgeom)->point, pt);
2150  case TRIANGLETYPE:
2151  case CIRCSTRINGTYPE:
2152  case LINETYPE:
2153  return ptarray_startpoint(((LWLINE*)lwgeom)->points, pt);
2154  case POLYGONTYPE:
2155  return lwpoly_startpoint((LWPOLY*)lwgeom, pt);
2156  case TINTYPE:
2157  case CURVEPOLYTYPE:
2158  case COMPOUNDTYPE:
2159  case MULTIPOINTTYPE:
2160  case MULTILINETYPE:
2161  case MULTIPOLYGONTYPE:
2162  case COLLECTIONTYPE:
2163  case POLYHEDRALSURFACETYPE:
2164  return lwcollection_startpoint((LWCOLLECTION*)lwgeom, pt);
2165  default:
2166  lwerror("lwgeom_startpoint: unsupported geometry type: %s", lwtype_name(lwgeom->type));
2167  return LW_FAILURE;
2168  }
2169 }
2170 
2171 void
2173 {
2174  if (!geom) return;
2175  if (lwgeom_is_empty(geom)) return;
2176  switch ( geom->type )
2177  {
2178  case POINTTYPE:
2179  {
2180  LWPOINT *pt = (LWPOINT*)(geom);
2181  ptarray_grid_in_place(pt->point, grid);
2182  return;
2183  }
2184  case CIRCSTRINGTYPE:
2185  case TRIANGLETYPE:
2186  case LINETYPE:
2187  {
2188  LWLINE *ln = (LWLINE*)(geom);
2189  ptarray_grid_in_place(ln->points, grid);
2190  /* For invalid line, return an EMPTY */
2191  if (ln->points->npoints < 2)
2192  ln->points->npoints = 0;
2193  return;
2194  }
2195  case POLYGONTYPE:
2196  {
2197  LWPOLY *ply = (LWPOLY*)(geom);
2198  if (!ply->rings) return;
2199 
2200  /* Check first the external ring */
2201  uint32_t i = 0;
2202  POINTARRAY *pa = ply->rings[0];
2203  ptarray_grid_in_place(pa, grid);
2204  if (pa->npoints < 4)
2205  {
2206  /* External ring collapsed: free everything */
2207  for (i = 0; i < ply->nrings; i++)
2208  {
2209  ptarray_free(ply->rings[i]);
2210  }
2211  ply->nrings = 0;
2212  return;
2213  }
2214 
2215  /* Check the other rings */
2216  uint32_t j = 1;
2217  for (i = 1; i < ply->nrings; i++)
2218  {
2219  POINTARRAY *pa = ply->rings[i];
2220  ptarray_grid_in_place(pa, grid);
2221 
2222  /* Skip bad rings */
2223  if (pa->npoints >= 4)
2224  {
2225  ply->rings[j++] = pa;
2226  }
2227  else
2228  {
2229  ptarray_free(pa);
2230  }
2231  }
2232  /* Adjust ring count appropriately */
2233  ply->nrings = j;
2234  return;
2235  }
2236  case MULTIPOINTTYPE:
2237  case MULTILINETYPE:
2238  case MULTIPOLYGONTYPE:
2239  case TINTYPE:
2240  case COLLECTIONTYPE:
2241  case COMPOUNDTYPE:
2242  {
2243  LWCOLLECTION *col = (LWCOLLECTION*)(geom);
2244  uint32_t i, j = 0;
2245  if (!col->geoms) return;
2246  for (i = 0; i < col->ngeoms; i++)
2247  {
2248  LWGEOM *g = col->geoms[i];
2249  lwgeom_grid_in_place(g, grid);
2250  /* Empty geoms need to be freed */
2251  /* before we move on */
2252  if (lwgeom_is_empty(g))
2253  {
2254  lwgeom_free(g);
2255  continue;
2256  }
2257  col->geoms[j++] = g;
2258  }
2259  col->ngeoms = j;
2260  return;
2261  }
2262  default:
2263  {
2264  lwerror("%s: Unsupported geometry type: %s", __func__,
2265  lwtype_name(geom->type));
2266  return;
2267  }
2268  }
2269 }
2270 
2271 
2272 LWGEOM *
2273 lwgeom_grid(const LWGEOM *lwgeom, const gridspec *grid)
2274 {
2275  LWGEOM *lwgeom_out = lwgeom_clone_deep(lwgeom);
2276  lwgeom_grid_in_place(lwgeom_out, grid);
2277  return lwgeom_out;
2278 }
2279 
2280 
2281 /* Prototype for recursion */
2282 static void lwgeom_subdivide_recursive(const LWGEOM *geom,
2283  uint8_t dimension,
2284  uint32_t maxvertices,
2285  uint32_t depth,
2286  LWCOLLECTION *col,
2287  double gridSize);
2288 
2289 static void
2291  uint8_t dimension,
2292  uint32_t maxvertices,
2293  uint32_t depth,
2294  LWCOLLECTION *col,
2295  double gridSize)
2296 {
2297  const uint32_t maxdepth = 50;
2298 
2299  if (!geom)
2300  return;
2301 
2302  const GBOX *box_in = lwgeom_get_bbox(geom);
2303  if (!box_in)
2304  return;
2305 
2306  LW_ON_INTERRUPT(return);
2307 
2308  GBOX clip;
2309  gbox_duplicate(box_in, &clip);
2310  double width = clip.xmax - clip.xmin;
2311  double height = clip.ymax - clip.ymin;
2312 
2313  if ( geom->type == POLYHEDRALSURFACETYPE || geom->type == TINTYPE )
2314  lwerror("%s: unsupported geometry type '%s'", __func__, lwtype_name(geom->type));
2315 
2316  if ( width == 0.0 && height == 0.0 )
2317  {
2318  if ( geom->type == POINTTYPE && dimension == 0)
2320  return;
2321  }
2322 
2323  if (width == 0.0)
2324  {
2325  clip.xmax += FP_TOLERANCE;
2326  clip.xmin -= FP_TOLERANCE;
2327  width = 2 * FP_TOLERANCE;
2328  }
2329  if (height == 0.0)
2330  {
2331  clip.ymax += FP_TOLERANCE;
2332  clip.ymin -= FP_TOLERANCE;
2333  height = 2 * FP_TOLERANCE;
2334  }
2335 
2336  /* Always just recurse into collections */
2337  if ( lwgeom_is_collection(geom) && geom->type != MULTIPOINTTYPE )
2338  {
2339  LWCOLLECTION *incol = (LWCOLLECTION*)geom;
2340  /* Don't increment depth yet, since we aren't actually
2341  * subdividing geometries yet */
2342  for (uint32_t i = 0; i < incol->ngeoms; i++ )
2343  lwgeom_subdivide_recursive(incol->geoms[i], dimension, maxvertices, depth, col, gridSize);
2344  return;
2345  }
2346 
2347  if (lwgeom_dimension(geom) < dimension)
2348  {
2349  /* We've hit a lower dimension object produced by clipping at
2350  * a shallower recursion level. Ignore it. */
2351  return;
2352  }
2353 
2354  /* But don't go too far. 2^50 ~= 10^15, that's enough subdivision */
2355  /* Just add what's left */
2356  if ( depth > maxdepth )
2357  {
2359  return;
2360  }
2361 
2362  uint32_t nvertices = lwgeom_count_vertices(geom);
2363 
2364  /* Skip empties entirely */
2365  if (nvertices == 0)
2366  return;
2367 
2368  /* If it is under the vertex tolerance, just add it, we're done */
2369  if (nvertices <= maxvertices)
2370  {
2372  return;
2373  }
2374 
2375  uint8_t split_ordinate = (width > height) ? 0 : 1;
2376  double center = (split_ordinate == 0) ? (clip.xmin + clip.xmax) / 2 : (clip.ymin + clip.ymax) / 2;
2377  double pivot = DBL_MAX;
2378  if (geom->type == POLYGONTYPE)
2379  {
2380  uint32_t ring_to_trim = 0;
2381  double ring_area = 0;
2382  double pivot_eps = DBL_MAX;
2383  double pt_eps = DBL_MAX;
2384  POINTARRAY *pa;
2385  LWPOLY *lwpoly = (LWPOLY *)geom;
2386 
2387  /* if there are more points in holes than in outer ring */
2388  if (nvertices >= 2 * lwpoly->rings[0]->npoints)
2389  {
2390  /* trim holes starting from biggest */
2391  for (uint32_t i = 1; i < lwpoly->nrings; i++)
2392  {
2393  double current_ring_area = fabs(ptarray_signed_area(lwpoly->rings[i]));
2394  if (current_ring_area >= ring_area)
2395  {
2396  ring_area = current_ring_area;
2397  ring_to_trim = i;
2398  }
2399  }
2400  }
2401 
2402  pa = lwpoly->rings[ring_to_trim];
2403 
2404  /* find most central point in chosen ring */
2405  for (uint32_t i = 0; i < pa->npoints; i++)
2406  {
2407  double pt;
2408  if (split_ordinate == 0)
2409  pt = getPoint2d_cp(pa, i)->x;
2410  else
2411  pt = getPoint2d_cp(pa, i)->y;
2412  pt_eps = fabs(pt - center);
2413  if (pivot_eps > pt_eps)
2414  {
2415  pivot = pt;
2416  pivot_eps = pt_eps;
2417  }
2418  }
2419  }
2420  GBOX subbox1, subbox2;
2421  gbox_duplicate(&clip, &subbox1);
2422  gbox_duplicate(&clip, &subbox2);
2423 
2424  if (pivot == DBL_MAX)
2425  pivot = center;
2426 
2427  if (split_ordinate == 0)
2428  {
2429  if (FP_NEQUALS(subbox1.xmax, pivot) && FP_NEQUALS(subbox1.xmin, pivot))
2430  subbox1.xmax = subbox2.xmin = pivot;
2431  else
2432  subbox1.xmax = subbox2.xmin = center;
2433  }
2434  else
2435  {
2436  if (FP_NEQUALS(subbox1.ymax, pivot) && FP_NEQUALS(subbox1.ymin, pivot))
2437  subbox1.ymax = subbox2.ymin = pivot;
2438  else
2439  subbox1.ymax = subbox2.ymin = center;
2440  }
2441 
2442  ++depth;
2443 
2444  {
2446  geom->srid, subbox1.xmin, subbox1.ymin, subbox1.xmax, subbox1.ymax);
2447  LWGEOM *clipped = lwgeom_intersection_prec(geom, subbox, gridSize);
2448  lwgeom_simplify_in_place(clipped, 0.0, LW_TRUE);
2449  lwgeom_free(subbox);
2450  if (clipped && !lwgeom_is_empty(clipped))
2451  {
2452  lwgeom_subdivide_recursive(clipped, dimension, maxvertices, depth, col, gridSize);
2453  lwgeom_free(clipped);
2454  }
2455  }
2456  {
2458  geom->srid, subbox2.xmin, subbox2.ymin, subbox2.xmax, subbox2.ymax);
2459  LWGEOM *clipped = lwgeom_intersection_prec(geom, subbox, gridSize);
2460  lwgeom_simplify_in_place(clipped, 0.0, LW_TRUE);
2461  lwgeom_free(subbox);
2462  if (clipped && !lwgeom_is_empty(clipped))
2463  {
2464  lwgeom_subdivide_recursive(clipped, dimension, maxvertices, depth, col, gridSize);
2465  lwgeom_free(clipped);
2466  }
2467  }
2468 }
2469 
2470 LWCOLLECTION *
2471 lwgeom_subdivide_prec(const LWGEOM *geom, uint32_t maxvertices, double gridSize)
2472 {
2473  static uint32_t startdepth = 0;
2474  static uint32_t minmaxvertices = 5;
2475  LWCOLLECTION *col;
2476 
2478 
2479  if ( lwgeom_is_empty(geom) )
2480  return col;
2481 
2482  if ( maxvertices < minmaxvertices )
2483  {
2484  lwcollection_free(col);
2485  lwerror("%s: cannot subdivide to fewer than %d vertices per output", __func__, minmaxvertices);
2486  }
2487 
2488  lwgeom_subdivide_recursive(geom, lwgeom_dimension(geom), maxvertices, startdepth, col, gridSize);
2489  lwgeom_set_srid((LWGEOM*)col, geom->srid);
2490  return col;
2491 }
2492 
2493 LWCOLLECTION *
2494 lwgeom_subdivide(const LWGEOM *geom, uint32_t maxvertices)
2495 {
2496  return lwgeom_subdivide_prec(geom, maxvertices, -1);
2497 }
2498 
2499 
2500 
2501 int
2503 {
2504  int type = geom->type;
2505 
2506  if( type != LINETYPE )
2507  {
2508  lwnotice("Geometry is not a LINESTRING");
2509  return LW_FALSE;
2510  }
2511  return lwline_is_trajectory((LWLINE*)geom);
2512 }
2513 
2514 static uint8_t
2515 bits_for_precision(int32_t significant_digits)
2516 {
2517  int32_t bits_needed = ceil(significant_digits / log10(2));
2518 
2519  if (bits_needed > 52)
2520  {
2521  return 52;
2522  }
2523  else if (bits_needed < 1)
2524  {
2525  return 1;
2526  }
2527 
2528  return bits_needed;
2529 }
2530 
2531 static double trim_preserve_decimal_digits(double d, int32_t decimal_digits)
2532 {
2533  if (d == 0)
2534  return 0;
2535 
2536  int digits_left_of_decimal = (int) (1 + log10(fabs(d)));
2537  uint8_t bits_needed = bits_for_precision(decimal_digits + digits_left_of_decimal);
2538  uint64_t mask = 0xffffffffffffffffULL << (52 - bits_needed);
2539  uint64_t dint = 0;
2540  size_t dsz = sizeof(d) < sizeof(dint) ? sizeof(d) : sizeof(dint);
2541 
2542  memcpy(&dint, &d, dsz);
2543  dint &= mask;
2544  memcpy(&d, &dint, dsz);
2545  return d;
2546 }
2547 
2548 void lwgeom_trim_bits_in_place(LWGEOM* geom, int32_t prec_x, int32_t prec_y, int32_t prec_z, int32_t prec_m)
2549 {
2551  POINT4D p;
2552 
2553  while (lwpointiterator_has_next(it))
2554  {
2555  lwpointiterator_peek(it, &p);
2556  p.x = trim_preserve_decimal_digits(p.x, prec_x);
2557  p.y = trim_preserve_decimal_digits(p.y, prec_y);
2558  if (lwgeom_has_z(geom))
2559  p.z = trim_preserve_decimal_digits(p.z, prec_z);
2560  if (lwgeom_has_m(geom))
2561  p.m = trim_preserve_decimal_digits(p.m, prec_m);
2563  }
2564 
2566 }
2567 
2568 LWGEOM *
2570 {
2571  int32_t srid = lwgeom_get_srid(lwgeom);
2572  uint8_t hasz = lwgeom_has_z(lwgeom);
2573  uint8_t hasm = lwgeom_has_m(lwgeom);
2574 
2575  switch (lwgeom->type)
2576  {
2577  case POINTTYPE:
2578  case MULTIPOINTTYPE: {
2579  return lwgeom_construct_empty(lwgeom->type, srid, hasz, hasm);
2580  }
2581  case LINETYPE:
2582  case CIRCSTRINGTYPE: {
2583  if (lwgeom_is_closed(lwgeom) || lwgeom_is_empty(lwgeom))
2584  return (LWGEOM *)lwmpoint_construct_empty(srid, hasz, hasm);
2585  else
2586  {
2587  LWLINE *lwline = (LWLINE *)lwgeom;
2588  LWMPOINT *lwmpoint = lwmpoint_construct_empty(srid, hasz, hasm);
2589  POINT4D pt;
2590  getPoint4d_p(lwline->points, 0, &pt);
2591  lwmpoint_add_lwpoint(lwmpoint, lwpoint_make(srid, hasz, hasm, &pt));
2592  getPoint4d_p(lwline->points, lwline->points->npoints - 1, &pt);
2593  lwmpoint_add_lwpoint(lwmpoint, lwpoint_make(srid, hasz, hasm, &pt));
2594 
2595  return (LWGEOM *)lwmpoint;
2596  }
2597  }
2598  case MULTILINETYPE:
2599  case MULTICURVETYPE: {
2600  LWMLINE *lwmline = (LWMLINE *)lwgeom;
2601  POINT4D *out = lwalloc(sizeof(POINT4D) * lwmline->ngeoms * 2);
2602  uint32_t n = 0;
2603 
2604  for (uint32_t i = 0; i < lwmline->ngeoms; i++)
2605  {
2606  LWMPOINT *points = lwgeom_as_lwmpoint(lwgeom_boundary((LWGEOM *)lwmline->geoms[i]));
2607  if (!points)
2608  continue;
2609 
2610  for (uint32_t k = 0; k < points->ngeoms; k++)
2611  {
2612  POINT4D pt = getPoint4d(points->geoms[k]->point, 0);
2613 
2614  uint8_t seen = LW_FALSE;
2615  for (uint32_t j = 0; j < n; j++)
2616  {
2617  if (memcmp(&(out[j]), &pt, sizeof(POINT4D)) == 0)
2618  {
2619  seen = LW_TRUE;
2620  out[j] = out[--n];
2621  break;
2622  }
2623  }
2624  if (!seen)
2625  out[n++] = pt;
2626  }
2627 
2628  lwgeom_free((LWGEOM *)points);
2629  }
2630 
2631  LWMPOINT *lwmpoint = lwmpoint_construct_empty(srid, hasz, hasm);
2632 
2633  for (uint32_t i = 0; i < n; i++)
2634  lwmpoint_add_lwpoint(lwmpoint, lwpoint_make(srid, hasz, hasm, &(out[i])));
2635 
2636  lwfree(out);
2637 
2638  return (LWGEOM *)lwmpoint;
2639  }
2640  case TRIANGLETYPE: {
2641  LWTRIANGLE *lwtriangle = (LWTRIANGLE *)lwgeom;
2642  POINTARRAY *points = ptarray_clone_deep(lwtriangle->points);
2643  return (LWGEOM *)lwline_construct(srid, 0, points);
2644  }
2645  case POLYGONTYPE: {
2646  LWPOLY *lwpoly = (LWPOLY *)lwgeom;
2647 
2648  LWMLINE *lwmline = lwmline_construct_empty(srid, hasz, hasm);
2649  for (uint32_t i = 0; i < lwpoly->nrings; i++)
2650  {
2651  POINTARRAY *ring = ptarray_clone_deep(lwpoly->rings[i]);
2652  lwmline_add_lwline(lwmline, lwline_construct(srid, 0, ring));
2653  }
2654 
2655  /* Homogenize the multilinestring to hopefully get a single LINESTRING */
2656  LWGEOM *lwout = lwgeom_homogenize((LWGEOM *)lwmline);
2657  lwgeom_free((LWGEOM *)lwmline);
2658  return lwout;
2659  }
2660  case CURVEPOLYTYPE: {
2661  LWCURVEPOLY *lwcurvepoly = (LWCURVEPOLY *)lwgeom;
2662  LWCOLLECTION *lwcol = lwcollection_construct_empty(MULTICURVETYPE, srid, hasz, hasm);
2663 
2664  for (uint32_t i = 0; i < lwcurvepoly->nrings; i++)
2665  lwcol = lwcollection_add_lwgeom(lwcol, lwgeom_clone_deep(lwcurvepoly->rings[i]));
2666 
2667  return (LWGEOM *)lwcol;
2668  }
2669  case MULTIPOLYGONTYPE:
2670  case COLLECTIONTYPE:
2671  case TINTYPE: {
2672  LWCOLLECTION *lwcol = (LWCOLLECTION *)lwgeom;
2673  LWCOLLECTION *lwcol_boundary = lwcollection_construct_empty(COLLECTIONTYPE, srid, hasz, hasm);
2674 
2675  for (uint32_t i = 0; i < lwcol->ngeoms; i++)
2676  lwcollection_add_lwgeom(lwcol_boundary, lwgeom_boundary(lwcol->geoms[i]));
2677 
2678  LWGEOM *lwout = lwgeom_homogenize((LWGEOM *)lwcol_boundary);
2679  lwgeom_free((LWGEOM *)lwcol_boundary);
2680 
2681  return lwout;
2682  }
2683  default:
2684  lwerror("%s: unsupported geometry type: %s", __func__, lwtype_name(lwgeom->type));
2685  return NULL;
2686  }
2687 }
2688 
2689 int
2690 lwgeom_isfinite(const LWGEOM *lwgeom)
2691 {
2693  int hasz = lwgeom_has_z(lwgeom);
2694  int hasm = lwgeom_has_m(lwgeom);
2695 
2696  while (lwpointiterator_has_next(it))
2697  {
2698  POINT4D p;
2699  lwpointiterator_next(it, &p);
2700  int finite = isfinite(p.x) &&
2701  isfinite(p.y) &&
2702  (hasz ? isfinite(p.z) : 1) &&
2703  (hasm ? isfinite(p.m) : 1);
2704 
2705  if (!finite)
2706  {
2708  return LW_FALSE;
2709  }
2710  }
2712  return LW_TRUE;
2713 }
char * r
Definition: cu_in_wkt.c:24
char result[OUT_DOUBLE_BUFFER_SIZE]
Definition: cu_print.c:267
int gbox_same(const GBOX *g1, const GBOX *g2)
Check if 2 given Gbox are the same.
Definition: gbox.c:164
void gbox_duplicate(const GBOX *original, GBOX *duplicate)
Copy the values of original GBOX into duplicate.
Definition: gbox.c:433
int lwgeom_calculate_gbox_cartesian(const LWGEOM *lwgeom, GBOX *gbox)
Calculate the 2-4D bounding box of a geometry.
Definition: gbox.c:740
GBOX * gbox_new(lwflags_t flags)
Create a new gbox with the dimensionality indicated by the flags.
Definition: gbox.c:32
GBOX * gbox_clone(const GBOX *gbox)
Definition: gbox.c:45
LWPOINT * lwpoint_construct_empty(int32_t srid, char hasz, char hasm)
Definition: lwpoint.c:151
POINT4D getPoint4d(const POINTARRAY *pa, uint32_t n)
Definition: lwgeom_api.c:108
#define LW_FALSE
Definition: liblwgeom.h:108
#define COLLECTIONTYPE
Definition: liblwgeom.h:122
#define COMPOUNDTYPE
Definition: liblwgeom.h:124
double distance2d_pt_pt(const POINT2D *p1, const POINT2D *p2)
Definition: measures.c:2397
void lwmpoint_free(LWMPOINT *mpt)
Definition: lwmpoint.c:72
LWPOINTITERATOR * lwpointiterator_create(const LWGEOM *g)
Create a new LWPOINTITERATOR over supplied LWGEOM*.
Definition: lwiterator.c:242
LWLINE * lwline_segmentize2d(const LWLINE *line, double dist)
Definition: lwline.c:132
void lwpoint_free(LWPOINT *pt)
Definition: lwpoint.c:213
#define LW_FAILURE
Definition: liblwgeom.h:110
LWCOLLECTION * lwcollection_segmentize2d(const LWCOLLECTION *coll, double dist)
Definition: lwcollection.c:251
#define CURVEPOLYTYPE
Definition: liblwgeom.h:125
void lwmpoly_free(LWMPOLY *mpoly)
Definition: lwmpoly.c:53
#define MULTILINETYPE
Definition: liblwgeom.h:120
LWCURVEPOLY * lwcurvepoly_construct_from_lwpoly(LWPOLY *lwpoly)
Construct an equivalent curve polygon from a polygon.
Definition: lwcurvepoly.c:52
int lwpointiterator_next(LWPOINTITERATOR *s, POINT4D *p)
Attempts to assign the next point in the iterator to p, and advances the iterator to the next point.
Definition: lwiterator.c:210
#define MULTISURFACETYPE
Definition: liblwgeom.h:127
#define LINETYPE
Definition: liblwgeom.h:117
void lwtin_free(LWTIN *tin)
Definition: lwtin.c:39
int lwgeom_calculate_gbox_geodetic(const LWGEOM *geom, GBOX *gbox)
Calculate the geodetic bounding box for an LWGEOM.
Definition: lwgeodetic.c:3002
#define WKT_EXTENDED
Definition: liblwgeom.h:2167
LWGEOM * lwgeom_stroke(const LWGEOM *geom, uint32_t perQuad)
Convert type with arcs into equivalent linearized type.
Definition: lwstroke.c:871
#define MULTIPOINTTYPE
Definition: liblwgeom.h:119
#define FLAGS_SET_BBOX(flags, value)
Definition: liblwgeom.h:188
int lwpointiterator_peek(LWPOINTITERATOR *s, POINT4D *p)
Attempts to assigns the next point in the iterator to p.
Definition: lwiterator.c:193
POINTARRAY * ptarray_clone_deep(const POINTARRAY *ptarray)
Deep clone a pointarray (also clones serialized pointlist)
Definition: ptarray.c:634
LWMLINE * lwmline_add_lwline(LWMLINE *mobj, const LWLINE *obj)
Definition: lwmline.c:46
void lwpointiterator_destroy(LWPOINTITERATOR *s)
Free all memory associated with the iterator.
Definition: lwiterator.c:267
LWMLINE * lwmline_construct_empty(int32_t srid, char hasz, char hasm)
Definition: lwmline.c:38
#define POINTTYPE
LWTYPE numbers, used internally by PostGIS.
Definition: liblwgeom.h:116
#define FLAGS_GET_Z(flags)
Definition: liblwgeom.h:179
LWLINE * lwline_construct(int32_t srid, GBOX *bbox, POINTARRAY *points)
Definition: lwline.c:42
LWTRIANGLE * lwtriangle_construct_empty(int32_t srid, char hasz, char hasm)
Definition: lwtriangle.c:58
LWMPOINT * lwmpoint_add_lwpoint(LWMPOINT *mobj, const LWPOINT *obj)
Definition: lwmpoint.c:45
#define TINTYPE
Definition: liblwgeom.h:130
LWPOLY * lwpoly_construct_envelope(int32_t srid, double x1, double y1, double x2, double y2)
Definition: lwpoly.c:98
#define MULTIPOLYGONTYPE
Definition: liblwgeom.h:121
int lwpointiterator_modify_next(LWPOINTITERATOR *s, const POINT4D *p)
Attempts to replace the next point int the iterator with p, and advances the iterator to the next poi...
Definition: lwiterator.c:224
LWGEOM * lwgeom_homogenize(const LWGEOM *geom)
Definition: lwhomogenize.c:208
void lwfree(void *mem)
Definition: lwutil.c:242
#define FLAGS_NDIMS(flags)
Definition: liblwgeom.h:193
#define POLYGONTYPE
Definition: liblwgeom.h:118
void lwcircstring_free(LWCIRCSTRING *curve)
Definition: lwcircstring.c:97
LWMPOINT * lwmpoint_construct_empty(int32_t srid, char hasz, char hasm)
Definition: lwmpoint.c:39
void lwtriangle_free(LWTRIANGLE *triangle)
Definition: lwtriangle.c:69
#define POLYHEDRALSURFACETYPE
Definition: liblwgeom.h:128
#define CIRCSTRINGTYPE
Definition: liblwgeom.h:123
LWPOINTITERATOR * lwpointiterator_create_rw(LWGEOM *g)
Create a new LWPOINTITERATOR over supplied LWGEOM* Supports modification of coordinates during iterat...
Definition: lwiterator.c:251
LWCOLLECTION * lwcollection_construct_empty(uint8_t type, int32_t srid, char hasz, char hasm)
Definition: lwcollection.c:92
const char * lwtype_name(uint8_t type)
Return the type name string associated with a type number (e.g.
Definition: lwutil.c:216
#define FLAGS_GET_M(flags)
Definition: liblwgeom.h:180
void lwcollection_free(LWCOLLECTION *col)
Definition: lwcollection.c:357
int getPoint4d_p(const POINTARRAY *pa, uint32_t n, POINT4D *point)
Definition: lwgeom_api.c:126
LWGEOM * lwgeom_intersection_prec(const LWGEOM *geom1, const LWGEOM *geom2, double gridSize)
void ptarray_free(POINTARRAY *pa)
Definition: ptarray.c:327
char * lwgeom_to_wkt(const LWGEOM *geom, uint8_t variant, int precision, size_t *size_out)
WKT emitter function.
Definition: lwout_wkt.c:704
#define FLAGS_GET_SOLID(flags)
Definition: liblwgeom.h:184
#define FLAGS_GET_ZM(flags)
Definition: liblwgeom.h:194
LWPOLY * lwpoly_segmentize2d(const LWPOLY *line, double dist)
Definition: lwpoly.c:312
int lwpointiterator_has_next(LWPOINTITERATOR *s)
Returns LW_TRUE if there is another point available in the iterator.
Definition: lwiterator.c:202
#define MULTICURVETYPE
Definition: liblwgeom.h:126
#define TRIANGLETYPE
Definition: liblwgeom.h:129
LWCOLLECTION * lwcollection_add_lwgeom(LWCOLLECTION *col, const LWGEOM *geom)
Appends geom to the collection managed by col.
Definition: lwcollection.c:188
#define FLAGS_SET_GEODETIC(flags, value)
Definition: liblwgeom.h:189
void * lwalloc(size_t size)
Definition: lwutil.c:227
LWCOLLECTION * lwcollection_construct(uint8_t type, int32_t srid, GBOX *bbox, uint32_t ngeoms, LWGEOM **geoms)
Definition: lwcollection.c:42
void lwpsurface_free(LWPSURFACE *psurf)
Definition: lwpsurface.c:39
void lwpoly_free(LWPOLY *poly)
Definition: lwpoly.c:175
void lwmline_free(LWMLINE *mline)
Definition: lwmline.c:112
#define LW_TRUE
Return types for functions with status returns.
Definition: liblwgeom.h:107
int lwline_is_trajectory(const LWLINE *geom)
Definition: lwline.c:454
LWCURVEPOLY * lwcurvepoly_construct_empty(int32_t srid, char hasz, char hasm)
Definition: lwcurvepoly.c:35
#define SRID_UNKNOWN
Unknown SRID value.
Definition: liblwgeom.h:229
LWPOLY * lwpoly_from_lwlines(const LWLINE *shell, uint32_t nholes, const LWLINE **holes)
Definition: lwpoly.c:360
#define NUMTYPES
Definition: liblwgeom.h:132
LWCOMPOUND * lwcompound_construct_from_lwline(const LWLINE *lwpoly)
Construct an equivalent compound curve from a linestring.
Definition: lwcompound.c:204
LWPOLY * lwpoly_construct_empty(int32_t srid, char hasz, char hasm)
Definition: lwpoly.c:161
LWCIRCSTRING * lwcircstring_construct_empty(int32_t srid, char hasz, char hasm)
Definition: lwcircstring.c:79
void lwline_free(LWLINE *line)
Definition: lwline.c:67
LWLINE * lwline_construct_empty(int32_t srid, char hasz, char hasm)
Definition: lwline.c:55
#define FLAGS_GET_GEODETIC(flags)
Definition: liblwgeom.h:182
LWPOINT * lwpoint_make(int32_t srid, int hasz, int hasm, const POINT4D *p)
Definition: lwpoint.c:206
enum LWORD_T LWORD
Ordinate names.
double lwline_length_2d(const LWLINE *line)
Definition: lwline.c:520
void ptarray_longitude_shift(POINTARRAY *pa)
Longitude shift for a pointarray.
Definition: ptarray.c:1485
LWLINE * lwline_clone_deep(const LWLINE *lwgeom)
Definition: lwline.c:109
double lwcircstring_length_2d(const LWCIRCSTRING *circ)
Definition: lwcircstring.c:274
LWPOINT * lwpoint_clone(const LWPOINT *lwgeom)
Definition: lwpoint.c:239
int lwcircstring_is_closed(const LWCIRCSTRING *curve)
Definition: lwcircstring.c:261
void ptarray_reverse_in_place(POINTARRAY *pa)
Definition: ptarray.c:339
LWLINE * lwline_clone(const LWLINE *lwgeom)
Definition: lwline.c:93
char lwcollection_same(const LWCOLLECTION *p1, const LWCOLLECTION *p2)
check for same geometry composition
Definition: lwcollection.c:279
double lwtriangle_area(const LWTRIANGLE *triangle)
Find the area of the outer ring.
Definition: lwtriangle.c:178
double lwcompound_length_2d(const LWCOMPOUND *comp)
Definition: lwcompound.c:74
uint32_t lwpoly_count_vertices(LWPOLY *poly)
Definition: lwpoly.c:418
double ptarray_signed_area(const POINTARRAY *pa)
Returns the area in cartesian units.
Definition: ptarray.c:1003
#define LW_ON_INTERRUPT(x)
int lwline_is_closed(const LWLINE *line)
Definition: lwline.c:445
char lwtriangle_same(const LWTRIANGLE *p1, const LWTRIANGLE *p2)
Definition: lwtriangle.c:126
int lwpoly_startpoint(const LWPOLY *lwpoly, POINT4D *pt)
Definition: lwpoly.c:524
LWPOINT * lwpoint_force_dims(const LWPOINT *lwpoint, int hasz, int hasm, double zval, double mval)
Definition: lwpoint.c:271
int lwcollection_startpoint(const LWCOLLECTION *col, POINT4D *pt)
Definition: lwcollection.c:550
int ptarray_startpoint(const POINTARRAY *pa, POINT4D *pt)
Definition: ptarray.c:2063
int lwcompound_is_closed(const LWCOMPOUND *curve)
Definition: lwcompound.c:35
double lwtriangle_perimeter_2d(const LWTRIANGLE *triangle)
Definition: lwtriangle.c:210
int lwpoly_is_clockwise(LWPOLY *poly)
Definition: lwpoly.c:288
void ptarray_grid_in_place(POINTARRAY *pa, const gridspec *grid)
Snap to grid.
Definition: ptarray.c:2078
double lwpoly_area(const LWPOLY *poly)
Find the area of the outer ring - sum (area of inner rings).
Definition: lwpoly.c:434
int lwtriangle_is_clockwise(LWTRIANGLE *triangle)
Definition: lwtriangle.c:113
char lwline_same(const LWLINE *p1, const LWLINE *p2)
Definition: lwline.c:141
void ptarray_remove_repeated_points_in_place(POINTARRAY *pa, double tolerance, uint32_t min_points)
Definition: ptarray.c:1522
double lwtriangle_perimeter(const LWTRIANGLE *triangle)
Definition: lwtriangle.c:201
int lwpoint_is_empty(const LWPOINT *point)
void ptarray_simplify_in_place(POINTARRAY *pa, double tolerance, uint32_t minpts)
Definition: ptarray.c:1700
LWLINE * lwline_force_dims(const LWLINE *lwline, int hasz, int hasm, double zval, double mval)
Definition: lwline.c:486
double lwcompound_length(const LWCOMPOUND *comp)
Definition: lwcompound.c:69
LWCOLLECTION * lwcollection_clone_deep(const LWCOLLECTION *lwgeom)
Deep clone LWCOLLECTION object.
Definition: lwcollection.c:150
LWCOLLECTION * lwcollection_force_dims(const LWCOLLECTION *lwcol, int hasz, int hasm, double zval, double mval)
Definition: lwcollection.c:476
void lwtriangle_force_clockwise(LWTRIANGLE *triangle)
Definition: lwtriangle.c:106
#define FP_TOLERANCE
Floating point comparators.
void ptarray_scale(POINTARRAY *pa, const POINT4D *factor)
WARNING, make sure you send in only 16-member double arrays or obviously things will go pear-shaped f...
Definition: ptarray.c:2045
void ptarray_affine(POINTARRAY *pa, const AFFINE *affine)
Affine transform a pointarray.
Definition: ptarray.c:1877
int lwpsurface_is_closed(const LWPSURFACE *psurface)
Definition: lwpsurface.c:99
double lwcurvepoly_perimeter(const LWCURVEPOLY *poly)
Definition: lwcurvepoly.c:147
char lwpoly_same(const LWPOLY *p1, const LWPOLY *p2)
Definition: lwpoly.c:339
double lwline_length(const LWLINE *line)
Definition: lwline.c:513
void lwpoly_force_clockwise(LWPOLY *poly)
Definition: lwpoly.c:268
void ptarray_swap_ordinates(POINTARRAY *pa, LWORD o1, LWORD o2)
Swap ordinate values o1 and o2 on a given POINTARRAY.
Definition: ptarray.c:387
double lwpoly_perimeter_2d(const LWPOLY *poly)
Compute the sum of polygon rings length (forcing 2d computation).
Definition: lwpoly.c:485
int lwtin_is_closed(const LWTIN *tin)
Definition: lwtin.c:93
char lwcircstring_same(const LWCIRCSTRING *p1, const LWCIRCSTRING *p2)
Definition: lwcircstring.c:131
LWPOLY * lwpoly_clone_deep(const LWPOLY *lwgeom)
Definition: lwpoly.c:228
double lwcircstring_length(const LWCIRCSTRING *circ)
Definition: lwcircstring.c:269
uint32_t lwcollection_count_vertices(LWCOLLECTION *col)
Definition: lwcollection.c:500
double lwcurvepoly_area(const LWCURVEPOLY *curvepoly)
This should be rewritten to make use of the curve itself.
Definition: lwcurvepoly.c:133
LWCOLLECTION * lwcollection_clone(const LWCOLLECTION *lwgeom)
Clone LWCOLLECTION object.
Definition: lwcollection.c:124
int lwpoly_is_closed(const LWPOLY *poly)
Definition: lwpoly.c:499
LWPOLY * lwpoly_clone(const LWPOLY *lwgeom)
Definition: lwpoly.c:213
void ptarray_copy_point(POINTARRAY *pa, uint32_t from, uint32_t to)
Definition: lwgeom_api.c:395
uint32_t lwline_count_vertices(LWLINE *line)
Definition: lwline.c:505
double lwcurvepoly_perimeter_2d(const LWCURVEPOLY *poly)
Definition: lwcurvepoly.c:159
double lwpoly_perimeter(const LWPOLY *poly)
Compute the sum of polygon rings length.
Definition: lwpoly.c:467
LWTRIANGLE * lwtriangle_clone(const LWTRIANGLE *lwgeom)
Definition: lwtriangle.c:99
#define FP_NEQUALS(A, B)
char lwpoint_same(const LWPOINT *p1, const LWPOINT *p2)
Definition: lwpoint.c:264
LWCIRCSTRING * lwcircstring_clone(const LWCIRCSTRING *curve)
Definition: lwcircstring.c:124
LWPOLY * lwpoly_force_dims(const LWPOLY *lwpoly, int hasz, int hasm, double zval, double mval)
Definition: lwpoly.c:394
int p2d_same(const POINT2D *p1, const POINT2D *p2)
Definition: lwalgorithm.c:50
int lwgeom_is_closed(const LWGEOM *geom)
Return true or false depending on whether a geometry is a linear feature that closes on itself.
Definition: lwgeom.c:1036
void lwgeom_refresh_bbox(LWGEOM *lwgeom)
Drop current bbox and calculate a fresh one.
Definition: lwgeom.c:690
LWGEOM * lwgeom_force_dims(const LWGEOM *geom, int hasz, int hasm, double zval, double mval)
Definition: lwgeom.c:800
static double trim_preserve_decimal_digits(double d, int32_t decimal_digits)
Definition: lwgeom.c:2531
LWLINE * lwgeom_as_lwline(const LWGEOM *lwgeom)
Definition: lwgeom.c:162
LWGEOM * lwcompound_as_lwgeom(const LWCOMPOUND *obj)
Definition: lwgeom.c:307
void lwgeom_set_geodetic(LWGEOM *geom, int value)
Set the FLAGS geodetic bit on geometry an all sub-geometries and pointlists.
Definition: lwgeom.c:947
LWGEOM * lwline_as_lwgeom(const LWLINE *obj)
Definition: lwgeom.c:322
LWGEOM * lwgeom_simplify(const LWGEOM *igeom, double dist, int preserve_collapsed)
Simplification.
Definition: lwgeom.c:1876
char lwgeom_same(const LWGEOM *lwgeom1, const LWGEOM *lwgeom2)
geom1 same as geom2 iff
Definition: lwgeom.c:574
LWGEOM * lwcollection_as_lwgeom(const LWCOLLECTION *obj)
Definition: lwgeom.c:292
int lwgeom_ndims(const LWGEOM *geom)
Return the number of dimensions (2, 3, 4) in a geometry.
Definition: lwgeom.c:938
uint32_t lwtype_get_collectiontype(uint8_t type)
Given an lwtype number, what homogeneous collection can hold it?
Definition: lwgeom.c:1114
int32_t lwgeom_get_srid(const LWGEOM *geom)
Return SRID number.
Definition: lwgeom.c:910
int lwgeom_startpoint(const LWGEOM *lwgeom, POINT4D *pt)
Definition: lwgeom.c:2141
int lwgeom_is_collection(const LWGEOM *geom)
Determine whether a LWGEOM can contain sub-geometries or not.
Definition: lwgeom.c:1080
LWMPOINT * lwgeom_as_lwmpoint(const LWGEOM *lwgeom)
Definition: lwgeom.c:225
void lwgeom_set_srid(LWGEOM *geom, int32_t srid)
Set the SRID on an LWGEOM For collections, only the parent gets an SRID, all the children get SRID_UN...
Definition: lwgeom.c:1530
void lwgeom_trim_bits_in_place(LWGEOM *geom, int32_t prec_x, int32_t prec_y, int32_t prec_z, int32_t prec_m)
Trim the bits of an LWGEOM in place, to optimize it for compression.
Definition: lwgeom.c:2548
void lwgeom_longitude_shift(LWGEOM *lwgeom)
Definition: lwgeom.c:991
LWGEOM * lwgeom_segmentize2d(const LWGEOM *lwgeom, double dist)
Definition: lwgeom.c:754
LWGEOM * lwgeom_as_multi(const LWGEOM *lwgeom)
Create a new LWGEOM of the appropriate MULTI* type.
Definition: lwgeom.c:363
double lwgeom_perimeter_2d(const LWGEOM *geom)
Definition: lwgeom.c:1936
int lwpoint_inside_circle(const LWPOINT *p, double cx, double cy, double rad)
Definition: lwgeom.c:645
static int cmp_point_x(const void *pa, const void *pb)
Definition: lwgeom.c:1553
LWGEOM * lwmline_as_lwgeom(const LWMLINE *obj)
Definition: lwgeom.c:282
static uint8_t bits_for_precision(int32_t significant_digits)
Definition: lwgeom.c:2515
LWGEOM * lwmpoint_as_lwgeom(const LWMPOINT *obj)
Definition: lwgeom.c:287
LWGEOM * lwpoly_as_lwgeom(const LWPOLY *obj)
Definition: lwgeom.c:312
int lwgeom_is_trajectory(const LWGEOM *geom)
Return LW_TRUE or LW_FALSE depending on whether or not a geometry is a linestring with measure value ...
Definition: lwgeom.c:2502
LWGEOM * lwgeom_force_sfs(LWGEOM *geom, int version)
Definition: lwgeom.c:832
LWGEOM * lwgeom_clone_deep(const LWGEOM *lwgeom)
Deep-clone an LWGEOM object.
Definition: lwgeom.c:512
LWMPOLY * lwgeom_as_lwmpoly(const LWGEOM *lwgeom)
Definition: lwgeom.c:243
LWGEOM * lwgeom_force_3dm(const LWGEOM *geom, double mval)
Definition: lwgeom.c:788
int lwgeom_has_srid(const LWGEOM *geom)
Return true or false depending on whether a geometry has a valid SRID set.
Definition: lwgeom.c:1387
LWGEOM * lwgeom_force_4d(const LWGEOM *geom, double zval, double mval)
Definition: lwgeom.c:794
double lwgeom_length(const LWGEOM *geom)
Definition: lwgeom.c:1958
LWGEOM * lwgeom_remove_repeated_points(const LWGEOM *in, double tolerance)
Definition: lwgeom.c:1454
int lwgeom_needs_bbox(const LWGEOM *geom)
Check whether or not a lwgeom is big enough to warrant a bounding box.
Definition: lwgeom.c:1191
int lwgeom_remove_repeated_points_in_place(LWGEOM *geom, double tolerance)
Definition: lwgeom.c:1587
LWGEOM * lwpsurface_as_lwgeom(const LWPSURFACE *obj)
Definition: lwgeom.c:272
LWGEOM * lwgeom_force_3dz(const LWGEOM *geom, double zval)
Definition: lwgeom.c:782
int lwgeom_has_z(const LWGEOM *geom)
Return LW_TRUE if geometry has Z ordinates.
Definition: lwgeom.c:917
LWGEOM * lwmpoly_as_lwgeom(const LWMPOLY *obj)
Definition: lwgeom.c:277
int lwtype_is_collection(uint8_t type)
Return TRUE if the geometry may contain sub-geometries, i.e.
Definition: lwgeom.c:1088
char * lwgeom_to_ewkt(const LWGEOM *lwgeom)
Return an alloced string.
Definition: lwgeom.c:548
void lwgeom_drop_bbox(LWGEOM *lwgeom)
Call this function to drop BBOX and SRID from LWGEOM.
Definition: lwgeom.c:665
int lwgeom_isfinite(const LWGEOM *lwgeom)
Check if a LWGEOM has any non-finite (NaN or Inf) coordinates.
Definition: lwgeom.c:2690
static int lwcollection_dimensionality(const LWCOLLECTION *col)
Definition: lwgeom.c:1396
LWTRIANGLE * lwgeom_as_lwtriangle(const LWGEOM *lwgeom)
Definition: lwgeom.c:207
double lwgeom_area(const LWGEOM *geom)
Definition: lwgeom.c:1891
LWMLINE * lwgeom_as_lwmline(const LWGEOM *lwgeom)
Definition: lwgeom.c:234
uint32_t lwgeom_count_vertices(const LWGEOM *geom)
Count points in an LWGEOM.
Definition: lwgeom.c:1229
LWGEOM * lwgeom_reverse(const LWGEOM *geom)
Definition: lwgeom.c:94
int lwgeom_dimension(const LWGEOM *geom)
For an LWGEOM, returns 0 for points, 1 for lines, 2 for polygons, 3 for volume, and the max dimension...
Definition: lwgeom.c:1281
LWCOLLECTION * lwgeom_subdivide(const LWGEOM *geom, uint32_t maxvertices)
Definition: lwgeom.c:2494
LWGEOM * lwpoint_as_lwgeom(const LWPOINT *obj)
Definition: lwgeom.c:327
LWGEOM * lwgeom_as_curve(const LWGEOM *lwgeom)
Create a new LWGEOM of the appropriate CURVE* type.
Definition: lwgeom.c:403
void lwgeom_swap_ordinates(LWGEOM *in, LWORD o1, LWORD o2)
Swap ordinate values in every vertex of the geometry.
Definition: lwgeom.c:1461
void lwgeom_affine(LWGEOM *geom, const AFFINE *affine)
Definition: lwgeom.c:2003
uint8_t lwtype_multitype(uint8_t type)
Definition: lwgeom.c:353
LWGEOM * lwgeom_boundary(LWGEOM *lwgeom)
Definition: lwgeom.c:2569
LWCURVEPOLY * lwgeom_as_lwcurvepoly(const LWGEOM *lwgeom)
Definition: lwgeom.c:189
const GBOX * lwgeom_get_bbox(const LWGEOM *lwg)
Get a non-empty geometry bounding box, computing and caching it if not already there.
Definition: lwgeom.c:726
static int cmp_point_y(const void *pa, const void *pb)
Definition: lwgeom.c:1570
double lwgeom_length_2d(const LWGEOM *geom)
Definition: lwgeom.c:1980
uint32_t lwgeom_count_rings(const LWGEOM *geom)
Count rings in an LWGEOM.
Definition: lwgeom.c:1339
void lwgeom_reverse_in_place(LWGEOM *geom)
Reverse vertex order of LWGEOM.
Definition: lwgeom.c:103
void lwgeom_force_clockwise(LWGEOM *lwgeom)
Force Right-hand-rule on LWGEOM polygons.
Definition: lwgeom.c:38
int lwgeom_is_solid(const LWGEOM *geom)
Return LW_TRUE if geometry has SOLID flag.
Definition: lwgeom.c:931
LWGEOM * lwtin_as_lwgeom(const LWTIN *obj)
Definition: lwgeom.c:267
LWGEOM * lwcurvepoly_as_lwgeom(const LWCURVEPOLY *obj)
Definition: lwgeom.c:302
uint8_t MULTITYPE[NUMTYPES]
Look-up for the correct MULTI* type promotion for singleton types.
Definition: lwgeom.c:337
LWTIN * lwgeom_as_lwtin(const LWGEOM *lwgeom)
Definition: lwgeom.c:260
LWGEOM * lwgeom_clone(const LWGEOM *lwgeom)
Clone LWGEOM object.
Definition: lwgeom.c:474
int lwgeom_calculate_gbox(const LWGEOM *lwgeom, GBOX *gbox)
Calculate the gbox for this geometry, a cartesian box or geodetic box, depending on how it is flagged...
Definition: lwgeom.c:738
LWGEOM * lwcircstring_as_lwgeom(const LWCIRCSTRING *obj)
Definition: lwgeom.c:297
LWPSURFACE * lwgeom_as_lwpsurface(const LWGEOM *lwgeom)
Definition: lwgeom.c:252
void lwgeom_add_bbox_deep(LWGEOM *lwgeom, GBOX *gbox)
Compute a box for geom and all sub-geometries, if not already computed.
Definition: lwgeom.c:697
int lwgeom_dimensionality(const LWGEOM *geom)
Return the dimensionality (relating to point/line/poly) of an lwgeom.
Definition: lwgeom.c:1409
LWCOLLECTION * lwgeom_as_lwcollection(const LWGEOM *lwgeom)
Definition: lwgeom.c:216
void lwgeom_free(LWGEOM *lwgeom)
Definition: lwgeom.c:1138
LWGEOM * lwtriangle_as_lwgeom(const LWTRIANGLE *obj)
Definition: lwgeom.c:317
LWCIRCSTRING * lwgeom_as_lwcircstring(const LWGEOM *lwgeom)
Definition: lwgeom.c:171
LWPOLY * lwgeom_as_lwpoly(const LWGEOM *lwgeom)
Definition: lwgeom.c:198
static void lwgeom_subdivide_recursive(const LWGEOM *geom, uint8_t dimension, uint32_t maxvertices, uint32_t depth, LWCOLLECTION *col, double gridSize)
Definition: lwgeom.c:2290
void lwgeom_release(LWGEOM *lwgeom)
Free the containing LWGEOM and the associated BOX.
Definition: lwgeom.c:451
double lwgeom_perimeter(const LWGEOM *geom)
Definition: lwgeom.c:1914
int lwgeom_has_m(const LWGEOM *geom)
Return LW_TRUE if geometry has M ordinates.
Definition: lwgeom.c:924
void lwgeom_scale(LWGEOM *geom, const POINT4D *factor)
Definition: lwgeom.c:2057
LWCOLLECTION * lwgeom_subdivide_prec(const LWGEOM *geom, uint32_t maxvertices, double gridSize)
Definition: lwgeom.c:2471
LWGEOM * lwgeom_construct_empty(uint8_t type, int32_t srid, char hasz, char hasm)
Definition: lwgeom.c:2111
LWGEOM * lwgeom_grid(const LWGEOM *lwgeom, const gridspec *grid)
Definition: lwgeom.c:2273
int lwgeom_is_clockwise(LWGEOM *lwgeom)
Check clockwise orientation on LWGEOM polygons.
Definition: lwgeom.c:66
int lwgeom_simplify_in_place(LWGEOM *geom, double epsilon, int preserve_collapsed)
Definition: lwgeom.c:1743
void lwgeom_add_bbox(LWGEOM *lwgeom)
Ensure there's a box in the LWGEOM.
Definition: lwgeom.c:678
LWGEOM * lwgeom_force_2d(const LWGEOM *geom)
Strip out the Z/M components of an LWGEOM.
Definition: lwgeom.c:776
void lwgeom_drop_srid(LWGEOM *lwgeom)
Definition: lwgeom.c:748
void lwgeom_grid_in_place(LWGEOM *geom, const gridspec *grid)
Definition: lwgeom.c:2172
LWCOMPOUND * lwgeom_as_lwcompound(const LWGEOM *lwgeom)
Definition: lwgeom.c:180
#define LWDEBUG(level, msg)
Definition: lwgeom_log.h:83
#define LWDEBUGF(level, msg,...)
Definition: lwgeom_log.h:88
void lwerror(const char *fmt,...)
Write a notice out to the error handler.
Definition: lwutil.c:190
void lwnotice(const char *fmt,...)
Write a notice out to the notice handler.
Definition: lwutil.c:177
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:101
static double distance2d_sqr_pt_pt(const POINT2D *p1, const POINT2D *p2)
Definition: lwinline.h:35
static int lwgeom_is_empty(const LWGEOM *geom)
Return true or false depending on whether a geometry is an "empty" geometry (no vertices members)
Definition: lwinline.h:203
static LWPOINT * lwgeom_as_lwpoint(const LWGEOM *lwgeom)
Definition: lwinline.h:131
int value
Definition: genraster.py:62
type
Definition: ovdump.py:42
static double pivot(double *left, double *right)
Definition: rt_statistics.c:40
double ymax
Definition: liblwgeom.h:371
double xmax
Definition: liblwgeom.h:369
double ymin
Definition: liblwgeom.h:370
double xmin
Definition: liblwgeom.h:368
lwflags_t flags
Definition: liblwgeom.h:367
uint32_t ngeoms
Definition: liblwgeom.h:594
uint8_t type
Definition: liblwgeom.h:592
LWGEOM ** geoms
Definition: liblwgeom.h:589
LWGEOM ** rings
Definition: liblwgeom.h:617
uint32_t nrings
Definition: liblwgeom.h:622
uint8_t type
Definition: liblwgeom.h:476
GBOX * bbox
Definition: liblwgeom.h:472
int32_t srid
Definition: liblwgeom.h:474
lwflags_t flags
Definition: liblwgeom.h:475
POINTARRAY * points
Definition: liblwgeom.h:497
uint8_t type
Definition: liblwgeom.h:500
LWLINE ** geoms
Definition: liblwgeom.h:561
uint32_t ngeoms
Definition: liblwgeom.h:566
uint32_t ngeoms
Definition: liblwgeom.h:552
LWPOINT ** geoms
Definition: liblwgeom.h:547
POINTARRAY * point
Definition: liblwgeom.h:485
uint8_t type
Definition: liblwgeom.h:488
POINTARRAY ** rings
Definition: liblwgeom.h:533
uint32_t nrings
Definition: liblwgeom.h:538
POINTARRAY * points
Definition: liblwgeom.h:509
double y
Definition: liblwgeom.h:404
double x
Definition: liblwgeom.h:404
double m
Definition: liblwgeom.h:428
double x
Definition: liblwgeom.h:428
double z
Definition: liblwgeom.h:428
double y
Definition: liblwgeom.h:428
lwflags_t flags
Definition: liblwgeom.h:445
uint32_t maxpoints
Definition: liblwgeom.h:442
uint32_t npoints
Definition: liblwgeom.h:441
Snap-to-grid.
Definition: liblwgeom.h:1382