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*
LWGEOM_GEOS_getPointN(const GEOSGeometry * g_in,uint32_t n)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
112 LWGEOM* lwcollection_make_geos_friendly(LWCOLLECTION* g);
113 LWGEOM* lwline_make_geos_friendly(LWLINE* line);
114 LWGEOM* lwpoly_make_geos_friendly(LWPOLY* poly);
115 POINTARRAY* ring_make_geos_friendly(POINTARRAY* ring);
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*
lwgeom_make_geos_friendly(LWGEOM * geom)124 lwgeom_make_geos_friendly(LWGEOM* geom)
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:
148 return lwcollection_make_geos_friendly((LWCOLLECTION*)geom);
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 */
171 POINTARRAY* ptarray_close2d(POINTARRAY* ring);
172 POINTARRAY*
ptarray_close2d(POINTARRAY * ring)173 ptarray_close2d(POINTARRAY* ring)
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*
ring_make_geos_friendly(POINTARRAY * ring)189 ring_make_geos_friendly(POINTARRAY* ring)
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*
lwpoly_make_geos_friendly(LWPOLY * poly)216 lwpoly_make_geos_friendly(LWPOLY* poly)
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*
lwline_make_geos_friendly(LWLINE * line)256 lwline_make_geos_friendly(LWLINE* line)
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*
lwcollection_make_geos_friendly(LWCOLLECTION * g)283 lwcollection_make_geos_friendly(LWCOLLECTION* g)
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*
LWGEOM_GEOS_nodeLines(const GEOSGeometry * lines)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*
LWGEOM_GEOS_makeValidPolygon(const GEOSGeometry * gin)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)),
458 lwgeom_geos_errmsg);
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*
LWGEOM_GEOS_makeValidLine(const GEOSGeometry * gin)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*
LWGEOM_GEOS_makeValidMultiLine(const GEOSGeometry * gin)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;
575 vg = LWGEOM_GEOS_makeValidLine(g);
576 /* Drop any invalid or empty geometry */
577 if (!vg)
578 continue;
579 if (GEOSisEmpty(vg))
580 {
581 GEOSGeom_destroy(vg);
582 continue;
583 }
584
585 if (GEOSGeomTypeId(vg) == GEOS_POINT)
586 points[npoints++] = vg;
587 else if (GEOSGeomTypeId(vg) == GEOS_LINESTRING)
588 lines[nlines++] = vg;
589 else if (GEOSGeomTypeId(vg) == GEOS_MULTILINESTRING)
590 {
591 nsubgeoms = GEOSGetNumGeometries(vg);
592 nlines_alloc += nsubgeoms;
593 lines = lwrealloc(lines, sizeof(GEOSGeometry*) * nlines_alloc);
594 for (j = 0; j < nsubgeoms; ++j)
595 {
596 const GEOSGeometry* gc = GEOSGetGeometryN(vg, j);
597 /* NOTE: ownership of the cloned geoms will be
598 * taken by final collection */
599 lines[nlines++] = GEOSGeom_clone(gc);
600 }
601 }
602 else
603 {
604 /* NOTE: return from GEOSGeomType will leak
605 * but we really don't expect this to happen */
606 lwerror("unexpected geom type returned by LWGEOM_GEOS_makeValid: %s", GEOSGeomType(vg));
607 }
608 }
609
610 if (npoints)
611 {
612 if (npoints > 1)
613 mpoint_out = GEOSGeom_createCollection(GEOS_MULTIPOINT, points, npoints);
614 else
615 mpoint_out = points[0];
616 }
617
618 if (nlines)
619 {
620 if (nlines > 1)
621 mline_out = GEOSGeom_createCollection(GEOS_MULTILINESTRING, lines, nlines);
622 else
623 mline_out = lines[0];
624 }
625
626 lwfree(lines);
627
628 if (mline_out && mpoint_out)
629 {
630 points[0] = mline_out;
631 points[1] = mpoint_out;
632 gout = GEOSGeom_createCollection(GEOS_GEOMETRYCOLLECTION, points, 2);
633 }
634 else if (mline_out)
635 gout = mline_out;
636
637 else if (mpoint_out)
638 gout = mpoint_out;
639
640 lwfree(points);
641
642 return gout;
643 }
644
645 /*
646 * We expect initGEOS being called already.
647 * Will return NULL on error (expect error handler being called by then)
648 */
649 static GEOSGeometry*
LWGEOM_GEOS_makeValidCollection(const GEOSGeometry * gin)650 LWGEOM_GEOS_makeValidCollection(const GEOSGeometry* gin)
651 {
652 int nvgeoms;
653 GEOSGeometry** vgeoms;
654 GEOSGeom gout;
655 int i;
656
657 nvgeoms = GEOSGetNumGeometries(gin);
658 if (nvgeoms == -1)
659 {
660 lwerror("GEOSGetNumGeometries: %s", lwgeom_geos_errmsg);
661 return 0;
662 }
663
664 vgeoms = lwalloc(sizeof(GEOSGeometry*) * nvgeoms);
665 if (!vgeoms)
666 {
667 lwerror("LWGEOM_GEOS_makeValidCollection: out of memory");
668 return 0;
669 }
670
671 for (i = 0; i < nvgeoms; ++i)
672 {
673 vgeoms[i] = LWGEOM_GEOS_makeValid(GEOSGetGeometryN(gin, i));
674 if (!vgeoms[i])
675 {
676 int j;
677 for (j = 0; j < i - 1; j++)
678 GEOSGeom_destroy(vgeoms[j]);
679 lwfree(vgeoms);
680 /* we expect lwerror being called already by makeValid */
681 return NULL;
682 }
683 }
684
685 /* Collect areas and lines (if any line) */
686 gout = GEOSGeom_createCollection(GEOS_GEOMETRYCOLLECTION, vgeoms, nvgeoms);
687 if (!gout) /* an exception again */
688 {
689 /* cleanup and throw */
690 for (i = 0; i < nvgeoms; ++i)
691 GEOSGeom_destroy(vgeoms[i]);
692 lwfree(vgeoms);
693 lwerror("GEOSGeom_createCollection() threw an error: %s", lwgeom_geos_errmsg);
694 return NULL;
695 }
696 lwfree(vgeoms);
697
698 return gout;
699 }
700
701 GEOSGeometry*
LWGEOM_GEOS_makeValid(const GEOSGeometry * gin)702 LWGEOM_GEOS_makeValid(const GEOSGeometry* gin)
703 {
704 GEOSGeometry* gout;
705 char ret_char;
706 #if POSTGIS_DEBUG_LEVEL >= 3
707 LWGEOM *geos_geom;
708 char *geom_ewkt;
709 #endif
710
711 /*
712 * Step 2: return what we got so far if already valid
713 */
714
715 ret_char = GEOSisValid(gin);
716 if (ret_char == 2)
717 {
718 /* I don't think should ever happen */
719 lwerror("GEOSisValid(): %s", lwgeom_geos_errmsg);
720 return NULL;
721 }
722 else if (ret_char)
723 {
724 #if POSTGIS_DEBUG_LEVEL >= 3
725 geos_geom = GEOS2LWGEOM(gin, 0);
726 geom_ewkt = lwgeom_to_ewkt(geos_geom);
727 LWDEBUGF(3, "Geometry [%s] is valid. ", geom_ewkt);
728 lwgeom_free(geos_geom);
729 lwfree(geom_ewkt);
730 #endif
731
732 /* It's valid at this step, return what we have */
733 return GEOSGeom_clone(gin);
734 }
735
736 #if POSTGIS_DEBUG_LEVEL >= 3
737 geos_geom = GEOS2LWGEOM(gin, 0);
738 geom_ewkt = lwgeom_to_ewkt(geos_geom);
739 LWDEBUGF(3,
740 "Geometry [%s] is still not valid: %s. Will try to clean up further.",
741 geom_ewkt,
742 lwgeom_geos_errmsg);
743 lwgeom_free(geos_geom);
744 lwfree(geom_ewkt);
745 #endif
746
747 /*
748 * Step 3 : make what we got valid
749 */
750
751 switch (GEOSGeomTypeId(gin))
752 {
753 case GEOS_MULTIPOINT:
754 case GEOS_POINT:
755 /* points are always valid, but we might have invalid ordinate values */
756 lwnotice("PUNTUAL geometry resulted invalid to GEOS -- dunno how to clean that up");
757 return NULL;
758 break;
759
760 case GEOS_LINESTRING:
761 gout = LWGEOM_GEOS_makeValidLine(gin);
762 if (!gout) /* an exception or something */
763 {
764 /* cleanup and throw */
765 lwerror("%s", lwgeom_geos_errmsg);
766 return NULL;
767 }
768 break; /* we've done */
769
770 case GEOS_MULTILINESTRING:
771 gout = LWGEOM_GEOS_makeValidMultiLine(gin);
772 if (!gout) /* an exception or something */
773 {
774 /* cleanup and throw */
775 lwerror("%s", lwgeom_geos_errmsg);
776 return NULL;
777 }
778 break; /* we've done */
779
780 case GEOS_POLYGON:
781 case GEOS_MULTIPOLYGON:
782 {
783 gout = LWGEOM_GEOS_makeValidPolygon(gin);
784 if (!gout) /* an exception or something */
785 {
786 /* cleanup and throw */
787 lwerror("%s", lwgeom_geos_errmsg);
788 return NULL;
789 }
790 break; /* we've done */
791 }
792
793 case GEOS_GEOMETRYCOLLECTION:
794 {
795 gout = LWGEOM_GEOS_makeValidCollection(gin);
796 if (!gout) /* an exception or something */
797 {
798 /* cleanup and throw */
799 lwerror("%s", lwgeom_geos_errmsg);
800 return NULL;
801 }
802 break; /* we've done */
803 }
804
805 default:
806 {
807 char* typname = GEOSGeomType(gin);
808 lwnotice("ST_MakeValid: doesn't support geometry type: %s", typname);
809 GEOSFree(typname);
810 return NULL;
811 break;
812 }
813 }
814
815 #if PARANOIA_LEVEL > 1
816 /*
817 * Now check if every point of input is also found in output, or abort by returning NULL
818 *
819 * Input geometry was lwgeom_in
820 */
821 {
822 int loss;
823 GEOSGeometry *pi, *po, *pd;
824
825 /* TODO: handle some errors here...
826 * Lack of exceptions is annoying indeed,
827 * I'm getting old --strk;
828 */
829 pi = GEOSGeom_extractUniquePoints(gin);
830 po = GEOSGeom_extractUniquePoints(gout);
831 pd = GEOSDifference(pi, po); /* input points - output points */
832 GEOSGeom_destroy(pi);
833 GEOSGeom_destroy(po);
834 loss = pd && !GEOSisEmpty(pd);
835 GEOSGeom_destroy(pd);
836 if (loss)
837 {
838 lwnotice("%s [%d] Vertices lost in LWGEOM_GEOS_makeValid", __FILE__, __LINE__);
839 /* return NULL */
840 }
841 }
842 #endif /* PARANOIA_LEVEL > 1 */
843
844 return gout;
845 }
846
847 /* Exported. Uses GEOS internally */
848 LWGEOM*
lwgeom_make_valid(LWGEOM * lwgeom_in)849 lwgeom_make_valid(LWGEOM* lwgeom_in)
850 {
851 int is3d;
852 GEOSGeom geosgeom;
853 GEOSGeometry* geosout;
854 LWGEOM* lwgeom_out;
855
856 is3d = FLAGS_GET_Z(lwgeom_in->flags);
857
858 /*
859 * Step 1 : try to convert to GEOS, if impossible, clean that up first
860 * otherwise (adding only duplicates of existing points)
861 */
862
863 initGEOS(lwgeom_geos_error, lwgeom_geos_error);
864
865 lwgeom_out = lwgeom_in;
866 geosgeom = LWGEOM2GEOS(lwgeom_out, 1);
867 if (!geosgeom)
868 {
869 LWDEBUGF(4,
870 "Original geom can't be converted to GEOS (%s)"
871 " - will try cleaning that up first",
872 lwgeom_geos_errmsg);
873
874 lwgeom_out = lwgeom_make_geos_friendly(lwgeom_out);
875 if (!lwgeom_out) lwerror("Could not make a valid geometry out of input");
876
877 /* try again as we did cleanup now */
878 /* TODO: invoke LWGEOM2GEOS directly with autoclean ? */
879 geosgeom = LWGEOM2GEOS(lwgeom_out, 0);
880 if (!geosgeom)
881 {
882 lwerror("Couldn't convert POSTGIS geom to GEOS: %s", lwgeom_geos_errmsg);
883 return NULL;
884 }
885 }
886 else
887 {
888 LWDEBUG(4, "original geom converted to GEOS");
889 lwgeom_out = lwgeom_in;
890 }
891
892 geosout = LWGEOM_GEOS_makeValid(geosgeom);
893 GEOSGeom_destroy(geosgeom);
894 if (!geosout) return NULL;
895
896 lwgeom_out = GEOS2LWGEOM(geosout, is3d);
897 GEOSGeom_destroy(geosout);
898
899 if (lwgeom_is_collection(lwgeom_in) && !lwgeom_is_collection(lwgeom_out))
900 {
901 LWGEOM** ogeoms = lwalloc(sizeof(LWGEOM*));
902 LWGEOM* ogeom;
903 LWDEBUG(3, "lwgeom_make_valid: forcing multi");
904 /* NOTE: this is safe because lwgeom_out is surely not lwgeom_in or
905 * otherwise we couldn't have a collection and a non-collection */
906 assert(lwgeom_in != lwgeom_out);
907 ogeoms[0] = lwgeom_out;
908 ogeom = (LWGEOM*)lwcollection_construct(
909 MULTITYPE[lwgeom_out->type], lwgeom_out->srid, lwgeom_out->bbox, 1, ogeoms);
910 lwgeom_out->bbox = NULL;
911 lwgeom_out = ogeom;
912 }
913
914 lwgeom_out->srid = lwgeom_in->srid;
915 return lwgeom_out;
916 }
917