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