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