Fix performance issue when splitting polygons (fixes #34326)

Using test data with a HUGE polygon with several hundreds of thousands of vertices and nearly a thousand of holes, split geometry was taking really long time. On my laptop ~67 seconds for a simple split line. With this fix, the time went down to ~10 seconds for the test polygon. The issue was that split tool was doing an expensive exact intersection calculation for each polygon part returned from polygonize, and then calculating areas to see whether the original polygon part and intersected part were more or less the same. I think that can be replaced by a much simpler point in polygon test to figure out whether the particular polygon part from polygonize is falls inside the original geometry or whether it falls into a hole or exterior. In my tests the time could be further improved to ~3 seconds by disabling mergeGeometriesMultiTypeSplit() function - it is a bit of a mistery to me what problem does it solve, however it has N^2 running time with N being number of rings - but I did not feel adventurous enough to poke into it. (mainly because there are no unit tests that would verify its behavior)
qgis · Jun 19, 2020 · 00a102c · 00a102c
1 parent 0cb8408
commit 00a102c
Showing 1 changed file with 10 additions and 18 deletions.
diff --git a/src/core/geometry/qgsgeos.cpp b/src/core/geometry/qgsgeos.cpp
@@ -896,32 +896,24 @@ QgsGeometryEngine::EngineOperationResult QgsGeos::splitPolygonGeometry( GEOSGeom
     return InvalidBaseGeometry;
   }
 
+  // we will need prepared geometry for intersection tests
+  const_cast<QgsGeos *>( this )->prepareGeometry();
+  if ( !mGeosPrepared )
+    return EngineError;
+
   //test every polygon if contained in original geometry
   //include in result if yes
   QVector<GEOSGeometry *> testedGeometries;
-  geos::unique_ptr intersectGeometry;
-
-  //ratio intersect geometry / geometry. This should be close to 1
-  //if the polygon belongs to the input geometry
 
+  // test whether the polygon parts returned from polygonize algorithm actually
+  // belong to the source polygon geometry (if the source polygon contains some holes,
+  // those would be also returned by polygonize and we need to skip them)
   for ( int i = 0; i < numberOfGeometries( polygons.get() ); i++ )
   {
     const GEOSGeometry *polygon = GEOSGetGeometryN_r( geosinit()->ctxt, polygons.get(), i );
-    intersectGeometry.reset( GEOSIntersection_r( geosinit()->ctxt, mGeos.get(), polygon ) );
-    if ( !intersectGeometry )
-    {
-      QgsDebugMsg( QStringLiteral( "intersectGeometry is nullptr" ) );
-      continue;
-    }
-
-    double intersectionArea;
-    GEOSArea_r( geosinit()->ctxt, intersectGeometry.get(), &intersectionArea );
-
-    double polygonArea;
-    GEOSArea_r( geosinit()->ctxt, polygon, &polygonArea );
 
-    const double areaRatio = intersectionArea / polygonArea;
-    if ( areaRatio > 0.99 && areaRatio < 1.01 )
+    geos::unique_ptr pointOnSurface( GEOSPointOnSurface_r( geosinit()->ctxt, polygon ) );
+    if ( pointOnSurface && GEOSPreparedIntersects_r( geosinit()->ctxt, mGeosPrepared.get(), pointOnSurface.get() ) )
       testedGeometries << GEOSGeom_clone_r( geosinit()->ctxt, polygon );
   }