[FEATURE] Subdivide algorithm for QgsGeometry

Subdivides the geometry. The returned geometry will be a collection containing subdivided parts from the original geometry, where no part has more then the specified maximum number of nodes.
qgis · Jun 14, 2017 · e74395d · e74395d
1 parent bde0c72
commit e74395d
Show file tree

Hide file tree

Showing 6 changed files with 210 additions and 0 deletions.
diff --git a/python/core/geometry/qgsgeometry.sip b/python/core/geometry/qgsgeometry.sip
@@ -782,6 +782,23 @@ Returns the smallest convex polygon that contains all the points in the geometry
  If ``edgesOnly`` is true than line string boundary geometries will be returned
  instead of polygons.
  An empty geometry will be returned if the diagram could not be calculated.
+.. versionadded:: 3.0
+ :rtype: QgsGeometry
+%End
+
+    QgsGeometry subdivide( int maxNodes = 256 ) const;
+%Docstring
+ Subdivides the geometry. The returned geometry will be a collection containing subdivided parts
+ from the original geometry, where no part has more then the specified maximum number of nodes (``maxNodes``).
+
+ This is useful for dividing a complex geometry into less complex parts, which are better able to be spatially
+ indexed and faster to perform further operations such as intersects on. The returned geometry parts may
+ not be valid and may contain self-intersections.
+
+ The minimum allowed value for ``maxNodes`` is 8.
+
+ Curved geometries will be segmentized before subdivision.
+
 .. versionadded:: 3.0
  :rtype: QgsGeometry
 %End

diff --git a/src/core/geometry/qgsgeometry.cpp b/src/core/geometry/qgsgeometry.cpp
@@ -1655,6 +1655,30 @@ QgsGeometry QgsGeometry::delaunayTriangulation( double tolerance, bool edgesOnly
   return geos.delaunayTriangulation( tolerance, edgesOnly );
 }
 
+QgsGeometry QgsGeometry::subdivide( int maxNodes ) const
+{
+  if ( !d->geometry )
+  {
+    return QgsGeometry();
+  }
+
+  const QgsAbstractGeometry *geom = d->geometry;
+  std::unique_ptr< QgsAbstractGeometry > segmentizedCopy;
+  if ( QgsWkbTypes::isCurvedType( d->geometry->wkbType() ) )
+  {
+    segmentizedCopy.reset( d->geometry->segmentize() );
+    geom = segmentizedCopy.get();
+  }
+
+  QgsGeos geos( geom );
+  QgsAbstractGeometry *result = geos.subdivide( maxNodes );
+  if ( !result )
+  {
+    return QgsGeometry();
+  }
+  return QgsGeometry( result );
+}
+
 QgsGeometry QgsGeometry::interpolate( double distance ) const
 {
   if ( !d->geometry )

diff --git a/src/core/geometry/qgsgeometry.h b/src/core/geometry/qgsgeometry.h
@@ -733,6 +733,22 @@ class CORE_EXPORT QgsGeometry
      */
     QgsGeometry delaunayTriangulation( double tolerance = 0.0, bool edgesOnly = false ) const;
 
+    /**
+     * Subdivides the geometry. The returned geometry will be a collection containing subdivided parts
+     * from the original geometry, where no part has more then the specified maximum number of nodes (\a maxNodes).
+     *
+     * This is useful for dividing a complex geometry into less complex parts, which are better able to be spatially
+     * indexed and faster to perform further operations such as intersects on. The returned geometry parts may
+     * not be valid and may contain self-intersections.
+     *
+     * The minimum allowed value for \a maxNodes is 8.
+     *
+     * Curved geometries will be segmentized before subdivision.
+     *
+     * \since QGIS 3.0
+     */
+    QgsGeometry subdivide( int maxNodes = 256 ) const;
+
     /**
      * Return interpolated point on line at distance
      * \since QGIS 1.9

diff --git a/src/core/geometry/qgsgeos.cpp b/src/core/geometry/qgsgeos.cpp
@@ -213,6 +213,122 @@ QgsAbstractGeometry *QgsGeos::clip( const QgsRectangle &rect, QString *errorMsg
   }
 }
 
+
+
+
+void QgsGeos::subdivideRecursive( const GEOSGeometry *currentPart, int maxNodes, int depth, QgsGeometryCollection *parts, const QgsRectangle &clipRect ) const
+{
+  int partType = GEOSGeomTypeId_r( geosinit.ctxt, currentPart );
+  if ( qgsDoubleNear( clipRect.width(), 0.0 ) && qgsDoubleNear( clipRect.height(), 0.0 ) )
+  {
+    if ( partType == GEOS_POINT )
+    {
+      parts->addGeometry( fromGeos( currentPart ) );
+      return;
+    }
+    else
+    {
+      return;
+    }
+  }
+
+  if ( partType == GEOS_MULTILINESTRING || partType == GEOS_MULTIPOLYGON || partType == GEOS_GEOMETRYCOLLECTION )
+  {
+    int partCount = GEOSGetNumGeometries_r( geosinit.ctxt, currentPart );
+    for ( int i = 0; i < partCount; ++i )
+    {
+      subdivideRecursive( GEOSGetGeometryN_r( geosinit.ctxt, currentPart, i ), maxNodes, depth, parts, clipRect );
+    }
+    return;
+  }
+
+  if ( depth > 50 )
+  {
+    parts->addGeometry( fromGeos( currentPart ) );
+    return;
+  }
+
+  int vertexCount = GEOSGetNumCoordinates_r( geosinit.ctxt, currentPart );
+  if ( vertexCount == 0 )
+  {
+    return;
+  }
+  else if ( vertexCount < maxNodes )
+  {
+    parts->addGeometry( fromGeos( currentPart ) );
+    return;
+  }
+
+  // chop clipping rect in half by longest side
+  double width = clipRect.width();
+  double height = clipRect.height();
+  QgsRectangle halfClipRect1 = clipRect;
+  QgsRectangle halfClipRect2 = clipRect;
+  if ( width > height )
+  {
+    halfClipRect1.setXMaximum( clipRect.xMinimum() + width / 2.0 );
+    halfClipRect2.setXMinimum( halfClipRect1.xMaximum() );
+  }
+  else
+  {
+    halfClipRect1.setYMaximum( clipRect.yMinimum() + height / 2.0 );
+    halfClipRect2.setYMinimum( halfClipRect1.yMaximum() );
+  }
+
+  if ( height <= 0 )
+  {
+    halfClipRect1.setYMinimum( halfClipRect1.yMinimum() - DBL_EPSILON );
+    halfClipRect2.setYMinimum( halfClipRect2.yMinimum() - DBL_EPSILON );
+    halfClipRect1.setYMaximum( halfClipRect1.yMaximum() + DBL_EPSILON );
+    halfClipRect2.setYMaximum( halfClipRect2.yMaximum() + DBL_EPSILON );
+  }
+  if ( width <= 0 )
+  {
+    halfClipRect1.setXMinimum( halfClipRect1.xMinimum() - DBL_EPSILON );
+    halfClipRect2.setXMinimum( halfClipRect2.xMinimum() - DBL_EPSILON );
+    halfClipRect1.setXMaximum( halfClipRect1.xMaximum() + DBL_EPSILON );
+    halfClipRect2.setXMaximum( halfClipRect2.xMaximum() + DBL_EPSILON );
+  }
+
+  GEOSGeomScopedPtr clipPart1;
+  clipPart1.reset( GEOSClipByRect_r( geosinit.ctxt, currentPart, halfClipRect1.xMinimum(), halfClipRect1.yMinimum(), halfClipRect1.xMaximum(), halfClipRect1.yMaximum() ) );
+  GEOSGeomScopedPtr clipPart2;
+  clipPart2.reset( GEOSClipByRect_r( geosinit.ctxt, currentPart, halfClipRect2.xMinimum(), halfClipRect2.yMinimum(), halfClipRect2.xMaximum(), halfClipRect2.yMaximum() ) );
+
+  ++depth;
+
+  if ( clipPart1 )
+  {
+    subdivideRecursive( clipPart1.get(), maxNodes, depth, parts, halfClipRect1 );
+  }
+  if ( clipPart2 )
+  {
+    subdivideRecursive( clipPart2.get(), maxNodes, depth, parts, halfClipRect2 );
+  }
+
+  return;
+}
+
+QgsAbstractGeometry *QgsGeos::subdivide( int maxNodes, QString *errorMsg ) const
+{
+  if ( !mGeos )
+  {
+    return nullptr;
+  }
+
+  // minimum allowed max is 8
+  maxNodes = qMax( maxNodes, 8 );
+
+  std::unique_ptr< QgsGeometryCollection > parts = QgsGeometryFactory::createCollectionOfType( mGeometry->wkbType() );
+  try
+  {
+    subdivideRecursive( mGeos, maxNodes, 0, parts.get(), mGeometry->boundingBox() );
+  }
+  CATCH_GEOS_WITH_ERRMSG( nullptr )
+
+  return parts.release();
+}
+
 QgsAbstractGeometry *QgsGeos::combine( const QgsAbstractGeometry &geom, QString *errorMsg ) const
 {
   return overlay( geom, UNION, errorMsg );

diff --git a/src/core/geometry/qgsgeos.h b/src/core/geometry/qgsgeos.h
@@ -23,6 +23,7 @@ email                : marco.hugentobler at sourcepole dot com
 class QgsLineString;
 class QgsPolygonV2;
 class QgsGeometry;
+class QgsGeometryCollection;
 
 /** \ingroup core
  * Does vector analysis using the geos library and handles import, export, exception handling*
@@ -52,6 +53,22 @@ class CORE_EXPORT QgsGeos: public QgsGeometryEngine
      */
     QgsAbstractGeometry *clip( const QgsRectangle &rectangle, QString *errorMsg = nullptr ) const;
 
+    /**
+     * Subdivides the geometry. The returned geometry will be a collection containing subdivided parts
+     * from the original geometry, where no part has more then the specified maximum number of nodes (\a maxNodes).
+     *
+     * This is useful for dividing a complex geometry into less complex parts, which are better able to be spatially
+     * indexed and faster to perform further operations such as intersects on. The returned geometry parts may
+     * not be valid and may contain self-intersections.
+     *
+     * The minimum allowed value for \a maxNodes is 8.
+     *
+     * Curved geometries are not supported.
+     *
+     * \since QGIS 3.0
+     */
+    QgsAbstractGeometry *subdivide( int maxNodes, QString *errorMsg = nullptr ) const;
+
     QgsAbstractGeometry *combine( const QgsAbstractGeometry &geom, QString *errorMsg = nullptr ) const override;
     QgsAbstractGeometry *combine( const QList< QgsAbstractGeometry *> &, QString *errorMsg = nullptr ) const override;
     QgsAbstractGeometry *symDifference( const QgsAbstractGeometry &geom, QString *errorMsg = nullptr ) const override;
@@ -253,6 +270,7 @@ class CORE_EXPORT QgsGeos: public QgsGeometryEngine
     static int lineContainedInLine( const GEOSGeometry *line1, const GEOSGeometry *line2 );
     static int pointContainedInLine( const GEOSGeometry *point, const GEOSGeometry *line );
     static int geomDigits( const GEOSGeometry *geom );
+    void subdivideRecursive( const GEOSGeometry *currentPart, int maxNodes, int depth, QgsGeometryCollection *parts, const QgsRectangle &clipRect ) const;
 };
 
 /// @cond PRIVATE

diff --git a/tests/src/python/test_qgsgeometry.py b/tests/src/python/test_qgsgeometry.py
@@ -4165,6 +4165,25 @@ def testPoint(self):
         self.assertEqual(point_zm.z(), 3)
         self.assertEqual(point_zm.m(), 4)
 
+    def testSubdivide(self):
+        tests = [["LINESTRING (1 1,1 9,9 9,9 1)", 8, "MULTILINESTRING ((1 1,1 9,9 9,9 1))"],
+                 ["Point (1 1)", 8, "MultiPoint ((1 1))"],
+                 ["GeometryCollection ()", 8, "GeometryCollection ()"],
+                 ["LINESTRING (1 1,1 2,1 3,1 4,1 5,1 6,1 7,1 8,1 9)", 8, "MultiLineString ((1 1, 1 2, 1 3, 1 4, 1 5),(1 5, 1 6, 1 7, 1 8, 1 9))"],
+                 ["LINESTRING (1 1,1 2,1 3,1 4,1 5,1 6,1 7,1 8,1 9)", 1,
+                  "MultiLineString ((1 1, 1 2, 1 3, 1 4, 1 5),(1 5, 1 6, 1 7, 1 8, 1 9))"],
+                 ["LINESTRING (1 1,1 2,1 3,1 4,1 5,1 6,1 7,1 8,1 9)", 16,
+                  "MultiLineString ((1 1, 1 2, 1 3, 1 4, 1 5, 1 6, 1 7, 1 8, 1 9))"],
+                 ["POLYGON ((0 0, 0 200, 200 200, 200 0, 0 0),(60 180, 20 180, 20 140, 60 140, 60 180),(180 60, 140 60, 140 20, 180 20, 180 60))", 8, "MultiPolygon (((0 0, 0 100, 100 100, 100 0, 0 0)),((100 0, 100 50, 140 50, 140 20, 150 20, 150 0, 100 0)),((150 0, 150 20, 180 20, 180 50, 200 50, 200 0, 150 0)),((100 50, 100 100, 150 100, 150 60, 140 60, 140 50, 100 50)),((150 60, 150 100, 200 100, 200 50, 180 50, 180 60, 150 60)),((0 100, 0 150, 20 150, 20 140, 50 140, 50 100, 0 100)),((50 100, 50 140, 60 140, 60 150, 100 150, 100 100, 50 100)),((0 150, 0 200, 50 200, 50 180, 20 180, 20 150, 0 150)),((50 180, 50 200, 100 200, 100 150, 60 150, 60 180, 50 180)),((100 100, 100 200, 200 200, 200 100, 100 100)))"]
+                 ]
+        for t in tests:
+            input = QgsGeometry.fromWkt(t[0])
+            o = input.subdivide(t[1])
+            exp = t[2]
+            result = o.exportToWkt()
+            self.assertTrue(compareWkt(result, exp, 0.00001),
+                            "clipped: mismatch Expected:\n{}\nGot:\n{}\n".format(exp, result))
+
     def testClipped(self):
         tests = [["LINESTRING (1 1,1 9,9 9,9 1)", QgsRectangle(0, 0, 10, 10), "LINESTRING (1 1,1 9,9 9,9 1)"],
                  ["LINESTRING (-1 -9,-1 11,9 11)", QgsRectangle(0, 0, 10, 10), "GEOMETRYCOLLECTION ()"],