PostGIS  2.5.2dev-r@@SVN_REVISION@@
liblwgeom/lwgeom_geos_clean.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 2009-2010 Sandro Santilli <strk@kbt.io>
22  *
23  **********************************************************************/
24 
25 #include "liblwgeom.h"
26 #include "lwgeom_geos.h"
27 #include "liblwgeom_internal.h"
28 #include "lwgeom_log.h"
29 
30 #include <string.h>
31 #include <stdlib.h>
32 #include <assert.h>
33 
34 /* #define POSTGIS_DEBUG_LEVEL 4 */
35 /* #define PARANOIA_LEVEL 2 */
36 #undef LWGEOM_PROFILE_MAKEVALID
37 
38 /*
39  * Return Nth vertex in GEOSGeometry as a POINT.
40  * May return NULL if the geometry has NO vertex.
41  */
42 GEOSGeometry* LWGEOM_GEOS_getPointN(const GEOSGeometry*, uint32_t);
43 GEOSGeometry*
44 LWGEOM_GEOS_getPointN(const GEOSGeometry* g_in, uint32_t n)
45 {
46  uint32_t dims = 0;
47  const GEOSCoordSequence* seq_in;
48  GEOSCoordSeq seq_out;
49  double val;
50  uint32_t sz = 0;
51  int gn;
52  GEOSGeometry* ret;
53 
54  switch (GEOSGeomTypeId(g_in))
55  {
56  case GEOS_MULTIPOINT:
57  case GEOS_MULTILINESTRING:
58  case GEOS_MULTIPOLYGON:
59  case GEOS_GEOMETRYCOLLECTION:
60  {
61  for (gn = 0; gn < GEOSGetNumGeometries(g_in); ++gn)
62  {
63  const GEOSGeometry* g = GEOSGetGeometryN(g_in, gn);
64  ret = LWGEOM_GEOS_getPointN(g, n);
65  if (ret) return ret;
66  }
67  break;
68  }
69 
70  case GEOS_POLYGON:
71  {
72  ret = LWGEOM_GEOS_getPointN(GEOSGetExteriorRing(g_in), n);
73  if (ret) return ret;
74  for (gn = 0; gn < GEOSGetNumInteriorRings(g_in); ++gn)
75  {
76  const GEOSGeometry* g = GEOSGetInteriorRingN(g_in, gn);
77  ret = LWGEOM_GEOS_getPointN(g, n);
78  if (ret) return ret;
79  }
80  break;
81  }
82 
83  case GEOS_POINT:
84  case GEOS_LINESTRING:
85  case GEOS_LINEARRING:
86  break;
87  }
88 
89  seq_in = GEOSGeom_getCoordSeq(g_in);
90  if (!seq_in) return NULL;
91  if (!GEOSCoordSeq_getSize(seq_in, &sz)) return NULL;
92  if (!sz) return NULL;
93 
94  if (!GEOSCoordSeq_getDimensions(seq_in, &dims)) return NULL;
95 
96  seq_out = GEOSCoordSeq_create(1, dims);
97  if (!seq_out) return NULL;
98 
99  if (!GEOSCoordSeq_getX(seq_in, n, &val)) return NULL;
100  if (!GEOSCoordSeq_setX(seq_out, n, val)) return NULL;
101  if (!GEOSCoordSeq_getY(seq_in, n, &val)) return NULL;
102  if (!GEOSCoordSeq_setY(seq_out, n, val)) return NULL;
103  if (dims > 2)
104  {
105  if (!GEOSCoordSeq_getZ(seq_in, n, &val)) return NULL;
106  if (!GEOSCoordSeq_setZ(seq_out, n, val)) return NULL;
107  }
108 
109  return GEOSGeom_createPoint(seq_out);
110 }
111 
116 
117 /*
118  * Ensure the geometry is "structurally" valid
119  * (enough for GEOS to accept it)
120  * May return the input untouched (if already valid).
121  * May return geometries of lower dimension (on collapses)
122  */
123 static LWGEOM*
125 {
126  LWDEBUGF(2, "lwgeom_make_geos_friendly enter (type %d)", geom->type);
127  switch (geom->type)
128  {
129  case POINTTYPE:
130  case MULTIPOINTTYPE:
131  /* a point is always valid */
132  return geom;
133  break;
134 
135  case LINETYPE:
136  /* lines need at least 2 points */
137  return lwline_make_geos_friendly((LWLINE*)geom);
138  break;
139 
140  case POLYGONTYPE:
141  /* polygons need all rings closed and with npoints > 3 */
142  return lwpoly_make_geos_friendly((LWPOLY*)geom);
143  break;
144 
145  case MULTILINETYPE:
146  case MULTIPOLYGONTYPE:
147  case COLLECTIONTYPE:
149  break;
150 
151  case CIRCSTRINGTYPE:
152  case COMPOUNDTYPE:
153  case CURVEPOLYTYPE:
154  case MULTISURFACETYPE:
155  case MULTICURVETYPE:
156  default:
157  lwerror("lwgeom_make_geos_friendly: unsupported input geometry type: %s (%d)",
158  lwtype_name(geom->type),
159  geom->type);
160  break;
161  }
162  return 0;
163 }
164 
165 /*
166  * Close the point array, if not already closed in 2d.
167  * Returns the input if already closed in 2d, or a newly
168  * constructed POINTARRAY.
169  * TODO: move in ptarray.c
170  */
172 POINTARRAY*
174 {
175  POINTARRAY* newring;
176 
177  /* close the ring if not already closed (2d only) */
178  if (!ptarray_is_closed_2d(ring))
179  {
180  /* close it up */
181  newring = ptarray_addPoint(ring, getPoint_internal(ring, 0), FLAGS_NDIMS(ring->flags), ring->npoints);
182  ring = newring;
183  }
184  return ring;
185 }
186 
187 /* May return the same input or a new one (never zero) */
188 POINTARRAY*
190 {
191  POINTARRAY* closedring;
192  POINTARRAY* ring_in = ring;
193 
194  /* close the ring if not already closed (2d only) */
195  closedring = ptarray_close2d(ring);
196  if (closedring != ring) ring = closedring;
197 
198  /* return 0 for collapsed ring (after closeup) */
199 
200  while (ring->npoints < 4)
201  {
202  POINTARRAY* oring = ring;
203  LWDEBUGF(4, "ring has %d points, adding another", ring->npoints);
204  /* let's add another... */
205  ring = ptarray_addPoint(ring, getPoint_internal(ring, 0), FLAGS_NDIMS(ring->flags), ring->npoints);
206  if (oring != ring_in) ptarray_free(oring);
207  }
208 
209  return ring;
210 }
211 
212 /* Make sure all rings are closed and have > 3 points.
213  * May return the input untouched.
214  */
215 LWGEOM*
217 {
218  LWGEOM* ret;
219  POINTARRAY** new_rings;
220  uint32_t i;
221 
222  /* If the polygon has no rings there's nothing to do */
223  if (!poly->nrings) return (LWGEOM*)poly;
224 
225  /* Allocate enough pointers for all rings */
226  new_rings = lwalloc(sizeof(POINTARRAY*) * poly->nrings);
227 
228  /* All rings must be closed and have > 3 points */
229  for (i = 0; i < poly->nrings; i++)
230  {
231  POINTARRAY* ring_in = poly->rings[i];
232  POINTARRAY* ring_out = ring_make_geos_friendly(ring_in);
233 
234  if (ring_in != ring_out)
235  {
236  LWDEBUGF(
237  3, "lwpoly_make_geos_friendly: ring %d cleaned, now has %d points", i, ring_out->npoints);
238  ptarray_free(ring_in);
239  }
240  else
241  LWDEBUGF(3, "lwpoly_make_geos_friendly: ring %d untouched", i);
242 
243  assert(ring_out);
244  new_rings[i] = ring_out;
245  }
246 
247  lwfree(poly->rings);
248  poly->rings = new_rings;
249  ret = (LWGEOM*)poly;
250 
251  return ret;
252 }
253 
254 /* Need NO or >1 points. Duplicate first if only one. */
255 LWGEOM*
257 {
258  LWGEOM* ret;
259 
260  if (line->points->npoints == 1) /* 0 is fine, 2 is fine */
261  {
262 #if 1
263  /* Duplicate point */
264  line->points = ptarray_addPoint(line->points,
265  getPoint_internal(line->points, 0),
266  FLAGS_NDIMS(line->points->flags),
267  line->points->npoints);
268  ret = (LWGEOM*)line;
269 #else
270  /* Turn into a point */
271  ret = (LWGEOM*)lwpoint_construct(line->srid, 0, line->points);
272 #endif
273  return ret;
274  }
275  else
276  {
277  return (LWGEOM*)line;
278  /* return lwline_clone(line); */
279  }
280 }
281 
282 LWGEOM*
284 {
285  LWGEOM** new_geoms;
286  uint32_t i, new_ngeoms = 0;
287  LWCOLLECTION* ret;
288 
289  /* enough space for all components */
290  new_geoms = lwalloc(sizeof(LWGEOM*) * g->ngeoms);
291 
292  ret = lwalloc(sizeof(LWCOLLECTION));
293  memcpy(ret, g, sizeof(LWCOLLECTION));
294  ret->maxgeoms = g->ngeoms;
295 
296  for (i = 0; i < g->ngeoms; i++)
297  {
298  LWGEOM* newg = lwgeom_make_geos_friendly(g->geoms[i]);
299  if (newg) new_geoms[new_ngeoms++] = newg;
300  }
301 
302  ret->bbox = NULL; /* recompute later... */
303 
304  ret->ngeoms = new_ngeoms;
305  if (new_ngeoms)
306  ret->geoms = new_geoms;
307  else
308  {
309  free(new_geoms);
310  ret->geoms = NULL;
311  ret->maxgeoms = 0;
312  }
313 
314  return (LWGEOM*)ret;
315 }
316 
317 /*
318  * Fully node given linework
319  */
320 static GEOSGeometry*
321 LWGEOM_GEOS_nodeLines(const GEOSGeometry* lines)
322 {
323  /* GEOS3.7 GEOSNode fails on regression tests */
324  /* GEOS3.7 GEOSUnaryUnion fails on regression tests */
325 
326  /* union of first point with geometry */
327  GEOSGeometry *unioned, *point;
328  point = LWGEOM_GEOS_getPointN(lines, 0);
329  if (!point) return NULL;
330  unioned = GEOSUnion(lines, point);
331  GEOSGeom_destroy(point);
332  return unioned;
333 }
334 
335 /*
336  * We expect initGEOS being called already.
337  * Will return NULL on error (expect error handler being called by then)
338  */
339 static GEOSGeometry*
340 LWGEOM_GEOS_makeValidPolygon(const GEOSGeometry* gin)
341 {
342  GEOSGeom gout;
343  GEOSGeom geos_bound;
344  GEOSGeom geos_cut_edges, geos_area, collapse_points;
345  GEOSGeometry* vgeoms[3]; /* One for area, one for cut-edges */
346  unsigned int nvgeoms = 0;
347 
348  assert(GEOSGeomTypeId(gin) == GEOS_POLYGON || GEOSGeomTypeId(gin) == GEOS_MULTIPOLYGON);
349 
350  geos_bound = GEOSBoundary(gin);
351  if (!geos_bound) return NULL;
352 
353  /* Use noded boundaries as initial "cut" edges */
354 
355  geos_cut_edges = LWGEOM_GEOS_nodeLines(geos_bound);
356  if (!geos_cut_edges)
357  {
358  GEOSGeom_destroy(geos_bound);
359  lwnotice("LWGEOM_GEOS_nodeLines(): %s", lwgeom_geos_errmsg);
360  return NULL;
361  }
362 
363  /* NOTE: the noding process may drop lines collapsing to points.
364  * We want to retrieve any of those */
365  {
366  GEOSGeometry* pi;
367  GEOSGeometry* po;
368 
369  pi = GEOSGeom_extractUniquePoints(geos_bound);
370  if (!pi)
371  {
372  GEOSGeom_destroy(geos_bound);
373  lwnotice("GEOSGeom_extractUniquePoints(): %s", lwgeom_geos_errmsg);
374  return NULL;
375  }
376 
377  po = GEOSGeom_extractUniquePoints(geos_cut_edges);
378  if (!po)
379  {
380  GEOSGeom_destroy(geos_bound);
381  GEOSGeom_destroy(pi);
382  lwnotice("GEOSGeom_extractUniquePoints(): %s", lwgeom_geos_errmsg);
383  return NULL;
384  }
385 
386  collapse_points = GEOSDifference(pi, po);
387  if (!collapse_points)
388  {
389  GEOSGeom_destroy(geos_bound);
390  GEOSGeom_destroy(pi);
391  GEOSGeom_destroy(po);
392  lwnotice("GEOSDifference(): %s", lwgeom_geos_errmsg);
393  return NULL;
394  }
395 
396  GEOSGeom_destroy(pi);
397  GEOSGeom_destroy(po);
398  }
399  GEOSGeom_destroy(geos_bound);
400 
401  /* And use an empty geometry as initial "area" */
402  geos_area = GEOSGeom_createEmptyPolygon();
403  if (!geos_area)
404  {
405  lwnotice("GEOSGeom_createEmptyPolygon(): %s", lwgeom_geos_errmsg);
406  GEOSGeom_destroy(geos_cut_edges);
407  return NULL;
408  }
409 
410  /*
411  * See if an area can be build with the remaining edges
412  * and if it can, symdifference with the original area.
413  * Iterate this until no more polygons can be created
414  * with left-over edges.
415  */
416  while (GEOSGetNumGeometries(geos_cut_edges))
417  {
418  GEOSGeometry* new_area = 0;
419  GEOSGeometry* new_area_bound = 0;
420  GEOSGeometry* symdif = 0;
421  GEOSGeometry* new_cut_edges = 0;
422 
423 #ifdef LWGEOM_PROFILE_MAKEVALID
424  lwnotice("ST_MakeValid: building area from %d edges", GEOSGetNumGeometries(geos_cut_edges));
425 #endif
426 
427  /*
428  * ASSUMPTION: cut_edges should already be fully noded
429  */
430 
431  new_area = LWGEOM_GEOS_buildArea(geos_cut_edges);
432  if (!new_area) /* must be an exception */
433  {
434  GEOSGeom_destroy(geos_cut_edges);
435  GEOSGeom_destroy(geos_area);
436  lwnotice("LWGEOM_GEOS_buildArea() threw an error: %s", lwgeom_geos_errmsg);
437  return NULL;
438  }
439 
440  if (GEOSisEmpty(new_area))
441  {
442  /* no more rings can be build with thes edges */
443  GEOSGeom_destroy(new_area);
444  break;
445  }
446 
447  /*
448  * We succeeded in building a ring!
449  * Save the new ring boundaries first (to compute
450  * further cut edges later)
451  */
452  new_area_bound = GEOSBoundary(new_area);
453  if (!new_area_bound)
454  {
455  /* We did check for empty area already so this must be some other error */
456  lwnotice("GEOSBoundary('%s') threw an error: %s",
457  lwgeom_to_ewkt(GEOS2LWGEOM(new_area, 0)),
459  GEOSGeom_destroy(new_area);
460  GEOSGeom_destroy(geos_area);
461  return NULL;
462  }
463 
464  /*
465  * Now symdif new and old area
466  */
467  symdif = GEOSSymDifference(geos_area, new_area);
468  if (!symdif) /* must be an exception */
469  {
470  GEOSGeom_destroy(geos_cut_edges);
471  GEOSGeom_destroy(new_area);
472  GEOSGeom_destroy(new_area_bound);
473  GEOSGeom_destroy(geos_area);
474  lwnotice("GEOSSymDifference() threw an error: %s", lwgeom_geos_errmsg);
475  return NULL;
476  }
477 
478  GEOSGeom_destroy(geos_area);
479  GEOSGeom_destroy(new_area);
480  geos_area = symdif;
481  symdif = 0;
482 
483  /*
484  * Now let's re-set geos_cut_edges with what's left
485  * from the original boundary.
486  * ASSUMPTION: only the previous cut-edges can be
487  * left, so we don't need to reconsider
488  * the whole original boundaries
489  *
490  * NOTE: this is an expensive operation.
491  *
492  */
493 
494  new_cut_edges = GEOSDifference(geos_cut_edges, new_area_bound);
495  GEOSGeom_destroy(new_area_bound);
496  if (!new_cut_edges) /* an exception ? */
497  {
498  /* cleanup and throw */
499  GEOSGeom_destroy(geos_cut_edges);
500  GEOSGeom_destroy(geos_area);
501  lwerror("GEOSDifference() threw an error: %s", lwgeom_geos_errmsg);
502  return NULL;
503  }
504  GEOSGeom_destroy(geos_cut_edges);
505  geos_cut_edges = new_cut_edges;
506  }
507 
508  if (!GEOSisEmpty(geos_area))
509  vgeoms[nvgeoms++] = geos_area;
510  else
511  GEOSGeom_destroy(geos_area);
512 
513  if (!GEOSisEmpty(geos_cut_edges))
514  vgeoms[nvgeoms++] = geos_cut_edges;
515  else
516  GEOSGeom_destroy(geos_cut_edges);
517 
518  if (!GEOSisEmpty(collapse_points))
519  vgeoms[nvgeoms++] = collapse_points;
520  else
521  GEOSGeom_destroy(collapse_points);
522 
523  if (1 == nvgeoms)
524  {
525  /* Return cut edges */
526  gout = vgeoms[0];
527  }
528  else
529  {
530  /* Collect areas and lines (if any line) */
531  gout = GEOSGeom_createCollection(GEOS_GEOMETRYCOLLECTION, vgeoms, nvgeoms);
532  if (!gout) /* an exception again */
533  {
534  /* cleanup and throw */
535  lwerror("GEOSGeom_createCollection() threw an error: %s", lwgeom_geos_errmsg);
536  /* TODO: cleanup! */
537  return NULL;
538  }
539  }
540 
541  return gout;
542 }
543 
544 static GEOSGeometry*
545 LWGEOM_GEOS_makeValidLine(const GEOSGeometry* gin)
546 {
547  GEOSGeometry* noded;
548  noded = LWGEOM_GEOS_nodeLines(gin);
549  return noded;
550 }
551 
552 static GEOSGeometry*
553 LWGEOM_GEOS_makeValidMultiLine(const GEOSGeometry* gin)
554 {
555  GEOSGeometry** lines;
556  GEOSGeometry** points;
557  GEOSGeometry* mline_out = 0;
558  GEOSGeometry* mpoint_out = 0;
559  GEOSGeometry* gout = 0;
560  uint32_t nlines = 0, nlines_alloc;
561  uint32_t npoints = 0;
562  uint32_t ngeoms = 0, nsubgeoms;
563  uint32_t i, j;
564 
565  ngeoms = GEOSGetNumGeometries(gin);
566 
567  nlines_alloc = ngeoms;
568  lines = lwalloc(sizeof(GEOSGeometry*) * nlines_alloc);
569  points = lwalloc(sizeof(GEOSGeometry*) * ngeoms);
570 
571  for (i = 0; i < ngeoms; ++i)
572  {
573  const GEOSGeometry* g = GEOSGetGeometryN(gin, i);
574  GEOSGeometry* vg;
576  if (GEOSisEmpty(vg))
577  {
578  /* we don't care about this one */
579  GEOSGeom_destroy(vg);
580  }
581  if (GEOSGeomTypeId(vg) == GEOS_POINT)
582  points[npoints++] = vg;
583  else if (GEOSGeomTypeId(vg) == GEOS_LINESTRING)
584  lines[nlines++] = vg;
585  else if (GEOSGeomTypeId(vg) == GEOS_MULTILINESTRING)
586  {
587  nsubgeoms = GEOSGetNumGeometries(vg);
588  nlines_alloc += nsubgeoms;
589  lines = lwrealloc(lines, sizeof(GEOSGeometry*) * nlines_alloc);
590  for (j = 0; j < nsubgeoms; ++j)
591  {
592  const GEOSGeometry* gc = GEOSGetGeometryN(vg, j);
593  /* NOTE: ownership of the cloned geoms will be
594  * taken by final collection */
595  lines[nlines++] = GEOSGeom_clone(gc);
596  }
597  }
598  else
599  {
600  /* NOTE: return from GEOSGeomType will leak
601  * but we really don't expect this to happen */
602  lwerror("unexpected geom type returned by LWGEOM_GEOS_makeValid: %s", GEOSGeomType(vg));
603  }
604  }
605 
606  if (npoints)
607  {
608  if (npoints > 1)
609  mpoint_out = GEOSGeom_createCollection(GEOS_MULTIPOINT, points, npoints);
610  else
611  mpoint_out = points[0];
612  }
613 
614  if (nlines)
615  {
616  if (nlines > 1)
617  mline_out = GEOSGeom_createCollection(GEOS_MULTILINESTRING, lines, nlines);
618  else
619  mline_out = lines[0];
620  }
621 
622  lwfree(lines);
623 
624  if (mline_out && mpoint_out)
625  {
626  points[0] = mline_out;
627  points[1] = mpoint_out;
628  gout = GEOSGeom_createCollection(GEOS_GEOMETRYCOLLECTION, points, 2);
629  }
630  else if (mline_out)
631  gout = mline_out;
632 
633  else if (mpoint_out)
634  gout = mpoint_out;
635 
636  lwfree(points);
637 
638  return gout;
639 }
640 
641 /*
642  * We expect initGEOS being called already.
643  * Will return NULL on error (expect error handler being called by then)
644  */
645 static GEOSGeometry*
646 LWGEOM_GEOS_makeValidCollection(const GEOSGeometry* gin)
647 {
648  int nvgeoms;
649  GEOSGeometry** vgeoms;
650  GEOSGeom gout;
651  int i;
652 
653  nvgeoms = GEOSGetNumGeometries(gin);
654  if (nvgeoms == -1)
655  {
656  lwerror("GEOSGetNumGeometries: %s", lwgeom_geos_errmsg);
657  return 0;
658  }
659 
660  vgeoms = lwalloc(sizeof(GEOSGeometry*) * nvgeoms);
661  if (!vgeoms)
662  {
663  lwerror("LWGEOM_GEOS_makeValidCollection: out of memory");
664  return 0;
665  }
666 
667  for (i = 0; i < nvgeoms; ++i)
668  {
669  vgeoms[i] = LWGEOM_GEOS_makeValid(GEOSGetGeometryN(gin, i));
670  if (!vgeoms[i])
671  {
672  int j;
673  for (j = 0; j < i - 1; j++)
674  GEOSGeom_destroy(vgeoms[j]);
675  lwfree(vgeoms);
676  /* we expect lwerror being called already by makeValid */
677  return NULL;
678  }
679  }
680 
681  /* Collect areas and lines (if any line) */
682  gout = GEOSGeom_createCollection(GEOS_GEOMETRYCOLLECTION, vgeoms, nvgeoms);
683  if (!gout) /* an exception again */
684  {
685  /* cleanup and throw */
686  for (i = 0; i < nvgeoms; ++i)
687  GEOSGeom_destroy(vgeoms[i]);
688  lwfree(vgeoms);
689  lwerror("GEOSGeom_createCollection() threw an error: %s", lwgeom_geos_errmsg);
690  return NULL;
691  }
692  lwfree(vgeoms);
693 
694  return gout;
695 }
696 
697 GEOSGeometry*
698 LWGEOM_GEOS_makeValid(const GEOSGeometry* gin)
699 {
700  GEOSGeometry* gout;
701  char ret_char;
702 #if POSTGIS_DEBUG_LEVEL >= 3
703  LWGEOM *geos_geom;
704  char *geom_ewkt;
705 #endif
706 
707  /*
708  * Step 2: return what we got so far if already valid
709  */
710 
711  ret_char = GEOSisValid(gin);
712  if (ret_char == 2)
713  {
714  /* I don't think should ever happen */
715  lwerror("GEOSisValid(): %s", lwgeom_geos_errmsg);
716  return NULL;
717  }
718  else if (ret_char)
719  {
720 #if POSTGIS_DEBUG_LEVEL >= 3
721  geos_geom = GEOS2LWGEOM(gin, 0);
722  geom_ewkt = lwgeom_to_ewkt(geos_geom);
723  LWDEBUGF(3, "Geometry [%s] is valid. ", geom_ewkt);
724  lwgeom_free(geos_geom);
725  lwfree(geom_ewkt);
726 #endif
727 
728  /* It's valid at this step, return what we have */
729  return GEOSGeom_clone(gin);
730  }
731 
732 #if POSTGIS_DEBUG_LEVEL >= 3
733  geos_geom = GEOS2LWGEOM(gin, 0);
734  geom_ewkt = lwgeom_to_ewkt(geos_geom);
735  LWDEBUGF(3,
736  "Geometry [%s] is still not valid: %s. Will try to clean up further.",
737  geom_ewkt,
739  lwgeom_free(geos_geom);
740  lwfree(geom_ewkt);
741 #endif
742 
743  /*
744  * Step 3 : make what we got valid
745  */
746 
747  switch (GEOSGeomTypeId(gin))
748  {
749  case GEOS_MULTIPOINT:
750  case GEOS_POINT:
751  /* points are always valid, but we might have invalid ordinate values */
752  lwnotice("PUNTUAL geometry resulted invalid to GEOS -- dunno how to clean that up");
753  return NULL;
754  break;
755 
756  case GEOS_LINESTRING:
757  gout = LWGEOM_GEOS_makeValidLine(gin);
758  if (!gout) /* an exception or something */
759  {
760  /* cleanup and throw */
762  return NULL;
763  }
764  break; /* we've done */
765 
766  case GEOS_MULTILINESTRING:
767  gout = LWGEOM_GEOS_makeValidMultiLine(gin);
768  if (!gout) /* an exception or something */
769  {
770  /* cleanup and throw */
772  return NULL;
773  }
774  break; /* we've done */
775 
776  case GEOS_POLYGON:
777  case GEOS_MULTIPOLYGON:
778  {
779  gout = LWGEOM_GEOS_makeValidPolygon(gin);
780  if (!gout) /* an exception or something */
781  {
782  /* cleanup and throw */
784  return NULL;
785  }
786  break; /* we've done */
787  }
788 
789  case GEOS_GEOMETRYCOLLECTION:
790  {
792  if (!gout) /* an exception or something */
793  {
794  /* cleanup and throw */
796  return NULL;
797  }
798  break; /* we've done */
799  }
800 
801  default:
802  {
803  char* typname = GEOSGeomType(gin);
804  lwnotice("ST_MakeValid: doesn't support geometry type: %s", typname);
805  GEOSFree(typname);
806  return NULL;
807  break;
808  }
809  }
810 
811 #if PARANOIA_LEVEL > 1
812  /*
813  * Now check if every point of input is also found in output, or abort by returning NULL
814  *
815  * Input geometry was lwgeom_in
816  */
817  {
818  int loss;
819  GEOSGeometry *pi, *po, *pd;
820 
821  /* TODO: handle some errors here...
822  * Lack of exceptions is annoying indeed,
823  * I'm getting old --strk;
824  */
825  pi = GEOSGeom_extractUniquePoints(gin);
826  po = GEOSGeom_extractUniquePoints(gout);
827  pd = GEOSDifference(pi, po); /* input points - output points */
828  GEOSGeom_destroy(pi);
829  GEOSGeom_destroy(po);
830  loss = !GEOSisEmpty(pd);
831  GEOSGeom_destroy(pd);
832  if (loss)
833  {
834  lwnotice("%s [%d] Vertices lost in LWGEOM_GEOS_makeValid", __FILE__, __LINE__);
835  /* return NULL */
836  }
837  }
838 #endif /* PARANOIA_LEVEL > 1 */
839 
840  return gout;
841 }
842 
843 /* Exported. Uses GEOS internally */
844 LWGEOM*
846 {
847  int is3d;
848  GEOSGeom geosgeom;
849  GEOSGeometry* geosout;
850  LWGEOM* lwgeom_out;
851 
852  is3d = FLAGS_GET_Z(lwgeom_in->flags);
853 
854  /*
855  * Step 1 : try to convert to GEOS, if impossible, clean that up first
856  * otherwise (adding only duplicates of existing points)
857  */
858 
860 
861  lwgeom_out = lwgeom_in;
862  geosgeom = LWGEOM2GEOS(lwgeom_out, 1);
863  if (!geosgeom)
864  {
865  LWDEBUGF(4,
866  "Original geom can't be converted to GEOS (%s)"
867  " - will try cleaning that up first",
869 
870  lwgeom_out = lwgeom_make_geos_friendly(lwgeom_out);
871  if (!lwgeom_out) lwerror("Could not make a valid geometry out of input");
872 
873  /* try again as we did cleanup now */
874  /* TODO: invoke LWGEOM2GEOS directly with autoclean ? */
875  geosgeom = LWGEOM2GEOS(lwgeom_out, 0);
876  if (!geosgeom)
877  {
878  lwerror("Couldn't convert POSTGIS geom to GEOS: %s", lwgeom_geos_errmsg);
879  return NULL;
880  }
881  }
882  else
883  {
884  LWDEBUG(4, "original geom converted to GEOS");
885  lwgeom_out = lwgeom_in;
886  }
887 
888  geosout = LWGEOM_GEOS_makeValid(geosgeom);
889  GEOSGeom_destroy(geosgeom);
890  if (!geosout) return NULL;
891 
892  lwgeom_out = GEOS2LWGEOM(geosout, is3d);
893  GEOSGeom_destroy(geosout);
894 
895  if (lwgeom_is_collection(lwgeom_in) && !lwgeom_is_collection(lwgeom_out))
896  {
897  LWGEOM** ogeoms = lwalloc(sizeof(LWGEOM*));
898  LWGEOM* ogeom;
899  LWDEBUG(3, "lwgeom_make_valid: forcing multi");
900  /* NOTE: this is safe because lwgeom_out is surely not lwgeom_in or
901  * otherwise we couldn't have a collection and a non-collection */
902  assert(lwgeom_in != lwgeom_out);
903  ogeoms[0] = lwgeom_out;
904  ogeom = (LWGEOM*)lwcollection_construct(
905  MULTITYPE[lwgeom_out->type], lwgeom_out->srid, lwgeom_out->bbox, 1, ogeoms);
906  lwgeom_out->bbox = NULL;
907  lwgeom_out = ogeom;
908  }
909 
910  lwgeom_out->srid = lwgeom_in->srid;
911  return lwgeom_out;
912 }
#define LINETYPE
Definition: liblwgeom.h:85
POINTARRAY * ptarray_close2d(POINTARRAY *ring)
GBOX * bbox
Definition: liblwgeom.h:400
LWGEOM * lwgeom_make_valid(LWGEOM *lwgeom_in)
Attempts to make an invalid geometries valid w/out losing points.
#define MULTICURVETYPE
Definition: liblwgeom.h:94
static GEOSGeometry * LWGEOM_GEOS_makeValidPolygon(const GEOSGeometry *gin)
int lwgeom_is_collection(const LWGEOM *lwgeom)
Determine whether a LWGEOM can contain sub-geometries or not.
Definition: lwgeom.c:1085
LWGEOM * GEOS2LWGEOM(const GEOSGeometry *geom, uint8_t want3d)
static LWGEOM * lwgeom_make_geos_friendly(LWGEOM *geom)
LWCOLLECTION * lwcollection_construct(uint8_t type, int srid, GBOX *bbox, uint32_t ngeoms, LWGEOM **geoms)
Definition: lwcollection.c:43
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
uint8_t type
Definition: liblwgeom.h:505
GEOSGeometry * LWGEOM_GEOS_makeValid(const GEOSGeometry *gin)
#define POLYGONTYPE
Definition: liblwgeom.h:86
char * lwgeom_to_ewkt(const LWGEOM *lwgeom)
Return an alloced string.
Definition: lwgeom.c:556
uint8_t flags
Definition: liblwgeom.h:399
#define CURVEPOLYTYPE
Definition: liblwgeom.h:93
void lwgeom_free(LWGEOM *geom)
Definition: lwgeom.c:1144
void ptarray_free(POINTARRAY *pa)
Definition: ptarray.c:328
#define COMPOUNDTYPE
Definition: liblwgeom.h:92
#define MULTIPOINTTYPE
Definition: liblwgeom.h:87
char lwgeom_geos_errmsg[LWGEOM_GEOS_ERRMSG_MAXSIZE]
GBOX * bbox
Definition: liblwgeom.h:507
#define LWDEBUG(level, msg)
Definition: lwgeom_log.h:83
LWGEOM * lwcollection_make_geos_friendly(LWCOLLECTION *g)
int32_t srid
Definition: liblwgeom.h:423
uint32_t maxgeoms
Definition: liblwgeom.h:510
uint32_t ngeoms
Definition: liblwgeom.h:509
int ptarray_is_closed_2d(const POINTARRAY *pa)
Definition: ptarray.c:695
int32_t srid
Definition: liblwgeom.h:401
uint32_t nrings
Definition: liblwgeom.h:457
unsigned int uint32_t
Definition: uthash.h:78
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 char * lwtype_name(uint8_t type)
Return the type name string associated with a type number (e.g.
Definition: lwutil.c:218
void lwgeom_geos_error(const char *fmt,...)
uint8_t flags
Definition: liblwgeom.h:371
LWGEOM ** geoms
Definition: liblwgeom.h:511
const GEOSGeometry * geom
POINTARRAY ** rings
Definition: liblwgeom.h:459
static GEOSGeometry * LWGEOM_GEOS_makeValidLine(const GEOSGeometry *gin)
#define FLAGS_GET_Z(flags)
Macros for manipulating the &#39;flags&#39; byte.
Definition: liblwgeom.h:139
GEOSGeometry * LWGEOM2GEOS(const LWGEOM *lwgeom, uint8_t autofix)
#define MULTIPOLYGONTYPE
Definition: liblwgeom.h:89
uint8_t MULTITYPE[NUMTYPES]
Look-up for the correct MULTI* type promotion for singleton types.
Definition: lwgeom.c:345
#define MULTISURFACETYPE
Definition: liblwgeom.h:95
POINTARRAY * ring_make_geos_friendly(POINTARRAY *ring)
GEOSGeometry * LWGEOM_GEOS_getPointN(const GEOSGeometry *, uint32_t)
uint8_t * getPoint_internal(const POINTARRAY *pa, uint32_t n)
Definition: ptarray.c:1750
void * lwrealloc(void *mem, size_t size)
Definition: lwutil.c:237
#define POINTTYPE
LWTYPE numbers, used internally by PostGIS.
Definition: liblwgeom.h:84
uint8_t type
Definition: liblwgeom.h:398
LWGEOM * lwpoly_make_geos_friendly(LWPOLY *poly)
void free(void *)
#define CIRCSTRINGTYPE
Definition: liblwgeom.h:91
LWPOINT * lwpoint_construct(int srid, GBOX *bbox, POINTARRAY *point)
Definition: lwpoint.c:129
void * lwalloc(size_t size)
Definition: lwutil.c:229
LWGEOM * lwline_make_geos_friendly(LWLINE *line)
#define MULTILINETYPE
Definition: liblwgeom.h:88
#define LWDEBUGF(level, msg,...)
Definition: lwgeom_log.h:88
#define FLAGS_NDIMS(flags)
Definition: liblwgeom.h:151
static GEOSGeometry * LWGEOM_GEOS_makeValidMultiLine(const GEOSGeometry *gin)
void lwerror(const char *fmt,...)
Write a notice out to the error handler.
Definition: lwutil.c:190
#define COLLECTIONTYPE
Definition: liblwgeom.h:90
This library is the generic geometry handling section of PostGIS.
static GEOSGeometry * LWGEOM_GEOS_nodeLines(const GEOSGeometry *lines)
GEOSGeometry * LWGEOM_GEOS_buildArea(const GEOSGeometry *geom_in)
static GEOSGeometry * LWGEOM_GEOS_makeValidCollection(const GEOSGeometry *gin)
POINTARRAY * points
Definition: liblwgeom.h:424
uint32_t npoints
Definition: liblwgeom.h:373