From 45e51c2a550843f086bacbeaa262cf1aba129cb7 Mon Sep 17 00:00:00 2001 From: Nyall Dawson Date: Wed, 7 Feb 2024 10:50:59 +1000 Subject: [PATCH] Expose GEOS constrainedDelaunayTriangulation to QgsGeometry Allows calling this method when QGIS is built against supported GEOS versions (>= 3.11) --- .../geometry/qgsgeometry.sip.in | 17 ++++++++++ .../geometry/qgsgeometry.sip.in | 17 ++++++++++ src/core/geometry/qgsgeometry.cpp | 14 ++++++++ src/core/geometry/qgsgeometry.h | 16 +++++++++ src/core/geometry/qgsgeos.cpp | 34 +++++++++++++++++++ src/core/geometry/qgsgeos.h | 17 +++++++++- tests/src/python/test_qgsgeometry.py | 25 ++++++++++++++ 7 files changed, 139 insertions(+), 1 deletion(-) diff --git a/python/PyQt6/core/auto_generated/geometry/qgsgeometry.sip.in b/python/PyQt6/core/auto_generated/geometry/qgsgeometry.sip.in index 3d6608153ee..dd78286c8f9 100644 --- a/python/PyQt6/core/auto_generated/geometry/qgsgeometry.sip.in +++ b/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 ); diff --git a/python/core/auto_generated/geometry/qgsgeometry.sip.in b/python/core/auto_generated/geometry/qgsgeometry.sip.in index 3d6608153ee..dd78286c8f9 100644 --- a/python/core/auto_generated/geometry/qgsgeometry.sip.in +++ b/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 ); diff --git a/src/core/geometry/qgsgeometry.cpp b/src/core/geometry/qgsgeometry.cpp index 4c4e95dcbaf..20fc41f4d0f 100644 --- a/src/core/geometry/qgsgeometry.cpp +++ b/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 ) diff --git a/src/core/geometry/qgsgeometry.h b/src/core/geometry/qgsgeometry.h index 9b37cc99406..e591a9dc57f 100644 --- a/src/core/geometry/qgsgeometry.h +++ b/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. diff --git a/src/core/geometry/qgsgeos.cpp b/src/core/geometry/qgsgeos.cpp index 02da31501f1..f05b834fedf 100644 --- a/src/core/geometry/qgsgeos.cpp +++ b/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 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 ) diff --git a/src/core/geometry/qgsgeos.h b/src/core/geometry/qgsgeos.h index 736e8ee09fb..fab48ca1c2b 100644 --- a/src/core/geometry/qgsgeos.h +++ b/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 diff --git a/tests/src/python/test_qgsgeometry.py b/tests/src/python/test_qgsgeometry.py index 1dcf2921793..ed6a9cccdfe 100644 --- a/tests/src/python/test_qgsgeometry.py +++ b/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()