PostGIS  2.3.8dev-r@@SVN_REVISION@@
liblwgeom/lwgeom_geos.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 2011-2014 Sandro Santilli <strk@kbt.io>
22  *
23  **********************************************************************/
24 
25 
26 #include "lwgeom_geos.h"
27 #include "liblwgeom.h"
28 #include "liblwgeom_internal.h"
29 #include "lwgeom_log.h"
30 
31 #include <stdlib.h>
32 #include <time.h>
33 
34 LWTIN *lwtin_from_geos(const GEOSGeometry *geom, int want3d);
35 
36 #undef LWGEOM_PROFILE_BUILDAREA
37 
38 #define LWGEOM_GEOS_ERRMSG_MAXSIZE 256
40 
41 extern void
42 lwgeom_geos_error(const char *fmt, ...)
43 {
44  va_list ap;
45 
46  va_start(ap, fmt);
47 
48  /* Call the supplied function */
50  {
52  }
53 
54  va_end(ap);
55 }
56 
57 
58 /*
59 ** GEOS <==> PostGIS conversion functions
60 **
61 ** Default conversion creates a GEOS point array, then iterates through the
62 ** PostGIS points, setting each value in the GEOS array one at a time.
63 **
64 */
65 
66 /* Return a POINTARRAY from a GEOSCoordSeq */
67 POINTARRAY *
68 ptarray_from_GEOSCoordSeq(const GEOSCoordSequence *cs, char want3d)
69 {
70  uint32_t dims=2;
71  uint32_t size, i;
72  POINTARRAY *pa;
73  POINT4D point;
74 
75  LWDEBUG(2, "ptarray_fromGEOSCoordSeq called");
76 
77  if ( ! GEOSCoordSeq_getSize(cs, &size) )
78  lwerror("Exception thrown");
79 
80  LWDEBUGF(4, " GEOSCoordSeq size: %d", size);
81 
82  if ( want3d )
83  {
84  if ( ! GEOSCoordSeq_getDimensions(cs, &dims) )
85  lwerror("Exception thrown");
86 
87  LWDEBUGF(4, " GEOSCoordSeq dimensions: %d", dims);
88 
89  /* forget higher dimensions (if any) */
90  if ( dims > 3 ) dims = 3;
91  }
92 
93  LWDEBUGF(4, " output dimensions: %d", dims);
94 
95  pa = ptarray_construct((dims==3), 0, size);
96 
97  for (i=0; i<size; i++)
98  {
99  GEOSCoordSeq_getX(cs, i, &(point.x));
100  GEOSCoordSeq_getY(cs, i, &(point.y));
101  if ( dims >= 3 ) GEOSCoordSeq_getZ(cs, i, &(point.z));
102  ptarray_set_point4d(pa,i,&point);
103  }
104 
105  return pa;
106 }
107 
108 /* Return an LWGEOM from a Geometry */
109 LWGEOM *
110 GEOS2LWGEOM(const GEOSGeometry *geom, char want3d)
111 {
112  int type = GEOSGeomTypeId(geom) ;
113  int hasZ;
114  int SRID = GEOSGetSRID(geom);
115 
116  /* GEOS's 0 is equivalent to our unknown as for SRID values */
117  if ( SRID == 0 ) SRID = SRID_UNKNOWN;
118 
119  if ( want3d )
120  {
121  hasZ = GEOSHasZ(geom);
122  if ( ! hasZ )
123  {
124  LWDEBUG(3, "Geometry has no Z, won't provide one");
125 
126  want3d = 0;
127  }
128  }
129 
130 /*
131  if ( GEOSisEmpty(geom) )
132  {
133  return (LWGEOM*)lwcollection_construct_empty(COLLECTIONTYPE, SRID, want3d, 0);
134  }
135 */
136 
137  switch (type)
138  {
139  const GEOSCoordSequence *cs;
140  POINTARRAY *pa, **ppaa;
141  const GEOSGeometry *g;
142  LWGEOM **geoms;
143  uint32_t i, ngeoms;
144 
145  case GEOS_POINT:
146  LWDEBUG(4, "lwgeom_from_geometry: it's a Point");
147  cs = GEOSGeom_getCoordSeq(geom);
148  if ( GEOSisEmpty(geom) )
149  return (LWGEOM*)lwpoint_construct_empty(SRID, want3d, 0);
150  pa = ptarray_from_GEOSCoordSeq(cs, want3d);
151  return (LWGEOM *)lwpoint_construct(SRID, NULL, pa);
152 
153  case GEOS_LINESTRING:
154  case GEOS_LINEARRING:
155  LWDEBUG(4, "lwgeom_from_geometry: it's a LineString or LinearRing");
156  if ( GEOSisEmpty(geom) )
157  return (LWGEOM*)lwline_construct_empty(SRID, want3d, 0);
158 
159  cs = GEOSGeom_getCoordSeq(geom);
160  pa = ptarray_from_GEOSCoordSeq(cs, want3d);
161  return (LWGEOM *)lwline_construct(SRID, NULL, pa);
162 
163  case GEOS_POLYGON:
164  LWDEBUG(4, "lwgeom_from_geometry: it's a Polygon");
165  if ( GEOSisEmpty(geom) )
166  return (LWGEOM*)lwpoly_construct_empty(SRID, want3d, 0);
167  ngeoms = GEOSGetNumInteriorRings(geom);
168  ppaa = lwalloc(sizeof(POINTARRAY *)*(ngeoms+1));
169  g = GEOSGetExteriorRing(geom);
170  cs = GEOSGeom_getCoordSeq(g);
171  ppaa[0] = ptarray_from_GEOSCoordSeq(cs, want3d);
172  for (i=0; i<ngeoms; i++)
173  {
174  g = GEOSGetInteriorRingN(geom, i);
175  cs = GEOSGeom_getCoordSeq(g);
176  ppaa[i+1] = ptarray_from_GEOSCoordSeq(cs,
177  want3d);
178  }
179  return (LWGEOM *)lwpoly_construct(SRID, NULL,
180  ngeoms+1, ppaa);
181 
182  case GEOS_MULTIPOINT:
183  case GEOS_MULTILINESTRING:
184  case GEOS_MULTIPOLYGON:
185  case GEOS_GEOMETRYCOLLECTION:
186  LWDEBUG(4, "lwgeom_from_geometry: it's a Collection or Multi");
187 
188  ngeoms = GEOSGetNumGeometries(geom);
189  geoms = NULL;
190  if ( ngeoms )
191  {
192  geoms = lwalloc(sizeof(LWGEOM *)*ngeoms);
193  for (i=0; i<ngeoms; i++)
194  {
195  g = GEOSGetGeometryN(geom, i);
196  geoms[i] = GEOS2LWGEOM(g, want3d);
197  }
198  }
199  return (LWGEOM *)lwcollection_construct(type,
200  SRID, NULL, ngeoms, geoms);
201 
202  default:
203  lwerror("GEOS2LWGEOM: unknown geometry type: %d", type);
204  return NULL;
205 
206  }
207 
208 }
209 
210 
211 
212 GEOSCoordSeq ptarray_to_GEOSCoordSeq(const POINTARRAY *);
213 
214 
215 GEOSCoordSeq
217 {
218  uint32_t dims = 2;
219  uint32_t i;
220  const POINT3DZ *p3d;
221  const POINT2D *p2d;
222  GEOSCoordSeq sq;
223 
224  if ( FLAGS_GET_Z(pa->flags) )
225  dims = 3;
226 
227  if ( ! (sq = GEOSCoordSeq_create(pa->npoints, dims)) )
228  lwerror("Error creating GEOS Coordinate Sequence");
229 
230  for ( i=0; i < pa->npoints; i++ )
231  {
232  if ( dims == 3 )
233  {
234  p3d = getPoint3dz_cp(pa, i);
235  p2d = (const POINT2D *)p3d;
236  LWDEBUGF(4, "Point: %g,%g,%g", p3d->x, p3d->y, p3d->z);
237  }
238  else
239  {
240  p2d = getPoint2d_cp(pa, i);
241  LWDEBUGF(4, "Point: %g,%g", p2d->x, p2d->y);
242  }
243 
244  GEOSCoordSeq_setX(sq, i, p2d->x);
245  GEOSCoordSeq_setY(sq, i, p2d->y);
246 
247  if ( dims == 3 )
248  GEOSCoordSeq_setZ(sq, i, p3d->z);
249  }
250  return sq;
251 }
252 
253 static GEOSGeometry *
254 ptarray_to_GEOSLinearRing(const POINTARRAY *pa, int autofix)
255 {
256  GEOSCoordSeq sq;
257  GEOSGeom g;
258  POINTARRAY *npa = 0;
259 
260  if ( autofix )
261  {
262  /* check ring for being closed and fix if not */
263  if ( ! ptarray_is_closed_2d(pa) )
264  {
265  npa = ptarray_addPoint(pa, getPoint_internal(pa, 0), FLAGS_NDIMS(pa->flags), pa->npoints);
266  pa = npa;
267  }
268  /* TODO: check ring for having at least 4 vertices */
269 #if 0
270  while ( pa->npoints < 4 )
271  {
272  npa = ptarray_addPoint(npa, getPoint_internal(pa, 0), FLAGS_NDIMS(pa->flags), pa->npoints);
273  }
274 #endif
275  }
276 
277  sq = ptarray_to_GEOSCoordSeq(pa);
278  if ( npa ) ptarray_free(npa);
279  g = GEOSGeom_createLinearRing(sq);
280  return g;
281 }
282 
283 GEOSGeometry *
284 GBOX2GEOS(const GBOX *box)
285 {
286  GEOSGeometry* envelope;
287  GEOSGeometry* ring;
288  GEOSCoordSequence* seq = GEOSCoordSeq_create(5, 2);
289  if (!seq)
290  {
291  return NULL;
292  }
293 
294  GEOSCoordSeq_setX(seq, 0, box->xmin);
295  GEOSCoordSeq_setY(seq, 0, box->ymin);
296 
297  GEOSCoordSeq_setX(seq, 1, box->xmax);
298  GEOSCoordSeq_setY(seq, 1, box->ymin);
299 
300  GEOSCoordSeq_setX(seq, 2, box->xmax);
301  GEOSCoordSeq_setY(seq, 2, box->ymax);
302 
303  GEOSCoordSeq_setX(seq, 3, box->xmin);
304  GEOSCoordSeq_setY(seq, 3, box->ymax);
305 
306  GEOSCoordSeq_setX(seq, 4, box->xmin);
307  GEOSCoordSeq_setY(seq, 4, box->ymin);
308 
309  ring = GEOSGeom_createLinearRing(seq);
310  if (!ring)
311  {
312  GEOSCoordSeq_destroy(seq);
313  return NULL;
314  }
315 
316  envelope = GEOSGeom_createPolygon(ring, NULL, 0);
317  if (!envelope)
318  {
319  GEOSGeom_destroy(ring);
320  return NULL;
321  }
322 
323  return envelope;
324 }
325 
326 GEOSGeometry *
327 LWGEOM2GEOS(const LWGEOM *lwgeom, int autofix)
328 {
329  GEOSCoordSeq sq;
330  GEOSGeom g, shell;
331  GEOSGeom *geoms = NULL;
332  /*
333  LWGEOM *tmp;
334  */
335  uint32_t ngeoms, i, j;
336  int geostype;
337 #if LWDEBUG_LEVEL >= 4
338  char *wkt;
339 #endif
340 
341  LWDEBUGF(4, "LWGEOM2GEOS got a %s", lwtype_name(lwgeom->type));
342 
343  if (lwgeom_has_arc(lwgeom))
344  {
345  LWGEOM *lwgeom_stroked = lwgeom_stroke(lwgeom, 32);
346  GEOSGeometry *g = LWGEOM2GEOS(lwgeom_stroked, autofix);
347  lwgeom_free(lwgeom_stroked);
348  return g;
349  }
350 
351  switch (lwgeom->type)
352  {
353  LWPOINT *lwp = NULL;
354  LWPOLY *lwpoly = NULL;
355  LWLINE *lwl = NULL;
356  LWCOLLECTION *lwc = NULL;
357 
358  case POINTTYPE:
359  lwp = (LWPOINT *)lwgeom;
360 
361  if ( lwgeom_is_empty(lwgeom) )
362  {
363  g = GEOSGeom_createEmptyPolygon();
364  }
365  else
366  {
367  sq = ptarray_to_GEOSCoordSeq(lwp->point);
368  g = GEOSGeom_createPoint(sq);
369  }
370  if ( ! g )
371  {
372  /* lwnotice("Exception in LWGEOM2GEOS"); */
373  return NULL;
374  }
375  break;
376  case LINETYPE:
377  lwl = (LWLINE *)lwgeom;
378  /* TODO: if (autofix) */
379  if ( lwl->points->npoints == 1 ) {
380  /* Duplicate point, to make geos-friendly */
381  lwl->points = ptarray_addPoint(lwl->points,
382  getPoint_internal(lwl->points, 0),
383  FLAGS_NDIMS(lwl->points->flags),
384  lwl->points->npoints);
385  }
386  sq = ptarray_to_GEOSCoordSeq(lwl->points);
387  g = GEOSGeom_createLineString(sq);
388  if ( ! g )
389  {
390  /* lwnotice("Exception in LWGEOM2GEOS"); */
391  return NULL;
392  }
393  break;
394 
395  case POLYGONTYPE:
396  lwpoly = (LWPOLY *)lwgeom;
397  if ( lwgeom_is_empty(lwgeom) )
398  {
399  g = GEOSGeom_createEmptyPolygon();
400  }
401  else
402  {
403  shell = ptarray_to_GEOSLinearRing(lwpoly->rings[0], autofix);
404  if ( ! shell ) return NULL;
405  /*lwerror("LWGEOM2GEOS: exception during polygon shell conversion"); */
406  ngeoms = lwpoly->nrings-1;
407  if ( ngeoms > 0 )
408  geoms = malloc(sizeof(GEOSGeom)*ngeoms);
409 
410  for (i=1; i<lwpoly->nrings; ++i)
411  {
412  geoms[i-1] = ptarray_to_GEOSLinearRing(lwpoly->rings[i], autofix);
413  if ( ! geoms[i-1] )
414  {
415  --i;
416  while (i) GEOSGeom_destroy(geoms[--i]);
417  free(geoms);
418  GEOSGeom_destroy(shell);
419  return NULL;
420  }
421  /*lwerror("LWGEOM2GEOS: exception during polygon hole conversion"); */
422  }
423  g = GEOSGeom_createPolygon(shell, geoms, ngeoms);
424  if (geoms) free(geoms);
425  }
426  if ( ! g ) return NULL;
427  break;
428  case MULTIPOINTTYPE:
429  case MULTILINETYPE:
430  case MULTIPOLYGONTYPE:
431  case COLLECTIONTYPE:
432  if ( lwgeom->type == MULTIPOINTTYPE )
433  geostype = GEOS_MULTIPOINT;
434  else if ( lwgeom->type == MULTILINETYPE )
435  geostype = GEOS_MULTILINESTRING;
436  else if ( lwgeom->type == MULTIPOLYGONTYPE )
437  geostype = GEOS_MULTIPOLYGON;
438  else
439  geostype = GEOS_GEOMETRYCOLLECTION;
440 
441  lwc = (LWCOLLECTION *)lwgeom;
442 
443  ngeoms = lwc->ngeoms;
444  if ( ngeoms > 0 )
445  geoms = malloc(sizeof(GEOSGeom)*ngeoms);
446 
447  j = 0;
448  for (i=0; i<ngeoms; ++i)
449  {
450  GEOSGeometry* g;
451 
452  if( lwgeom_is_empty(lwc->geoms[i]) )
453  continue;
454 
455  g = LWGEOM2GEOS(lwc->geoms[i], 0);
456  if ( ! g )
457  {
458  while (j) GEOSGeom_destroy(geoms[--j]);
459  free(geoms);
460  return NULL;
461  }
462  geoms[j++] = g;
463  }
464  g = GEOSGeom_createCollection(geostype, geoms, j);
465  if ( geoms ) free(geoms);
466  if ( ! g ) return NULL;
467  break;
468 
469  default:
470  lwerror("Unknown geometry type: %d - %s", lwgeom->type, lwtype_name(lwgeom->type));
471  return NULL;
472  }
473 
474  GEOSSetSRID(g, lwgeom->srid);
475 
476 #if LWDEBUG_LEVEL >= 4
477  wkt = GEOSGeomToWKT(g);
478  LWDEBUGF(4, "LWGEOM2GEOS: GEOSGeom: %s", wkt);
479  free(wkt);
480 #endif
481 
482  return g;
483 }
484 
485 GEOSGeometry*
486 make_geos_point(double x, double y)
487 {
488  GEOSCoordSequence* seq = GEOSCoordSeq_create(1, 2);
489  GEOSGeometry* geom = NULL;
490 
491  if (!seq)
492  return NULL;
493 
494  GEOSCoordSeq_setX(seq, 0, x);
495  GEOSCoordSeq_setY(seq, 0, y);
496 
497  geom = GEOSGeom_createPoint(seq);
498  if (!geom)
499  GEOSCoordSeq_destroy(seq);
500  return geom;
501 }
502 
503 GEOSGeometry*
504 make_geos_segment(double x1, double y1, double x2, double y2)
505 {
506  GEOSCoordSequence* seq = GEOSCoordSeq_create(2, 2);
507  GEOSGeometry* geom = NULL;
508 
509  if (!seq)
510  return NULL;
511 
512  GEOSCoordSeq_setX(seq, 0, x1);
513  GEOSCoordSeq_setY(seq, 0, y1);
514  GEOSCoordSeq_setX(seq, 1, x2);
515  GEOSCoordSeq_setY(seq, 1, y2);
516 
517  geom = GEOSGeom_createLineString(seq);
518  if (!geom)
519  GEOSCoordSeq_destroy(seq);
520  return geom;
521 }
522 
523 const char*
525 {
526  const char *ver = GEOSversion();
527  return ver;
528 }
529 
530 LWGEOM *
532 {
533  LWGEOM *result ;
534  GEOSGeometry *g1;
535  int is3d ;
536  int srid ;
537 
538  srid = (int)(geom1->srid);
539  is3d = FLAGS_GET_Z(geom1->flags);
540 
541  initGEOS(lwnotice, lwgeom_geos_error);
542 
543  g1 = LWGEOM2GEOS(geom1, 0);
544  if ( 0 == g1 ) /* exception thrown at construction */
545  {
546  lwerror("First argument geometry could not be converted to GEOS: %s", lwgeom_geos_errmsg);
547  return NULL ;
548  }
549 
550  if ( -1 == GEOSNormalize(g1) )
551  {
552  lwerror("Error in GEOSNormalize: %s", lwgeom_geos_errmsg);
553  return NULL; /* never get here */
554  }
555 
556  GEOSSetSRID(g1, srid); /* needed ? */
557  result = GEOS2LWGEOM(g1, is3d);
558  GEOSGeom_destroy(g1);
559 
560  if (result == NULL)
561  {
562  lwerror("Error performing intersection: GEOS2LWGEOM: %s",
564  return NULL ; /* never get here */
565  }
566 
567  return result ;
568 }
569 
570 LWGEOM *
571 lwgeom_intersection(const LWGEOM *geom1, const LWGEOM *geom2)
572 {
573  LWGEOM *result ;
574  GEOSGeometry *g1, *g2, *g3 ;
575  int is3d ;
576  int srid ;
577 
578  /* A.Intersection(Empty) == Empty */
579  if ( lwgeom_is_empty(geom2) )
580  return lwgeom_clone_deep(geom2);
581 
582  /* Empty.Intersection(A) == Empty */
583  if ( lwgeom_is_empty(geom1) )
584  return lwgeom_clone_deep(geom1);
585 
586  /* ensure srids are identical */
587  srid = (int)(geom1->srid);
588  error_if_srid_mismatch(srid, (int)(geom2->srid));
589 
590  is3d = (FLAGS_GET_Z(geom1->flags) || FLAGS_GET_Z(geom2->flags)) ;
591 
592  initGEOS(lwnotice, lwgeom_geos_error);
593 
594  LWDEBUG(3, "intersection() START");
595 
596  g1 = LWGEOM2GEOS(geom1, 0);
597  if ( 0 == g1 ) /* exception thrown at construction */
598  {
599  lwerror("First argument geometry could not be converted to GEOS: %s", lwgeom_geos_errmsg);
600  return NULL ;
601  }
602 
603  g2 = LWGEOM2GEOS(geom2, 0);
604  if ( 0 == g2 ) /* exception thrown at construction */
605  {
606  lwerror("Second argument geometry could not be converted to GEOS.");
607  GEOSGeom_destroy(g1);
608  return NULL ;
609  }
610 
611  LWDEBUG(3, " constructed geometrys - calling geos");
612  LWDEBUGF(3, " g1 = %s", GEOSGeomToWKT(g1));
613  LWDEBUGF(3, " g2 = %s", GEOSGeomToWKT(g2));
614  /*LWDEBUGF(3, "g2 is valid = %i",GEOSisvalid(g2)); */
615  /*LWDEBUGF(3, "g1 is valid = %i",GEOSisvalid(g1)); */
616 
617  g3 = GEOSIntersection(g1,g2);
618 
619  LWDEBUG(3, " intersection finished");
620 
621  if (g3 == NULL)
622  {
623  GEOSGeom_destroy(g1);
624  GEOSGeom_destroy(g2);
625  lwerror("Error performing intersection: %s",
627  return NULL; /* never get here */
628  }
629 
630  LWDEBUGF(3, "result: %s", GEOSGeomToWKT(g3) ) ;
631 
632  GEOSSetSRID(g3, srid);
633 
634  result = GEOS2LWGEOM(g3, is3d);
635 
636  if (result == NULL)
637  {
638  GEOSGeom_destroy(g1);
639  GEOSGeom_destroy(g2);
640  GEOSGeom_destroy(g3);
641  lwerror("Error performing intersection: GEOS2LWGEOM: %s",
643  return NULL ; /* never get here */
644  }
645 
646  GEOSGeom_destroy(g1);
647  GEOSGeom_destroy(g2);
648  GEOSGeom_destroy(g3);
649 
650  return result ;
651 }
652 
653 LWGEOM *
655 {
656  LWGEOM *result ;
657  GEOSGeometry *g1, *g3 ;
658  int is3d = FLAGS_GET_Z(geom1->flags);
659  int srid = geom1->srid;
660 
661  /* Empty.Linemerge() == Empty */
662  if ( lwgeom_is_empty(geom1) )
663  return (LWGEOM*)lwcollection_construct_empty( COLLECTIONTYPE, srid, is3d,
664  lwgeom_has_m(geom1) );
665 
666  initGEOS(lwnotice, lwgeom_geos_error);
667 
668  LWDEBUG(3, "linemerge() START");
669 
670  g1 = LWGEOM2GEOS(geom1, 0);
671  if ( 0 == g1 ) /* exception thrown at construction */
672  {
673  lwerror("First argument geometry could not be converted to GEOS: %s", lwgeom_geos_errmsg);
674  return NULL ;
675  }
676 
677  LWDEBUG(3, " constructed geometrys - calling geos");
678  LWDEBUGF(3, " g1 = %s", GEOSGeomToWKT(g1));
679  /*LWDEBUGF(3, "g1 is valid = %i",GEOSisvalid(g1)); */
680 
681  g3 = GEOSLineMerge(g1);
682 
683  LWDEBUG(3, " linemerge finished");
684 
685  if (g3 == NULL)
686  {
687  GEOSGeom_destroy(g1);
688  lwerror("Error performing linemerge: %s",
690  return NULL; /* never get here */
691  }
692 
693  LWDEBUGF(3, "result: %s", GEOSGeomToWKT(g3) ) ;
694 
695  GEOSSetSRID(g3, srid);
696 
697  result = GEOS2LWGEOM(g3, is3d);
698 
699  if (result == NULL)
700  {
701  GEOSGeom_destroy(g1);
702  GEOSGeom_destroy(g3);
703  lwerror("Error performing linemerge: GEOS2LWGEOM: %s",
705  return NULL ; /* never get here */
706  }
707 
708  GEOSGeom_destroy(g1);
709  GEOSGeom_destroy(g3);
710 
711  return result ;
712 }
713 
714 LWGEOM *
716 {
717  LWGEOM *result ;
718  GEOSGeometry *g1, *g3 ;
719  int is3d = FLAGS_GET_Z(geom1->flags);
720  int srid = geom1->srid;
721 
722  /* Empty.UnaryUnion() == Empty */
723  if ( lwgeom_is_empty(geom1) )
724  return lwgeom_clone_deep(geom1);
725 
726  initGEOS(lwnotice, lwgeom_geos_error);
727 
728  g1 = LWGEOM2GEOS(geom1, 0);
729  if ( 0 == g1 ) /* exception thrown at construction */
730  {
731  lwerror("First argument geometry could not be converted to GEOS: %s", lwgeom_geos_errmsg);
732  return NULL ;
733  }
734 
735  g3 = GEOSUnaryUnion(g1);
736 
737  if (g3 == NULL)
738  {
739  GEOSGeom_destroy(g1);
740  lwerror("Error performing unaryunion: %s",
742  return NULL; /* never get here */
743  }
744 
745  LWDEBUGF(3, "result: %s", GEOSGeomToWKT(g3) ) ;
746 
747  GEOSSetSRID(g3, srid);
748 
749  result = GEOS2LWGEOM(g3, is3d);
750 
751  if (result == NULL)
752  {
753  GEOSGeom_destroy(g1);
754  GEOSGeom_destroy(g3);
755  lwerror("Error performing unaryunion: GEOS2LWGEOM: %s",
757  return NULL ; /* never get here */
758  }
759 
760  GEOSGeom_destroy(g1);
761  GEOSGeom_destroy(g3);
762 
763  return result ;
764 }
765 
766 LWGEOM *
767 lwgeom_difference(const LWGEOM *geom1, const LWGEOM *geom2)
768 {
769  GEOSGeometry *g1, *g2, *g3;
770  LWGEOM *result;
771  int is3d;
772  int srid;
773 
774  /* A.Difference(Empty) == A */
775  if ( lwgeom_is_empty(geom2) )
776  return lwgeom_clone_deep(geom1);
777 
778  /* Empty.Intersection(A) == Empty */
779  if ( lwgeom_is_empty(geom1) )
780  return lwgeom_clone_deep(geom1);
781 
782  /* ensure srids are identical */
783  srid = (int)(geom1->srid);
784  error_if_srid_mismatch(srid, (int)(geom2->srid));
785 
786  is3d = (FLAGS_GET_Z(geom1->flags) || FLAGS_GET_Z(geom2->flags)) ;
787 
788  initGEOS(lwnotice, lwgeom_geos_error);
789 
790  g1 = LWGEOM2GEOS(geom1, 0);
791  if ( 0 == g1 ) /* exception thrown at construction */
792  {
793  lwerror("First argument geometry could not be converted to GEOS: %s", lwgeom_geos_errmsg);
794  return NULL;
795  }
796 
797  g2 = LWGEOM2GEOS(geom2, 0);
798  if ( 0 == g2 ) /* exception thrown at construction */
799  {
800  GEOSGeom_destroy(g1);
801  lwerror("Second argument geometry could not be converted to GEOS: %s", lwgeom_geos_errmsg);
802  return NULL;
803  }
804 
805  g3 = GEOSDifference(g1,g2);
806 
807  if (g3 == NULL)
808  {
809  GEOSGeom_destroy(g1);
810  GEOSGeom_destroy(g2);
811  lwerror("GEOSDifference: %s", lwgeom_geos_errmsg);
812  return NULL ; /* never get here */
813  }
814 
815  LWDEBUGF(3, "result: %s", GEOSGeomToWKT(g3) ) ;
816 
817  GEOSSetSRID(g3, srid);
818 
819  result = GEOS2LWGEOM(g3, is3d);
820 
821  if (result == NULL)
822  {
823  GEOSGeom_destroy(g1);
824  GEOSGeom_destroy(g2);
825  GEOSGeom_destroy(g3);
826  lwerror("Error performing difference: GEOS2LWGEOM: %s",
828  return NULL; /* never get here */
829  }
830 
831  GEOSGeom_destroy(g1);
832  GEOSGeom_destroy(g2);
833  GEOSGeom_destroy(g3);
834 
835  /* compressType(result); */
836 
837  return result;
838 }
839 
840 LWGEOM *
841 lwgeom_symdifference(const LWGEOM* geom1, const LWGEOM* geom2)
842 {
843  GEOSGeometry *g1, *g2, *g3;
844  LWGEOM *result;
845  int is3d;
846  int srid;
847 
848  /* A.SymDifference(Empty) == A */
849  if ( lwgeom_is_empty(geom2) )
850  return lwgeom_clone_deep(geom1);
851 
852  /* Empty.DymDifference(B) == B */
853  if ( lwgeom_is_empty(geom1) )
854  return lwgeom_clone_deep(geom2);
855 
856  /* ensure srids are identical */
857  srid = (int)(geom1->srid);
858  error_if_srid_mismatch(srid, (int)(geom2->srid));
859 
860  is3d = (FLAGS_GET_Z(geom1->flags) || FLAGS_GET_Z(geom2->flags)) ;
861 
862  initGEOS(lwnotice, lwgeom_geos_error);
863 
864  g1 = LWGEOM2GEOS(geom1, 0);
865 
866  if ( 0 == g1 ) /* exception thrown at construction */
867  {
868  lwerror("First argument geometry could not be converted to GEOS: %s", lwgeom_geos_errmsg);
869  return NULL;
870  }
871 
872  g2 = LWGEOM2GEOS(geom2, 0);
873 
874  if ( 0 == g2 ) /* exception thrown at construction */
875  {
876  lwerror("Second argument geometry could not be converted to GEOS: %s", lwgeom_geos_errmsg);
877  GEOSGeom_destroy(g1);
878  return NULL;
879  }
880 
881  g3 = GEOSSymDifference(g1,g2);
882 
883  if (g3 == NULL)
884  {
885  GEOSGeom_destroy(g1);
886  GEOSGeom_destroy(g2);
887  lwerror("GEOSSymDifference: %s", lwgeom_geos_errmsg);
888  return NULL; /*never get here */
889  }
890 
891  LWDEBUGF(3, "result: %s", GEOSGeomToWKT(g3));
892 
893  GEOSSetSRID(g3, srid);
894 
895  result = GEOS2LWGEOM(g3, is3d);
896 
897  if (result == NULL)
898  {
899  GEOSGeom_destroy(g1);
900  GEOSGeom_destroy(g2);
901  GEOSGeom_destroy(g3);
902  lwerror("GEOS symdifference() threw an error (result postgis geometry formation)!");
903  return NULL ; /*never get here */
904  }
905 
906  GEOSGeom_destroy(g1);
907  GEOSGeom_destroy(g2);
908  GEOSGeom_destroy(g3);
909 
910  return result;
911 }
912 
913 LWGEOM *
915 {
916  GEOSGeometry *g, *g_centroid;
917  LWGEOM *centroid;
918  int srid, is3d;
919 
920  if (lwgeom_is_empty(geom))
921  {
923  lwgeom_get_srid(geom),
924  lwgeom_has_z(geom),
925  lwgeom_has_m(geom));
926  return lwpoint_as_lwgeom(lwp);
927  }
928 
929  srid = lwgeom_get_srid(geom);
930  is3d = lwgeom_has_z(geom);
931 
932  initGEOS(lwnotice, lwgeom_geos_error);
933 
934  g = LWGEOM2GEOS(geom, 0);
935 
936  if (0 == g) /* exception thrown at construction */
937  {
938  lwerror("Geometry could not be converted to GEOS: %s", lwgeom_geos_errmsg);
939  return NULL;
940  }
941 
942  g_centroid = GEOSGetCentroid(g);
943  GEOSGeom_destroy(g);
944 
945  if (g_centroid == NULL)
946  {
947  lwerror("GEOSGetCentroid: %s", lwgeom_geos_errmsg);
948  return NULL; /*never get here */
949  }
950 
951  LWDEBUGF(3, "result: %s", GEOSGeomToWKT(g_centroid));
952 
953  GEOSSetSRID(g_centroid, srid);
954 
955  centroid = GEOS2LWGEOM(g_centroid, is3d);
956  GEOSGeom_destroy(g_centroid);
957 
958  if (centroid == NULL)
959  {
960  lwerror("GEOS GEOSGetCentroid() threw an error (result postgis geometry formation)!");
961  return NULL ; /*never get here */
962  }
963 
964  return centroid;
965 }
966 
967 LWGEOM*
968 lwgeom_union(const LWGEOM *geom1, const LWGEOM *geom2)
969 {
970  int is3d;
971  int srid;
972  GEOSGeometry *g1, *g2, *g3;
973  LWGEOM *result;
974 
975  LWDEBUG(2, "in geomunion");
976 
977  /* A.Union(empty) == A */
978  if ( lwgeom_is_empty(geom1) )
979  return lwgeom_clone_deep(geom2);
980 
981  /* B.Union(empty) == B */
982  if ( lwgeom_is_empty(geom2) )
983  return lwgeom_clone_deep(geom1);
984 
985 
986  /* ensure srids are identical */
987  srid = (int)(geom1->srid);
988  error_if_srid_mismatch(srid, (int)(geom2->srid));
989 
990  is3d = (FLAGS_GET_Z(geom1->flags) || FLAGS_GET_Z(geom2->flags)) ;
991 
992  initGEOS(lwnotice, lwgeom_geos_error);
993 
994  g1 = LWGEOM2GEOS(geom1, 0);
995 
996  if ( 0 == g1 ) /* exception thrown at construction */
997  {
998  lwerror("First argument geometry could not be converted to GEOS: %s", lwgeom_geos_errmsg);
999  return NULL;
1000  }
1001 
1002  g2 = LWGEOM2GEOS(geom2, 0);
1003 
1004  if ( 0 == g2 ) /* exception thrown at construction */
1005  {
1006  GEOSGeom_destroy(g1);
1007  lwerror("Second argument geometry could not be converted to GEOS: %s", lwgeom_geos_errmsg);
1008  return NULL;
1009  }
1010 
1011  LWDEBUGF(3, "g1=%s", GEOSGeomToWKT(g1));
1012  LWDEBUGF(3, "g2=%s", GEOSGeomToWKT(g2));
1013 
1014  g3 = GEOSUnion(g1,g2);
1015 
1016  LWDEBUGF(3, "g3=%s", GEOSGeomToWKT(g3));
1017 
1018  GEOSGeom_destroy(g1);
1019  GEOSGeom_destroy(g2);
1020 
1021  if (g3 == NULL)
1022  {
1023  lwerror("GEOSUnion: %s", lwgeom_geos_errmsg);
1024  return NULL; /* never get here */
1025  }
1026 
1027 
1028  GEOSSetSRID(g3, srid);
1029 
1030  result = GEOS2LWGEOM(g3, is3d);
1031 
1032  GEOSGeom_destroy(g3);
1033 
1034  if (result == NULL)
1035  {
1036  lwerror("Error performing union: GEOS2LWGEOM: %s",
1038  return NULL; /*never get here */
1039  }
1040 
1041  return result;
1042 }
1043 
1044 LWGEOM *
1045 lwgeom_clip_by_rect(const LWGEOM *geom1, double x0, double y0, double x1, double y1)
1046 {
1047 #if POSTGIS_GEOS_VERSION < 35
1048  lwerror("The GEOS version this postgis binary "
1049  "was compiled against (%d) doesn't support "
1050  "'GEOSClipByRect' function (3.5.0+ required)",
1052  return NULL;
1053 #else /* POSTGIS_GEOS_VERSION >= 35 */
1054  LWGEOM *result ;
1055  GEOSGeometry *g1, *g3 ;
1056  int is3d ;
1057 
1058  /* A.Intersection(Empty) == Empty */
1059  if ( lwgeom_is_empty(geom1) )
1060  return lwgeom_clone_deep(geom1);
1061 
1062  is3d = FLAGS_GET_Z(geom1->flags);
1063 
1064  initGEOS(lwnotice, lwgeom_geos_error);
1065 
1066  LWDEBUG(3, "clip_by_rect() START");
1067 
1068  g1 = LWGEOM2GEOS(geom1, 1); /* auto-fix structure */
1069  if ( 0 == g1 ) /* exception thrown at construction */
1070  {
1071  lwerror("First argument geometry could not be converted to GEOS: %s", lwgeom_geos_errmsg);
1072  return NULL ;
1073  }
1074 
1075  LWDEBUG(3, " constructed geometrys - calling geos");
1076  LWDEBUGF(3, " g1 = %s", GEOSGeomToWKT(g1));
1077  /*LWDEBUGF(3, "g1 is valid = %i",GEOSisvalid(g1)); */
1078 
1079  g3 = GEOSClipByRect(g1,x0,y0,x1,y1);
1080  GEOSGeom_destroy(g1);
1081 
1082  LWDEBUG(3, " clip_by_rect finished");
1083 
1084  if (g3 == NULL)
1085  {
1086  lwnotice("Error performing rectangular clipping: %s", lwgeom_geos_errmsg);
1087  return NULL;
1088  }
1089 
1090  LWDEBUGF(3, "result: %s", GEOSGeomToWKT(g3) ) ;
1091 
1092  result = GEOS2LWGEOM(g3, is3d);
1093  GEOSGeom_destroy(g3);
1094 
1095  if (result == NULL)
1096  {
1097  lwerror("Error performing intersection: GEOS2LWGEOM: %s", lwgeom_geos_errmsg);
1098  return NULL ; /* never get here */
1099  }
1100 
1101  result->srid = geom1->srid;
1102 
1103  return result ;
1104 #endif /* POSTGIS_GEOS_VERSION >= 35 */
1105 }
1106 
1107 
1108 /* ------------ BuildArea stuff ---------------------------------------------------------------------{ */
1109 
1110 typedef struct Face_t {
1111  const GEOSGeometry* geom;
1112  GEOSGeometry* env;
1113  double envarea;
1114  struct Face_t* parent; /* if this face is an hole of another one, or NULL */
1115 } Face;
1116 
1117 static Face* newFace(const GEOSGeometry* g);
1118 static void delFace(Face* f);
1119 static unsigned int countParens(const Face* f);
1120 static void findFaceHoles(Face** faces, int nfaces);
1121 
1122 static Face*
1123 newFace(const GEOSGeometry* g)
1124 {
1125  Face* f = lwalloc(sizeof(Face));
1126  f->geom = g;
1127  f->env = GEOSEnvelope(f->geom);
1128  GEOSArea(f->env, &f->envarea);
1129  f->parent = NULL;
1130  /* lwnotice("Built Face with area %g and %d holes", f->envarea, GEOSGetNumInteriorRings(f->geom)); */
1131  return f;
1132 }
1133 
1134 static unsigned int
1136 {
1137  unsigned int pcount = 0;
1138  while ( f->parent ) {
1139  ++pcount;
1140  f = f->parent;
1141  }
1142  return pcount;
1143 }
1144 
1145 /* Destroy the face and release memory associated with it */
1146 static void
1148 {
1149  GEOSGeom_destroy(f->env);
1150  lwfree(f);
1151 }
1152 
1153 
1154 static int
1155 compare_by_envarea(const void* g1, const void* g2)
1156 {
1157  Face* f1 = *(Face**)g1;
1158  Face* f2 = *(Face**)g2;
1159  double n1 = f1->envarea;
1160  double n2 = f2->envarea;
1161 
1162  if ( n1 < n2 ) return 1;
1163  if ( n1 > n2 ) return -1;
1164  return 0;
1165 }
1166 
1167 /* Find holes of each face */
1168 static void
1169 findFaceHoles(Face** faces, int nfaces)
1170 {
1171  int i, j, h;
1172 
1173  /* We sort by envelope area so that we know holes are only
1174  * after their shells */
1175  qsort(faces, nfaces, sizeof(Face*), compare_by_envarea);
1176  for (i=0; i<nfaces; ++i) {
1177  Face* f = faces[i];
1178  int nholes = GEOSGetNumInteriorRings(f->geom);
1179  LWDEBUGF(2, "Scanning face %d with env area %g and %d holes", i, f->envarea, nholes);
1180  for (h=0; h<nholes; ++h) {
1181  const GEOSGeometry *hole = GEOSGetInteriorRingN(f->geom, h);
1182  LWDEBUGF(2, "Looking for hole %d/%d of face %d among %d other faces", h+1, nholes, i, nfaces-i-1);
1183  for (j=i+1; j<nfaces; ++j) {
1184  const GEOSGeometry *f2er;
1185  Face* f2 = faces[j];
1186  if ( f2->parent ) continue; /* hole already assigned */
1187  f2er = GEOSGetExteriorRing(f2->geom);
1188  /* TODO: can be optimized as the ring would have the
1189  * same vertices, possibly in different order.
1190  * maybe comparing number of points could already be
1191  * useful.
1192  */
1193  if ( GEOSEquals(f2er, hole) ) {
1194  LWDEBUGF(2, "Hole %d/%d of face %d is face %d", h+1, nholes, i, j);
1195  f2->parent = f;
1196  break;
1197  }
1198  }
1199  }
1200  }
1201 }
1202 
1203 static GEOSGeometry*
1205 {
1206  GEOSGeometry **geoms = lwalloc(sizeof(GEOSGeometry*)*nfaces);
1207  GEOSGeometry *ret;
1208  unsigned int ngeoms = 0;
1209  int i;
1210 
1211  for (i=0; i<nfaces; ++i) {
1212  Face *f = faces[i];
1213  if ( countParens(f) % 2 ) continue; /* we skip odd parents geoms */
1214  geoms[ngeoms++] = GEOSGeom_clone(f->geom);
1215  }
1216 
1217  ret = GEOSGeom_createCollection(GEOS_MULTIPOLYGON, geoms, ngeoms);
1218  lwfree(geoms);
1219  return ret;
1220 }
1221 
1222 GEOSGeometry*
1223 LWGEOM_GEOS_buildArea(const GEOSGeometry* geom_in)
1224 {
1225  GEOSGeometry *tmp;
1226  GEOSGeometry *geos_result, *shp;
1227  GEOSGeometry const *vgeoms[1];
1228  uint32_t i, ngeoms;
1229  int srid = GEOSGetSRID(geom_in);
1230  Face ** geoms;
1231 
1232  vgeoms[0] = geom_in;
1233 #ifdef LWGEOM_PROFILE_BUILDAREA
1234  lwnotice("Polygonizing");
1235 #endif
1236  geos_result = GEOSPolygonize(vgeoms, 1);
1237 
1238  LWDEBUGF(3, "GEOSpolygonize returned @ %p", geos_result);
1239 
1240  /* Null return from GEOSpolygonize (an exception) */
1241  if ( ! geos_result ) return 0;
1242 
1243  /*
1244  * We should now have a collection
1245  */
1246 #if PARANOIA_LEVEL > 0
1247  if ( GEOSGeometryTypeId(geos_result) != COLLECTIONTYPE )
1248  {
1249  GEOSGeom_destroy(geos_result);
1250  lwerror("Unexpected return from GEOSpolygonize");
1251  return 0;
1252  }
1253 #endif
1254 
1255  ngeoms = GEOSGetNumGeometries(geos_result);
1256 #ifdef LWGEOM_PROFILE_BUILDAREA
1257  lwnotice("Num geometries from polygonizer: %d", ngeoms);
1258 #endif
1259 
1260 
1261  LWDEBUGF(3, "GEOSpolygonize: ngeoms in polygonize output: %d", ngeoms);
1262  LWDEBUGF(3, "GEOSpolygonize: polygonized:%s",
1263  lwgeom_to_ewkt(GEOS2LWGEOM(geos_result, 0)));
1264 
1265  /*
1266  * No geometries in collection, early out
1267  */
1268  if ( ngeoms == 0 )
1269  {
1270  GEOSSetSRID(geos_result, srid);
1271  return geos_result;
1272  }
1273 
1274  /*
1275  * Return first geometry if we only have one in collection,
1276  * to avoid the unnecessary Geometry clone below.
1277  */
1278  if ( ngeoms == 1 )
1279  {
1280  tmp = (GEOSGeometry *)GEOSGetGeometryN(geos_result, 0);
1281  if ( ! tmp )
1282  {
1283  GEOSGeom_destroy(geos_result);
1284  return 0; /* exception */
1285  }
1286  shp = GEOSGeom_clone(tmp);
1287  GEOSGeom_destroy(geos_result); /* only safe after the clone above */
1288  GEOSSetSRID(shp, srid);
1289  return shp;
1290  }
1291 
1292  LWDEBUGF(2, "Polygonize returned %d geoms", ngeoms);
1293 
1294  /*
1295  * Polygonizer returns a polygon for each face in the built topology.
1296  *
1297  * This means that for any face with holes we'll have other faces
1298  * representing each hole. We can imagine a parent-child relationship
1299  * between these faces.
1300  *
1301  * In order to maximize the number of visible rings in output we
1302  * only use those faces which have an even number of parents.
1303  *
1304  * Example:
1305  *
1306  * +---------------+
1307  * | L0 | L0 has no parents
1308  * | +---------+ |
1309  * | | L1 | | L1 is an hole of L0
1310  * | | +---+ | |
1311  * | | |L2 | | | L2 is an hole of L1 (which is an hole of L0)
1312  * | | | | | |
1313  * | | +---+ | |
1314  * | +---------+ |
1315  * | |
1316  * +---------------+
1317  *
1318  * See http://trac.osgeo.org/postgis/ticket/1806
1319  *
1320  */
1321 
1322 #ifdef LWGEOM_PROFILE_BUILDAREA
1323  lwnotice("Preparing face structures");
1324 #endif
1325 
1326  /* Prepare face structures for later analysis */
1327  geoms = lwalloc(sizeof(Face**)*ngeoms);
1328  for (i=0; i<ngeoms; ++i)
1329  geoms[i] = newFace(GEOSGetGeometryN(geos_result, i));
1330 
1331 #ifdef LWGEOM_PROFILE_BUILDAREA
1332  lwnotice("Finding face holes");
1333 #endif
1334 
1335  /* Find faces representing other faces holes */
1336  findFaceHoles(geoms, ngeoms);
1337 
1338 #ifdef LWGEOM_PROFILE_BUILDAREA
1339  lwnotice("Colletting even ancestor faces");
1340 #endif
1341 
1342  /* Build a MultiPolygon composed only by faces with an
1343  * even number of ancestors */
1344  tmp = collectFacesWithEvenAncestors(geoms, ngeoms);
1345 
1346 #ifdef LWGEOM_PROFILE_BUILDAREA
1347  lwnotice("Cleaning up");
1348 #endif
1349 
1350  /* Cleanup face structures */
1351  for (i=0; i<ngeoms; ++i) delFace(geoms[i]);
1352  lwfree(geoms);
1353 
1354  /* Faces referenced memory owned by geos_result.
1355  * It is safe to destroy geos_result after deleting them. */
1356  GEOSGeom_destroy(geos_result);
1357 
1358 #ifdef LWGEOM_PROFILE_BUILDAREA
1359  lwnotice("Self-unioning");
1360 #endif
1361 
1362  /* Run a single overlay operation to dissolve shared edges */
1363  shp = GEOSUnionCascaded(tmp);
1364  if ( ! shp )
1365  {
1366  GEOSGeom_destroy(tmp);
1367  return 0; /* exception */
1368  }
1369 
1370 #ifdef LWGEOM_PROFILE_BUILDAREA
1371  lwnotice("Final cleanup");
1372 #endif
1373 
1374  GEOSGeom_destroy(tmp);
1375 
1376  GEOSSetSRID(shp, srid);
1377 
1378  return shp;
1379 }
1380 
1381 LWGEOM*
1383 {
1384  GEOSGeometry* geos_in;
1385  GEOSGeometry* geos_out;
1386  LWGEOM* geom_out;
1387  int SRID = (int)(geom->srid);
1388  int is3d = FLAGS_GET_Z(geom->flags);
1389 
1390  /* Can't build an area from an empty! */
1391  if ( lwgeom_is_empty(geom) )
1392  {
1393  return (LWGEOM*)lwpoly_construct_empty(SRID, is3d, 0);
1394  }
1395 
1396  LWDEBUG(3, "buildarea called");
1397 
1398  LWDEBUGF(3, "ST_BuildArea got geom @ %p", geom);
1399 
1400  initGEOS(lwnotice, lwgeom_geos_error);
1401 
1402  geos_in = LWGEOM2GEOS(geom, 0);
1403 
1404  if ( 0 == geos_in ) /* exception thrown at construction */
1405  {
1406  lwerror("First argument geometry could not be converted to GEOS: %s", lwgeom_geos_errmsg);
1407  return NULL;
1408  }
1409  geos_out = LWGEOM_GEOS_buildArea(geos_in);
1410  GEOSGeom_destroy(geos_in);
1411 
1412  if ( ! geos_out ) /* exception thrown.. */
1413  {
1414  lwerror("LWGEOM_GEOS_buildArea: %s", lwgeom_geos_errmsg);
1415  return NULL;
1416  }
1417 
1418  /* If no geometries are in result collection, return NULL */
1419  if ( GEOSGetNumGeometries(geos_out) == 0 )
1420  {
1421  GEOSGeom_destroy(geos_out);
1422  return NULL;
1423  }
1424 
1425  geom_out = GEOS2LWGEOM(geos_out, is3d);
1426  GEOSGeom_destroy(geos_out);
1427 
1428 #if PARANOIA_LEVEL > 0
1429  if ( geom_out == NULL )
1430  {
1431  lwerror("serialization error");
1432  return NULL;
1433  }
1434 
1435 #endif
1436 
1437  return geom_out;
1438 }
1439 
1440 int
1442 {
1443  GEOSGeometry* geos_in;
1444  int simple;
1445 
1446  /* Empty is always simple */
1447  if ( lwgeom_is_empty(geom) )
1448  {
1449  return 1;
1450  }
1451 
1452  initGEOS(lwnotice, lwgeom_geos_error);
1453 
1454  geos_in = LWGEOM2GEOS(geom, 0);
1455  if ( 0 == geos_in ) /* exception thrown at construction */
1456  {
1457  lwerror("First argument geometry could not be converted to GEOS: %s", lwgeom_geos_errmsg);
1458  return -1;
1459  }
1460  simple = GEOSisSimple(geos_in);
1461  GEOSGeom_destroy(geos_in);
1462 
1463  if ( simple == 2 ) /* exception thrown */
1464  {
1465  lwerror("lwgeom_is_simple: %s", lwgeom_geos_errmsg);
1466  return -1;
1467  }
1468 
1469  return simple ? 1 : 0;
1470 }
1471 
1472 /* ------------ end of BuildArea stuff ---------------------------------------------------------------------} */
1473 
1474 LWGEOM*
1475 lwgeom_geos_noop(const LWGEOM* geom_in)
1476 {
1477  GEOSGeometry *geosgeom;
1478  LWGEOM* geom_out;
1479 
1480  int is3d = FLAGS_GET_Z(geom_in->flags);
1481 
1482  initGEOS(lwnotice, lwgeom_geos_error);
1483  geosgeom = LWGEOM2GEOS(geom_in, 0);
1484  if ( ! geosgeom ) {
1485  lwerror("Geometry could not be converted to GEOS: %s",
1487  return NULL;
1488  }
1489  geom_out = GEOS2LWGEOM(geosgeom, is3d);
1490  GEOSGeom_destroy(geosgeom);
1491  if ( ! geom_out ) {
1492  lwerror("GEOS Geometry could not be converted to LWGEOM: %s",
1494  }
1495  return geom_out;
1496 
1497 }
1498 
1499 LWGEOM*
1500 lwgeom_snap(const LWGEOM* geom1, const LWGEOM* geom2, double tolerance)
1501 {
1502  int srid, is3d;
1503  GEOSGeometry *g1, *g2, *g3;
1504  LWGEOM* out;
1505 
1506  srid = geom1->srid;
1507  error_if_srid_mismatch(srid, (int)(geom2->srid));
1508 
1509  is3d = (FLAGS_GET_Z(geom1->flags) || FLAGS_GET_Z(geom2->flags)) ;
1510 
1511  initGEOS(lwnotice, lwgeom_geos_error);
1512 
1513  g1 = (GEOSGeometry *)LWGEOM2GEOS(geom1, 0);
1514  if ( 0 == g1 ) /* exception thrown at construction */
1515  {
1516  lwerror("First argument geometry could not be converted to GEOS: %s", lwgeom_geos_errmsg);
1517  return NULL;
1518  }
1519 
1520  g2 = (GEOSGeometry *)LWGEOM2GEOS(geom2, 0);
1521  if ( 0 == g2 ) /* exception thrown at construction */
1522  {
1523  lwerror("Second argument geometry could not be converted to GEOS: %s", lwgeom_geos_errmsg);
1524  GEOSGeom_destroy(g1);
1525  return NULL;
1526  }
1527 
1528  g3 = GEOSSnap(g1, g2, tolerance);
1529  if (g3 == NULL)
1530  {
1531  GEOSGeom_destroy(g1);
1532  GEOSGeom_destroy(g2);
1533  lwerror("GEOSSnap: %s", lwgeom_geos_errmsg);
1534  return NULL;
1535  }
1536 
1537  GEOSGeom_destroy(g1);
1538  GEOSGeom_destroy(g2);
1539 
1540  GEOSSetSRID(g3, srid);
1541  out = GEOS2LWGEOM(g3, is3d);
1542  if (out == NULL)
1543  {
1544  GEOSGeom_destroy(g3);
1545  lwerror("GEOSSnap() threw an error (result LWGEOM geometry formation)!");
1546  return NULL;
1547  }
1548  GEOSGeom_destroy(g3);
1549 
1550  return out;
1551 }
1552 
1553 LWGEOM*
1554 lwgeom_sharedpaths(const LWGEOM* geom1, const LWGEOM* geom2)
1555 {
1556  GEOSGeometry *g1, *g2, *g3;
1557  LWGEOM *out;
1558  int is3d, srid;
1559 
1560  srid = geom1->srid;
1561  error_if_srid_mismatch(srid, (int)(geom2->srid));
1562 
1563  is3d = (FLAGS_GET_Z(geom1->flags) || FLAGS_GET_Z(geom2->flags)) ;
1564 
1565  initGEOS(lwnotice, lwgeom_geos_error);
1566 
1567  g1 = (GEOSGeometry *)LWGEOM2GEOS(geom1, 0);
1568  if ( 0 == g1 ) /* exception thrown at construction */
1569  {
1570  lwerror("First argument geometry could not be converted to GEOS: %s", lwgeom_geos_errmsg);
1571  return NULL;
1572  }
1573 
1574  g2 = (GEOSGeometry *)LWGEOM2GEOS(geom2, 0);
1575  if ( 0 == g2 ) /* exception thrown at construction */
1576  {
1577  lwerror("Second argument geometry could not be converted to GEOS: %s", lwgeom_geos_errmsg);
1578  GEOSGeom_destroy(g1);
1579  return NULL;
1580  }
1581 
1582  g3 = GEOSSharedPaths(g1,g2);
1583 
1584  GEOSGeom_destroy(g1);
1585  GEOSGeom_destroy(g2);
1586 
1587  if (g3 == NULL)
1588  {
1589  lwerror("GEOSSharedPaths: %s", lwgeom_geos_errmsg);
1590  return NULL;
1591  }
1592 
1593  GEOSSetSRID(g3, srid);
1594  out = GEOS2LWGEOM(g3, is3d);
1595  GEOSGeom_destroy(g3);
1596 
1597  if (out == NULL)
1598  {
1599  lwerror("GEOS2LWGEOM threw an error");
1600  return NULL;
1601  }
1602 
1603  return out;
1604 }
1605 
1606 LWGEOM*
1607 lwgeom_offsetcurve(const LWLINE *lwline, double size, int quadsegs, int joinStyle, double mitreLimit)
1608 {
1609  GEOSGeometry *g1, *g3;
1610  LWGEOM *lwgeom_result;
1611  LWGEOM *lwgeom_in = lwline_as_lwgeom(lwline);
1612 
1613  initGEOS(lwnotice, lwgeom_geos_error);
1614 
1615  g1 = (GEOSGeometry *)LWGEOM2GEOS(lwgeom_in, 0);
1616  if ( ! g1 )
1617  {
1618  lwerror("lwgeom_offsetcurve: Geometry could not be converted to GEOS: %s", lwgeom_geos_errmsg);
1619  return NULL;
1620  }
1621 
1622  g3 = GEOSOffsetCurve(g1, size, quadsegs, joinStyle, mitreLimit);
1623 
1624  /* Don't need input geometry anymore */
1625  GEOSGeom_destroy(g1);
1626 
1627  if (g3 == NULL)
1628  {
1629  lwerror("GEOSOffsetCurve: %s", lwgeom_geos_errmsg);
1630  return NULL;
1631  }
1632 
1633  LWDEBUGF(3, "result: %s", GEOSGeomToWKT(g3));
1634 
1635  GEOSSetSRID(g3, lwgeom_get_srid(lwgeom_in));
1636  lwgeom_result = GEOS2LWGEOM(g3, lwgeom_has_z(lwgeom_in));
1637  GEOSGeom_destroy(g3);
1638 
1639  if (lwgeom_result == NULL)
1640  {
1641  lwerror("lwgeom_offsetcurve: GEOS2LWGEOM returned null");
1642  return NULL;
1643  }
1644 
1645  return lwgeom_result;
1646 }
1647 
1648 
1649 static void shuffle(void *array, size_t n, size_t size) {
1650  char tmp[size];
1651  char *arr = array;
1652  size_t stride = size;
1653 
1654  if (n > 1) {
1655  size_t i;
1656  for (i = 0; i < n - 1; ++i) {
1657  size_t rnd = (size_t) rand();
1658  size_t j = i + rnd / (RAND_MAX / (n - i) + 1);
1659 
1660  memcpy(tmp, arr + j * stride, size);
1661  memcpy(arr + j * stride, arr + i * stride, size);
1662  memcpy(arr + i * stride, tmp, size);
1663  }
1664  }
1665 }
1666 
1667 LWMPOINT*
1668 lwpoly_to_points(const LWPOLY *lwpoly, int npoints)
1669 {
1670  double area, bbox_area, bbox_width, bbox_height;
1671  GBOX bbox;
1672  const LWGEOM *lwgeom = (LWGEOM*)lwpoly;
1673  int sample_npoints, sample_sqrt, sample_width, sample_height;
1674  double sample_cell_size;
1675  int i, j;
1676  int iterations = 0;
1677  int npoints_generated = 0;
1678  int npoints_tested = 0;
1679  GEOSGeometry *g;
1680  const GEOSPreparedGeometry *gprep;
1681  GEOSGeometry *gpt;
1682  GEOSCoordSequence *gseq;
1683  LWMPOINT *mpt;
1684  int srid = lwgeom_get_srid(lwgeom);
1685  int done = 0;
1686  int *cells;
1687 
1688  if (lwgeom_get_type(lwgeom) != POLYGONTYPE)
1689  {
1690  lwerror("%s: only polygons supported", __func__);
1691  return NULL;
1692  }
1693 
1694  if (npoints == 0 || lwgeom_is_empty(lwgeom))
1695  {
1696  return NULL;
1697  // return lwmpoint_construct_empty(lwgeom_get_srid(poly), lwgeom_has_z(poly), lwgeom_has_m(poly));
1698  }
1699 
1700  if (!lwpoly->bbox)
1701  {
1702  lwgeom_calculate_gbox(lwgeom, &bbox);
1703  }
1704  else
1705  {
1706  bbox = *(lwpoly->bbox);
1707  }
1708  area = lwpoly_area(lwpoly);
1709  bbox_width = bbox.xmax - bbox.xmin;
1710  bbox_height = bbox.ymax - bbox.ymin;
1711  bbox_area = bbox_width * bbox_height;
1712 
1713  if (area == 0.0 || bbox_area == 0.0)
1714  {
1715  lwerror("%s: zero area input polygon, TBD", __func__);
1716  return NULL;
1717  }
1718 
1719  /* Gross up our test set a bit to increase odds of getting */
1720  /* coverage in one pass */
1721  sample_npoints = npoints * bbox_area / area;
1722 
1723  /* We're going to generate points using a sample grid */
1724  /* as described http://lin-ear-th-inking.blogspot.ca/2010/05/more-random-points-in-jts.html */
1725  /* to try and get a more uniform "random" set of points */
1726  /* So we have to figure out how to stick a grid into our box */
1727  sample_sqrt = lround(sqrt(sample_npoints));
1728  if (sample_sqrt == 0)
1729  sample_sqrt = 1;
1730 
1731  /* Calculate the grids we're going to randomize within */
1732  if (bbox_width > bbox_height)
1733  {
1734  sample_width = sample_sqrt;
1735  sample_height = ceil((double)sample_npoints / (double)sample_width);
1736  sample_cell_size = bbox_width / sample_width;
1737  }
1738  else
1739  {
1740  sample_height = sample_sqrt;
1741  sample_width = ceil((double)sample_npoints / (double)sample_height);
1742  sample_cell_size = bbox_height / sample_height;
1743  }
1744 
1745  /* Prepare the polygon for fast true/false testing */
1746  initGEOS(lwnotice, lwgeom_geos_error);
1747  g = (GEOSGeometry *)LWGEOM2GEOS(lwgeom, 0);
1748  if (!g)
1749  {
1750  lwerror("%s: Geometry could not be converted to GEOS: %s", __func__, lwgeom_geos_errmsg);
1751  return NULL;
1752  }
1753  gprep = GEOSPrepare(g);
1754 
1755  /* Get an empty multi-point ready to return */
1756  mpt = lwmpoint_construct_empty(srid, 0, 0);
1757 
1758  /* Init random number generator */
1759  srand(time(NULL));
1760 
1761  /* Now we fill in an array of cells, and then shuffle that array, */
1762  /* so we can visit the cells in random order to avoid visual ugliness */
1763  /* caused by visiting them sequentially */
1764  cells = lwalloc(2*sizeof(int)*sample_height*sample_width);
1765  for (i = 0; i < sample_width; i++)
1766  {
1767  for (j = 0; j < sample_height; j++)
1768  {
1769  cells[2*(i*sample_height+j)] = i;
1770  cells[2*(i*sample_height+j)+1] = j;
1771  }
1772  }
1773  shuffle(cells, sample_height*sample_width, 2*sizeof(int));
1774 
1775  /* Start testing points */
1776  while (npoints_generated < npoints)
1777  {
1778  iterations++;
1779  for (i = 0; i < sample_width*sample_height; i++)
1780  {
1781  int contains = 0;
1782  double y = bbox.ymin + cells[2*i] * sample_cell_size;
1783  double x = bbox.xmin + cells[2*i+1] * sample_cell_size;
1784  x += rand() * sample_cell_size / RAND_MAX;
1785  y += rand() * sample_cell_size / RAND_MAX;
1786  if (x >= bbox.xmax || y >= bbox.ymax)
1787  continue;
1788 
1789  gseq = GEOSCoordSeq_create(1, 2);
1790  GEOSCoordSeq_setX(gseq, 0, x);
1791  GEOSCoordSeq_setY(gseq, 0, y);
1792  gpt = GEOSGeom_createPoint(gseq);
1793 
1794  contains = GEOSPreparedIntersects(gprep, gpt);
1795 
1796  GEOSGeom_destroy(gpt);
1797 
1798  if (contains == 2)
1799  {
1800  GEOSPreparedGeom_destroy(gprep);
1801  GEOSGeom_destroy(g);
1802  lwerror("%s: GEOS exception on PreparedContains: %s", __func__, lwgeom_geos_errmsg);
1803  return NULL;
1804  }
1805  if (contains)
1806  {
1807  npoints_generated++;
1808  mpt = lwmpoint_add_lwpoint(mpt, lwpoint_make2d(srid, x, y));
1809  if (npoints_generated == npoints)
1810  {
1811  done = 1;
1812  break;
1813  }
1814  }
1815 
1816  /* Short-circuit check for ctrl-c occasionally */
1817  npoints_tested++;
1818  if (npoints_tested % 10000 == 0)
1819  {
1820  LW_ON_INTERRUPT(GEOSPreparedGeom_destroy(gprep); GEOSGeom_destroy(g); return NULL);
1821  }
1822 
1823  if (done) break;
1824  }
1825  if (done || iterations > 100) break;
1826  }
1827 
1828  GEOSPreparedGeom_destroy(gprep);
1829  GEOSGeom_destroy(g);
1830  lwfree(cells);
1831 
1832  return mpt;
1833 }
1834 
1835 
1836 /*
1837 * Allocate points to sub-geometries by area, then call lwgeom_poly_to_points
1838 * and bundle up final result in a single multipoint.
1839 */
1840 LWMPOINT*
1841 lwmpoly_to_points(const LWMPOLY *lwmpoly, int npoints)
1842 {
1843  const LWGEOM *lwgeom = (LWGEOM*)lwmpoly;
1844  double area;
1845  int i;
1846  LWMPOINT *mpt = NULL;
1847 
1848  if (lwgeom_get_type(lwgeom) != MULTIPOLYGONTYPE)
1849  {
1850  lwerror("%s: only multipolygons supported", __func__);
1851  return NULL;
1852  }
1853  if (npoints == 0 || lwgeom_is_empty(lwgeom))
1854  {
1855  return NULL;
1856  }
1857 
1858  area = lwgeom_area(lwgeom);
1859 
1860  for (i = 0; i < lwmpoly->ngeoms; i++)
1861  {
1862  double sub_area = lwpoly_area(lwmpoly->geoms[i]);
1863  int sub_npoints = lround(npoints * sub_area / area);
1864  if(sub_npoints > 0)
1865  {
1866  LWMPOINT *sub_mpt = lwpoly_to_points(lwmpoly->geoms[i], sub_npoints);
1867  if (!mpt)
1868  {
1869  mpt = sub_mpt;
1870  }
1871  else
1872  {
1873  int j;
1874  for (j = 0; j < sub_mpt->ngeoms; j++)
1875  {
1876  mpt = lwmpoint_add_lwpoint(mpt, sub_mpt->geoms[j]);
1877  }
1878  /* Just free the shell, leave the underlying lwpoints alone, as they
1879  are now owed by the returning multipoint */
1880  lwgeom_release((LWGEOM*)sub_mpt);
1881  }
1882  }
1883  }
1884 
1885  return mpt;
1886 }
1887 
1888 
1889 LWMPOINT*
1890 lwgeom_to_points(const LWGEOM *lwgeom, int npoints)
1891 {
1892  switch(lwgeom_get_type(lwgeom))
1893  {
1894  case MULTIPOLYGONTYPE:
1895  return lwmpoly_to_points((LWMPOLY*)lwgeom, npoints);
1896  case POLYGONTYPE:
1897  return lwpoly_to_points((LWPOLY*)lwgeom, npoints);
1898  default:
1899  lwerror("%s: unsupported geometry type '%s'", __func__, lwtype_name(lwgeom_get_type(lwgeom)));
1900  return NULL;
1901  }
1902 }
1903 
1904 
1905 
1906 
1907 LWTIN *lwtin_from_geos(const GEOSGeometry *geom, int want3d) {
1908  int type = GEOSGeomTypeId(geom);
1909  int hasZ;
1910  int SRID = GEOSGetSRID(geom);
1911 
1912  /* GEOS's 0 is equivalent to our unknown as for SRID values */
1913  if ( SRID == 0 ) SRID = SRID_UNKNOWN;
1914 
1915  if ( want3d ) {
1916  hasZ = GEOSHasZ(geom);
1917  if ( ! hasZ ) {
1918  LWDEBUG(3, "Geometry has no Z, won't provide one");
1919  want3d = 0;
1920  }
1921  }
1922 
1923  switch (type) {
1924  LWTRIANGLE **geoms;
1925  uint32_t i, ngeoms;
1926  case GEOS_GEOMETRYCOLLECTION:
1927  LWDEBUG(4, "lwgeom_from_geometry: it's a Collection or Multi");
1928 
1929  ngeoms = GEOSGetNumGeometries(geom);
1930  geoms = NULL;
1931  if ( ngeoms ) {
1932  geoms = lwalloc(ngeoms * sizeof *geoms);
1933  if (!geoms) {
1934  lwerror("lwtin_from_geos: can't allocate geoms");
1935  return NULL;
1936  }
1937  for (i=0; i<ngeoms; i++) {
1938  const GEOSGeometry *poly, *ring;
1939  const GEOSCoordSequence *cs;
1940  POINTARRAY *pa;
1941 
1942  poly = GEOSGetGeometryN(geom, i);
1943  ring = GEOSGetExteriorRing(poly);
1944  cs = GEOSGeom_getCoordSeq(ring);
1945  pa = ptarray_from_GEOSCoordSeq(cs, want3d);
1946 
1947  geoms[i] = lwtriangle_construct(SRID, NULL, pa);
1948  }
1949  }
1950  return (LWTIN *)lwcollection_construct(TINTYPE, SRID, NULL, ngeoms, (LWGEOM **)geoms);
1951  case GEOS_POLYGON:
1952  case GEOS_MULTIPOINT:
1953  case GEOS_MULTILINESTRING:
1954  case GEOS_MULTIPOLYGON:
1955  case GEOS_LINESTRING:
1956  case GEOS_LINEARRING:
1957  case GEOS_POINT:
1958  lwerror("lwtin_from_geos: invalid geometry type for tin: %d", type);
1959  break;
1960 
1961  default:
1962  lwerror("GEOS2LWGEOM: unknown geometry type: %d", type);
1963  return NULL;
1964  }
1965 
1966  /* shouldn't get here */
1967  return NULL;
1968 }
1969 /*
1970  * output = 1 for edges, 2 for TIN, 0 for polygons
1971  */
1972 LWGEOM* lwgeom_delaunay_triangulation(const LWGEOM *lwgeom_in, double tolerance, int output) {
1973 #if POSTGIS_GEOS_VERSION < 34
1974  lwerror("lwgeom_delaunay_triangulation: GEOS 3.4 or higher required");
1975  return NULL;
1976 #else
1977  GEOSGeometry *g1, *g3;
1978  LWGEOM *lwgeom_result;
1979 
1980  if (output < 0 || output > 2) {
1981  lwerror("lwgeom_delaunay_triangulation: invalid output type specified %d", output);
1982  return NULL;
1983  }
1984 
1985  initGEOS(lwnotice, lwgeom_geos_error);
1986 
1987  g1 = (GEOSGeometry *)LWGEOM2GEOS(lwgeom_in, 0);
1988  if ( ! g1 )
1989  {
1990  lwerror("lwgeom_delaunay_triangulation: Geometry could not be converted to GEOS: %s", lwgeom_geos_errmsg);
1991  return NULL;
1992  }
1993 
1994  /* if output != 1 we want polys */
1995  g3 = GEOSDelaunayTriangulation(g1, tolerance, output == 1);
1996 
1997  /* Don't need input geometry anymore */
1998  GEOSGeom_destroy(g1);
1999 
2000  if (g3 == NULL)
2001  {
2002  lwerror("GEOSDelaunayTriangulation: %s", lwgeom_geos_errmsg);
2003  return NULL;
2004  }
2005 
2006  /* LWDEBUGF(3, "result: %s", GEOSGeomToWKT(g3)); */
2007 
2008  GEOSSetSRID(g3, lwgeom_get_srid(lwgeom_in));
2009 
2010  if (output == 2) {
2011  lwgeom_result = (LWGEOM *)lwtin_from_geos(g3, lwgeom_has_z(lwgeom_in));
2012  } else {
2013  lwgeom_result = GEOS2LWGEOM(g3, lwgeom_has_z(lwgeom_in));
2014  }
2015 
2016  GEOSGeom_destroy(g3);
2017 
2018  if (lwgeom_result == NULL) {
2019  if (output != 2) {
2020  lwerror("lwgeom_delaunay_triangulation: GEOS2LWGEOM returned null");
2021  } else {
2022  lwerror("lwgeom_delaunay_triangulation: lwtin_from_geos returned null");
2023  }
2024  return NULL;
2025  }
2026 
2027  return lwgeom_result;
2028 
2029 #endif /* POSTGIS_GEOS_VERSION < 34 */
2030 }
2031 
2032 static
2033 GEOSCoordSequence* lwgeom_get_geos_coordseq_2d(const LWGEOM* g, uint32_t num_points)
2034 {
2035  uint32_t i = 0;
2036  uint8_t num_dims = 2;
2037  LWPOINTITERATOR* it;
2038  GEOSCoordSequence* coords;
2039  POINT4D tmp;
2040 
2041  coords = GEOSCoordSeq_create(num_points, num_dims);
2042  if (!coords)
2043  return NULL;
2044 
2045  it = lwpointiterator_create(g);
2046  while(lwpointiterator_next(it, &tmp))
2047  {
2048  if(i >= num_points)
2049  {
2050  lwerror("Incorrect num_points provided to lwgeom_get_geos_coordseq_2d");
2051  GEOSCoordSeq_destroy(coords);
2053  return NULL;
2054  }
2055 
2056  if(!GEOSCoordSeq_setX(coords, i, tmp.x) || !GEOSCoordSeq_setY(coords, i, tmp.y))
2057  {
2058  GEOSCoordSeq_destroy(coords);
2060  return NULL;
2061  }
2062  i++;
2063  }
2065 
2066  return coords;
2067 }
2068 
2069 LWGEOM* lwgeom_voronoi_diagram(const LWGEOM* g, const GBOX* env, double tolerance, int output_edges) {
2070 #if POSTGIS_GEOS_VERSION < 35
2071  lwerror("lwgeom_voronoi_diagram: GEOS 3.5 or higher required");
2072  return NULL;
2073 #else
2074  uint32_t num_points = lwgeom_count_vertices(g);
2075  LWGEOM *lwgeom_result;
2076  char is_3d = LW_FALSE;
2077  int srid = lwgeom_get_srid(g);
2078  GEOSCoordSequence* coords;
2079  GEOSGeometry* geos_geom;
2080  GEOSGeometry* geos_env = NULL;
2081  GEOSGeometry* geos_result;
2082 
2083  if (num_points < 2)
2084  {
2086  return lwcollection_as_lwgeom(empty);
2087  }
2088 
2089  initGEOS(lwnotice, lwgeom_geos_error);
2090 
2091  /* Instead of using the standard LWGEOM2GEOS transformer, we read the vertices of the
2092  * LWGEOM directly and put them into a single GEOS CoordinateSeq that can be used to
2093  * define a LineString. This allows us to process geometry types that may not be
2094  * supported by GEOS, and reduces the memory requirements in cases of many geometries
2095  * with few points (such as LWMPOINT).
2096  */
2097  coords = lwgeom_get_geos_coordseq_2d(g, num_points);
2098  if (!coords)
2099  return NULL;
2100 
2101  geos_geom = GEOSGeom_createLineString(coords);
2102  if (!geos_geom)
2103  {
2104  GEOSCoordSeq_destroy(coords);
2105  return NULL;
2106  }
2107 
2108  if (env)
2109  geos_env = GBOX2GEOS(env);
2110 
2111  geos_result = GEOSVoronoiDiagram(geos_geom, geos_env, tolerance, output_edges);
2112 
2113  GEOSGeom_destroy(geos_geom);
2114  if (env)
2115  GEOSGeom_destroy(geos_env);
2116 
2117  if (!geos_result)
2118  {
2119  lwerror("GEOSVoronoiDiagram: %s", lwgeom_geos_errmsg);
2120  return NULL;
2121  }
2122 
2123  lwgeom_result = GEOS2LWGEOM(geos_result, is_3d);
2124  GEOSGeom_destroy(geos_result);
2125 
2126  lwgeom_set_srid(lwgeom_result, srid);
2127 
2128  return lwgeom_result;
2129 #endif /* POSTGIS_GEOS_VERSION < 35 */
2130 }
2131 
void ptarray_set_point4d(POINTARRAY *pa, int n, const POINT4D *p4d)
Definition: lwgeom_api.c:549
double x
Definition: liblwgeom.h:351
#define LINETYPE
Definition: liblwgeom.h:85
LWGEOM * lwgeom_stroke(const LWGEOM *geom, uint32_t perQuad)
Definition: lwstroke.c:462
LWGEOM * lwgeom_voronoi_diagram(const LWGEOM *g, const GBOX *env, double tolerance, int output_edges)
Take vertices of a geometry and build the Voronoi diagram.
double z
Definition: liblwgeom.h:333
#define POSTGIS_GEOS_VERSION
Definition: sqldefines.h:10
double y
Definition: liblwgeom.h:333
POINTARRAY * ptarray_construct(char hasz, char hasm, uint32_t npoints)
Construct an empty pointarray, allocating storage and setting the npoints, but not filling in any inf...
Definition: ptarray.c:62
int lwgeom_is_simple(const LWGEOM *geom)
LWLINE * lwline_construct_empty(int srid, char hasz, char hasm)
Definition: lwline.c:64
GEOSGeometry * make_geos_segment(double x1, double y1, double x2, double y2)
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:212
static void shuffle(void *array, size_t n, size_t size)
LWCOLLECTION * lwcollection_construct(uint8_t type, int srid, GBOX *bbox, uint32_t ngeoms, LWGEOM **geoms)
Definition: lwcollection.c:43
double x
Definition: liblwgeom.h:333
static void delFace(Face *f)
void lwnotice(const char *fmt,...)
Write a notice out to the notice handler.
Definition: lwutil.c:89
void lwfree(void *mem)
Definition: lwutil.c:242
LWMPOINT * lwpoly_to_points(const LWPOLY *lwpoly, int npoints)
int npoints
Definition: liblwgeom.h:370
uint32_t lwgeom_get_type(const LWGEOM *geom)
Return LWTYPE number.
Definition: lwgeom.c:842
static int compare_by_envarea(const void *g1, const void *g2)
#define POLYGONTYPE
Definition: liblwgeom.h:86
LWPOINT * lwpoint_make2d(int srid, double x, double y)
Definition: lwpoint.c:145
Datum area(PG_FUNCTION_ARGS)
char * lwgeom_to_ewkt(const LWGEOM *lwgeom)
Return an alloced string.
Definition: lwgeom.c:482
uint8_t flags
Definition: liblwgeom.h:396
double xmax
Definition: liblwgeom.h:292
void lwgeom_free(LWGEOM *geom)
Definition: lwgeom.c:1063
void ptarray_free(POINTARRAY *pa)
Definition: ptarray.c:330
#define MULTIPOINTTYPE
Definition: liblwgeom.h:87
def fmt
Definition: pixval.py:92
LWGEOM * lwgeom_symdifference(const LWGEOM *geom1, const LWGEOM *geom2)
Datum contains(PG_FUNCTION_ARGS)
char lwgeom_geos_errmsg[LWGEOM_GEOS_ERRMSG_MAXSIZE]
Datum centroid(PG_FUNCTION_ARGS)
GEOSGeometry * GBOX2GEOS(const GBOX *box)
#define LWDEBUG(level, msg)
Definition: lwgeom_log.h:83
int32_t lwgeom_get_srid(const LWGEOM *geom)
Return SRID number.
Definition: lwgeom.c:835
#define LW_ON_INTERRUPT(x)
struct Face_t Face
LWTIN * lwtin_from_geos(const GEOSGeometry *geom, int want3d)
static void findFaceHoles(Face **faces, int nfaces)
LWGEOM * lwgeom_clone_deep(const LWGEOM *lwgeom)
Deep clone an LWGEOM, everything is copied.
Definition: lwgeom.c:446
GBOX * bbox
Definition: liblwgeom.h:452
void error_if_srid_mismatch(int srid1, int srid2)
Definition: lwutil.c:369
LWPOINT * lwpoint_construct_empty(int srid, char hasz, char hasm)
Definition: lwpoint.c:133
LWGEOM * lwgeom_centroid(const LWGEOM *geom)
POINTARRAY * point
Definition: liblwgeom.h:410
int ptarray_is_closed_2d(const POINTARRAY *pa)
Definition: ptarray.c:694
int lwgeom_has_z(const LWGEOM *geom)
Return LW_TRUE if geometry has Z ordinates.
Definition: lwgeom.c:849
LWPOINTITERATOR * lwpointiterator_create(const LWGEOM *g)
Create a new LWPOINTITERATOR over supplied LWGEOM*.
Definition: lwiterator.c:244
int32_t srid
Definition: liblwgeom.h:398
GEOSGeometry * make_geos_point(double x, double y)
LWGEOM * lwgeom_intersection(const LWGEOM *geom1, const LWGEOM *geom2)
static GEOSCoordSequence * lwgeom_get_geos_coordseq_2d(const LWGEOM *g, uint32_t num_points)
int lwgeom_calculate_gbox(const LWGEOM *lwgeom, GBOX *gbox)
Calculate bounding box of a geometry, automatically taking into account whether it is cartesian or ge...
Definition: lwgeom.c:665
double x
Definition: liblwgeom.h:327
GEOSCoordSeq ptarray_to_GEOSCoordSeq(const POINTARRAY *)
LWGEOM * lwgeom_linemerge(const LWGEOM *geom1)
LWGEOM * lwgeom_buildarea(const LWGEOM *geom)
Take a geometry and return an areal geometry (Polygon or MultiPolygon).
POINTARRAY * ptarray_addPoint(const POINTARRAY *pa, uint8_t *p, size_t pdims, uint32_t where)
Add a point in a pointarray.
Definition: ptarray.c:509
double ymin
Definition: liblwgeom.h:293
void lwpointiterator_destroy(LWPOINTITERATOR *s)
Free all memory associated with the iterator.
Definition: lwiterator.c:269
const char * lwtype_name(uint8_t type)
Return the type name string associated with a type number (e.g.
Definition: lwutil.c:216
LWGEOM * geom
LWGEOM * lwline_as_lwgeom(const LWLINE *obj)
Definition: lwgeom.c:262
double xmin
Definition: liblwgeom.h:291
void lwgeom_geos_error(const char *fmt,...)
#define LW_FALSE
Definition: liblwgeom.h:76
LWGEOM * lwgeom_sharedpaths(const LWGEOM *geom1, const LWGEOM *geom2)
const POINT2D * getPoint2d_cp(const POINTARRAY *pa, int n)
Returns a POINT2D pointer into the POINTARRAY serialized_ptlist, suitable for reading from...
Definition: lwgeom_api.c:485
uint8_t flags
Definition: liblwgeom.h:368
LWPOLY * lwpoly_construct(int srid, GBOX *bbox, uint32_t nrings, POINTARRAY **points)
Definition: lwpoly.c:43
LWPOLY ** geoms
Definition: liblwgeom.h:495
LWLINE * lwline_construct(int srid, GBOX *bbox, POINTARRAY *points)
Definition: lwline.c:42
LWGEOM * lwgeom_geos_noop(const LWGEOM *geom_in)
Convert an LWGEOM to a GEOS Geometry and convert back – for debug only.
LWGEOM ** geoms
Definition: liblwgeom.h:508
static GEOSGeometry * collectFacesWithEvenAncestors(Face **faces, int nfaces)
#define SRID_UNKNOWN
Unknown SRID value.
Definition: liblwgeom.h:187
#define TINTYPE
Definition: liblwgeom.h:98
const GEOSGeometry * geom
LWGEOM * lwgeom_union(const LWGEOM *geom1, const LWGEOM *geom2)
static unsigned int countParens(const Face *f)
uint8_t * getPoint_internal(const POINTARRAY *pa, int n)
Definition: ptarray.c:1706
LWGEOM * lwgeom_normalize(const LWGEOM *geom1)
LWGEOM * lwgeom_difference(const LWGEOM *geom1, const LWGEOM *geom2)
POINTARRAY ** rings
Definition: liblwgeom.h:456
LWPOINT ** geoms
Definition: liblwgeom.h:469
double lwpoly_area(const LWPOLY *poly)
Find the area of the outer ring - sum (area of inner rings).
Definition: lwpoly.c:484
LWGEOM * lwgeom_delaunay_triangulation(const LWGEOM *lwgeom_in, double tolerance, int output)
Take vertices of a geometry and build a delaunay triangulation on them.
LWGEOM * lwgeom_clip_by_rect(const LWGEOM *geom1, double x0, double y0, double x1, double y1)
int nrings
Definition: liblwgeom.h:454
LWMPOINT * lwmpoint_construct_empty(int srid, char hasz, char hasm)
Definition: lwmpoint.c:39
double ymax
Definition: liblwgeom.h:294
int lwgeom_has_arc(const LWGEOM *geom)
Definition: lwstroke.c:55
static Face * newFace(const GEOSGeometry *g)
double y
Definition: liblwgeom.h:327
#define LWGEOM_GEOS_ERRMSG_MAXSIZE
LWMPOINT * lwmpoly_to_points(const LWMPOLY *lwmpoly, int npoints)
#define FLAGS_GET_Z(flags)
Macros for manipulating the &#39;flags&#39; byte.
Definition: liblwgeom.h:139
GEOSGeometry * LWGEOM2GEOS(const LWGEOM *lwgeom, int autofix)
double z
Definition: liblwgeom.h:351
LWTRIANGLE * lwtriangle_construct(int srid, GBOX *bbox, POINTARRAY *points)
Definition: lwtriangle.c:40
int ngeoms
Definition: liblwgeom.h:493
GEOSGeometry * env
#define MULTIPOLYGONTYPE
Definition: liblwgeom.h:89
LWGEOM * lwgeom_unaryunion(const LWGEOM *geom1)
const POINT3DZ * getPoint3dz_cp(const POINTARRAY *pa, int n)
Returns a POINT3DZ pointer into the POINTARRAY serialized_ptlist, suitable for reading from...
Definition: lwgeom_api.c:499
LWGEOM * lwgeom_snap(const LWGEOM *geom1, const LWGEOM *geom2, double tolerance)
Snap vertices and segments of a geometry to another using a given tolerance.
double lwgeom_area(const LWGEOM *geom)
Definition: lwgeom.c:1592
LWMPOINT * lwmpoint_add_lwpoint(LWMPOINT *mobj, const LWPOINT *obj)
Definition: lwmpoint.c:45
static GEOSGeometry * ptarray_to_GEOSLinearRing(const POINTARRAY *pa, int autofix)
#define POINTTYPE
LWTYPE numbers, used internally by PostGIS.
Definition: liblwgeom.h:84
LWPOLY * lwpoly_construct_empty(int srid, char hasz, char hasm)
Definition: lwpoly.c:139
LWGEOM * GEOS2LWGEOM(const GEOSGeometry *geom, char want3d)
void lwgeom_release(LWGEOM *lwgeom)
Free the containing LWGEOM and the associated BOX.
Definition: lwgeom.c:385
uint8_t type
Definition: liblwgeom.h:395
type
Definition: ovdump.py:41
struct Face_t * parent
void free(void *)
LWGEOM * lwpoint_as_lwgeom(const LWPOINT *obj)
Definition: lwgeom.c:267
const char * lwgeom_geos_version()
Return GEOS version string (not to be freed)
void * malloc(YYSIZE_T)
void lwgeom_set_srid(LWGEOM *geom, int srid)
Set the SRID on an LWGEOM For collections, only the parent gets an SRID, all the children get SRID_UN...
LWGEOM * lwgeom_offsetcurve(const LWLINE *lwline, double size, int quadsegs, int joinStyle, double mitreLimit)
LWPOINT * lwpoint_construct(int srid, GBOX *bbox, POINTARRAY *point)
Definition: lwpoint.c:111
void * lwalloc(size_t size)
Definition: lwutil.c:227
int lwgeom_is_empty(const LWGEOM *geom)
Return true or false depending on whether a geometry is an "empty" geometry (no vertices members) ...
Definition: lwgeom.c:1310
int lwgeom_count_vertices(const LWGEOM *geom)
Count the total number of vertices in any LWGEOM.
Definition: lwgeom.c:1153
double y
Definition: liblwgeom.h:351
#define MULTILINETYPE
Definition: liblwgeom.h:88
LWCOLLECTION * lwcollection_construct_empty(uint8_t type, int srid, char hasz, char hasm)
Definition: lwcollection.c:94
POINTARRAY * ptarray_from_GEOSCoordSeq(const GEOSCoordSequence *cs, char want3d)
int ngeoms
Definition: liblwgeom.h:467
#define LWDEBUGF(level, msg,...)
Definition: lwgeom_log.h:88
#define FLAGS_NDIMS(flags)
Definition: liblwgeom.h:151
int lwgeom_has_m(const LWGEOM *geom)
Return LW_TRUE if geometry has M ordinates.
Definition: lwgeom.c:856
void lwerror(const char *fmt,...)
Write a notice out to the error handler.
Definition: lwutil.c:102
LWMPOINT * lwgeom_to_points(const LWGEOM *lwgeom, int npoints)
#define COLLECTIONTYPE
Definition: liblwgeom.h:90
This library is the generic geometry handling section of PostGIS.
GEOSGeometry * LWGEOM_GEOS_buildArea(const GEOSGeometry *geom_in)
POINTARRAY * points
Definition: liblwgeom.h:421
LWGEOM * lwcollection_as_lwgeom(const LWCOLLECTION *obj)
Definition: lwgeom.c:232