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