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