Expose GEOS constrainedDelaunayTriangulation to QgsGeometry

Allows calling this method when QGIS is built against supported
GEOS versions (>= 3.11)

python/PyQt6/core/auto_generated/geometry/qgsgeometry.sip.in

@@ -1807,7 +1807,24 @@ If ``edgesOnly`` is ``True`` than line string boundary geometries will be return
 instead of polygons.
 An empty geometry will be returned if the diagram could not be calculated.
 
+.. seealso:: :py:func:`constrainedDelaunayTriangulation`
+
 .. versionadded:: 3.0
+%End
+
+    QgsGeometry constrainedDelaunayTriangulation() const throw( QgsNotSupportedException );
+%Docstring
+Returns a constrained Delaunay triangulation for the vertices of the geometry.
+
+An empty geometry will be returned if the triangulation could not be calculated.
+
+This method requires a QGIS build based on GEOS 3.11 or later.
+
+:raises QgsNotSupportedException: on QGIS builds based on GEOS 3.10 or earlier.
+
+.. seealso:: :py:func:`delaunayTriangulation`
+
+.. versionadded:: 3.36
 %End
 
     Qgis::CoverageValidityResult validateCoverage( double gapWidth, QgsGeometry *invalidEdges /Out/ = 0 ) const throw( QgsNotSupportedException );

python/core/auto_generated/geometry/qgsgeometry.sip.in

@@ -1807,7 +1807,24 @@ If ``edgesOnly`` is ``True`` than line string boundary geometries will be return
 instead of polygons.
 An empty geometry will be returned if the diagram could not be calculated.
 
+.. seealso:: :py:func:`constrainedDelaunayTriangulation`
+
 .. versionadded:: 3.0
+%End
+
+    QgsGeometry constrainedDelaunayTriangulation() const throw( QgsNotSupportedException );
+%Docstring
+Returns a constrained Delaunay triangulation for the vertices of the geometry.
+
+An empty geometry will be returned if the triangulation could not be calculated.
+
+This method requires a QGIS build based on GEOS 3.11 or later.
+
+:raises QgsNotSupportedException: on QGIS builds based on GEOS 3.10 or earlier.
+
+.. seealso:: :py:func:`delaunayTriangulation`
+
+.. versionadded:: 3.36
 %End
 
     Qgis::CoverageValidityResult validateCoverage( double gapWidth, QgsGeometry *invalidEdges /Out/ = 0 ) const throw( QgsNotSupportedException );

src/core/geometry/qgsgeometry.cpp

@@ -2641,6 +2641,20 @@ QgsGeometry QgsGeometry::delaunayTriangulation( double tolerance, bool edgesOnly
   return result;
 }
 
+QgsGeometry QgsGeometry::constrainedDelaunayTriangulation() const
+{
+  if ( !d->geometry )
+  {
+    return QgsGeometry();
+  }
+
+  QgsGeos geos( d->geometry.get() );
+  mLastError.clear();
+  QgsGeometry result( geos.constrainedDelaunayTriangulation() );
+  result.mLastError = mLastError;
+  return result;
+}
+
 QgsGeometry QgsGeometry::unionCoverage() const
 {
   if ( !d->geometry )

src/core/geometry/qgsgeometry.h

@@ -1829,10 +1829,26 @@ class CORE_EXPORT QgsGeometry
      * If \a 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.
+     *
+     * \see constrainedDelaunayTriangulation()
      * \since QGIS 3.0
      */
     QgsGeometry delaunayTriangulation( double tolerance = 0.0, bool edgesOnly = false ) const;
 
+    /**
+     * Returns a constrained Delaunay triangulation for the vertices of the geometry.
+     *
+     * An empty geometry will be returned if the triangulation could not be calculated.
+     *
+     * This method requires a QGIS build based on GEOS 3.11 or later.
+     *
+     * \throws QgsNotSupportedException on QGIS builds based on GEOS 3.10 or earlier.
+     * \see delaunayTriangulation()
+     *
+     * \since QGIS 3.36
+     */
+    QgsGeometry constrainedDelaunayTriangulation() const SIP_THROW( QgsNotSupportedException );
+
     /**
      * Analyze a coverage (represented as a collection of polygonal geometry with exactly matching edge
      * geometry) to find places where the assumption of exactly matching edges is not met.

src/core/geometry/qgsgeos.cpp

@@ -3066,6 +3066,40 @@ QgsGeometry QgsGeos::delaunayTriangulation( double tolerance, bool edgesOnly, QS
   CATCH_GEOS_WITH_ERRMSG( QgsGeometry() )
 }
 
+std::unique_ptr<QgsAbstractGeometry> QgsGeos::constrainedDelaunayTriangulation( QString *errorMsg ) const
+{
+#if GEOS_VERSION_MAJOR==3 && GEOS_VERSION_MINOR<11
+  ( void )errorMsg;
+  throw QgsNotSupportedException( QObject::tr( "Calculating constrainedDelaunayTriangulation requires a QGIS build based on GEOS 3.11 or later" ) );
+#else
+  if ( !mGeos )
+  {
+    return nullptr;
+  }
+
+  geos::unique_ptr geos;
+  try
+  {
+    geos.reset( GEOSConstrainedDelaunayTriangulation_r( geosinit()->ctxt, mGeos.get() ) );
+
+    if ( !geos || GEOSisEmpty_r( geosinit()->ctxt, geos.get() ) != 0 )
+    {
+      return nullptr;
+    }
+
+    std::unique_ptr< QgsAbstractGeometry > res = fromGeos( geos.get() );
+    if ( const QgsGeometryCollection *collection = qgsgeometry_cast< const QgsGeometryCollection * >( res.get() ) )
+    {
+      return std::unique_ptr< QgsAbstractGeometry >( collection->extractPartsByType( Qgis::WkbType::Polygon, true ) );
+    }
+    else
+    {
+      return res;
+    }
+  }
+  CATCH_GEOS_WITH_ERRMSG( nullptr )
+#endif
+}
 
 //! Extract coordinates of linestring's endpoints. Returns false on error.
 static bool _linestringEndpoints( const GEOSGeometry *linestring, double &x1, double &y1, double &x2, double &y2 )

src/core/geometry/qgsgeos.h

@@ -547,10 +547,25 @@ class CORE_EXPORT QgsGeos: public QgsGeometryEngine
      * If \a 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.
+     * \see constrainedDelaunayTriangulation()
      * \since QGIS 3.0
      */
     QgsGeometry delaunayTriangulation( double tolerance = 0.0, bool edgesOnly = false, QString *errorMsg = nullptr ) const;
 
+    /**
+     * Returns a constrained Delaunay triangulation for the vertices of the geometry.
+     *
+     * An empty geometry will be returned if the triangulation could not be calculated.
+     *
+     * This method requires a QGIS build based on GEOS 3.11 or later.
+     *
+     * \throws QgsNotSupportedException on QGIS builds based on GEOS 3.10 or earlier.
+     * \see delaunayTriangulation()
+     *
+     * \since QGIS 3.36
+     */
+    std::unique_ptr< QgsAbstractGeometry > constrainedDelaunayTriangulation( QString *errorMsg = nullptr ) const;
+
     /**
      * Returns a possibly concave geometry that encloses the input geometry.
      *
@@ -572,7 +587,7 @@ class CORE_EXPORT QgsGeos: public QgsGeometryEngine
      * \see convexHull()
      * \since QGIS 3.28
      */
-    QgsAbstractGeometry  *concaveHull( double targetPercent, bool allowHoles = false, QString *errorMsg = nullptr ) const;
+    QgsAbstractGeometry *concaveHull( double targetPercent, bool allowHoles = false, QString *errorMsg = nullptr ) const;
 
     /**
      * Analyze a coverage (represented as a collection of polygonal geometry with exactly matching edge

tests/src/python/test_qgsgeometry.py

@@ -7545,6 +7545,31 @@ class TestQgsGeometry(QgisTestCase):
         self.assertEqual(res_isClockwise, False)
         self.assertEqual(res_isCounterClockwise, True)
 
+    @unittest.skipIf(Qgis.geosVersionInt() < 31100, "GEOS 3.11 required")
+    def testConstrainedDelaunayTriangulation(self):
+        """
+        Test QgsGeometry.constrainedDelaunayTriangulation
+        """
+        empty = QgsGeometry()
+        o = empty.constrainedDelaunayTriangulation()
+        self.assertFalse(o)
+        line = QgsGeometry.fromWkt('LineString EMPTY')
+        o = line.constrainedDelaunayTriangulation()
+        self.assertFalse(o)
+
+        input = QgsGeometry.fromWkt("MULTIPOINT ((10 10), (10 20), (20 20))")
+        o = input.constrainedDelaunayTriangulation()
+        self.assertTrue(o.isNull())
+
+        input = QgsGeometry.fromWkt(
+            "POLYGON ((42 30, 41.96 29.61, 41.85 29.23, 41.66 28.89, 41.41 28.59, 41.11 28.34, 40.77 28.15, 40.39 28.04, 40 28, 39.61 28.04, 39.23 28.15, 38.89 28.34, 38.59 28.59, 38.34 28.89, 38.15 29.23, 38.04 29.61, 38 30, 38.04 30.39, 38.15 30.77, 38.34 31.11, 38.59 31.41, 38.89 31.66, 39.23 31.85, 39.61 31.96, 40 32, 40.39 31.96, 40.77 31.85, 41.11 31.66, 41.41 31.41, 41.66 31.11, 41.85 30.77, 41.96 30.39, 42 30))")
+        o = input.constrainedDelaunayTriangulation()
+        o.normalize()
+        self.assertEqual(
+            o.asWkt(2),
+            "MultiPolygon (((41.96 29.61, 41.96 30.39, 42 30, 41.96 29.61)),((41.66 31.11, 41.85 30.77, 41.96 30.39, 41.66 31.11)),((41.66 28.89, 41.96 29.61, 41.85 29.23, 41.66 28.89)),((41.41 31.41, 41.96 30.39, 41.96 29.61, 41.41 31.41)),((41.41 31.41, 41.66 31.11, 41.96 30.39, 41.41 31.41)),((41.41 28.59, 41.96 29.61, 41.66 28.89, 41.41 28.59)),((41.41 28.59, 41.41 31.41, 41.96 29.61, 41.41 28.59)),((40.39 31.96, 41.11 31.66, 41.41 31.41, 40.39 31.96)),((40.39 31.96, 40.77 31.85, 41.11 31.66, 40.39 31.96)),((40.39 28.04, 41.41 28.59, 41.11 28.34, 40.39 28.04)),((40.39 28.04, 41.11 28.34, 40.77 28.15, 40.39 28.04)),((39.61 31.96, 40.39 31.96, 41.41 31.41, 39.61 31.96)),((39.61 31.96, 40 32, 40.39 31.96, 39.61 31.96)),((39.61 28.04, 40.39 28.04, 40 28, 39.61 28.04)),((38.89 31.66, 39.23 31.85, 39.61 31.96, 38.89 31.66)),((38.89 28.34, 39.61 28.04, 39.23 28.15, 38.89 28.34)),((38.59 31.41, 41.41 31.41, 41.41 28.59, 38.59 31.41)),((38.59 31.41, 39.61 31.96, 41.41 31.41, 38.59 31.41)),((38.59 31.41, 38.89 31.66, 39.61 31.96, 38.59 31.41)),((38.59 28.59, 41.41 28.59, 40.39 28.04, 38.59 28.59)),((38.59 28.59, 40.39 28.04, 39.61 28.04, 38.59 28.59)),((38.59 28.59, 39.61 28.04, 38.89 28.34, 38.59 28.59)),((38.59 28.59, 38.59 31.41, 41.41 28.59, 38.59 28.59)),((38.04 30.39, 38.59 31.41, 38.59 28.59, 38.04 30.39)),((38.04 30.39, 38.34 31.11, 38.59 31.41, 38.04 30.39)),((38.04 30.39, 38.15 30.77, 38.34 31.11, 38.04 30.39)),((38.04 29.61, 38.59 28.59, 38.34 28.89, 38.04 29.61)),((38.04 29.61, 38.34 28.89, 38.15 29.23, 38.04 29.61)),((38.04 29.61, 38.04 30.39, 38.59 28.59, 38.04 29.61)),((38 30, 38.04 30.39, 38.04 29.61, 38 30)))"
+        )
+
 
 if __name__ == '__main__':
     unittest.main()