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)

(cherry picked from commit 00a102c)
(cherry picked from commit b75ba12)

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 );
   }