Branch data Line data Source code
1 : : /***************************************************************************
2 : : qgsgeometrymakevalid.cpp
3 : : --------------------------------------
4 : : Date : January 2017
5 : : Copyright : (C) 2017 by Martin Dobias
6 : : Email : wonder dot sk at gmail dot com
7 : : ***************************************************************************
8 : : * *
9 : : * This program is free software; you can redistribute it and/or modify *
10 : : * it under the terms of the GNU General Public License as published by *
11 : : * the Free Software Foundation; either version 2 of the License, or *
12 : : * (at your option) any later version. *
13 : : * *
14 : : ***************************************************************************/
15 : :
16 : : // The code in this file has been ported from lwgeom library (part of PostGIS).
17 : : // See lwgeom_geos_clean.c for the original code by Sandro Santilli.
18 : : //
19 : : // Ideally one day the implementation will go to GEOS library...
20 : :
21 : : #include "qgsgeometry.h"
22 : : #include "qgsgeos.h"
23 : : #include "qgslogger.h"
24 : :
25 : : #include "qgsgeometrycollection.h"
26 : : #include "qgsmultipoint.h"
27 : : #include "qgsmultilinestring.h"
28 : : #include "qgsmultipolygon.h"
29 : : #include "qgspolygon.h"
30 : :
31 : : #include <memory>
32 : :
33 : :
34 : : // ------------ BuildArea stuff ---------------------------------------------------------------------
35 : :
36 : 4 : typedef struct Face_t
37 : : {
38 : 4 : const GEOSGeometry *geom = nullptr;
39 : 4 : GEOSGeometry *env = nullptr;
40 : : double envarea;
41 : : struct Face_t *parent; // if this face is an hole of another one, or NULL
42 : : } Face;
43 : :
44 : : static Face *newFace( const GEOSGeometry *g );
45 : : static void delFace( Face *f );
46 : : static unsigned int countParens( const Face *f );
47 : : static void findFaceHoles( Face **faces, int nfaces );
48 : :
49 : 4 : static Face *newFace( const GEOSGeometry *g )
50 : : {
51 : 4 : GEOSContextHandle_t handle = QgsGeos::getGEOSHandler();
52 : :
53 : 4 : Face *f = new Face;
54 : 4 : f->geom = g;
55 : 4 : f->env = GEOSEnvelope_r( handle, f->geom );
56 : 4 : GEOSArea_r( handle, f->env, &f->envarea );
57 : 4 : f->parent = nullptr;
58 : 4 : return f;
59 : : }
60 : :
61 : 4 : static unsigned int countParens( const Face *f )
62 : : {
63 : 4 : unsigned int pcount = 0;
64 : 4 : while ( f->parent )
65 : : {
66 : 0 : ++pcount;
67 : 0 : f = f->parent;
68 : : }
69 : 4 : return pcount;
70 : : }
71 : :
72 : : // Destroy the face and release memory associated with it
73 : 4 : static void delFace( Face *f )
74 : : {
75 : 4 : GEOSContextHandle_t handle = QgsGeos::getGEOSHandler();
76 : 4 : GEOSGeom_destroy_r( handle, f->env );
77 : 4 : delete f;
78 : 4 : }
79 : :
80 : :
81 : 2 : static int compare_by_envarea( const void *g1, const void *g2 )
82 : : {
83 : 2 : Face *f1 = *( Face ** )g1;
84 : 2 : Face *f2 = *( Face ** )g2;
85 : 2 : double n1 = f1->envarea;
86 : 2 : double n2 = f2->envarea;
87 : :
88 : 2 : if ( n1 < n2 ) return 1;
89 : 2 : if ( n1 > n2 ) return -1;
90 : 2 : return 0;
91 : 2 : }
92 : :
93 : : // Find holes of each face
94 : 2 : static void findFaceHoles( Face **faces, int nfaces )
95 : : {
96 : 2 : GEOSContextHandle_t handle = QgsGeos::getGEOSHandler();
97 : :
98 : : // We sort by envelope area so that we know holes are only
99 : : // after their shells
100 : 2 : qsort( faces, nfaces, sizeof( Face * ), compare_by_envarea );
101 : 6 : for ( int i = 0; i < nfaces; ++i )
102 : : {
103 : 4 : Face *f = faces[i];
104 : 4 : int nholes = GEOSGetNumInteriorRings_r( handle, f->geom );
105 : 4 : for ( int h = 0; h < nholes; ++h )
106 : : {
107 : 0 : const GEOSGeometry *hole = GEOSGetInteriorRingN_r( handle, f->geom, h );
108 : 0 : for ( int j = i + 1; j < nfaces; ++j )
109 : : {
110 : 0 : Face *f2 = faces[j];
111 : 0 : if ( f2->parent )
112 : 0 : continue; // hole already assigned
113 : : /* TODO: can be optimized as the ring would have the
114 : : * same vertices, possibly in different order.
115 : : * maybe comparing number of points could already be
116 : : * useful.
117 : : */
118 : 0 : if ( GEOSEquals_r( handle, GEOSGetExteriorRing_r( handle, f2->geom ), hole ) )
119 : : {
120 : 0 : f2->parent = f;
121 : 0 : break;
122 : : }
123 : 0 : }
124 : 0 : }
125 : 4 : }
126 : 2 : }
127 : :
128 : 2 : static GEOSGeometry *collectFacesWithEvenAncestors( Face **faces, int nfaces )
129 : : {
130 : 2 : GEOSContextHandle_t handle = QgsGeos::getGEOSHandler();
131 : :
132 : 2 : unsigned int ngeoms = 0;
133 : 2 : GEOSGeometry **geoms = new GEOSGeometry*[nfaces];
134 : 6 : for ( int i = 0; i < nfaces; ++i )
135 : : {
136 : 4 : Face *f = faces[i];
137 : 4 : if ( countParens( f ) % 2 )
138 : 0 : continue; // we skip odd parents geoms
139 : 4 : geoms[ngeoms++] = GEOSGeom_clone_r( handle, f->geom );
140 : 4 : }
141 : :
142 : 2 : GEOSGeometry *ret = GEOSGeom_createCollection_r( handle, GEOS_MULTIPOLYGON, geoms, ngeoms );
143 : 2 : delete [] geoms;
144 : 2 : return ret;
145 : : }
146 : :
147 : : Q_NOWARN_UNREACHABLE_PUSH
148 : 6 : static GEOSGeometry *LWGEOM_GEOS_buildArea( const GEOSGeometry *geom_in, QString &errorMessage )
149 : : {
150 : 6 : GEOSContextHandle_t handle = QgsGeos::getGEOSHandler();
151 : :
152 : 6 : GEOSGeometry *tmp = nullptr;
153 : 6 : GEOSGeometry *shp = nullptr;
154 : 6 : GEOSGeometry *geos_result = nullptr;
155 : 6 : int srid = GEOSGetSRID_r( handle, geom_in );
156 : :
157 : : GEOSGeometry const *vgeoms[1];
158 : 6 : vgeoms[0] = geom_in;
159 : : try
160 : : {
161 : 6 : geos_result = GEOSPolygonize_r( handle, vgeoms, 1 );
162 : 6 : }
163 : : catch ( GEOSException &e )
164 : : {
165 : 0 : errorMessage = QStringLiteral( "GEOSPolygonize(): %1" ).arg( e.what() );
166 : 0 : return nullptr;
167 : 0 : }
168 : :
169 : : // We should now have a collection
170 : : #if PARANOIA_LEVEL > 0
171 : : if ( GEOSGeometryTypeId_r( handle, geos_result ) != COLLECTIONTYPE )
172 : : {
173 : : GEOSGeom_destroy_r( handle, geos_result );
174 : : errorMessage = "Unexpected return from GEOSpolygonize";
175 : : return nullptr;
176 : : }
177 : : #endif
178 : :
179 : 6 : int ngeoms = GEOSGetNumGeometries_r( handle, geos_result );
180 : :
181 : : // No geometries in collection, early out
182 : 6 : if ( ngeoms == 0 )
183 : : {
184 : 3 : GEOSSetSRID_r( handle, geos_result, srid );
185 : 3 : return geos_result;
186 : : }
187 : :
188 : : // Return first geometry if we only have one in collection,
189 : : // to avoid the unnecessary Geometry clone below.
190 : 3 : if ( ngeoms == 1 )
191 : : {
192 : 1 : tmp = ( GEOSGeometry * )GEOSGetGeometryN_r( handle, geos_result, 0 );
193 : 1 : if ( ! tmp )
194 : : {
195 : 0 : GEOSGeom_destroy_r( handle, geos_result );
196 : 0 : return nullptr;
197 : : }
198 : 1 : shp = GEOSGeom_clone_r( handle, tmp );
199 : 1 : GEOSGeom_destroy_r( handle, geos_result ); // only safe after the clone above
200 : 1 : GEOSSetSRID_r( handle, shp, srid );
201 : 1 : return shp;
202 : : }
203 : :
204 : : /*
205 : : * Polygonizer returns a polygon for each face in the built topology.
206 : : *
207 : : * This means that for any face with holes we'll have other faces
208 : : * representing each hole. We can imagine a parent-child relationship
209 : : * between these faces.
210 : : *
211 : : * In order to maximize the number of visible rings in output we
212 : : * only use those faces which have an even number of parents.
213 : : *
214 : : * Example:
215 : : *
216 : : * +---------------+
217 : : * | L0 | L0 has no parents
218 : : * | +---------+ |
219 : : * | | L1 | | L1 is an hole of L0
220 : : * | | +---+ | |
221 : : * | | |L2 | | | L2 is an hole of L1 (which is an hole of L0)
222 : : * | | | | | |
223 : : * | | +---+ | |
224 : : * | +---------+ |
225 : : * | |
226 : : * +---------------+
227 : : *
228 : : * See http://trac.osgeo.org/postgis/ticket/1806
229 : : *
230 : : */
231 : :
232 : : // Prepare face structures for later analysis
233 : 2 : Face **geoms = new Face*[ngeoms];
234 : 6 : for ( int i = 0; i < ngeoms; ++i )
235 : 4 : geoms[i] = newFace( GEOSGetGeometryN_r( handle, geos_result, i ) );
236 : :
237 : : // Find faces representing other faces holes
238 : 2 : findFaceHoles( geoms, ngeoms );
239 : :
240 : : // Build a MultiPolygon composed only by faces with an
241 : : // even number of ancestors
242 : 2 : tmp = collectFacesWithEvenAncestors( geoms, ngeoms );
243 : :
244 : : // Cleanup face structures
245 : 6 : for ( int i = 0; i < ngeoms; ++i )
246 : 4 : delFace( geoms[i] );
247 : 2 : delete [] geoms;
248 : :
249 : : // Faces referenced memory owned by geos_result.
250 : : // It is safe to destroy geos_result after deleting them.
251 : 2 : GEOSGeom_destroy_r( handle, geos_result );
252 : :
253 : : // Run a single overlay operation to dissolve shared edges
254 : 2 : shp = GEOSUnionCascaded_r( handle, tmp );
255 : 2 : if ( !shp )
256 : : {
257 : 0 : GEOSGeom_destroy_r( handle, tmp );
258 : 0 : return nullptr;
259 : : }
260 : :
261 : 2 : GEOSGeom_destroy_r( handle, tmp );
262 : :
263 : 2 : GEOSSetSRID_r( handle, shp, srid );
264 : :
265 : 2 : return shp;
266 : 6 : }
267 : : Q_NOWARN_UNREACHABLE_POP
268 : :
269 : : // Return Nth vertex in GEOSGeometry as a POINT.
270 : : // May return NULL if the geometry has NO vertexex.
271 : 6 : static GEOSGeometry *LWGEOM_GEOS_getPointN( const GEOSGeometry *g_in, uint32_t n )
272 : : {
273 : 6 : GEOSContextHandle_t handle = QgsGeos::getGEOSHandler();
274 : :
275 : : uint32_t dims;
276 : 6 : const GEOSCoordSequence *seq_in = nullptr;
277 : : GEOSCoordSeq seq_out;
278 : : double val;
279 : : uint32_t sz;
280 : : int gn;
281 : 6 : GEOSGeometry *ret = nullptr;
282 : :
283 : 6 : switch ( GEOSGeomTypeId_r( handle, g_in ) )
284 : : {
285 : : case GEOS_MULTIPOINT:
286 : : case GEOS_MULTILINESTRING:
287 : : case GEOS_MULTIPOLYGON:
288 : : case GEOS_GEOMETRYCOLLECTION:
289 : : {
290 : 0 : for ( gn = 0; gn < GEOSGetNumGeometries_r( handle, g_in ); ++gn )
291 : : {
292 : 0 : const GEOSGeometry *g = GEOSGetGeometryN_r( handle, g_in, gn );
293 : 0 : ret = LWGEOM_GEOS_getPointN( g, n );
294 : 0 : if ( ret ) return ret;
295 : 0 : }
296 : 0 : break;
297 : : }
298 : :
299 : : case GEOS_POLYGON:
300 : : {
301 : 0 : ret = LWGEOM_GEOS_getPointN( GEOSGetExteriorRing_r( handle, g_in ), n );
302 : 0 : if ( ret ) return ret;
303 : 0 : for ( gn = 0; gn < GEOSGetNumInteriorRings_r( handle, g_in ); ++gn )
304 : : {
305 : 0 : const GEOSGeometry *g = GEOSGetInteriorRingN_r( handle, g_in, gn );
306 : 0 : ret = LWGEOM_GEOS_getPointN( g, n );
307 : 0 : if ( ret ) return ret;
308 : 0 : }
309 : 0 : break;
310 : : }
311 : :
312 : : case GEOS_POINT:
313 : : case GEOS_LINESTRING:
314 : : case GEOS_LINEARRING:
315 : 6 : break;
316 : :
317 : : }
318 : :
319 : 6 : seq_in = GEOSGeom_getCoordSeq_r( handle, g_in );
320 : 6 : if ( ! seq_in ) return nullptr;
321 : 6 : if ( ! GEOSCoordSeq_getSize_r( handle, seq_in, &sz ) ) return nullptr;
322 : 6 : if ( ! sz ) return nullptr;
323 : :
324 : 6 : if ( ! GEOSCoordSeq_getDimensions_r( handle, seq_in, &dims ) ) return nullptr;
325 : :
326 : 6 : seq_out = GEOSCoordSeq_create_r( handle, 1, dims );
327 : 6 : if ( ! seq_out ) return nullptr;
328 : :
329 : 6 : if ( ! GEOSCoordSeq_getX_r( handle, seq_in, n, &val ) ) return nullptr;
330 : 6 : if ( ! GEOSCoordSeq_setX_r( handle, seq_out, n, val ) ) return nullptr;
331 : 6 : if ( ! GEOSCoordSeq_getY_r( handle, seq_in, n, &val ) ) return nullptr;
332 : 6 : if ( ! GEOSCoordSeq_setY_r( handle, seq_out, n, val ) ) return nullptr;
333 : 6 : if ( dims > 2 )
334 : : {
335 : 0 : if ( ! GEOSCoordSeq_getZ_r( handle, seq_in, n, &val ) ) return nullptr;
336 : 0 : if ( ! GEOSCoordSeq_setZ_r( handle, seq_out, n, val ) ) return nullptr;
337 : 0 : }
338 : :
339 : 6 : return GEOSGeom_createPoint_r( handle, seq_out );
340 : 6 : }
341 : :
342 : :
343 : : static bool lwline_make_geos_friendly( QgsLineString *line );
344 : : static bool lwpoly_make_geos_friendly( QgsPolygon *poly );
345 : : static bool lwcollection_make_geos_friendly( QgsGeometryCollection *g );
346 : :
347 : :
348 : : // Ensure the geometry is "structurally" valid (enough for GEOS to accept it)
349 : 1 : static bool lwgeom_make_geos_friendly( QgsAbstractGeometry *geom )
350 : : {
351 : 1 : QgsDebugMsgLevel( QStringLiteral( "lwgeom_make_geos_friendly enter (type %1)" ).arg( geom->wkbType() ), 3 );
352 : 1 : switch ( QgsWkbTypes::flatType( geom->wkbType() ) )
353 : : {
354 : : case QgsWkbTypes::Point:
355 : : case QgsWkbTypes::MultiPoint:
356 : : // a point is always valid
357 : 0 : return true;
358 : : break;
359 : :
360 : : case QgsWkbTypes::LineString:
361 : : // lines need at least 2 points
362 : 1 : return lwline_make_geos_friendly( qgsgeometry_cast<QgsLineString *>( geom ) );
363 : : break;
364 : :
365 : : case QgsWkbTypes::Polygon:
366 : : // polygons need all rings closed and with npoints > 3
367 : 0 : return lwpoly_make_geos_friendly( qgsgeometry_cast<QgsPolygon *>( geom ) );
368 : : break;
369 : :
370 : : case QgsWkbTypes::MultiLineString:
371 : : case QgsWkbTypes::MultiPolygon:
372 : : case QgsWkbTypes::GeometryCollection:
373 : 0 : return lwcollection_make_geos_friendly( qgsgeometry_cast<QgsGeometryCollection *>( geom ) );
374 : : break;
375 : :
376 : : default:
377 : 0 : QgsDebugMsg( QStringLiteral( "lwgeom_make_geos_friendly: unsupported input geometry type: %1" ).arg( geom->wkbType() ) );
378 : 0 : break;
379 : : }
380 : 0 : return false;
381 : 1 : }
382 : :
383 : :
384 : 0 : static bool ring_make_geos_friendly( QgsCurve *ring )
385 : : {
386 : 0 : if ( ring->nCoordinates() == 0 )
387 : 0 : return false;
388 : :
389 : : // earlier we allowed in only geometries with straight segments
390 : 0 : QgsLineString *linestring = qgsgeometry_cast<QgsLineString *>( ring );
391 : : Q_ASSERT_X( linestring, "ring_make_geos_friendly", "ring was not a linestring" );
392 : :
393 : : // close the ring if not already closed (2d only)
394 : :
395 : 0 : QgsPoint p1 = linestring->startPoint(), p2 = linestring->endPoint();
396 : 0 : if ( p1.x() != p2.x() || p1.y() != p2.y() )
397 : 0 : linestring->addVertex( p1 );
398 : :
399 : : // must have at least 4 coordinates to be accepted by GEOS
400 : :
401 : 0 : while ( linestring->nCoordinates() < 4 )
402 : 0 : linestring->addVertex( p1 );
403 : :
404 : 0 : return true;
405 : 0 : }
406 : :
407 : : // Make sure all rings are closed and have > 3 points.
408 : 0 : static bool lwpoly_make_geos_friendly( QgsPolygon *poly )
409 : : {
410 : : // If the polygon has no rings there's nothing to do
411 : : // TODO: in qgis representation there always is exterior ring
412 : : //if ( ! poly->nrings ) return true;
413 : :
414 : : // All rings must be closed and have > 3 points
415 : 0 : for ( int i = 0; i < poly->numInteriorRings(); i++ )
416 : : {
417 : 0 : if ( !ring_make_geos_friendly( qgsgeometry_cast<QgsCurve *>( poly->interiorRing( i ) ) ) )
418 : 0 : return false;
419 : 0 : }
420 : :
421 : 0 : return true;
422 : 0 : }
423 : :
424 : : // Need NO or >1 points. Duplicate first if only one.
425 : 1 : static bool lwline_make_geos_friendly( QgsLineString *line )
426 : : {
427 : 1 : if ( line->numPoints() == 1 ) // 0 is fine, 2 is fine
428 : : {
429 : 1 : line->addVertex( line->startPoint() );
430 : 1 : }
431 : 1 : return true;
432 : 0 : }
433 : :
434 : 0 : static bool lwcollection_make_geos_friendly( QgsGeometryCollection *g )
435 : : {
436 : 0 : for ( int i = 0; i < g->numGeometries(); i++ )
437 : : {
438 : 0 : if ( !lwgeom_make_geos_friendly( g->geometryN( i ) ) )
439 : 0 : return false;
440 : 0 : }
441 : :
442 : 0 : return true;
443 : 0 : }
444 : :
445 : :
446 : : // Fully node given linework
447 : 6 : static GEOSGeometry *LWGEOM_GEOS_nodeLines( const GEOSGeometry *lines )
448 : : {
449 : 6 : GEOSContextHandle_t handle = QgsGeos::getGEOSHandler();
450 : :
451 : : // Union with first geometry point, obtaining full noding
452 : : // and dissolving of duplicated repeated points
453 : : //
454 : : // TODO: substitute this with UnaryUnion?
455 : :
456 : 6 : GEOSGeometry *point = LWGEOM_GEOS_getPointN( lines, 0 );
457 : 6 : if ( ! point )
458 : 0 : return nullptr;
459 : :
460 : 6 : GEOSGeometry *noded = nullptr;
461 : : try
462 : : {
463 : 6 : noded = GEOSUnion_r( handle, lines, point );
464 : 6 : }
465 : : catch ( GEOSException & )
466 : : {
467 : : // no need to do anything here - we'll return nullptr anyway
468 : 0 : }
469 : 6 : GEOSGeom_destroy_r( handle, point );
470 : 6 : return noded;
471 : 6 : }
472 : :
473 : : // Will return NULL on error (expect error handler being called by then)
474 : : Q_NOWARN_UNREACHABLE_PUSH
475 : 3 : static GEOSGeometry *LWGEOM_GEOS_makeValidPolygon( const GEOSGeometry *gin, QString &errorMessage )
476 : : {
477 : 3 : GEOSContextHandle_t handle = QgsGeos::getGEOSHandler();
478 : :
479 : : GEOSGeom gout;
480 : : GEOSGeom geos_bound;
481 : : GEOSGeom geos_cut_edges, geos_area, collapse_points;
482 : : GEOSGeometry *vgeoms[3]; // One for area, one for cut-edges
483 : 3 : unsigned int nvgeoms = 0;
484 : :
485 : : Q_ASSERT( GEOSGeomTypeId_r( handle, gin ) == GEOS_POLYGON ||
486 : : GEOSGeomTypeId_r( handle, gin ) == GEOS_MULTIPOLYGON );
487 : :
488 : 3 : geos_bound = GEOSBoundary_r( handle, gin );
489 : 3 : if ( !geos_bound )
490 : 0 : return nullptr;
491 : :
492 : : // Use noded boundaries as initial "cut" edges
493 : :
494 : 3 : geos_cut_edges = LWGEOM_GEOS_nodeLines( geos_bound );
495 : 3 : if ( !geos_cut_edges )
496 : : {
497 : 0 : GEOSGeom_destroy_r( handle, geos_bound );
498 : 0 : errorMessage = QStringLiteral( "LWGEOM_GEOS_nodeLines() failed" );
499 : 0 : return nullptr;
500 : : }
501 : :
502 : : // NOTE: the noding process may drop lines collapsing to points.
503 : : // We want to retrieve any of those
504 : : {
505 : 3 : GEOSGeometry *pi = nullptr;
506 : 3 : GEOSGeometry *po = nullptr;
507 : :
508 : : try
509 : : {
510 : 3 : pi = GEOSGeom_extractUniquePoints_r( handle, geos_bound );
511 : 3 : }
512 : : catch ( GEOSException &e )
513 : : {
514 : 0 : GEOSGeom_destroy_r( handle, geos_bound );
515 : 0 : errorMessage = QStringLiteral( "GEOSGeom_extractUniquePoints(): %1" ).arg( e.what() );
516 : 0 : return nullptr;
517 : 0 : }
518 : :
519 : : try
520 : : {
521 : 3 : po = GEOSGeom_extractUniquePoints_r( handle, geos_cut_edges );
522 : 3 : }
523 : : catch ( GEOSException &e )
524 : : {
525 : 0 : GEOSGeom_destroy_r( handle, geos_bound );
526 : 0 : GEOSGeom_destroy_r( handle, pi );
527 : 0 : errorMessage = QStringLiteral( "GEOSGeom_extractUniquePoints(): %1" ).arg( e.what() );
528 : 0 : return nullptr;
529 : 0 : }
530 : :
531 : : try
532 : : {
533 : 3 : collapse_points = GEOSDifference_r( handle, pi, po );
534 : 3 : }
535 : : catch ( GEOSException &e )
536 : : {
537 : 0 : GEOSGeom_destroy_r( handle, geos_bound );
538 : 0 : GEOSGeom_destroy_r( handle, pi );
539 : 0 : GEOSGeom_destroy_r( handle, po );
540 : 0 : errorMessage = QStringLiteral( "GEOSDifference(): %1" ).arg( e.what() );
541 : 0 : return nullptr;
542 : 0 : }
543 : :
544 : 3 : GEOSGeom_destroy_r( handle, pi );
545 : 3 : GEOSGeom_destroy_r( handle, po );
546 : : }
547 : 3 : GEOSGeom_destroy_r( handle, geos_bound );
548 : :
549 : : // And use an empty geometry as initial "area"
550 : : try
551 : : {
552 : 3 : geos_area = GEOSGeom_createEmptyPolygon_r( handle );
553 : 3 : }
554 : : catch ( GEOSException &e )
555 : : {
556 : 0 : errorMessage = QStringLiteral( "GEOSGeom_createEmptyPolygon(): %1" ).arg( e.what() );
557 : 0 : GEOSGeom_destroy_r( handle, geos_cut_edges );
558 : 0 : return nullptr;
559 : 0 : }
560 : :
561 : : // See if an area can be build with the remaining edges
562 : : // and if it can, symdifference with the original area.
563 : : // Iterate this until no more polygons can be created
564 : : // with left-over edges.
565 : 6 : while ( GEOSGetNumGeometries_r( handle, geos_cut_edges ) )
566 : : {
567 : 6 : GEOSGeometry *new_area = nullptr;
568 : 6 : GEOSGeometry *new_area_bound = nullptr;
569 : 6 : GEOSGeometry *symdif = nullptr;
570 : 6 : GEOSGeometry *new_cut_edges = nullptr;
571 : :
572 : : // ASSUMPTION: cut_edges should already be fully noded
573 : :
574 : : try
575 : : {
576 : 6 : new_area = LWGEOM_GEOS_buildArea( geos_cut_edges, errorMessage );
577 : 6 : }
578 : : catch ( GEOSException &e )
579 : : {
580 : 0 : GEOSGeom_destroy_r( handle, geos_cut_edges );
581 : 0 : GEOSGeom_destroy_r( handle, geos_area );
582 : 0 : errorMessage = QStringLiteral( "LWGEOM_GEOS_buildArea() threw an error: %1" ).arg( e.what() );
583 : 0 : return nullptr;
584 : 0 : }
585 : :
586 : 6 : if ( GEOSisEmpty_r( handle, new_area ) )
587 : : {
588 : : // no more rings can be build with the edges
589 : 3 : GEOSGeom_destroy_r( handle, new_area );
590 : 3 : break;
591 : : }
592 : :
593 : : // We succeeded in building a ring !
594 : :
595 : : // Save the new ring boundaries first (to compute
596 : : // further cut edges later)
597 : : try
598 : : {
599 : 3 : new_area_bound = GEOSBoundary_r( handle, new_area );
600 : 3 : }
601 : : catch ( GEOSException &e )
602 : : {
603 : : // We did check for empty area already so
604 : : // this must be some other error
605 : 0 : errorMessage = QStringLiteral( "GEOSBoundary() threw an error: %1" ).arg( e.what() );
606 : 0 : GEOSGeom_destroy_r( handle, new_area );
607 : 0 : GEOSGeom_destroy_r( handle, geos_area );
608 : 0 : return nullptr;
609 : 0 : }
610 : :
611 : : // Now symdiff new and old area
612 : : try
613 : : {
614 : 3 : symdif = GEOSSymDifference_r( handle, geos_area, new_area );
615 : 3 : }
616 : : catch ( GEOSException &e )
617 : : {
618 : 0 : GEOSGeom_destroy_r( handle, geos_cut_edges );
619 : 0 : GEOSGeom_destroy_r( handle, new_area );
620 : 0 : GEOSGeom_destroy_r( handle, new_area_bound );
621 : 0 : GEOSGeom_destroy_r( handle, geos_area );
622 : 0 : errorMessage = QStringLiteral( "GEOSSymDifference() threw an error: %1" ).arg( e.what() );
623 : 0 : return nullptr;
624 : 0 : }
625 : :
626 : 3 : GEOSGeom_destroy_r( handle, geos_area );
627 : 3 : GEOSGeom_destroy_r( handle, new_area );
628 : 3 : geos_area = symdif;
629 : 3 : symdif = nullptr;
630 : :
631 : : // Now let's re-set geos_cut_edges with what's left
632 : : // from the original boundary.
633 : : // ASSUMPTION: only the previous cut-edges can be
634 : : // left, so we don't need to reconsider
635 : : // the whole original boundaries
636 : : //
637 : : // NOTE: this is an expensive operation.
638 : :
639 : : try
640 : : {
641 : 3 : new_cut_edges = GEOSDifference_r( handle, geos_cut_edges, new_area_bound );
642 : 3 : }
643 : : catch ( GEOSException &e )
644 : : {
645 : 0 : GEOSGeom_destroy_r( handle, geos_cut_edges );
646 : 0 : GEOSGeom_destroy_r( handle, new_area_bound );
647 : 0 : GEOSGeom_destroy_r( handle, geos_area );
648 : 0 : errorMessage = QStringLiteral( "GEOSDifference() threw an error: %1" ).arg( e.what() );
649 : 0 : return nullptr;
650 : 0 : }
651 : 3 : GEOSGeom_destroy_r( handle, geos_cut_edges );
652 : 3 : GEOSGeom_destroy_r( handle, new_area_bound );
653 : 3 : geos_cut_edges = new_cut_edges;
654 : : }
655 : :
656 : 3 : if ( !GEOSisEmpty_r( handle, geos_area ) )
657 : : {
658 : 3 : vgeoms[nvgeoms++] = geos_area;
659 : 3 : }
660 : : else
661 : : {
662 : 0 : GEOSGeom_destroy_r( handle, geos_area );
663 : : }
664 : :
665 : 3 : if ( !GEOSisEmpty_r( handle, geos_cut_edges ) )
666 : : {
667 : 1 : vgeoms[nvgeoms++] = geos_cut_edges;
668 : 1 : }
669 : : else
670 : : {
671 : 2 : GEOSGeom_destroy_r( handle, geos_cut_edges );
672 : : }
673 : :
674 : 3 : if ( !GEOSisEmpty_r( handle, collapse_points ) )
675 : : {
676 : 0 : vgeoms[nvgeoms++] = collapse_points;
677 : 0 : }
678 : : else
679 : : {
680 : 3 : GEOSGeom_destroy_r( handle, collapse_points );
681 : : }
682 : :
683 : 3 : if ( 1 == nvgeoms )
684 : : {
685 : : // Return cut edges
686 : 2 : gout = vgeoms[0];
687 : 2 : }
688 : : else
689 : : {
690 : : // Collect areas and lines (if any line)
691 : : try
692 : : {
693 : 1 : gout = GEOSGeom_createCollection_r( handle, GEOS_GEOMETRYCOLLECTION, vgeoms, nvgeoms );
694 : 1 : }
695 : : catch ( GEOSException &e )
696 : : {
697 : 0 : errorMessage = QStringLiteral( "GEOSGeom_createCollection() threw an error: %1" ).arg( e.what() );
698 : : // TODO: cleanup!
699 : 0 : return nullptr;
700 : 0 : }
701 : : }
702 : :
703 : 3 : return gout;
704 : :
705 : 3 : }
706 : : Q_NOWARN_UNREACHABLE_PUSH
707 : :
708 : 3 : static GEOSGeometry *LWGEOM_GEOS_makeValidLine( const GEOSGeometry *gin, QString &errorMessage )
709 : : {
710 : 3 : Q_UNUSED( errorMessage )
711 : 3 : return LWGEOM_GEOS_nodeLines( gin );
712 : : }
713 : :
714 : 1 : static GEOSGeometry *LWGEOM_GEOS_makeValidMultiLine( const GEOSGeometry *gin, QString &errorMessage )
715 : : {
716 : 1 : GEOSContextHandle_t handle = QgsGeos::getGEOSHandler();
717 : :
718 : 1 : int ngeoms = GEOSGetNumGeometries_r( handle, gin );
719 : 1 : uint32_t nlines_alloc = ngeoms;
720 : 1 : QVector<GEOSGeometry *> lines;
721 : 1 : QVector<GEOSGeometry *> points;
722 : 1 : lines.reserve( nlines_alloc );
723 : 1 : points.reserve( ngeoms );
724 : :
725 : 2 : for ( int i = 0; i < ngeoms; ++i )
726 : : {
727 : 1 : const GEOSGeometry *g = GEOSGetGeometryN_r( handle, gin, i );
728 : 1 : GEOSGeometry *vg = LWGEOM_GEOS_makeValidLine( g, errorMessage );
729 : 1 : if ( GEOSisEmpty_r( handle, vg ) )
730 : : {
731 : : // we don't care about this one
732 : 0 : GEOSGeom_destroy_r( handle, vg );
733 : 0 : }
734 : 1 : if ( GEOSGeomTypeId_r( handle, vg ) == GEOS_POINT )
735 : : {
736 : 1 : points.append( vg );
737 : 1 : }
738 : 0 : else if ( GEOSGeomTypeId_r( handle, vg ) == GEOS_LINESTRING )
739 : : {
740 : 0 : lines.append( vg );
741 : 0 : }
742 : 0 : else if ( GEOSGeomTypeId_r( handle, vg ) == GEOS_MULTILINESTRING )
743 : : {
744 : 0 : int nsubgeoms = GEOSGetNumGeometries_r( handle, vg );
745 : 0 : nlines_alloc += nsubgeoms;
746 : 0 : lines.reserve( nlines_alloc );
747 : 0 : for ( int j = 0; j < nsubgeoms; ++j )
748 : : {
749 : : // NOTE: ownership of the cloned geoms will be taken by final collection
750 : 0 : lines.append( GEOSGeom_clone_r( handle, GEOSGetGeometryN_r( handle, vg, j ) ) );
751 : 0 : }
752 : 0 : }
753 : : else
754 : : {
755 : : // NOTE: return from GEOSGeomType will leak but we really don't expect this to happen
756 : 0 : errorMessage = QStringLiteral( "unexpected geom type returned by LWGEOM_GEOS_makeValid: %1" ).arg( GEOSGeomTypeId_r( handle, vg ) );
757 : : }
758 : 1 : }
759 : :
760 : 1 : GEOSGeometry *mpoint_out = nullptr;
761 : 1 : if ( !points.isEmpty() )
762 : : {
763 : 1 : if ( points.count() > 1 )
764 : : {
765 : 0 : mpoint_out = GEOSGeom_createCollection_r( handle, GEOS_MULTIPOINT, points.data(), points.count() );
766 : 0 : }
767 : : else
768 : : {
769 : 1 : mpoint_out = points[0];
770 : : }
771 : 1 : }
772 : :
773 : 1 : GEOSGeometry *mline_out = nullptr;
774 : 1 : if ( !lines.isEmpty() )
775 : : {
776 : 0 : if ( lines.count() > 1 )
777 : : {
778 : 0 : mline_out = GEOSGeom_createCollection_r( handle, GEOS_MULTILINESTRING, lines.data(), lines.count() );
779 : 0 : }
780 : : else
781 : : {
782 : 0 : mline_out = lines[0];
783 : : }
784 : 0 : }
785 : :
786 : 1 : if ( mline_out && mpoint_out )
787 : : {
788 : : GEOSGeometry *collection[2];
789 : 0 : collection[0] = mline_out;
790 : 0 : collection[1] = mpoint_out;
791 : 0 : return GEOSGeom_createCollection_r( handle, GEOS_GEOMETRYCOLLECTION, collection, 2 );
792 : : }
793 : 1 : else if ( mline_out )
794 : : {
795 : 0 : return mline_out;
796 : : }
797 : 1 : else if ( mpoint_out )
798 : : {
799 : 1 : return mpoint_out;
800 : : }
801 : : else
802 : : {
803 : 0 : return nullptr;
804 : : }
805 : 1 : }
806 : :
807 : :
808 : : static GEOSGeometry *LWGEOM_GEOS_makeValid( const GEOSGeometry *gin, QString &errorMessage );
809 : :
810 : : // Will return NULL on error (expect error handler being called by then)
811 : 1 : static GEOSGeometry *LWGEOM_GEOS_makeValidCollection( const GEOSGeometry *gin, QString &errorMessage )
812 : : {
813 : 1 : GEOSContextHandle_t handle = QgsGeos::getGEOSHandler();
814 : :
815 : 1 : int nvgeoms = GEOSGetNumGeometries_r( handle, gin );
816 : 1 : if ( nvgeoms == -1 )
817 : : {
818 : 0 : errorMessage = QStringLiteral( "GEOSGetNumGeometries: %1" ).arg( QLatin1String( "?" ) );
819 : 0 : return nullptr;
820 : : }
821 : :
822 : 1 : QVector<GEOSGeometry *> vgeoms( nvgeoms );
823 : 4 : for ( int i = 0; i < nvgeoms; ++i )
824 : : {
825 : 3 : vgeoms[i] = LWGEOM_GEOS_makeValid( GEOSGetGeometryN_r( handle, gin, i ), errorMessage );
826 : 3 : if ( ! vgeoms[i] )
827 : : {
828 : 0 : while ( i-- )
829 : 0 : GEOSGeom_destroy_r( handle, vgeoms[i] );
830 : : // we expect lwerror being called already by makeValid
831 : 0 : return nullptr;
832 : : }
833 : 3 : }
834 : :
835 : : // Collect areas and lines (if any line)
836 : : try
837 : : {
838 : 1 : return GEOSGeom_createCollection_r( handle, GEOS_GEOMETRYCOLLECTION, vgeoms.data(), nvgeoms );
839 : 0 : }
840 : : catch ( GEOSException &e )
841 : : {
842 : : // cleanup and throw
843 : 0 : for ( int i = 0; i < nvgeoms; ++i )
844 : 0 : GEOSGeom_destroy_r( handle, vgeoms[i] );
845 : 0 : errorMessage = QStringLiteral( "GEOSGeom_createCollection() threw an error: %1" ).arg( e.what() );
846 : 0 : return nullptr;
847 : 0 : }
848 : 1 : }
849 : :
850 : :
851 : 9 : static GEOSGeometry *LWGEOM_GEOS_makeValid( const GEOSGeometry *gin, QString &errorMessage )
852 : : {
853 : 9 : GEOSContextHandle_t handle = QgsGeos::getGEOSHandler();
854 : :
855 : : // return what we got so far if already valid
856 : : Q_NOWARN_UNREACHABLE_PUSH
857 : : try
858 : : {
859 : 9 : if ( GEOSisValid_r( handle, gin ) )
860 : : {
861 : : // It's valid at this step, return what we have
862 : 2 : return GEOSGeom_clone_r( handle, gin );
863 : : }
864 : 7 : }
865 : : catch ( GEOSException &e )
866 : : {
867 : : // I don't think should ever happen
868 : 0 : errorMessage = QStringLiteral( "GEOSisValid(): %1" ).arg( e.what() );
869 : 0 : return nullptr;
870 : 0 : }
871 : : Q_NOWARN_UNREACHABLE_POP
872 : :
873 : : // make what we got valid
874 : :
875 : 7 : switch ( GEOSGeomTypeId_r( handle, gin ) )
876 : : {
877 : : case GEOS_MULTIPOINT:
878 : : case GEOS_POINT:
879 : : // points are always valid, but we might have invalid ordinate values
880 : 0 : QgsDebugMsg( QStringLiteral( "PUNTUAL geometry resulted invalid to GEOS -- dunno how to clean that up" ) );
881 : 0 : return nullptr;
882 : :
883 : : case GEOS_LINESTRING:
884 : 2 : return LWGEOM_GEOS_makeValidLine( gin, errorMessage );
885 : :
886 : : case GEOS_MULTILINESTRING:
887 : 1 : return LWGEOM_GEOS_makeValidMultiLine( gin, errorMessage );
888 : :
889 : : case GEOS_POLYGON:
890 : : case GEOS_MULTIPOLYGON:
891 : 3 : return LWGEOM_GEOS_makeValidPolygon( gin, errorMessage );
892 : :
893 : : case GEOS_GEOMETRYCOLLECTION:
894 : 1 : return LWGEOM_GEOS_makeValidCollection( gin, errorMessage );
895 : :
896 : : default:
897 : 0 : errorMessage = QStringLiteral( "ST_MakeValid: doesn't support geometry type: %1" ).arg( GEOSGeomTypeId_r( handle, gin ) );
898 : 0 : return nullptr;
899 : : }
900 : 9 : }
901 : :
902 : :
903 : 6 : std::unique_ptr< QgsAbstractGeometry > _qgis_lwgeom_make_valid( const QgsAbstractGeometry *lwgeom_in, QString &errorMessage )
904 : : {
905 : : //bool is3d = FLAGS_GET_Z(lwgeom_in->flags);
906 : :
907 : : // try to convert to GEOS, if impossible, clean that up first
908 : : // otherwise (adding only duplicates of existing points)
909 : 6 : GEOSContextHandle_t handle = QgsGeos::getGEOSHandler();
910 : :
911 : 6 : geos::unique_ptr geosgeom = QgsGeos::asGeos( lwgeom_in );
912 : 6 : if ( !geosgeom )
913 : : {
914 : 1 : QgsDebugMsgLevel( QStringLiteral( "Original geom can't be converted to GEOS - will try cleaning that up first" ), 3 );
915 : :
916 : 1 : std::unique_ptr<QgsAbstractGeometry> lwgeom_in_clone( lwgeom_in->clone() );
917 : 1 : if ( !lwgeom_make_geos_friendly( lwgeom_in_clone.get() ) )
918 : : {
919 : 0 : QgsDebugMsg( QStringLiteral( "Could not make a valid geometry out of input" ) );
920 : 0 : }
921 : :
922 : : // try again as we did cleanup now
923 : : // TODO: invoke LWGEOM2GEOS directly with autoclean ?
924 : 1 : geosgeom = QgsGeos::asGeos( lwgeom_in_clone.get() );
925 : :
926 : 1 : if ( ! geosgeom )
927 : : {
928 : 0 : errorMessage = QStringLiteral( "Could not convert QGIS geom to GEOS" );
929 : 0 : return nullptr;
930 : : }
931 : 1 : }
932 : : else
933 : : {
934 : 5 : QgsDebugMsgLevel( QStringLiteral( "original geom converted to GEOS" ), 4 );
935 : : }
936 : :
937 : 6 : GEOSGeometry *geosout = LWGEOM_GEOS_makeValid( geosgeom.get(), errorMessage );
938 : 6 : if ( !geosout )
939 : 0 : return nullptr;
940 : :
941 : 6 : std::unique_ptr< QgsAbstractGeometry > lwgeom_out = QgsGeos::fromGeos( geosout );
942 : 6 : GEOSGeom_destroy_r( handle, geosout );
943 : 6 : if ( !lwgeom_out )
944 : 0 : return nullptr;
945 : :
946 : : // force multi-type if we had a multi-type before
947 : 6 : if ( QgsWkbTypes::isMultiType( lwgeom_in->wkbType() ) && !QgsWkbTypes::isMultiType( lwgeom_out->wkbType() ) )
948 : : {
949 : 1 : QgsGeometryCollection *collection = nullptr;
950 : 1 : switch ( QgsWkbTypes::multiType( lwgeom_out->wkbType() ) )
951 : : {
952 : : case QgsWkbTypes::MultiPoint:
953 : 1 : collection = new QgsMultiPoint();
954 : 1 : break;
955 : : case QgsWkbTypes::MultiLineString:
956 : 0 : collection = new QgsMultiLineString();
957 : 0 : break;
958 : : case QgsWkbTypes::MultiPolygon:
959 : 0 : collection = new QgsMultiPolygon();
960 : 0 : break;
961 : : default:
962 : 0 : collection = new QgsGeometryCollection();
963 : 0 : break;
964 : : }
965 : 1 : collection->addGeometry( lwgeom_out.release() ); // takes ownership
966 : 1 : lwgeom_out.reset( collection );
967 : 1 : }
968 : :
969 : 6 : return lwgeom_out;
970 : 6 : }
|
