qgsgeometrymakevalid.cpp
Go to the documentation of this file.
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 *
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 typedef struct Face_t
37 {
38  const GEOSGeometry *geom = nullptr;
39  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 static Face *newFace( const GEOSGeometry *g )
50 {
51  GEOSContextHandle_t handle = QgsGeos::getGEOSHandler();
52
53  Face *f = new Face;
54  f->geom = g;
55  f->env = GEOSEnvelope_r( handle, f->geom );
56  GEOSArea_r( handle, f->env, &f->envarea );
57  f->parent = nullptr;
58  return f;
59 }
60
61 static unsigned int countParens( const Face *f )
62 {
63  unsigned int pcount = 0;
64  while ( f->parent )
65  {
66  ++pcount;
67  f = f->parent;
68  }
69  return pcount;
70 }
71
72 // Destroy the face and release memory associated with it
73 static void delFace( Face *f )
74 {
75  GEOSContextHandle_t handle = QgsGeos::getGEOSHandler();
76  GEOSGeom_destroy_r( handle, f->env );
77  delete f;
78 }
79
80
81 static int compare_by_envarea( const void *g1, const void *g2 )
82 {
83  Face *f1 = *( Face ** )g1;
84  Face *f2 = *( Face ** )g2;
85  double n1 = f1->envarea;
86  double n2 = f2->envarea;
87
88  if ( n1 < n2 ) return 1;
89  if ( n1 > n2 ) return -1;
90  return 0;
91 }
92
93 // Find holes of each face
94 static void findFaceHoles( Face **faces, int nfaces )
95 {
96  GEOSContextHandle_t handle = QgsGeos::getGEOSHandler();
97
98  // We sort by envelope area so that we know holes are only
99  // after their shells
100  qsort( faces, nfaces, sizeof( Face * ), compare_by_envarea );
101  for ( int i = 0; i < nfaces; ++i )
102  {
103  Face *f = faces[i];
104  int nholes = GEOSGetNumInteriorRings_r( handle, f->geom );
105  for ( int h = 0; h < nholes; ++h )
106  {
107  const GEOSGeometry *hole = GEOSGetInteriorRingN_r( handle, f->geom, h );
108  for ( int j = i + 1; j < nfaces; ++j )
109  {
110  Face *f2 = faces[j];
111  if ( f2->parent )
112  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  if ( GEOSEquals_r( handle, GEOSGetExteriorRing_r( handle, f2->geom ), hole ) )
119  {
120  f2->parent = f;
121  break;
122  }
123  }
124  }
125  }
126 }
127
128 static GEOSGeometry *collectFacesWithEvenAncestors( Face **faces, int nfaces )
129 {
130  GEOSContextHandle_t handle = QgsGeos::getGEOSHandler();
131
132  unsigned int ngeoms = 0;
133  GEOSGeometry **geoms = new GEOSGeometry*[nfaces];
134  for ( int i = 0; i < nfaces; ++i )
135  {
136  Face *f = faces[i];
137  if ( countParens( f ) % 2 )
138  continue; // we skip odd parents geoms
139  geoms[ngeoms++] = GEOSGeom_clone_r( handle, f->geom );
140  }
141
142  GEOSGeometry *ret = GEOSGeom_createCollection_r( handle, GEOS_MULTIPOLYGON, geoms, ngeoms );
143  delete [] geoms;
144  return ret;
145 }
146
148 static GEOSGeometry *LWGEOM_GEOS_buildArea( const GEOSGeometry *geom_in, QString &errorMessage )
149 {
150  GEOSContextHandle_t handle = QgsGeos::getGEOSHandler();
151
152  GEOSGeometry *tmp = nullptr;
153  GEOSGeometry *shp = nullptr;
154  GEOSGeometry *geos_result = nullptr;
155  int srid = GEOSGetSRID_r( handle, geom_in );
156
157  GEOSGeometry const *vgeoms[1];
158  vgeoms[0] = geom_in;
159  try
160  {
161  geos_result = GEOSPolygonize_r( handle, vgeoms, 1 );
162  }
163  catch ( GEOSException &e )
164  {
165  errorMessage = QStringLiteral( "GEOSPolygonize(): %1" ).arg( e.what() );
166  return nullptr;
167  }
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  int ngeoms = GEOSGetNumGeometries_r( handle, geos_result );
180
181  // No geometries in collection, early out
182  if ( ngeoms == 0 )
183  {
184  GEOSSetSRID_r( handle, geos_result, srid );
185  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  if ( ngeoms == 1 )
191  {
192  tmp = ( GEOSGeometry * )GEOSGetGeometryN_r( handle, geos_result, 0 );
193  if ( ! tmp )
194  {
195  GEOSGeom_destroy_r( handle, geos_result );
196  return nullptr;
197  }
198  shp = GEOSGeom_clone_r( handle, tmp );
199  GEOSGeom_destroy_r( handle, geos_result ); // only safe after the clone above
200  GEOSSetSRID_r( handle, shp, srid );
201  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  Face **geoms = new Face*[ngeoms];
234  for ( int i = 0; i < ngeoms; ++i )
235  geoms[i] = newFace( GEOSGetGeometryN_r( handle, geos_result, i ) );
236
237  // Find faces representing other faces holes
238  findFaceHoles( geoms, ngeoms );
239
240  // Build a MultiPolygon composed only by faces with an
241  // even number of ancestors
242  tmp = collectFacesWithEvenAncestors( geoms, ngeoms );
243
244  // Cleanup face structures
245  for ( int i = 0; i < ngeoms; ++i )
246  delFace( geoms[i] );
247  delete [] geoms;
248
249  // Faces referenced memory owned by geos_result.
250  // It is safe to destroy geos_result after deleting them.
251  GEOSGeom_destroy_r( handle, geos_result );
252
253  // Run a single overlay operation to dissolve shared edges
254  shp = GEOSUnionCascaded_r( handle, tmp );
255  if ( !shp )
256  {
257  GEOSGeom_destroy_r( handle, tmp );
258  return nullptr;
259  }
260
261  GEOSGeom_destroy_r( handle, tmp );
262
263  GEOSSetSRID_r( handle, shp, srid );
264
265  return shp;
266 }
268
269 // Return Nth vertex in GEOSGeometry as a POINT.
270 // May return NULL if the geometry has NO vertexex.
271 static GEOSGeometry *LWGEOM_GEOS_getPointN( const GEOSGeometry *g_in, uint32_t n )
272 {
273  GEOSContextHandle_t handle = QgsGeos::getGEOSHandler();
274
275  uint32_t dims;
276  const GEOSCoordSequence *seq_in = nullptr;
277  GEOSCoordSeq seq_out;
278  double val;
279  uint32_t sz;
280  int gn;
281  GEOSGeometry *ret = nullptr;
282
283  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  for ( gn = 0; gn < GEOSGetNumGeometries_r( handle, g_in ); ++gn )
291  {
292  const GEOSGeometry *g = GEOSGetGeometryN_r( handle, g_in, gn );
293  ret = LWGEOM_GEOS_getPointN( g, n );
294  if ( ret ) return ret;
295  }
296  break;
297  }
298
299  case GEOS_POLYGON:
300  {
301  ret = LWGEOM_GEOS_getPointN( GEOSGetExteriorRing_r( handle, g_in ), n );
302  if ( ret ) return ret;
303  for ( gn = 0; gn < GEOSGetNumInteriorRings_r( handle, g_in ); ++gn )
304  {
305  const GEOSGeometry *g = GEOSGetInteriorRingN_r( handle, g_in, gn );
306  ret = LWGEOM_GEOS_getPointN( g, n );
307  if ( ret ) return ret;
308  }
309  break;
310  }
311
312  case GEOS_POINT:
313  case GEOS_LINESTRING:
314  case GEOS_LINEARRING:
315  break;
316
317  }
318
319  seq_in = GEOSGeom_getCoordSeq_r( handle, g_in );
320  if ( ! seq_in ) return nullptr;
321  if ( ! GEOSCoordSeq_getSize_r( handle, seq_in, &sz ) ) return nullptr;
322  if ( ! sz ) return nullptr;
323
324  if ( ! GEOSCoordSeq_getDimensions_r( handle, seq_in, &dims ) ) return nullptr;
325
326  seq_out = GEOSCoordSeq_create_r( handle, 1, dims );
327  if ( ! seq_out ) return nullptr;
328
329  if ( ! GEOSCoordSeq_getX_r( handle, seq_in, n, &val ) ) return nullptr;
330  if ( ! GEOSCoordSeq_setX_r( handle, seq_out, n, val ) ) return nullptr;
331  if ( ! GEOSCoordSeq_getY_r( handle, seq_in, n, &val ) ) return nullptr;
332  if ( ! GEOSCoordSeq_setY_r( handle, seq_out, n, val ) ) return nullptr;
333  if ( dims > 2 )
334  {
335  if ( ! GEOSCoordSeq_getZ_r( handle, seq_in, n, &val ) ) return nullptr;
336  if ( ! GEOSCoordSeq_setZ_r( handle, seq_out, n, val ) ) return nullptr;
337  }
338
339  return GEOSGeom_createPoint_r( handle, seq_out );
340 }
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 static bool lwgeom_make_geos_friendly( QgsAbstractGeometry *geom )
350 {
351  QgsDebugMsg( QStringLiteral( "lwgeom_make_geos_friendly enter (type %1)" ).arg( geom->wkbType() ) );
352  switch ( QgsWkbTypes::flatType( geom->wkbType() ) )
353  {
354  case QgsWkbTypes::Point:
356  // a point is always valid
357  return true;
358  break;
359
361  // lines need at least 2 points
362  return lwline_make_geos_friendly( qgsgeometry_cast<QgsLineString *>( geom ) );
363  break;
364
366  // polygons need all rings closed and with npoints > 3
367  return lwpoly_make_geos_friendly( qgsgeometry_cast<QgsPolygon *>( geom ) );
368  break;
369
373  return lwcollection_make_geos_friendly( qgsgeometry_cast<QgsGeometryCollection *>( geom ) );
374  break;
375
376  default:
377  QgsDebugMsg( QStringLiteral( "lwgeom_make_geos_friendly: unsupported input geometry type: %1" ).arg( geom->wkbType() ) );
378  break;
379  }
380  return false;
381 }
382
383
384 static bool ring_make_geos_friendly( QgsCurve *ring )
385 {
386  if ( ring->nCoordinates() == 0 )
387  return false;
388
389  // earlier we allowed in only geometries with straight segments
390  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  QgsPoint p1 = linestring->startPoint(), p2 = linestring->endPoint();
396  if ( p1.x() != p2.x() || p1.y() != p2.y() )
398
399  // must have at least 4 coordinates to be accepted by GEOS
400
401  while ( linestring->nCoordinates() < 4 )
403
404  return true;
405 }
406
407 // Make sure all rings are closed and have > 3 points.
408 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  for ( int i = 0; i < poly->numInteriorRings(); i++ )
416  {
417  if ( !ring_make_geos_friendly( qgsgeometry_cast<QgsCurve *>( poly->interiorRing( i ) ) ) )
418  return false;
419  }
420
421  return true;
422 }
423
424 // Need NO or >1 points. Duplicate first if only one.
425 static bool lwline_make_geos_friendly( QgsLineString *line )
426 {
427  if ( line->numPoints() == 1 ) // 0 is fine, 2 is fine
428  {
430  }
431  return true;
432 }
433
434 static bool lwcollection_make_geos_friendly( QgsGeometryCollection *g )
435 {
436  for ( int i = 0; i < g->numGeometries(); i++ )
437  {
438  if ( !lwgeom_make_geos_friendly( g->geometryN( i ) ) )
439  return false;
440  }
441
442  return true;
443 }
444
445
446 // Fully node given linework
447 static GEOSGeometry *LWGEOM_GEOS_nodeLines( const GEOSGeometry *lines )
448 {
449  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  GEOSGeometry *point = LWGEOM_GEOS_getPointN( lines, 0 );
457  if ( ! point )
458  return nullptr;
459
460  GEOSGeometry *noded = nullptr;
461  try
462  {
463  noded = GEOSUnion_r( handle, lines, point );
464  }
465  catch ( GEOSException & )
466  {
467  // no need to do anything here - we'll return nullptr anyway
468  }
469  GEOSGeom_destroy_r( handle, point );
470  return noded;
471 }
472
473 // Will return NULL on error (expect error handler being called by then)
475 static GEOSGeometry *LWGEOM_GEOS_makeValidPolygon( const GEOSGeometry *gin, QString &errorMessage )
476 {
477  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  unsigned int nvgeoms = 0;
484
485  Q_ASSERT( GEOSGeomTypeId_r( handle, gin ) == GEOS_POLYGON ||
486  GEOSGeomTypeId_r( handle, gin ) == GEOS_MULTIPOLYGON );
487
488  geos_bound = GEOSBoundary_r( handle, gin );
489  if ( !geos_bound )
490  return nullptr;
491
492  // Use noded boundaries as initial "cut" edges
493
494  geos_cut_edges = LWGEOM_GEOS_nodeLines( geos_bound );
495  if ( !geos_cut_edges )
496  {
497  GEOSGeom_destroy_r( handle, geos_bound );
498  errorMessage = QStringLiteral( "LWGEOM_GEOS_nodeLines() failed" );
499  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  GEOSGeometry *pi = nullptr;
506  GEOSGeometry *po = nullptr;
507
508  try
509  {
510  pi = GEOSGeom_extractUniquePoints_r( handle, geos_bound );
511  }
512  catch ( GEOSException &e )
513  {
514  GEOSGeom_destroy_r( handle, geos_bound );
515  errorMessage = QStringLiteral( "GEOSGeom_extractUniquePoints(): %1" ).arg( e.what() );
516  return nullptr;
517  }
518
519  try
520  {
521  po = GEOSGeom_extractUniquePoints_r( handle, geos_cut_edges );
522  }
523  catch ( GEOSException &e )
524  {
525  GEOSGeom_destroy_r( handle, geos_bound );
526  GEOSGeom_destroy_r( handle, pi );
527  errorMessage = QStringLiteral( "GEOSGeom_extractUniquePoints(): %1" ).arg( e.what() );
528  return nullptr;
529  }
530
531  try
532  {
533  collapse_points = GEOSDifference_r( handle, pi, po );
534  }
535  catch ( GEOSException &e )
536  {
537  GEOSGeom_destroy_r( handle, geos_bound );
538  GEOSGeom_destroy_r( handle, pi );
539  GEOSGeom_destroy_r( handle, po );
540  errorMessage = QStringLiteral( "GEOSDifference(): %1" ).arg( e.what() );
541  return nullptr;
542  }
543
544  GEOSGeom_destroy_r( handle, pi );
545  GEOSGeom_destroy_r( handle, po );
546  }
547  GEOSGeom_destroy_r( handle, geos_bound );
548
549  // And use an empty geometry as initial "area"
550  try
551  {
552  geos_area = GEOSGeom_createEmptyPolygon_r( handle );
553  }
554  catch ( GEOSException &e )
555  {
556  errorMessage = QStringLiteral( "GEOSGeom_createEmptyPolygon(): %1" ).arg( e.what() );
557  GEOSGeom_destroy_r( handle, geos_cut_edges );
558  return nullptr;
559  }
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  while ( GEOSGetNumGeometries_r( handle, geos_cut_edges ) )
566  {
567  GEOSGeometry *new_area = nullptr;
568  GEOSGeometry *new_area_bound = nullptr;
569  GEOSGeometry *symdif = nullptr;
570  GEOSGeometry *new_cut_edges = nullptr;
571
572  // ASSUMPTION: cut_edges should already be fully noded
573
574  try
575  {
576  new_area = LWGEOM_GEOS_buildArea( geos_cut_edges, errorMessage );
577  }
578  catch ( GEOSException &e )
579  {
580  GEOSGeom_destroy_r( handle, geos_cut_edges );
581  GEOSGeom_destroy_r( handle, geos_area );
582  errorMessage = QStringLiteral( "LWGEOM_GEOS_buildArea() threw an error: %1" ).arg( e.what() );
583  return nullptr;
584  }
585
586  if ( GEOSisEmpty_r( handle, new_area ) )
587  {
588  // no more rings can be build with thes edges
589  GEOSGeom_destroy_r( handle, new_area );
590  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  new_area_bound = GEOSBoundary_r( handle, new_area );
600  }
601  catch ( GEOSException &e )
602  {
603  // We did check for empty area already so
604  // this must be some other error
605  errorMessage = QStringLiteral( "GEOSBoundary() threw an error: %1" ).arg( e.what() );
606  GEOSGeom_destroy_r( handle, new_area );
607  GEOSGeom_destroy_r( handle, geos_area );
608  return nullptr;
609  }
610
611  // Now symdiff new and old area
612  try
613  {
614  symdif = GEOSSymDifference_r( handle, geos_area, new_area );
615  }
616  catch ( GEOSException &e )
617  {
618  GEOSGeom_destroy_r( handle, geos_cut_edges );
619  GEOSGeom_destroy_r( handle, new_area );
620  GEOSGeom_destroy_r( handle, new_area_bound );
621  GEOSGeom_destroy_r( handle, geos_area );
622  errorMessage = QStringLiteral( "GEOSSymDifference() threw an error: %1" ).arg( e.what() );
623  return nullptr;
624  }
625
626  GEOSGeom_destroy_r( handle, geos_area );
627  GEOSGeom_destroy_r( handle, new_area );
628  geos_area = symdif;
629  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  new_cut_edges = GEOSDifference_r( handle, geos_cut_edges, new_area_bound );
642  }
643  catch ( GEOSException &e )
644  {
645  GEOSGeom_destroy_r( handle, geos_cut_edges );
646  GEOSGeom_destroy_r( handle, new_area_bound );
647  GEOSGeom_destroy_r( handle, geos_area );
648  errorMessage = QStringLiteral( "GEOSDifference() threw an error: %1" ).arg( e.what() );
649  return nullptr;
650  }
651  GEOSGeom_destroy_r( handle, geos_cut_edges );
652  GEOSGeom_destroy_r( handle, new_area_bound );
653  geos_cut_edges = new_cut_edges;
654  }
655
656  if ( !GEOSisEmpty_r( handle, geos_area ) )
657  {
658  vgeoms[nvgeoms++] = geos_area;
659  }
660  else
661  {
662  GEOSGeom_destroy_r( handle, geos_area );
663  }
664
665  if ( !GEOSisEmpty_r( handle, geos_cut_edges ) )
666  {
667  vgeoms[nvgeoms++] = geos_cut_edges;
668  }
669  else
670  {
671  GEOSGeom_destroy_r( handle, geos_cut_edges );
672  }
673
674  if ( !GEOSisEmpty_r( handle, collapse_points ) )
675  {
676  vgeoms[nvgeoms++] = collapse_points;
677  }
678  else
679  {
680  GEOSGeom_destroy_r( handle, collapse_points );
681  }
682
683  if ( 1 == nvgeoms )
684  {
685  // Return cut edges
686  gout = vgeoms[0];
687  }
688  else
689  {
690  // Collect areas and lines (if any line)
691  try
692  {
693  gout = GEOSGeom_createCollection_r( handle, GEOS_GEOMETRYCOLLECTION, vgeoms, nvgeoms );
694  }
695  catch ( GEOSException &e )
696  {
697  errorMessage = QStringLiteral( "GEOSGeom_createCollection() threw an error: %1" ).arg( e.what() );
698  // TODO: cleanup!
699  return nullptr;
700  }
701  }
702
703  return gout;
704
705 }
707
708 static GEOSGeometry *LWGEOM_GEOS_makeValidLine( const GEOSGeometry *gin, QString &errorMessage )
709 {
710  Q_UNUSED( errorMessage );
711  return LWGEOM_GEOS_nodeLines( gin );
712 }
713
714 static GEOSGeometry *LWGEOM_GEOS_makeValidMultiLine( const GEOSGeometry *gin, QString &errorMessage )
715 {
716  GEOSContextHandle_t handle = QgsGeos::getGEOSHandler();
717
718  int ngeoms = GEOSGetNumGeometries_r( handle, gin );
719  uint32_t nlines_alloc = ngeoms;
720  QVector<GEOSGeometry *> lines;
721  QVector<GEOSGeometry *> points;
722  lines.reserve( nlines_alloc );
723  points.reserve( ngeoms );
724
725  for ( int i = 0; i < ngeoms; ++i )
726  {
727  const GEOSGeometry *g = GEOSGetGeometryN_r( handle, gin, i );
728  GEOSGeometry *vg = LWGEOM_GEOS_makeValidLine( g, errorMessage );
729  if ( GEOSisEmpty_r( handle, vg ) )
730  {
732  GEOSGeom_destroy_r( handle, vg );
733  }
734  if ( GEOSGeomTypeId_r( handle, vg ) == GEOS_POINT )
735  {
736  points.append( vg );
737  }
738  else if ( GEOSGeomTypeId_r( handle, vg ) == GEOS_LINESTRING )
739  {
740  lines.append( vg );
741  }
742  else if ( GEOSGeomTypeId_r( handle, vg ) == GEOS_MULTILINESTRING )
743  {
744  int nsubgeoms = GEOSGetNumGeometries_r( handle, vg );
745  nlines_alloc += nsubgeoms;
746  lines.reserve( nlines_alloc );
747  for ( int j = 0; j < nsubgeoms; ++j )
748  {
749  // NOTE: ownership of the cloned geoms will be taken by final collection
750  lines.append( GEOSGeom_clone_r( handle, GEOSGetGeometryN_r( handle, vg, j ) ) );
751  }
752  }
753  else
754  {
755  // NOTE: return from GEOSGeomType will leak but we really don't expect this to happen
756  errorMessage = QStringLiteral( "unexpected geom type returned by LWGEOM_GEOS_makeValid: %1" ).arg( GEOSGeomTypeId_r( handle, vg ) );
757  }
758  }
759
760  GEOSGeometry *mpoint_out = nullptr;
761  if ( points.count() )
762  {
763  if ( points.count() > 1 )
764  {
765  mpoint_out = GEOSGeom_createCollection_r( handle, GEOS_MULTIPOINT, points.data(), points.count() );
766  }
767  else
768  {
769  mpoint_out = points[0];
770  }
771  }
772
773  GEOSGeometry *mline_out = nullptr;
774  if ( lines.count() )
775  {
776  if ( lines.count() > 1 )
777  {
778  mline_out = GEOSGeom_createCollection_r( handle, GEOS_MULTILINESTRING, lines.data(), lines.count() );
779  }
780  else
781  {
782  mline_out = lines[0];
783  }
784  }
785
786  if ( mline_out && mpoint_out )
787  {
788  GEOSGeometry *collection[2];
789  collection[0] = mline_out;
790  collection[1] = mpoint_out;
791  return GEOSGeom_createCollection_r( handle, GEOS_GEOMETRYCOLLECTION, collection, 2 );
792  }
793  else if ( mline_out )
794  {
795  return mline_out;
796  }
797  else if ( mpoint_out )
798  {
799  return mpoint_out;
800  }
801  else
802  {
803  return nullptr;
804  }
805 }
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 static GEOSGeometry *LWGEOM_GEOS_makeValidCollection( const GEOSGeometry *gin, QString &errorMessage )
812 {
813  GEOSContextHandle_t handle = QgsGeos::getGEOSHandler();
814
815  int nvgeoms = GEOSGetNumGeometries_r( handle, gin );
816  if ( nvgeoms == -1 )
817  {
818  errorMessage = QStringLiteral( "GEOSGetNumGeometries: %1" ).arg( QStringLiteral( "?" ) );
819  return nullptr;
820  }
821
822  QVector<GEOSGeometry *> vgeoms( nvgeoms );
823  for ( int i = 0; i < nvgeoms; ++i )
824  {
825  vgeoms[i] = LWGEOM_GEOS_makeValid( GEOSGetGeometryN_r( handle, gin, i ), errorMessage );
826  if ( ! vgeoms[i] )
827  {
828  while ( i-- )
829  GEOSGeom_destroy_r( handle, vgeoms[i] );
830  // we expect lwerror being called already by makeValid
831  return nullptr;
832  }
833  }
834
835  // Collect areas and lines (if any line)
836  try
837  {
838  return GEOSGeom_createCollection_r( handle, GEOS_GEOMETRYCOLLECTION, vgeoms.data(), nvgeoms );
839  }
840  catch ( GEOSException &e )
841  {
842  // cleanup and throw
843  for ( int i = 0; i < nvgeoms; ++i )
844  GEOSGeom_destroy_r( handle, vgeoms[i] );
845  errorMessage = QStringLiteral( "GEOSGeom_createCollection() threw an error: %1" ).arg( e.what() );
846  return nullptr;
847  }
848 }
849
850
851 static GEOSGeometry *LWGEOM_GEOS_makeValid( const GEOSGeometry *gin, QString &errorMessage )
852 {
853  GEOSContextHandle_t handle = QgsGeos::getGEOSHandler();
854
855  // return what we got so far if already valid
857  try
858  {
859  if ( GEOSisValid_r( handle, gin ) )
860  {
861  // It's valid at this step, return what we have
862  return GEOSGeom_clone_r( handle, gin );
863  }
864  }
865  catch ( GEOSException &e )
866  {
867  // I don't think should ever happen
868  errorMessage = QStringLiteral( "GEOSisValid(): %1" ).arg( e.what() );
869  return nullptr;
870  }
872
873  // make what we got valid
874
875  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  QgsDebugMsg( QStringLiteral( "PUNTUAL geometry resulted invalid to GEOS -- dunno how to clean that up" ) );
881  return nullptr;
882
883  case GEOS_LINESTRING:
884  return LWGEOM_GEOS_makeValidLine( gin, errorMessage );
885
886  case GEOS_MULTILINESTRING:
887  return LWGEOM_GEOS_makeValidMultiLine( gin, errorMessage );
888
889  case GEOS_POLYGON:
890  case GEOS_MULTIPOLYGON:
891  return LWGEOM_GEOS_makeValidPolygon( gin, errorMessage );
892
893  case GEOS_GEOMETRYCOLLECTION:
894  return LWGEOM_GEOS_makeValidCollection( gin, errorMessage );
895
896  default:
897  errorMessage = QStringLiteral( "ST_MakeValid: doesn't support geometry type: %1" ).arg( GEOSGeomTypeId_r( handle, gin ) );
898  return nullptr;
899  }
900 }
901
902
903 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  GEOSContextHandle_t handle = QgsGeos::getGEOSHandler();
910
911  geos::unique_ptr geosgeom = QgsGeos::asGeos( lwgeom_in );
912  if ( !geosgeom )
913  {
914  QgsDebugMsg( QStringLiteral( "Original geom can't be converted to GEOS - will try cleaning that up first" ) );
915
916  std::unique_ptr<QgsAbstractGeometry> lwgeom_in_clone( lwgeom_in->clone() );
917  if ( !lwgeom_make_geos_friendly( lwgeom_in_clone.get() ) )
918  {
919  QgsDebugMsg( QStringLiteral( "Could not make a valid geometry out of input" ) );
920  }
921
922  // try again as we did cleanup now
923  // TODO: invoke LWGEOM2GEOS directly with autoclean ?
924  geosgeom = QgsGeos::asGeos( lwgeom_in_clone.get() );
925
926  if ( ! geosgeom )
927  {
928  errorMessage = QStringLiteral( "Could not convert QGIS geom to GEOS" );
929  return nullptr;
930  }
931  }
932  else
933  {
934  QgsDebugMsgLevel( QStringLiteral( "original geom converted to GEOS" ), 4 );
935  }
936
937  GEOSGeometry *geosout = LWGEOM_GEOS_makeValid( geosgeom.get(), errorMessage );
938  if ( !geosout )
939  return nullptr;
940
941  std::unique_ptr< QgsAbstractGeometry > lwgeom_out = QgsGeos::fromGeos( geosout );
942  GEOSGeom_destroy_r( handle, geosout );
943  if ( !lwgeom_out )
944  return nullptr;
945
946  // force multi-type if we had a multi-type before
947  if ( QgsWkbTypes::isMultiType( lwgeom_in->wkbType() ) && !QgsWkbTypes::isMultiType( lwgeom_out->wkbType() ) )
948  {
949  QgsGeometryCollection *collection = nullptr;
950  switch ( QgsWkbTypes::multiType( lwgeom_out->wkbType() ) )
951  {
953  collection = new QgsMultiPoint();
954  break;
956  collection = new QgsMultiLineString();
957  break;
959  collection = new QgsMultiPolygon();
960  break;
961  default:
962  collection = new QgsGeometryCollection();
963  break;
964  }
965  collection->addGeometry( lwgeom_out.release() ); // takes ownership
966  lwgeom_out.reset( collection );
967  }
968
969  return lwgeom_out;
970 }
double y
Definition: qgspoint.h:42
static Type multiType(Type type)
Returns the multi type for a WKB type.
Definition: qgswkbtypes.h:298
Multi point geometry collection.
Definition: qgsmultipoint.h:29
static bool isMultiType(Type type)
Returns true if the WKB type is a multi type.
Definition: qgswkbtypes.h:695
#define QgsDebugMsg(str)
Definition: qgslogger.h:38
QgsWkbTypes::Type wkbType() const
Returns the WKB type of the geometry.
Multi line string geometry collection.
const QgsAbstractGeometry * geometryN(int n) const
Returns a const reference to a geometry from within the collection.
virtual QgsAbstractGeometry * clone() const =0
Clones the geometry by performing a deep copy.
QgsPoint endPoint() const override
Returns the end point of the curve.
static GEOSContextHandle_t getGEOSHandler()
Definition: qgsgeos.cpp:2821
int numInteriorRings() const
Returns the number of interior rings contained with the curve polygon.
int numPoints() const override
Returns the number of points in the curve.
virtual int nCoordinates() const
Returns the number of nodes contained in the geometry.
static std::unique_ptr< QgsAbstractGeometry > fromGeos(const GEOSGeometry *geos)
Create a geometry from a GEOSGeometry.
Definition: qgsgeos.cpp:1081
Adds a new vertex to the end of the line string.
#define QgsDebugMsgLevel(str, level)
Definition: qgslogger.h:39
Geometry collection.
GEOSGeometry * env
std::unique_ptr< GEOSGeometry, GeosDeleter > unique_ptr
Scoped GEOS pointer.
Definition: qgsgeos.h:79
struct Face_t Face
T qgsgeometry_cast(const QgsAbstractGeometry *geom)
Abstract base class for curved geometry type.
Definition: qgscurve.h:35
int numGeometries() const
Returns the number of geometries within the collection.
Abstract base class for all geometries.
Point geometry type, with support for z-dimension and m-values.
Definition: qgspoint.h:37
QgsPoint startPoint() const override
Returns the starting point of the curve.
struct Face_t * parent
#define Q_NOWARN_UNREACHABLE_POP
Definition: qgis.h:617
static geos::unique_ptr asGeos(const QgsGeometry &geometry, double precision=0)
Returns a geos geometry - caller takes ownership of the object (should be deleted with GEOSGeom_destr...
Definition: qgsgeos.cpp:166
Multi polygon geometry collection.
Line string geometry type, with support for z-dimension and m-values.
Definition: qgslinestring.h:43
const QgsCurve * interiorRing(int i) const
Retrieves an interior ring from the curve polygon.
std::unique_ptr< QgsAbstractGeometry > _qgis_lwgeom_make_valid(const QgsAbstractGeometry *lwgeom_in, QString &errorMessage)
Implementation of QgsGeometry::makeValid(). Not a public API.
const GEOSGeometry * geom
#define Q_NOWARN_UNREACHABLE_PUSH
Definition: qgis.h:616
Polygon geometry type.
Definition: qgspolygon.h:31
int nCoordinates() const override
Returns the number of nodes contained in the geometry.
static Type flatType(Type type)
Returns the flat type for a WKB type.
Definition: qgswkbtypes.h:565