PostGIS  2.1.10dev-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  * Copyright 2011-2012 Sandro Santilli <strk@keybit.net>
7  *
8  * This is free software; you can redistribute and/or modify it under
9  * the terms of the GNU General Public Licence. See the COPYING file.
10  *
11  **********************************************************************/
12 
13 #include "lwgeom_geos.h"
14 #include "liblwgeom.h"
15 #include "liblwgeom_internal.h"
16 #include "lwgeom_log.h"
17 
18 #include <stdlib.h>
19 
20 LWTIN *lwtin_from_geos(const GEOSGeometry *geom, int want3d);
21 
22 #undef LWGEOM_PROFILE_BUILDAREA
23 
24 #define LWGEOM_GEOS_ERRMSG_MAXSIZE 256
26 
27 extern void
28 lwgeom_geos_error(const char *fmt, ...)
29 {
30  va_list ap;
31 
32  va_start(ap, fmt);
33 
34  /* Call the supplied function */
36  {
38  }
39 
40  va_end(ap);
41 }
42 
43 
44 /*
45 ** GEOS <==> PostGIS conversion functions
46 **
47 ** Default conversion creates a GEOS point array, then iterates through the
48 ** PostGIS points, setting each value in the GEOS array one at a time.
49 **
50 */
51 
52 /* Return a POINTARRAY from a GEOSCoordSeq */
53 POINTARRAY *
54 ptarray_from_GEOSCoordSeq(const GEOSCoordSequence *cs, char want3d)
55 {
56  uint32_t dims=2;
57  uint32_t size, i;
58  POINTARRAY *pa;
59  POINT4D point;
60 
61  LWDEBUG(2, "ptarray_fromGEOSCoordSeq called");
62 
63  if ( ! GEOSCoordSeq_getSize(cs, &size) )
64  lwerror("Exception thrown");
65 
66  LWDEBUGF(4, " GEOSCoordSeq size: %d", size);
67 
68  if ( want3d )
69  {
70  if ( ! GEOSCoordSeq_getDimensions(cs, &dims) )
71  lwerror("Exception thrown");
72 
73  LWDEBUGF(4, " GEOSCoordSeq dimensions: %d", dims);
74 
75  /* forget higher dimensions (if any) */
76  if ( dims > 3 ) dims = 3;
77  }
78 
79  LWDEBUGF(4, " output dimensions: %d", dims);
80 
81  pa = ptarray_construct((dims==3), 0, size);
82 
83  for (i=0; i<size; i++)
84  {
85  GEOSCoordSeq_getX(cs, i, &(point.x));
86  GEOSCoordSeq_getY(cs, i, &(point.y));
87  if ( dims >= 3 ) GEOSCoordSeq_getZ(cs, i, &(point.z));
88  ptarray_set_point4d(pa,i,&point);
89  }
90 
91  return pa;
92 }
93 
94 /* Return an LWGEOM from a Geometry */
95 LWGEOM *
96 GEOS2LWGEOM(const GEOSGeometry *geom, char want3d)
97 {
98  int type = GEOSGeomTypeId(geom) ;
99  int hasZ;
100  int SRID = GEOSGetSRID(geom);
101 
102  /* GEOS's 0 is equivalent to our unknown as for SRID values */
103  if ( SRID == 0 ) SRID = SRID_UNKNOWN;
104 
105  if ( want3d )
106  {
107  hasZ = GEOSHasZ(geom);
108  if ( ! hasZ )
109  {
110  LWDEBUG(3, "Geometry has no Z, won't provide one");
111 
112  want3d = 0;
113  }
114  }
115 
116 /*
117  if ( GEOSisEmpty(geom) )
118  {
119  return (LWGEOM*)lwcollection_construct_empty(COLLECTIONTYPE, SRID, want3d, 0);
120  }
121 */
122 
123  switch (type)
124  {
125  const GEOSCoordSequence *cs;
126  POINTARRAY *pa, **ppaa;
127  const GEOSGeometry *g;
128  LWGEOM **geoms;
129  uint32_t i, ngeoms;
130 
131  case GEOS_POINT:
132  LWDEBUG(4, "lwgeom_from_geometry: it's a Point");
133  cs = GEOSGeom_getCoordSeq(geom);
134  if ( GEOSisEmpty(geom) )
135  return (LWGEOM*)lwpoint_construct_empty(SRID, want3d, 0);
136  pa = ptarray_from_GEOSCoordSeq(cs, want3d);
137  return (LWGEOM *)lwpoint_construct(SRID, NULL, pa);
138 
139  case GEOS_LINESTRING:
140  case GEOS_LINEARRING:
141  LWDEBUG(4, "lwgeom_from_geometry: it's a LineString or LinearRing");
142  if ( GEOSisEmpty(geom) )
143  return (LWGEOM*)lwline_construct_empty(SRID, want3d, 0);
144 
145  cs = GEOSGeom_getCoordSeq(geom);
146  pa = ptarray_from_GEOSCoordSeq(cs, want3d);
147  return (LWGEOM *)lwline_construct(SRID, NULL, pa);
148 
149  case GEOS_POLYGON:
150  LWDEBUG(4, "lwgeom_from_geometry: it's a Polygon");
151  if ( GEOSisEmpty(geom) )
152  return (LWGEOM*)lwpoly_construct_empty(SRID, want3d, 0);
153  ngeoms = GEOSGetNumInteriorRings(geom);
154  ppaa = lwalloc(sizeof(POINTARRAY *)*(ngeoms+1));
155  g = GEOSGetExteriorRing(geom);
156  cs = GEOSGeom_getCoordSeq(g);
157  ppaa[0] = ptarray_from_GEOSCoordSeq(cs, want3d);
158  for (i=0; i<ngeoms; i++)
159  {
160  g = GEOSGetInteriorRingN(geom, i);
161  cs = GEOSGeom_getCoordSeq(g);
162  ppaa[i+1] = ptarray_from_GEOSCoordSeq(cs,
163  want3d);
164  }
165  return (LWGEOM *)lwpoly_construct(SRID, NULL,
166  ngeoms+1, ppaa);
167 
168  case GEOS_MULTIPOINT:
169  case GEOS_MULTILINESTRING:
170  case GEOS_MULTIPOLYGON:
171  case GEOS_GEOMETRYCOLLECTION:
172  LWDEBUG(4, "lwgeom_from_geometry: it's a Collection or Multi");
173 
174  ngeoms = GEOSGetNumGeometries(geom);
175  geoms = NULL;
176  if ( ngeoms )
177  {
178  geoms = lwalloc(sizeof(LWGEOM *)*ngeoms);
179  for (i=0; i<ngeoms; i++)
180  {
181  g = GEOSGetGeometryN(geom, i);
182  geoms[i] = GEOS2LWGEOM(g, want3d);
183  }
184  }
185  return (LWGEOM *)lwcollection_construct(type,
186  SRID, NULL, ngeoms, geoms);
187 
188  default:
189  lwerror("GEOS2LWGEOM: unknown geometry type: %d", type);
190  return NULL;
191 
192  }
193 
194 }
195 
196 
197 
198 GEOSCoordSeq ptarray_to_GEOSCoordSeq(const POINTARRAY *);
199 
200 
201 GEOSCoordSeq
203 {
204  uint32_t dims = 2;
205  uint32_t i;
206  const POINT3DZ *p3d;
207  const POINT2D *p2d;
208  GEOSCoordSeq sq;
209 
210  if ( FLAGS_GET_Z(pa->flags) )
211  dims = 3;
212 
213  if ( ! (sq = GEOSCoordSeq_create(pa->npoints, dims)) )
214  lwerror("Error creating GEOS Coordinate Sequence");
215 
216  for ( i=0; i < pa->npoints; i++ )
217  {
218  if ( dims == 3 )
219  {
220  p3d = getPoint3dz_cp(pa, i);
221  p2d = (const POINT2D *)p3d;
222  LWDEBUGF(4, "Point: %g,%g,%g", p3d->x, p3d->y, p3d->z);
223  }
224  else
225  {
226  p2d = getPoint2d_cp(pa, i);
227  LWDEBUGF(4, "Point: %g,%g", p2d->x, p2d->y);
228  }
229 
230 #if POSTGIS_GEOS_VERSION < 33
231  /* Make sure we don't pass any infinite values down into GEOS */
232  /* GEOS 3.3+ is supposed to handle this stuff OK */
233  if ( isinf(p2d->x) || isinf(p2d->y) || (dims == 3 && isinf(p3d->z)) )
234  lwerror("Infinite coordinate value found in geometry.");
235  if ( isnan(p2d->x) || isnan(p2d->y) || (dims == 3 && isnan(p3d->z)) )
236  lwerror("NaN coordinate value found in geometry.");
237 #endif
238 
239  GEOSCoordSeq_setX(sq, i, p2d->x);
240  GEOSCoordSeq_setY(sq, i, p2d->y);
241 
242  if ( dims == 3 )
243  GEOSCoordSeq_setZ(sq, i, p3d->z);
244  }
245  return sq;
246 }
247 
248 
249 
250 
251 GEOSGeometry *
252 LWGEOM2GEOS(const LWGEOM *lwgeom)
253 {
254  GEOSCoordSeq sq;
255  GEOSGeom g, shell;
256  GEOSGeom *geoms = NULL;
257  /*
258  LWGEOM *tmp;
259  */
260  uint32_t ngeoms, i, j;
261  int geostype;
262 #if LWDEBUG_LEVEL >= 4
263  char *wkt;
264 #endif
265 
266  LWDEBUGF(4, "LWGEOM2GEOS got a %s", lwtype_name(lwgeom->type));
267 
268  if (lwgeom_has_arc(lwgeom))
269  {
270  LWDEBUG(3, "LWGEOM2GEOS: arced geometry found.");
271 
272  lwerror("Exception in LWGEOM2GEOS: curved geometry not supported.");
273  return NULL;
274  }
275 
276  switch (lwgeom->type)
277  {
278  LWPOINT *lwp = NULL;
279  LWPOLY *lwpoly = NULL;
280  LWLINE *lwl = NULL;
281  LWCOLLECTION *lwc = NULL;
282 #if POSTGIS_GEOS_VERSION < 33
283  POINTARRAY *pa = NULL;
284 #endif
285 
286  case POINTTYPE:
287  lwp = (LWPOINT *)lwgeom;
288 
289  if ( lwgeom_is_empty(lwgeom) )
290  {
291 #if POSTGIS_GEOS_VERSION < 33
292  pa = ptarray_construct_empty(lwgeom_has_z(lwgeom), lwgeom_has_m(lwgeom), 2);
293  sq = ptarray_to_GEOSCoordSeq(pa);
294  shell = GEOSGeom_createLinearRing(sq);
295  g = GEOSGeom_createPolygon(shell, NULL, 0);
296 #else
297  g = GEOSGeom_createEmptyPolygon();
298 #endif
299  }
300  else
301  {
302  sq = ptarray_to_GEOSCoordSeq(lwp->point);
303  g = GEOSGeom_createPoint(sq);
304  }
305  if ( ! g )
306  {
307  /* lwnotice("Exception in LWGEOM2GEOS"); */
308  return NULL;
309  }
310  break;
311  case LINETYPE:
312  lwl = (LWLINE *)lwgeom;
313  if ( lwl->points->npoints == 1 ) {
314  /* Duplicate point, to make geos-friendly */
315  lwl->points = ptarray_addPoint(lwl->points,
316  getPoint_internal(lwl->points, 0),
317  FLAGS_NDIMS(lwl->points->flags),
318  lwl->points->npoints);
319  }
320  sq = ptarray_to_GEOSCoordSeq(lwl->points);
321  g = GEOSGeom_createLineString(sq);
322  if ( ! g )
323  {
324  /* lwnotice("Exception in LWGEOM2GEOS"); */
325  return NULL;
326  }
327  break;
328 
329  case POLYGONTYPE:
330  lwpoly = (LWPOLY *)lwgeom;
331  if ( lwgeom_is_empty(lwgeom) )
332  {
333 #if POSTGIS_GEOS_VERSION < 33
334  POINTARRAY *pa = ptarray_construct_empty(lwgeom_has_z(lwgeom), lwgeom_has_m(lwgeom), 2);
335  sq = ptarray_to_GEOSCoordSeq(pa);
336  shell = GEOSGeom_createLinearRing(sq);
337  g = GEOSGeom_createPolygon(shell, NULL, 0);
338 #else
339  g = GEOSGeom_createEmptyPolygon();
340 #endif
341  }
342  else
343  {
344  sq = ptarray_to_GEOSCoordSeq(lwpoly->rings[0]);
345  /* TODO: check ring for being closed and fix if not */
346  shell = GEOSGeom_createLinearRing(sq);
347  if ( ! shell ) return NULL;
348  /*lwerror("LWGEOM2GEOS: exception during polygon shell conversion"); */
349  ngeoms = lwpoly->nrings-1;
350  if ( ngeoms > 0 )
351  geoms = malloc(sizeof(GEOSGeom)*ngeoms);
352 
353  for (i=1; i<lwpoly->nrings; ++i)
354  {
355  sq = ptarray_to_GEOSCoordSeq(lwpoly->rings[i]);
356  geoms[i-1] = GEOSGeom_createLinearRing(sq);
357  if ( ! geoms[i-1] )
358  {
359  --i;
360  while (i) GEOSGeom_destroy(geoms[--i]);
361  free(geoms);
362  GEOSGeom_destroy(shell);
363  return NULL;
364  }
365  /*lwerror("LWGEOM2GEOS: exception during polygon hole conversion"); */
366  }
367  g = GEOSGeom_createPolygon(shell, geoms, ngeoms);
368  if (geoms) free(geoms);
369  }
370  if ( ! g ) return NULL;
371  break;
372  case MULTIPOINTTYPE:
373  case MULTILINETYPE:
374  case MULTIPOLYGONTYPE:
375  case COLLECTIONTYPE:
376  if ( lwgeom->type == MULTIPOINTTYPE )
377  geostype = GEOS_MULTIPOINT;
378  else if ( lwgeom->type == MULTILINETYPE )
379  geostype = GEOS_MULTILINESTRING;
380  else if ( lwgeom->type == MULTIPOLYGONTYPE )
381  geostype = GEOS_MULTIPOLYGON;
382  else
383  geostype = GEOS_GEOMETRYCOLLECTION;
384 
385  lwc = (LWCOLLECTION *)lwgeom;
386 
387  ngeoms = lwc->ngeoms;
388 
389  if ( ngeoms > 0 )
390  geoms = malloc(sizeof(GEOSGeom)*ngeoms);
391 
392  for (i=0, j=0; i<ngeoms; ++i)
393  {
394  GEOSGeometry* g;
395 
396  if( lwgeom_is_empty(lwc->geoms[i]) )
397  continue;
398 
399  g = LWGEOM2GEOS(lwc->geoms[i]);
400  if ( ! g )
401  {
402  while (j) GEOSGeom_destroy(geoms[--j]);
403  free(geoms);
404  return NULL;
405  }
406  geoms[j++] = g;
407  }
408  g = GEOSGeom_createCollection(geostype, geoms, j);
409  if ( geoms ) free(geoms);
410  if ( ! g ) return NULL;
411  break;
412 
413  default:
414  lwerror("Unknown geometry type: %d - %s", lwgeom->type, lwtype_name(lwgeom->type));
415  return NULL;
416  }
417 
418  GEOSSetSRID(g, lwgeom->srid);
419 
420 #if LWDEBUG_LEVEL >= 4
421  wkt = GEOSGeomToWKT(g);
422  LWDEBUGF(4, "LWGEOM2GEOS: GEOSGeom: %s", wkt);
423  free(wkt);
424 #endif
425 
426  return g;
427 }
428 
429 const char*
431 {
432  const char *ver = GEOSversion();
433  return ver;
434 }
435 
436 LWGEOM *
438 {
439  LWGEOM *result ;
440  GEOSGeometry *g1;
441  int is3d ;
442  int srid ;
443 
444  srid = (int)(geom1->srid);
445  is3d = FLAGS_GET_Z(geom1->flags);
446 
447  initGEOS(lwnotice, lwgeom_geos_error);
448 
449  g1 = LWGEOM2GEOS(geom1);
450  if ( 0 == g1 ) /* exception thrown at construction */
451  {
452  lwerror("First argument geometry could not be converted to GEOS: %s", lwgeom_geos_errmsg);
453  return NULL ;
454  }
455 
456  if ( -1 == GEOSNormalize(g1) )
457  {
458  lwerror("Error in GEOSNormalize: %s", lwgeom_geos_errmsg);
459  return NULL; /* never get here */
460  }
461 
462  GEOSSetSRID(g1, srid); /* needed ? */
463  result = GEOS2LWGEOM(g1, is3d);
464  GEOSGeom_destroy(g1);
465 
466  if (result == NULL)
467  {
468  lwerror("Error performing intersection: GEOS2LWGEOM: %s",
470  return NULL ; /* never get here */
471  }
472 
473  return result ;
474 }
475 
476 LWGEOM *
477 lwgeom_intersection(const LWGEOM *geom1, const LWGEOM *geom2)
478 {
479  LWGEOM *result ;
480  GEOSGeometry *g1, *g2, *g3 ;
481  int is3d ;
482  int srid ;
483 
484  /* A.Intersection(Empty) == Empty */
485  if ( lwgeom_is_empty(geom2) )
486  return lwgeom_clone_deep(geom2);
487 
488  /* Empty.Intersection(A) == Empty */
489  if ( lwgeom_is_empty(geom1) )
490  return lwgeom_clone_deep(geom1);
491 
492  /* ensure srids are identical */
493  srid = (int)(geom1->srid);
494  error_if_srid_mismatch(srid, (int)(geom2->srid));
495 
496  is3d = (FLAGS_GET_Z(geom1->flags) || FLAGS_GET_Z(geom2->flags)) ;
497 
498  initGEOS(lwnotice, lwgeom_geos_error);
499 
500  LWDEBUG(3, "intersection() START");
501 
502  g1 = LWGEOM2GEOS(geom1);
503  if ( 0 == g1 ) /* exception thrown at construction */
504  {
505  lwerror("First argument geometry could not be converted to GEOS: %s", lwgeom_geos_errmsg);
506  return NULL ;
507  }
508 
509  g2 = LWGEOM2GEOS(geom2);
510  if ( 0 == g2 ) /* exception thrown at construction */
511  {
512  lwerror("Second argument geometry could not be converted to GEOS.");
513  GEOSGeom_destroy(g1);
514  return NULL ;
515  }
516 
517  LWDEBUG(3, " constructed geometrys - calling geos");
518  LWDEBUGF(3, " g1 = %s", GEOSGeomToWKT(g1));
519  LWDEBUGF(3, " g2 = %s", GEOSGeomToWKT(g2));
520  /*LWDEBUGF(3, "g2 is valid = %i",GEOSisvalid(g2)); */
521  /*LWDEBUGF(3, "g1 is valid = %i",GEOSisvalid(g1)); */
522 
523  g3 = GEOSIntersection(g1,g2);
524 
525  LWDEBUG(3, " intersection finished");
526 
527  if (g3 == NULL)
528  {
529  GEOSGeom_destroy(g1);
530  GEOSGeom_destroy(g2);
531  lwerror("Error performing intersection: %s",
533  return NULL; /* never get here */
534  }
535 
536  LWDEBUGF(3, "result: %s", GEOSGeomToWKT(g3) ) ;
537 
538  GEOSSetSRID(g3, srid);
539 
540  result = GEOS2LWGEOM(g3, is3d);
541 
542  if (result == NULL)
543  {
544  GEOSGeom_destroy(g1);
545  GEOSGeom_destroy(g2);
546  GEOSGeom_destroy(g3);
547  lwerror("Error performing intersection: GEOS2LWGEOM: %s",
549  return NULL ; /* never get here */
550  }
551 
552  GEOSGeom_destroy(g1);
553  GEOSGeom_destroy(g2);
554  GEOSGeom_destroy(g3);
555 
556  return result ;
557 }
558 
559 LWGEOM *
560 lwgeom_difference(const LWGEOM *geom1, const LWGEOM *geom2)
561 {
562  GEOSGeometry *g1, *g2, *g3;
563  LWGEOM *result;
564  int is3d;
565  int srid;
566 
567  /* A.Difference(Empty) == A */
568  if ( lwgeom_is_empty(geom2) )
569  return lwgeom_clone_deep(geom1);
570 
571  /* Empty.Intersection(A) == Empty */
572  if ( lwgeom_is_empty(geom1) )
573  return lwgeom_clone_deep(geom1);
574 
575  /* ensure srids are identical */
576  srid = (int)(geom1->srid);
577  error_if_srid_mismatch(srid, (int)(geom2->srid));
578 
579  is3d = (FLAGS_GET_Z(geom1->flags) || FLAGS_GET_Z(geom2->flags)) ;
580 
581  initGEOS(lwnotice, lwgeom_geos_error);
582 
583  g1 = LWGEOM2GEOS(geom1);
584  if ( 0 == g1 ) /* exception thrown at construction */
585  {
586  lwerror("First argument geometry could not be converted to GEOS: %s", lwgeom_geos_errmsg);
587  return NULL;
588  }
589 
590  g2 = LWGEOM2GEOS(geom2);
591  if ( 0 == g2 ) /* exception thrown at construction */
592  {
593  GEOSGeom_destroy(g1);
594  lwerror("Second argument geometry could not be converted to GEOS: %s", lwgeom_geos_errmsg);
595  return NULL;
596  }
597 
598  g3 = GEOSDifference(g1,g2);
599 
600  if (g3 == NULL)
601  {
602  GEOSGeom_destroy(g1);
603  GEOSGeom_destroy(g2);
604  lwerror("GEOSDifference: %s", lwgeom_geos_errmsg);
605  return NULL ; /* never get here */
606  }
607 
608  LWDEBUGF(3, "result: %s", GEOSGeomToWKT(g3) ) ;
609 
610  GEOSSetSRID(g3, srid);
611 
612  result = GEOS2LWGEOM(g3, is3d);
613 
614  if (result == NULL)
615  {
616  GEOSGeom_destroy(g1);
617  GEOSGeom_destroy(g2);
618  GEOSGeom_destroy(g3);
619  lwerror("Error performing difference: GEOS2LWGEOM: %s",
621  return NULL; /* never get here */
622  }
623 
624  GEOSGeom_destroy(g1);
625  GEOSGeom_destroy(g2);
626  GEOSGeom_destroy(g3);
627 
628  /* compressType(result); */
629 
630  return result;
631 }
632 
633 LWGEOM *
634 lwgeom_symdifference(const LWGEOM* geom1, const LWGEOM* geom2)
635 {
636  GEOSGeometry *g1, *g2, *g3;
637  LWGEOM *result;
638  int is3d;
639  int srid;
640 
641  /* A.SymDifference(Empty) == A */
642  if ( lwgeom_is_empty(geom2) )
643  return lwgeom_clone_deep(geom1);
644 
645  /* Empty.DymDifference(B) == B */
646  if ( lwgeom_is_empty(geom1) )
647  return lwgeom_clone_deep(geom2);
648 
649  /* ensure srids are identical */
650  srid = (int)(geom1->srid);
651  error_if_srid_mismatch(srid, (int)(geom2->srid));
652 
653  is3d = (FLAGS_GET_Z(geom1->flags) || FLAGS_GET_Z(geom2->flags)) ;
654 
655  initGEOS(lwnotice, lwgeom_geos_error);
656 
657  g1 = LWGEOM2GEOS(geom1);
658 
659  if ( 0 == g1 ) /* exception thrown at construction */
660  {
661  lwerror("First argument geometry could not be converted to GEOS: %s", lwgeom_geos_errmsg);
662  return NULL;
663  }
664 
665  g2 = LWGEOM2GEOS(geom2);
666 
667  if ( 0 == g2 ) /* exception thrown at construction */
668  {
669  lwerror("Second argument geometry could not be converted to GEOS: %s", lwgeom_geos_errmsg);
670  GEOSGeom_destroy(g1);
671  return NULL;
672  }
673 
674  g3 = GEOSSymDifference(g1,g2);
675 
676  if (g3 == NULL)
677  {
678  GEOSGeom_destroy(g1);
679  GEOSGeom_destroy(g2);
680  lwerror("GEOSSymDifference: %s", lwgeom_geos_errmsg);
681  return NULL; /*never get here */
682  }
683 
684  LWDEBUGF(3, "result: %s", GEOSGeomToWKT(g3));
685 
686  GEOSSetSRID(g3, srid);
687 
688  result = GEOS2LWGEOM(g3, is3d);
689 
690  if (result == NULL)
691  {
692  GEOSGeom_destroy(g1);
693  GEOSGeom_destroy(g2);
694  GEOSGeom_destroy(g3);
695  lwerror("GEOS symdifference() threw an error (result postgis geometry formation)!");
696  return NULL ; /*never get here */
697  }
698 
699  GEOSGeom_destroy(g1);
700  GEOSGeom_destroy(g2);
701  GEOSGeom_destroy(g3);
702 
703  return result;
704 }
705 
706 LWGEOM*
707 lwgeom_union(const LWGEOM *geom1, const LWGEOM *geom2)
708 {
709  int is3d;
710  int srid;
711  GEOSGeometry *g1, *g2, *g3;
712  LWGEOM *result;
713 
714  LWDEBUG(2, "in geomunion");
715 
716  /* A.Union(empty) == A */
717  if ( lwgeom_is_empty(geom1) )
718  return lwgeom_clone_deep(geom2);
719 
720  /* B.Union(empty) == B */
721  if ( lwgeom_is_empty(geom2) )
722  return lwgeom_clone_deep(geom1);
723 
724 
725  /* ensure srids are identical */
726  srid = (int)(geom1->srid);
727  error_if_srid_mismatch(srid, (int)(geom2->srid));
728 
729  is3d = (FLAGS_GET_Z(geom1->flags) || FLAGS_GET_Z(geom2->flags)) ;
730 
731  initGEOS(lwnotice, lwgeom_geos_error);
732 
733  g1 = LWGEOM2GEOS(geom1);
734 
735  if ( 0 == g1 ) /* exception thrown at construction */
736  {
737  lwerror("First argument geometry could not be converted to GEOS: %s", lwgeom_geos_errmsg);
738  return NULL;
739  }
740 
741  g2 = LWGEOM2GEOS(geom2);
742 
743  if ( 0 == g2 ) /* exception thrown at construction */
744  {
745  GEOSGeom_destroy(g1);
746  lwerror("Second argument geometry could not be converted to GEOS: %s", lwgeom_geos_errmsg);
747  return NULL;
748  }
749 
750  LWDEBUGF(3, "g1=%s", GEOSGeomToWKT(g1));
751  LWDEBUGF(3, "g2=%s", GEOSGeomToWKT(g2));
752 
753  g3 = GEOSUnion(g1,g2);
754 
755  LWDEBUGF(3, "g3=%s", GEOSGeomToWKT(g3));
756 
757  GEOSGeom_destroy(g1);
758  GEOSGeom_destroy(g2);
759 
760  if (g3 == NULL)
761  {
762  lwerror("GEOSUnion: %s", lwgeom_geos_errmsg);
763  return NULL; /* never get here */
764  }
765 
766 
767  GEOSSetSRID(g3, srid);
768 
769  result = GEOS2LWGEOM(g3, is3d);
770 
771  GEOSGeom_destroy(g3);
772 
773  if (result == NULL)
774  {
775  lwerror("Error performing union: GEOS2LWGEOM: %s",
777  return NULL; /*never get here */
778  }
779 
780  return result;
781 }
782 
783 /* ------------ BuildArea stuff ---------------------------------------------------------------------{ */
784 
785 typedef struct Face_t {
786  const GEOSGeometry* geom;
787  GEOSGeometry* env;
788  double envarea;
789  struct Face_t* parent; /* if this face is an hole of another one, or NULL */
790 } Face;
791 
792 static Face* newFace(const GEOSGeometry* g);
793 static void delFace(Face* f);
794 static unsigned int countParens(const Face* f);
795 static void findFaceHoles(Face** faces, int nfaces);
796 
797 static Face*
798 newFace(const GEOSGeometry* g)
799 {
800  Face* f = lwalloc(sizeof(Face));
801  f->geom = g;
802  f->env = GEOSEnvelope(f->geom);
803  GEOSArea(f->env, &f->envarea);
804  f->parent = NULL;
805  /* lwnotice("Built Face with area %g and %d holes", f->envarea, GEOSGetNumInteriorRings(f->geom)); */
806  return f;
807 }
808 
809 static unsigned int
810 countParens(const Face* f)
811 {
812  unsigned int pcount = 0;
813  while ( f->parent ) {
814  ++pcount;
815  f = f->parent;
816  }
817  return pcount;
818 }
819 
820 /* Destroy the face and release memory associated with it */
821 static void
823 {
824  GEOSGeom_destroy(f->env);
825  lwfree(f);
826 }
827 
828 
829 static int
830 compare_by_envarea(const void* g1, const void* g2)
831 {
832  Face* f1 = *(Face**)g1;
833  Face* f2 = *(Face**)g2;
834  double n1 = f1->envarea;
835  double n2 = f2->envarea;
836 
837  if ( n1 < n2 ) return 1;
838  if ( n1 > n2 ) return -1;
839  return 0;
840 }
841 
842 /* Find holes of each face */
843 static void
844 findFaceHoles(Face** faces, int nfaces)
845 {
846  int i, j, h;
847 
848  /* We sort by envelope area so that we know holes are only
849  * after their shells */
850  qsort(faces, nfaces, sizeof(Face*), compare_by_envarea);
851  for (i=0; i<nfaces; ++i) {
852  Face* f = faces[i];
853  int nholes = GEOSGetNumInteriorRings(f->geom);
854  LWDEBUGF(2, "Scanning face %d with env area %g and %d holes", i, f->envarea, nholes);
855  for (h=0; h<nholes; ++h) {
856  const GEOSGeometry *hole = GEOSGetInteriorRingN(f->geom, h);
857  LWDEBUGF(2, "Looking for hole %d/%d of face %d among %d other faces", h+1, nholes, i, nfaces-i-1);
858  for (j=i+1; j<nfaces; ++j) {
859  const GEOSGeometry *f2er;
860  Face* f2 = faces[j];
861  if ( f2->parent ) continue; /* hole already assigned */
862  f2er = GEOSGetExteriorRing(f2->geom);
863  /* TODO: can be optimized as the ring would have the
864  * same vertices, possibly in different order.
865  * maybe comparing number of points could already be
866  * useful.
867  */
868  if ( GEOSEquals(f2er, hole) ) {
869  LWDEBUGF(2, "Hole %d/%d of face %d is face %d", h+1, nholes, i, j);
870  f2->parent = f;
871  break;
872  }
873  }
874  }
875  }
876 }
877 
878 static GEOSGeometry*
880 {
881  GEOSGeometry **geoms = lwalloc(sizeof(GEOSGeometry*)*nfaces);
882  GEOSGeometry *ret;
883  unsigned int ngeoms = 0;
884  int i;
885 
886  for (i=0; i<nfaces; ++i) {
887  Face *f = faces[i];
888  if ( countParens(f) % 2 ) continue; /* we skip odd parents geoms */
889  geoms[ngeoms++] = GEOSGeom_clone(f->geom);
890  }
891 
892  ret = GEOSGeom_createCollection(GEOS_MULTIPOLYGON, geoms, ngeoms);
893  lwfree(geoms);
894  return ret;
895 }
896 
897 GEOSGeometry*
898 LWGEOM_GEOS_buildArea(const GEOSGeometry* geom_in)
899 {
900  GEOSGeometry *tmp;
901  GEOSGeometry *geos_result, *shp;
902  GEOSGeometry const *vgeoms[1];
903  uint32_t i, ngeoms;
904  int srid = GEOSGetSRID(geom_in);
905  Face ** geoms;
906 
907  vgeoms[0] = geom_in;
908 #ifdef LWGEOM_PROFILE_BUILDAREA
909  lwnotice("Polygonizing");
910 #endif
911  geos_result = GEOSPolygonize(vgeoms, 1);
912 
913  LWDEBUGF(3, "GEOSpolygonize returned @ %p", geos_result);
914 
915  /* Null return from GEOSpolygonize (an exception) */
916  if ( ! geos_result ) return 0;
917 
918  /*
919  * We should now have a collection
920  */
921 #if PARANOIA_LEVEL > 0
922  if ( GEOSGeometryTypeId(geos_result) != COLLECTIONTYPE )
923  {
924  GEOSGeom_destroy(geos_result);
925  lwerror("Unexpected return from GEOSpolygonize");
926  return 0;
927  }
928 #endif
929 
930  ngeoms = GEOSGetNumGeometries(geos_result);
931 #ifdef LWGEOM_PROFILE_BUILDAREA
932  lwnotice("Num geometries from polygonizer: %d", ngeoms);
933 #endif
934 
935 
936  LWDEBUGF(3, "GEOSpolygonize: ngeoms in polygonize output: %d", ngeoms);
937  LWDEBUGF(3, "GEOSpolygonize: polygonized:%s",
938  lwgeom_to_ewkt(GEOS2LWGEOM(geos_result, 0)));
939 
940  /*
941  * No geometries in collection, early out
942  */
943  if ( ngeoms == 0 )
944  {
945  GEOSSetSRID(geos_result, srid);
946  return geos_result;
947  }
948 
949  /*
950  * Return first geometry if we only have one in collection,
951  * to avoid the unnecessary Geometry clone below.
952  */
953  if ( ngeoms == 1 )
954  {
955  tmp = (GEOSGeometry *)GEOSGetGeometryN(geos_result, 0);
956  if ( ! tmp )
957  {
958  GEOSGeom_destroy(geos_result);
959  return 0; /* exception */
960  }
961  shp = GEOSGeom_clone(tmp);
962  GEOSGeom_destroy(geos_result); /* only safe after the clone above */
963  GEOSSetSRID(shp, srid);
964  return shp;
965  }
966 
967  LWDEBUGF(2, "Polygonize returned %d geoms", ngeoms);
968 
969  /*
970  * Polygonizer returns a polygon for each face in the built topology.
971  *
972  * This means that for any face with holes we'll have other faces
973  * representing each hole. We can imagine a parent-child relationship
974  * between these faces.
975  *
976  * In order to maximize the number of visible rings in output we
977  * only use those faces which have an even number of parents.
978  *
979  * Example:
980  *
981  * +---------------+
982  * | L0 | L0 has no parents
983  * | +---------+ |
984  * | | L1 | | L1 is an hole of L0
985  * | | +---+ | |
986  * | | |L2 | | | L2 is an hole of L1 (which is an hole of L0)
987  * | | | | | |
988  * | | +---+ | |
989  * | +---------+ |
990  * | |
991  * +---------------+
992  *
993  * See http://trac.osgeo.org/postgis/ticket/1806
994  *
995  */
996 
997 #ifdef LWGEOM_PROFILE_BUILDAREA
998  lwnotice("Preparing face structures");
999 #endif
1000 
1001  /* Prepare face structures for later analysis */
1002  geoms = lwalloc(sizeof(Face**)*ngeoms);
1003  for (i=0; i<ngeoms; ++i)
1004  geoms[i] = newFace(GEOSGetGeometryN(geos_result, i));
1005 
1006 #ifdef LWGEOM_PROFILE_BUILDAREA
1007  lwnotice("Finding face holes");
1008 #endif
1009 
1010  /* Find faces representing other faces holes */
1011  findFaceHoles(geoms, ngeoms);
1012 
1013 #ifdef LWGEOM_PROFILE_BUILDAREA
1014  lwnotice("Colletting even ancestor faces");
1015 #endif
1016 
1017  /* Build a MultiPolygon composed only by faces with an
1018  * even number of ancestors */
1019  tmp = collectFacesWithEvenAncestors(geoms, ngeoms);
1020 
1021 #ifdef LWGEOM_PROFILE_BUILDAREA
1022  lwnotice("Cleaning up");
1023 #endif
1024 
1025  /* Cleanup face structures */
1026  for (i=0; i<ngeoms; ++i) delFace(geoms[i]);
1027  lwfree(geoms);
1028 
1029  /* Faces referenced memory owned by geos_result.
1030  * It is safe to destroy geos_result after deleting them. */
1031  GEOSGeom_destroy(geos_result);
1032 
1033 #ifdef LWGEOM_PROFILE_BUILDAREA
1034  lwnotice("Self-unioning");
1035 #endif
1036 
1037  /* Run a single overlay operation to dissolve shared edges */
1038  shp = GEOSUnionCascaded(tmp);
1039  if ( ! shp )
1040  {
1041  GEOSGeom_destroy(tmp);
1042  return 0; /* exception */
1043  }
1044 
1045 #ifdef LWGEOM_PROFILE_BUILDAREA
1046  lwnotice("Final cleanup");
1047 #endif
1048 
1049  GEOSGeom_destroy(tmp);
1050 
1051  GEOSSetSRID(shp, srid);
1052 
1053  return shp;
1054 }
1055 
1056 LWGEOM*
1058 {
1059  GEOSGeometry* geos_in;
1060  GEOSGeometry* geos_out;
1061  LWGEOM* geom_out;
1062  int SRID = (int)(geom->srid);
1063  int is3d = FLAGS_GET_Z(geom->flags);
1064 
1065  /* Can't build an area from an empty! */
1066  if ( lwgeom_is_empty(geom) )
1067  {
1068  return (LWGEOM*)lwpoly_construct_empty(SRID, is3d, 0);
1069  }
1070 
1071  LWDEBUG(3, "buildarea called");
1072 
1073  LWDEBUGF(3, "ST_BuildArea got geom @ %p", geom);
1074 
1075  initGEOS(lwnotice, lwgeom_geos_error);
1076 
1077  geos_in = LWGEOM2GEOS(geom);
1078 
1079  if ( 0 == geos_in ) /* exception thrown at construction */
1080  {
1081  lwerror("First argument geometry could not be converted to GEOS: %s", lwgeom_geos_errmsg);
1082  return NULL;
1083  }
1084  geos_out = LWGEOM_GEOS_buildArea(geos_in);
1085  GEOSGeom_destroy(geos_in);
1086 
1087  if ( ! geos_out ) /* exception thrown.. */
1088  {
1089  lwerror("LWGEOM_GEOS_buildArea: %s", lwgeom_geos_errmsg);
1090  return NULL;
1091  }
1092 
1093  /* If no geometries are in result collection, return NULL */
1094  if ( GEOSGetNumGeometries(geos_out) == 0 )
1095  {
1096  GEOSGeom_destroy(geos_out);
1097  return NULL;
1098  }
1099 
1100  geom_out = GEOS2LWGEOM(geos_out, is3d);
1101  GEOSGeom_destroy(geos_out);
1102 
1103 #if PARANOIA_LEVEL > 0
1104  if ( geom_out == NULL )
1105  {
1106  lwerror("serialization error");
1107  return NULL;
1108  }
1109 
1110 #endif
1111 
1112  return geom_out;
1113 }
1114 
1115 /* ------------ end of BuildArea stuff ---------------------------------------------------------------------} */
1116 
1117 LWGEOM*
1118 lwgeom_geos_noop(const LWGEOM* geom_in)
1119 {
1120  GEOSGeometry *geosgeom;
1121  LWGEOM* geom_out;
1122 
1123  int is3d = FLAGS_GET_Z(geom_in->flags);
1124 
1125  initGEOS(lwnotice, lwgeom_geos_error);
1126  geosgeom = LWGEOM2GEOS(geom_in);
1127  if ( ! geosgeom ) {
1128  lwerror("Geometry could not be converted to GEOS: %s",
1130  return NULL;
1131  }
1132  geom_out = GEOS2LWGEOM(geosgeom, is3d);
1133  GEOSGeom_destroy(geosgeom);
1134  if ( ! geom_out ) {
1135  lwerror("GEOS Geometry could not be converted to LWGEOM: %s",
1137  }
1138  return geom_out;
1139 
1140 }
1141 
1142 LWGEOM*
1143 lwgeom_snap(const LWGEOM* geom1, const LWGEOM* geom2, double tolerance)
1144 {
1145 #if POSTGIS_GEOS_VERSION < 33
1146  lwerror("The GEOS version this lwgeom library "
1147  "was compiled against (%d) doesn't support "
1148  "'Snap' function (3.3.0+ required)",
1150  return NULL;
1151 #else /* POSTGIS_GEOS_VERSION >= 33 */
1152 
1153  int srid, is3d;
1154  GEOSGeometry *g1, *g2, *g3;
1155  LWGEOM* out;
1156 
1157  srid = geom1->srid;
1158  error_if_srid_mismatch(srid, (int)(geom2->srid));
1159 
1160  is3d = (FLAGS_GET_Z(geom1->flags) || FLAGS_GET_Z(geom2->flags)) ;
1161 
1162  initGEOS(lwnotice, lwgeom_geos_error);
1163 
1164  g1 = (GEOSGeometry *)LWGEOM2GEOS(geom1);
1165  if ( 0 == g1 ) /* exception thrown at construction */
1166  {
1167  lwerror("First argument geometry could not be converted to GEOS: %s", lwgeom_geos_errmsg);
1168  return NULL;
1169  }
1170 
1171  g2 = (GEOSGeometry *)LWGEOM2GEOS(geom2);
1172  if ( 0 == g2 ) /* exception thrown at construction */
1173  {
1174  lwerror("Second argument geometry could not be converted to GEOS: %s", lwgeom_geos_errmsg);
1175  GEOSGeom_destroy(g1);
1176  return NULL;
1177  }
1178 
1179  g3 = GEOSSnap(g1, g2, tolerance);
1180  if (g3 == NULL)
1181  {
1182  GEOSGeom_destroy(g1);
1183  GEOSGeom_destroy(g2);
1184  lwerror("GEOSSnap: %s", lwgeom_geos_errmsg);
1185  return NULL;
1186  }
1187 
1188  GEOSGeom_destroy(g1);
1189  GEOSGeom_destroy(g2);
1190 
1191  GEOSSetSRID(g3, srid);
1192  out = GEOS2LWGEOM(g3, is3d);
1193  if (out == NULL)
1194  {
1195  GEOSGeom_destroy(g3);
1196  lwerror("GEOSSnap() threw an error (result LWGEOM geometry formation)!");
1197  return NULL;
1198  }
1199  GEOSGeom_destroy(g3);
1200 
1201  return out;
1202 
1203 #endif /* POSTGIS_GEOS_VERSION >= 33 */
1204 }
1205 
1206 LWGEOM*
1207 lwgeom_sharedpaths(const LWGEOM* geom1, const LWGEOM* geom2)
1208 {
1209 #if POSTGIS_GEOS_VERSION < 33
1210  lwerror("The GEOS version this postgis binary "
1211  "was compiled against (%d) doesn't support "
1212  "'SharedPaths' function (3.3.0+ required)",
1214  return NULL;
1215 #else /* POSTGIS_GEOS_VERSION >= 33 */
1216  GEOSGeometry *g1, *g2, *g3;
1217  LWGEOM *out;
1218  int is3d, srid;
1219 
1220  srid = geom1->srid;
1221  error_if_srid_mismatch(srid, (int)(geom2->srid));
1222 
1223  is3d = (FLAGS_GET_Z(geom1->flags) || FLAGS_GET_Z(geom2->flags)) ;
1224 
1225  initGEOS(lwnotice, lwgeom_geos_error);
1226 
1227  g1 = (GEOSGeometry *)LWGEOM2GEOS(geom1);
1228  if ( 0 == g1 ) /* exception thrown at construction */
1229  {
1230  lwerror("First argument geometry could not be converted to GEOS: %s", lwgeom_geos_errmsg);
1231  return NULL;
1232  }
1233 
1234  g2 = (GEOSGeometry *)LWGEOM2GEOS(geom2);
1235  if ( 0 == g2 ) /* exception thrown at construction */
1236  {
1237  lwerror("Second argument geometry could not be converted to GEOS: %s", lwgeom_geos_errmsg);
1238  GEOSGeom_destroy(g1);
1239  return NULL;
1240  }
1241 
1242  g3 = GEOSSharedPaths(g1,g2);
1243 
1244  GEOSGeom_destroy(g1);
1245  GEOSGeom_destroy(g2);
1246 
1247  if (g3 == NULL)
1248  {
1249  lwerror("GEOSSharedPaths: %s", lwgeom_geos_errmsg);
1250  return NULL;
1251  }
1252 
1253  GEOSSetSRID(g3, srid);
1254  out = GEOS2LWGEOM(g3, is3d);
1255  GEOSGeom_destroy(g3);
1256 
1257  if (out == NULL)
1258  {
1259  lwerror("GEOS2LWGEOM threw an error");
1260  return NULL;
1261  }
1262 
1263  return out;
1264 #endif /* POSTGIS_GEOS_VERSION >= 33 */
1265 }
1266 
1267 LWGEOM*
1268 lwgeom_offsetcurve(const LWLINE *lwline, double size, int quadsegs, int joinStyle, double mitreLimit)
1269 {
1270 #if POSTGIS_GEOS_VERSION < 32
1271  lwerror("lwgeom_offsetcurve: GEOS 3.2 or higher required");
1272 #else
1273  GEOSGeometry *g1, *g3;
1274  LWGEOM *lwgeom_result;
1275  LWGEOM *lwgeom_in = lwline_as_lwgeom(lwline);
1276 
1277  initGEOS(lwnotice, lwgeom_geos_error);
1278 
1279  g1 = (GEOSGeometry *)LWGEOM2GEOS(lwgeom_in);
1280  if ( ! g1 )
1281  {
1282  lwerror("lwgeom_offsetcurve: Geometry could not be converted to GEOS: %s", lwgeom_geos_errmsg);
1283  return NULL;
1284  }
1285 
1286 #if POSTGIS_GEOS_VERSION < 33
1287  /* Size is always positive for GEOSSingleSidedBuffer, and a flag determines left/right */
1288  g3 = GEOSSingleSidedBuffer(g1, size < 0 ? -size : size,
1289  quadsegs, joinStyle, mitreLimit,
1290  size < 0 ? 0 : 1);
1291 #else
1292  g3 = GEOSOffsetCurve(g1, size, quadsegs, joinStyle, mitreLimit);
1293 #endif
1294  /* Don't need input geometry anymore */
1295  GEOSGeom_destroy(g1);
1296 
1297  if (g3 == NULL)
1298  {
1299  lwerror("GEOSOffsetCurve: %s", lwgeom_geos_errmsg);
1300  return NULL;
1301  }
1302 
1303  LWDEBUGF(3, "result: %s", GEOSGeomToWKT(g3));
1304 
1305  GEOSSetSRID(g3, lwgeom_get_srid(lwgeom_in));
1306  lwgeom_result = GEOS2LWGEOM(g3, lwgeom_has_z(lwgeom_in));
1307  GEOSGeom_destroy(g3);
1308 
1309  if (lwgeom_result == NULL)
1310  {
1311  lwerror("lwgeom_offsetcurve: GEOS2LWGEOM returned null");
1312  return NULL;
1313  }
1314 
1315  return lwgeom_result;
1316 
1317 #endif /* POSTGIS_GEOS_VERSION < 32 */
1318 }
1319 
1320 LWTIN *lwtin_from_geos(const GEOSGeometry *geom, int want3d) {
1321  int type = GEOSGeomTypeId(geom);
1322  int hasZ;
1323  int SRID = GEOSGetSRID(geom);
1324 
1325  /* GEOS's 0 is equivalent to our unknown as for SRID values */
1326  if ( SRID == 0 ) SRID = SRID_UNKNOWN;
1327 
1328  if ( want3d ) {
1329  hasZ = GEOSHasZ(geom);
1330  if ( ! hasZ ) {
1331  LWDEBUG(3, "Geometry has no Z, won't provide one");
1332  want3d = 0;
1333  }
1334  }
1335 
1336  switch (type) {
1337  LWTRIANGLE **geoms;
1338  uint32_t i, ngeoms;
1339  case GEOS_GEOMETRYCOLLECTION:
1340  LWDEBUG(4, "lwgeom_from_geometry: it's a Collection or Multi");
1341 
1342  ngeoms = GEOSGetNumGeometries(geom);
1343  geoms = NULL;
1344  if ( ngeoms ) {
1345  geoms = lwalloc(ngeoms * sizeof *geoms);
1346  if (!geoms) {
1347  lwerror("lwtin_from_geos: can't allocate geoms");
1348  return NULL;
1349  }
1350  for (i=0; i<ngeoms; i++) {
1351  const GEOSGeometry *poly, *ring;
1352  const GEOSCoordSequence *cs;
1353  POINTARRAY *pa;
1354 
1355  poly = GEOSGetGeometryN(geom, i);
1356  ring = GEOSGetExteriorRing(poly);
1357  cs = GEOSGeom_getCoordSeq(ring);
1358  pa = ptarray_from_GEOSCoordSeq(cs, want3d);
1359 
1360  geoms[i] = lwtriangle_construct(SRID, NULL, pa);
1361  }
1362  }
1363  return (LWTIN *)lwcollection_construct(TINTYPE, SRID, NULL, ngeoms, (LWGEOM **)geoms);
1364  case GEOS_POLYGON:
1365  case GEOS_MULTIPOINT:
1366  case GEOS_MULTILINESTRING:
1367  case GEOS_MULTIPOLYGON:
1368  case GEOS_LINESTRING:
1369  case GEOS_LINEARRING:
1370  case GEOS_POINT:
1371  lwerror("lwtin_from_geos: invalid geometry type for tin: %d", type);
1372  break;
1373 
1374  default:
1375  lwerror("GEOS2LWGEOM: unknown geometry type: %d", type);
1376  return NULL;
1377  }
1378 
1379  /* shouldn't get here */
1380  return NULL;
1381 }
1382 /*
1383  * output = 1 for edges, 2 for TIN, 0 for polygons
1384  */
1385 LWGEOM* lwgeom_delaunay_triangulation(const LWGEOM *lwgeom_in, double tolerance, int output) {
1386 #if POSTGIS_GEOS_VERSION < 34
1387  lwerror("lwgeom_delaunay_triangulation: GEOS 3.4 or higher required");
1388  return NULL;
1389 #else
1390  GEOSGeometry *g1, *g3;
1391  LWGEOM *lwgeom_result;
1392 
1393  if (output < 0 || output > 2) {
1394  lwerror("lwgeom_delaunay_triangulation: invalid output type specified %d", output);
1395  return NULL;
1396  }
1397 
1398  initGEOS(lwnotice, lwgeom_geos_error);
1399 
1400  g1 = (GEOSGeometry *)LWGEOM2GEOS(lwgeom_in);
1401  if ( ! g1 )
1402  {
1403  lwerror("lwgeom_delaunay_triangulation: Geometry could not be converted to GEOS: %s", lwgeom_geos_errmsg);
1404  return NULL;
1405  }
1406 
1407  /* if output != 1 we want polys */
1408  g3 = GEOSDelaunayTriangulation(g1, tolerance, output == 1);
1409 
1410  /* Don't need input geometry anymore */
1411  GEOSGeom_destroy(g1);
1412 
1413  if (g3 == NULL)
1414  {
1415  lwerror("GEOSDelaunayTriangulation: %s", lwgeom_geos_errmsg);
1416  return NULL;
1417  }
1418 
1419  /* LWDEBUGF(3, "result: %s", GEOSGeomToWKT(g3)); */
1420 
1421  GEOSSetSRID(g3, lwgeom_get_srid(lwgeom_in));
1422 
1423  if (output == 2) {
1424  lwgeom_result = (LWGEOM *)lwtin_from_geos(g3, lwgeom_has_z(lwgeom_in));
1425  } else {
1426  lwgeom_result = GEOS2LWGEOM(g3, lwgeom_has_z(lwgeom_in));
1427  }
1428 
1429  GEOSGeom_destroy(g3);
1430 
1431  if (lwgeom_result == NULL) {
1432  if (output != 2) {
1433  lwerror("lwgeom_delaunay_triangulation: GEOS2LWGEOM returned null");
1434  } else {
1435  lwerror("lwgeom_delaunay_triangulation: lwtin_from_geos returned null");
1436  }
1437  return NULL;
1438  }
1439 
1440  return lwgeom_result;
1441 
1442 #endif /* POSTGIS_GEOS_VERSION < 34 */
1443 }
void ptarray_set_point4d(POINTARRAY *pa, int n, const POINT4D *p4d)
Definition: lwgeom_api.c:501
double x
Definition: liblwgeom.h:308
#define LINETYPE
Definition: liblwgeom.h:61
double z
Definition: liblwgeom.h:290
#define POSTGIS_GEOS_VERSION
Definition: sqldefines.h:10
double y
Definition: liblwgeom.h:290
POINTARRAY * ptarray_construct(char hasz, char hasm, uint32_t npoints)
Construct an empty pointarray, allocating storage and setting the npoints, but not filling in any inf...
Definition: ptarray.c:49
LWLINE * lwline_construct_empty(int srid, char hasz, char hasm)
Definition: lwline.c:51
LWCOLLECTION * lwcollection_construct(uint8_t type, int srid, GBOX *bbox, uint32_t ngeoms, LWGEOM **geoms)
Definition: lwcollection.c:30
double x
Definition: liblwgeom.h:290
static void delFace(Face *f)
void lwfree(void *mem)
Definition: lwutil.c:190
int npoints
Definition: liblwgeom.h:327
void error_if_srid_mismatch(int srid1, int srid2)
Definition: lwutil.c:317
static int compare_by_envarea(const void *g1, const void *g2)
#define POLYGONTYPE
Definition: liblwgeom.h:62
char * lwgeom_to_ewkt(const LWGEOM *lwgeom)
Return an alloced string.
Definition: lwgeom.c:425
uint8_t flags
Definition: liblwgeom.h:353
POINTARRAY * ptarray_construct_empty(char hasz, char hasm, uint32_t maxpoints)
Create a new POINTARRAY with no points.
Definition: ptarray.c:57
#define MULTIPOINTTYPE
Definition: liblwgeom.h:63
LWGEOM * lwgeom_symdifference(const LWGEOM *geom1, const LWGEOM *geom2)
tuple fmt
Definition: pixval.py:92
char lwgeom_geos_errmsg[LWGEOM_GEOS_ERRMSG_MAXSIZE]
#define LWDEBUG(level, msg)
Definition: lwgeom_log.h:50
int32_t lwgeom_get_srid(const LWGEOM *geom)
Return SRID number.
Definition: lwgeom.c:778
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:389
GEOSGeometry * LWGEOM2GEOS(const LWGEOM *lwgeom)
LWPOINT * lwpoint_construct_empty(int srid, char hasz, char hasm)
Definition: lwpoint.c:118
POINTARRAY * point
Definition: liblwgeom.h:367
int lwgeom_has_z(const LWGEOM *geom)
Return LW_TRUE if geometry has Z ordinates.
Definition: lwgeom.c:792
char ** result
Definition: liblwgeom.h:218
int32_t srid
Definition: liblwgeom.h:355
LWGEOM * lwgeom_intersection(const LWGEOM *geom1, const LWGEOM *geom2)
void lwerror(const char *fmt,...)
Write a notice out to the error handler.
Definition: lwutil.c:67
double x
Definition: liblwgeom.h:284
GEOSCoordSeq ptarray_to_GEOSCoordSeq(const POINTARRAY *)
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:481
void lwnotice(const char *fmt,...)
Write a notice out to the notice handler.
Definition: lwutil.c:54
const char * lwtype_name(uint8_t type)
Return the type name string associated with a type number (e.g.
Definition: lwutil.c:164
LWGEOM * lwline_as_lwgeom(const LWLINE *obj)
Definition: lwgeom.c:249
void lwgeom_geos_error(const char *fmt,...)
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:458
uint8_t flags
Definition: liblwgeom.h:325
LWPOLY * lwpoly_construct(int srid, GBOX *bbox, uint32_t nrings, POINTARRAY **points)
Definition: lwpoly.c:29
LWLINE * lwline_construct(int srid, GBOX *bbox, POINTARRAY *points)
Definition: lwline.c:29
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:465
static GEOSGeometry * collectFacesWithEvenAncestors(Face **faces, int nfaces)
#define SRID_UNKNOWN
Unknown SRID value.
Definition: liblwgeom.h:154
#define TINTYPE
Definition: liblwgeom.h:74
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:1645
LWGEOM * lwgeom_normalize(const LWGEOM *geom1)
LWGEOM * lwgeom_difference(const LWGEOM *geom1, const LWGEOM *geom2)
POINTARRAY ** rings
Definition: liblwgeom.h:413
LWGEOM * lwgeom_delaunay_triangulation(const LWGEOM *lwgeom_in, double tolerance, int output)
Take vertices of a geometry and build a delaunay triangulation on them.
int nrings
Definition: liblwgeom.h:411
int lwgeom_has_arc(const LWGEOM *geom)
Definition: lwsegmentize.c:42
static Face * newFace(const GEOSGeometry *g)
double y
Definition: liblwgeom.h:284
#define LWGEOM_GEOS_ERRMSG_MAXSIZE
#define FLAGS_GET_Z(flags)
Macros for manipulating the 'flags' byte.
Definition: liblwgeom.h:106
double z
Definition: liblwgeom.h:308
LWTRIANGLE * lwtriangle_construct(int srid, GBOX *bbox, POINTARRAY *points)
Definition: lwtriangle.c:27
GEOSGeometry * env
#define MULTIPOLYGONTYPE
Definition: liblwgeom.h:65
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:472
LWGEOM * lwgeom_snap(const LWGEOM *geom1, const LWGEOM *geom2, double tolerance)
Snap vertices and segments of a geometry to another using a given tolerance.
#define POINTTYPE
LWTYPE numbers, used internally by PostGIS.
Definition: liblwgeom.h:60
LWPOLY * lwpoly_construct_empty(int srid, char hasz, char hasm)
Definition: lwpoly.c:66
LWGEOM * GEOS2LWGEOM(const GEOSGeometry *geom, char want3d)
uint8_t type
Definition: liblwgeom.h:352
struct Face_t * parent
const char * lwgeom_geos_version()
Return GEOS version string (not to be freed)
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:96
void * lwalloc(size_t size)
Definition: lwutil.c:175
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:1229
double y
Definition: liblwgeom.h:308
#define MULTILINETYPE
Definition: liblwgeom.h:64
POINTARRAY * ptarray_from_GEOSCoordSeq(const GEOSCoordSequence *cs, char want3d)
#define LWDEBUGF(level, msg,...)
Definition: lwgeom_log.h:55
#define FLAGS_NDIMS(flags)
Definition: liblwgeom.h:118
int lwgeom_has_m(const LWGEOM *geom)
Return LW_TRUE if geometry has M ordinates.
Definition: lwgeom.c:799
#define COLLECTIONTYPE
Definition: liblwgeom.h:66
This library is the generic geometry handling section of PostGIS.
GEOSGeometry * LWGEOM_GEOS_buildArea(const GEOSGeometry *geom_in)
POINTARRAY * points
Definition: liblwgeom.h:378