Add constrainedDelaunayTriangulation (#114)

Add constrainedDelaunayTriangulation

Author: jaakkor2

Project.toml

@@ -14,7 +14,7 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 [compat]
 julia = "1.6"
 CEnum = "0.2, 0.3, 0.4"
-GEOS_jll = "3.9"
+GEOS_jll = "3.10"
 GeoInterface = "0.4, 0.5"
 
 [targets]

src/LibGEOS.jl

@@ -43,6 +43,7 @@ export Point,
     snap,
     delaunayTriangulation,
     delaunayTriangulationEdges,
+    constrainedDelaunayTriangulation,
     disjoint,
     touches,
     intersects,

src/geos_functions.jl

@@ -642,6 +642,17 @@ function delaunayTriangulation(
     result
 end
 
+function constrainedDelaunayTriangulation(
+    ptr::GEOSGeom,
+    context::GEOSContext = _context,
+)
+    result = GEOSConstrainedDelaunayTriangulation_r(context.ptr, ptr)
+    if result == C_NULL
+        error("LibGEOS: Error in GEOSConstrainedDelaunayTriangulation")
+    end
+    result
+end
+
 # -----
 # Binary predicates - return 2 on exception, 1 on true, 0 on false
 # -----

src/geos_operations.jl

@@ -182,6 +182,8 @@ for geom in (
         MultiLineString(delaunayTriangulation(obj.ptr, tol, true))
     @eval delaunayTriangulation(obj::$geom, tol::Real = 0.0) =
         GeometryCollection(delaunayTriangulation(obj.ptr, tol, false))
+    @eval constrainedDelaunayTriangulation(obj::$geom) =
+        GeometryCollection(constrainedDelaunayTriangulation(obj.ptr))
 end
 
 sharedPaths(obj1::LineString, obj2::LineString) =

test/test_geos_operations.jl

@@ -70,6 +70,22 @@ end
     g2 = delaunayTriangulationEdges(g1, 2.0)
     equivalent_to_wkt(g2, "MULTILINESTRING ((0 0, 10 10), (0 0, 10 0), (10 0, 10 10))")
 
+    # LibGEOS.constrainedDelaunayTriangulationTest
+    g1 = readgeom("POLYGON EMPTY")
+    g2 = constrainedDelaunayTriangulation(g1)
+    @test isEmpty(g1)
+    @test isEmpty(g2)
+    @test GeoInterface.geotype(g2) == :GeometryCollection
+
+    g1 = readgeom("POINT(0 0)")
+    g2 = constrainedDelaunayTriangulation(g1)
+    @test isEmpty(g2)
+    @test GeoInterface.geotype(g2) == :GeometryCollection
+
+    g1 = readgeom("POLYGON ((10 10, 20 40, 90 90, 90 10, 10 10))")
+    g2 = constrainedDelaunayTriangulation(g1)
+    equivalent_to_wkt(g2, "GEOMETRYCOLLECTION (POLYGON ((10 10, 20 40, 90 10, 10 10)), POLYGON ((90 90, 20 40, 90 10, 90 90)))")
+
     # GEOSDistanceTest
     g1 = readgeom("POINT(10 10)")
     g2 = readgeom("POINT(3 6)")