Backport GH-953 to 3.11 branch (#1035)

* Fix IncrementalDelaunayTriangulator to ensure triangulation boundary is convex (#953)
* Fix compile errors
* Add to NEWS.md

---------

Co-authored-by: Martin Davis <mtnclimb@gmail.com>

NEWS.md

@@ -2,6 +2,7 @@
 xxxx-xx-xx
 
 - Fixes/Improvements:
+  - Fix IncrementalDelaunayTriangulator to ensure triangulation boundary is convex (GH-953, Martin Davis)
   - Fix build on Illumus (GH-971)
   - Fix DiscreteHausdorffDistance for LinearRing (GH-1000, Martin Davis)
   - PointOnSurface crashes with a collection containing a empty linestring (GH-1002, Paul Ramsey)

include/geos/triangulate/IncrementalDelaunayTriangulator.h

@@ -42,6 +42,7 @@ class GEOS_DLL IncrementalDelaunayTriangulator {
 private:
     quadedge::QuadEdgeSubdivision* subdiv;
     bool isUsingTolerance;
+    bool m_isForceConvex;
 
 public:
     /**
@@ -55,6 +56,19 @@ class GEOS_DLL IncrementalDelaunayTriangulator {
 
     typedef std::vector<quadedge::Vertex> VertexList;
 
+    /**
+     * Sets whether the triangulation is forced to have a convex boundary. Because
+     * of the use of a finite-size frame, this condition requires special logic to
+     * enforce. The default is true, since this is a requirement for some uses of
+     * Delaunay Triangulations (such as Concave Hull generation). However, forcing
+     * the triangulation boundary to be convex may cause the overall frame
+     * triangulation to be non-Delaunay. This can cause a problem for Voronoi
+     * generation, so the logic can be disabled via this method.
+     * 
+     * @param isForceConvex true if the triangulation boundary is forced to be convex
+     */
+    void forceConvex(bool isForceConvex);
+
     /**
      * Inserts all sites in a collection. The inserted vertices <b>MUST</b> be
      * unique up to the provided tolerance value. (i.e. no two vertices should be
@@ -77,6 +91,14 @@ class GEOS_DLL IncrementalDelaunayTriangulator {
      * @return a quadedge containing the inserted vertex
      */
     quadedge::QuadEdge& insertSite(const quadedge::Vertex& v);
+
+private:
+
+    bool isConcaveBoundary(const quadedge::QuadEdge& e);
+
+    bool isConcaveAtOrigin(const quadedge::QuadEdge& e);
+
+    bool isBetweenFrameAndInserted(const quadedge::QuadEdge& e, const quadedge::Vertex& vInsert);
 };
 
 } //namespace geos.triangulate

include/geos/triangulate/quadedge/QuadEdgeSubdivision.h

@@ -50,11 +50,6 @@ namespace quadedge { //geos.triangulate.quadedge
 
 class TriangleVisitor;
 
-const double EDGE_COINCIDENCE_TOL_FACTOR = 1000;
-
-//-- Frame size factor for initializing subdivision frame.  Larger is more robust
-const double FRAME_SIZE_FACTOR = 100;
-
 /** \brief
  * A class that contains the [QuadEdges](@ref QuadEdge) representing a planar
  * subdivision that models a triangulation.

src/triangulate/IncrementalDelaunayTriangulator.cpp

@@ -21,18 +21,27 @@
 #include <geos/triangulate/quadedge/QuadEdge.h>
 #include <geos/triangulate/quadedge/QuadEdgeSubdivision.h>
 #include <geos/triangulate/quadedge/LocateFailureException.h>
+#include <geos/algorithm/Orientation.h>
 
 namespace geos {
 namespace triangulate { //geos.triangulate
 
+using namespace algorithm;
 using namespace quadedge;
 
 IncrementalDelaunayTriangulator::IncrementalDelaunayTriangulator(
     QuadEdgeSubdivision* p_subdiv) :
-    subdiv(p_subdiv), isUsingTolerance(p_subdiv->getTolerance() > 0.0)
+    subdiv(p_subdiv), isUsingTolerance(p_subdiv->getTolerance() > 0.0),
+    m_isForceConvex(true)
 {
 }
 
+void 
+IncrementalDelaunayTriangulator::forceConvex(bool isForceConvex) 
+{
+    m_isForceConvex = isForceConvex;
+}
+
 void
 IncrementalDelaunayTriangulator::insertSites(const VertexList& vertices)
 {
@@ -82,25 +91,76 @@ IncrementalDelaunayTriangulator::insertSite(const Vertex& v)
     }
     while(&e->lNext() != startEdge);
 
-
-    // Examine suspect edges to ensure that the Delaunay condition
-    // is satisfied.
+    /**
+     * Examine suspect edges to ensure that the Delaunay condition is satisfied.
+     * If it is not, flip the edge and continue scanning.
+     * 
+     * Since the frame is not infinitely far away,
+     * edges which touch the frame or are adjacent to it require special logic
+     * to ensure the inner triangulation maintains a convex boundary.
+     */
     for(;;) {
+        //-- general case - flip if vertex is in circumcircle
         QuadEdge* t = &e->oPrev();
-        if(t->dest().rightOf(*e) &&
-                v.isInCircle(e->orig(), t->dest(), e->dest())) {
+        bool doFlip = t->dest().rightOf(*e) &&
+                v.isInCircle(e->orig(), t->dest(), e->dest());
+
+        if (m_isForceConvex) {
+            //-- special cases to ensure triangulation boundary is convex
+            if (isConcaveBoundary(*e)) {
+            //-- flip if the triangulation boundary is concave
+                doFlip = true;
+            }
+            else if (isBetweenFrameAndInserted(*e, v)) {
+            //-- don't flip if edge lies between the inserted vertex and a frame vertex
+                doFlip = false;
+            }
+        }
+
+        if (doFlip) {
             QuadEdge::swap(*e);
             e = &e->oPrev();
+            continue;
         }
-        else if(&e->oNext() == startEdge) {
+        if (&e->oNext() == startEdge) {
             return *base; // no more suspect edges.
         }
-        else {
-            e = &e->oNext().lPrev();
-        }
+        //-- check next edge
+        e = &e->oNext().lPrev();
     }
 }
 
+bool 
+IncrementalDelaunayTriangulator::isConcaveBoundary(const QuadEdge& e) 
+{
+    if (subdiv->isFrameVertex(e.dest())) {
+        return isConcaveAtOrigin(e);
+    }
+    if (subdiv->isFrameVertex(e.orig())) {
+        return isConcaveAtOrigin(e.sym());
+    }
+    return false;
+}
+
+bool 
+IncrementalDelaunayTriangulator::isConcaveAtOrigin(const QuadEdge& e) 
+{
+    Coordinate p = e.orig().getCoordinate();
+    Coordinate pp = e.oPrev().dest().getCoordinate();
+    Coordinate pn = e.oNext().dest().getCoordinate();
+    bool isConcave = Orientation::COUNTERCLOCKWISE == Orientation::index(pp, pn, p);
+    return isConcave;
+}
+
+bool 
+IncrementalDelaunayTriangulator::isBetweenFrameAndInserted(const QuadEdge& e, const Vertex& vInsert) 
+{
+    const Vertex v1 = e.oNext().dest();
+    const Vertex v2 = e.oPrev().dest();
+    return (v1.getCoordinate() == vInsert.getCoordinate() && subdiv->isFrameVertex(v2))
+        || (v2.getCoordinate() == vInsert.getCoordinate() && subdiv->isFrameVertex(v1));
+}
+
 } //namespace geos.triangulate
-} //namespace goes
+} //namespace geos
 

src/triangulate/VoronoiDiagramBuilder.cpp

@@ -94,6 +94,11 @@ VoronoiDiagramBuilder::create()
 
     subdiv.reset(new quadedge::QuadEdgeSubdivision(diagramEnv, tolerance));
     IncrementalDelaunayTriangulator triangulator(subdiv.get());
+    /**
+     * Avoid creating very narrow triangles along triangulation boundary.
+     * These otherwise can cause malformed Voronoi cells.
+     */
+    triangulator.forceConvex(false);
     triangulator.insertSites(vertices);
 }
 

src/triangulate/quadedge/QuadEdgeSubdivision.cpp

@@ -49,6 +49,11 @@ namespace geos {
 namespace triangulate { //geos.triangulate
 namespace quadedge { //geos.triangulate.quadedge
 
+const double EDGE_COINCIDENCE_TOL_FACTOR = 1000;
+
+//-- Frame size factor for initializing subdivision frame.  Larger is more robust
+const double FRAME_SIZE_FACTOR = 10;
+
 void
 QuadEdgeSubdivision::getTriangleEdges(const QuadEdge& startQE,
                                       const QuadEdge* triEdge[3])

tests/unit/algorithm/hull/ConcaveHullTest.cpp

@@ -243,6 +243,28 @@ void object::test<15>()
        0, "POLYGON ((20 90, 40 96, 56 95, 70 80, 80 90, 90 70, 80 60, 95 45, 80 40, 70 20, 90 20, 80 10, 60 15, 45 5, 40 20, 40 80, 15 45, 21 30, 20 10, 10 20, 5 40, 11 60, 20 70, 20 90))" );
 }
 
+//------------------------------------------------
+
+// These tests test that the computed Delaunay triangulation is correct
+// See https://github.com/locationtech/jts/pull/1004
+
+// testRobust_GEOS946
+template<>
+template<>
+void object::test<20>()
+{
+    checkHullByLengthRatio("MULTIPOINT ((113.56577197798602 22.80081530883069),(113.565723279387 22.800815316487014),(113.56571548761124 22.80081531771092),(113.56571548780202 22.800815317674463),(113.56577197817877 22.8008153088047),(113.56577197798602 22.80081530883069))", 
+       0.75, "POLYGON ((113.56571548761124 22.80081531771092, 113.565723279387 22.800815316487014, 113.56577197798602 22.80081530883069, 113.56577197817877 22.8008153088047, 113.56571548780202 22.800815317674463, 113.56571548761124 22.80081531771092))" );
+}
+
+// testRobust_GEOS946_2
+template<>
+template<>
+void object::test<21>()
+{
+    checkHullByLengthRatio("MULTIPOINT ((584245.72096874 7549593.72686167), (584251.71398371 7549594.01629478), (584242.72446125 7549593.58214511), (584230.73978847 7549592.9760418), (584233.73581213 7549593.13045099), (584236.7318358 7549593.28486019), (584239.72795377 7549593.43742855), (584227.74314188 7549592.83423486))", 
+       0.75, "POLYGON ((584227.74314188 7549592.83423486, 584239.72795377 7549593.43742855, 584242.72446125 7549593.58214511, 584245.72096874 7549593.72686167, 584251.71398371 7549594.01629478, 584230.73978847 7549592.9760418, 584227.74314188 7549592.83423486))" );
+}
 
 
 } // namespace tut

tests/unit/triangulate/DelaunayTest.cpp

@@ -5,6 +5,9 @@
 // tut
 #include <tut/tut.hpp>
 // geos
+#include <geos/algorithm/ConvexHull.h>
+// #include <geos/coverage/CoverageUnion.h>
+#include <geos/operation/overlayng/CoverageUnion.h>
 #include <geos/triangulate/quadedge/QuadEdge.h>
 #include <geos/triangulate/quadedge/QuadEdgeSubdivision.h>
 #include <geos/triangulate/IncrementalDelaunayTriangulator.h>
@@ -39,6 +42,28 @@ typedef group::object object;
 group test_incdelaunaytri_group("geos::triangulate::Delaunay");
 
 //helper function for running triangulation
+void
+checkDelaunayHull(const char* sitesWkt)
+{
+    WKTReader reader;
+    auto sites = reader.read(sitesWkt);
+
+    DelaunayTriangulationBuilder builder;
+    const GeometryFactory& geomFact(*GeometryFactory::getDefaultInstance());
+    builder.setSites(*sites);
+    std::unique_ptr<Geometry> tris = builder.getTriangles(geomFact);
+
+    // std::unique_ptr<Geometry> hullTris = geos::coverage::CoverageUnion::Union(tris.get());
+    std::unique_ptr<Geometry> hullTris = geos::operation::overlayng::CoverageUnion::geomunion(tris.get());
+    std::unique_ptr<Geometry> hullSites = sites->convexHull();
+
+    //std::cout << "hullTris: " << hullTris->toString() << std::endl;
+    //std::cout << "hullSites: " << hullSites->toString() << std::endl;
+
+    //-- use topological equality, because there may be collinear vertices in the union
+    ensure(hullTris->equals(hullSites.get()));
+}
+
 void
 runDelaunay(const char* sitesWkt, bool computeTriangles, const char* expectedWkt, double tolerance = 0.0)
 {
@@ -271,4 +296,78 @@ void object::test<13>
     runDelaunay(wkt, true, expected);
 }
 
+// see https://github.com/libgeos/geos/issues/719
+// testNarrow_GEOS_719()
+template<>
+template<>
+void object::test<14>()
+{
+    const char* wkt = "MULTIPOINT ((1139294.6389832513 8201313.534695469), (1139360.8549531854 8201271.189805277), (1139497.5995843115 8201199.995542546), (1139567.7837303514 8201163.348533507), (1139635.3942210067 8201119.902527407))";
+    const char* expected = "GEOMETRYCOLLECTION (POLYGON ((1139294.6389832513 8201313.534695469, 1139360.8549531854 8201271.189805277, 1139497.5995843115 8201199.995542546, 1139294.6389832513 8201313.534695469)), POLYGON ((1139294.6389832513 8201313.534695469, 1139497.5995843115 8201199.995542546, 1139567.7837303514 8201163.348533507, 1139294.6389832513 8201313.534695469)), POLYGON ((1139567.7837303514 8201163.348533507, 1139497.5995843115 8201199.995542546, 1139635.3942210067 8201119.902527407, 1139567.7837303514 8201163.348533507)), POLYGON ((1139635.3942210067 8201119.902527407, 1139497.5995843115 8201199.995542546, 1139360.8549531854 8201271.189805277, 1139635.3942210067 8201119.902527407)))";
+    runDelaunay(wkt, true, expected);
+}
+
+// testNarrowTriangle()
+template<>
+template<>
+void object::test<15>()
+  {
+    const char*  wkt = "MULTIPOINT ((100 200), (200 190), (300 200))";
+    const char*  expected = "GEOMETRYCOLLECTION (POLYGON ((100 200, 300 200, 200 190, 100 200)))";
+    runDelaunay(wkt, true, expected);
+  }
+
+// seee https://github.com/locationtech/jts/issues/477
+// testNarrow_GH477_1()
+template<>
+template<>
+void object::test<16>()
+{
+    const char*  wkt = "MULTIPOINT ((0 0), (1 0), (-1 0.05), (0 0))";
+    const char*  expected = "GEOMETRYCOLLECTION (POLYGON ((-1 0.05, 1 0, 0 0, -1 0.05)))";
+    runDelaunay(wkt, true, expected);
+}
+
+// see https://github.com/locationtech/jts/issues/477
+// testNarrow_GH477_2()
+template<>
+template<>
+void object::test<17>()
+{
+    const char*  wkt = "MULTIPOINT ((0 0), (0 486), (1 486), (1 22), (2 22), (2 0))";
+    const char*  expected = "GEOMETRYCOLLECTION (POLYGON ((0 0, 0 486, 1 22, 0 0)), POLYGON ((0 0, 1 22, 2 0, 0 0)), POLYGON ((0 486, 1 486, 1 22, 0 486)), POLYGON ((1 22, 1 486, 2 22, 1 22)), POLYGON ((1 22, 2 22, 2 0, 1 22)))";
+    runDelaunay(wkt, true, expected);
+}
+
+// see https://github.com/libgeos/geos/issues/946
+// testNarrow_GEOS_946()
+template<>
+template<>
+void object::test<18>()
+{
+    const char*  wkt = "MULTIPOINT ((113.56577197798602 22.80081530883069),(113.565723279387 22.800815316487014),(113.56571548761124 22.80081531771092),(113.56571548780202 22.800815317674463),(113.56577197817877 22.8008153088047),(113.56577197798602 22.80081530883069))";
+    const char*  expected = "GEOMETRYCOLLECTION (POLYGON ((113.56571548761124 22.80081531771092, 113.565723279387 22.800815316487014, 113.56571548780202 22.800815317674463, 113.56571548761124 22.80081531771092)), POLYGON ((113.56571548780202 22.800815317674463, 113.565723279387 22.800815316487014, 113.56577197817877 22.8008153088047, 113.56571548780202 22.800815317674463)), POLYGON ((113.565723279387 22.800815316487014, 113.56577197798602 22.80081530883069, 113.56577197817877 22.8008153088047, 113.565723279387 22.800815316487014)))";
+    runDelaunay(wkt, true, expected);
+}
+
+// see https://github.com/shapely/shapely/issues/1873
+// testNarrow_Shapely_1873()
+template<>
+template<>
+void object::test<19>()
+{
+    const char*  wkt = "MULTIPOINT ((584245.72096874 7549593.72686167), (584251.71398371 7549594.01629478), (584242.72446125 7549593.58214511), (584230.73978847 7549592.9760418), (584233.73581213 7549593.13045099), (584236.7318358 7549593.28486019), (584239.72795377 7549593.43742855), (584227.74314188 7549592.83423486))";
+    const char*  expected = "GEOMETRYCOLLECTION (POLYGON ((584227.74314188 7549592.83423486, 584233.73581213 7549593.13045099, 584230.73978847 7549592.9760418, 584227.74314188 7549592.83423486)), POLYGON ((584227.74314188 7549592.83423486, 584236.7318358 7549593.28486019, 584233.73581213 7549593.13045099, 584227.74314188 7549592.83423486)), POLYGON ((584227.74314188 7549592.83423486, 584239.72795377 7549593.43742855, 584236.7318358 7549593.28486019, 584227.74314188 7549592.83423486)), POLYGON ((584230.73978847 7549592.9760418, 584233.73581213 7549593.13045099, 584245.72096874 7549593.72686167, 584230.73978847 7549592.9760418)), POLYGON ((584230.73978847 7549592.9760418, 584245.72096874 7549593.72686167, 584251.71398371 7549594.01629478, 584230.73978847 7549592.9760418)), POLYGON ((584233.73581213 7549593.13045099, 584236.7318358 7549593.28486019, 584242.72446125 7549593.58214511, 584233.73581213 7549593.13045099)), POLYGON ((584233.73581213 7549593.13045099, 584242.72446125 7549593.58214511, 584245.72096874 7549593.72686167, 584233.73581213 7549593.13045099)), POLYGON ((584236.7318358 7549593.28486019, 584239.72795377 7549593.43742855, 584242.72446125 7549593.58214511, 584236.7318358 7549593.28486019)))";
+    runDelaunay(wkt, true, expected);
+}
+
+// testNarrowPoints()
+template<>
+template<>
+void object::test<20>()
+{
+    const char*  wkt = "MULTIPOINT ((2 204), (3 66), (1 96), (0 236), (3 173), (2 114), (3 201), (0 46), (1 181))";
+    checkDelaunayHull(wkt);
+}
+
 } // namespace tut

tests/unit/triangulate/quadedge/QuadEdgeSubdivisionTest.cpp

@@ -95,6 +95,11 @@ void object::test<2>
     IncrementalDelaunayTriangulator::VertexList vertices = DelaunayTriangulationBuilder::toVertices(*siteCoords);
     std::unique_ptr<QuadEdgeSubdivision> subdiv(new quadedge::QuadEdgeSubdivision(Env, 0));
     IncrementalDelaunayTriangulator triangulator(subdiv.get());
+        /**
+     * Avoid creating very narrow triangles along triangulation boundary.
+     * These otherwise can cause malformed Voronoi cells.
+     */
+    triangulator.forceConvex(false);
     triangulator.insertSites(vertices);
 
     //Test for getVoronoiDiagram::
@@ -105,11 +110,11 @@ void object::test<2>
 
     // return value depends on subdivision frame vertices
     auto expected = reader.read(
-                "GEOMETRYCOLLECTION (POLYGON ((-45175 15275, -30075 15250, 150 137.5, 150 -30050, -45175 15275)), POLYGON ((-30075 15250, 30375 15250, 150 137.5, -30075 15250)), POLYGON ((30375 15250, 45475 15275, 150 -30050, 150 137.5, 30375 15250)))"
+                "GEOMETRYCOLLECTION (POLYGON ((150 -3050, 150 137.5, 3375 1750, 4975 1775, 150 -3050)), POLYGON ((-4675 1775, -3075 1750, 150 137.5, 150 -3050, -4675 1775)), POLYGON ((-3075 1750, 3375 1750, 150 137.5, -3075 1750)))"
     );
     polys->normalize();
     expected->normalize();
-    ensure(polys->equalsExact(expected.get(), 1e-7));
+    ensure(polys->toString(), polys->equalsExact(expected.get(), 1e-7));
 //		ensure(polys->getCoordinateDimension() == expected->getCoordinateDimension());
 }