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