[geos-commits] [SCM] GEOS branch 3.12 updated. 4af2a56bc7fb93bdc6567b8eb3bf184c04f73cb9
git at osgeo.org
git at osgeo.org
Tue Sep 5 18:42:36 PDT 2023
This is an automated email from the git hooks/post-receive script. It was
generated because a ref change was pushed to the repository containing
the project "GEOS".
The branch, 3.12 has been updated
via 4af2a56bc7fb93bdc6567b8eb3bf184c04f73cb9 (commit)
from 87bf94a3227626bba6a94fe65abd3b6bd2344fa6 (commit)
Those revisions listed above that are new to this repository have
not appeared on any other notification email; so we list those
revisions in full, below.
- Log -----------------------------------------------------------------
commit 4af2a56bc7fb93bdc6567b8eb3bf184c04f73cb9
Author: Martin Davis <mtnclimb at gmail.com>
Date: Tue Sep 5 16:27:34 2023 -0700
Fix IncrementalDelaunayTriangulator to ensure triangulation boundary is convex (#953)
diff --git a/include/geos/triangulate/IncrementalDelaunayTriangulator.h b/include/geos/triangulate/IncrementalDelaunayTriangulator.h
index 95f14d243..942d62d33 100644
--- a/include/geos/triangulate/IncrementalDelaunayTriangulator.h
+++ b/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 @@ public:
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 @@ public:
* @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
diff --git a/include/geos/triangulate/quadedge/QuadEdgeSubdivision.h b/include/geos/triangulate/quadedge/QuadEdgeSubdivision.h
index edd2aa3ca..24492d51c 100644
--- a/include/geos/triangulate/quadedge/QuadEdgeSubdivision.h
+++ b/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.
diff --git a/src/triangulate/IncrementalDelaunayTriangulator.cpp b/src/triangulate/IncrementalDelaunayTriangulator.cpp
index d4cf79d1a..a742b7ee1 100644
--- a/src/triangulate/IncrementalDelaunayTriangulator.cpp
+++ b/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();
}
}
-} //namespace geos.triangulate
-} //namespace goes
+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 geos
diff --git a/src/triangulate/VoronoiDiagramBuilder.cpp b/src/triangulate/VoronoiDiagramBuilder.cpp
index f1244acc6..d4b09c8c5 100644
--- a/src/triangulate/VoronoiDiagramBuilder.cpp
+++ b/src/triangulate/VoronoiDiagramBuilder.cpp
@@ -104,6 +104,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);
}
diff --git a/src/triangulate/quadedge/QuadEdgeSubdivision.cpp b/src/triangulate/quadedge/QuadEdgeSubdivision.cpp
index 511ff87e4..ba92a5e54 100644
--- a/src/triangulate/quadedge/QuadEdgeSubdivision.cpp
+++ b/src/triangulate/quadedge/QuadEdgeSubdivision.cpp
@@ -47,6 +47,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])
diff --git a/tests/unit/algorithm/hull/ConcaveHullTest.cpp b/tests/unit/algorithm/hull/ConcaveHullTest.cpp
index f805a3d79..3d47bbc93 100644
--- a/tests/unit/algorithm/hull/ConcaveHullTest.cpp
+++ b/tests/unit/algorithm/hull/ConcaveHullTest.cpp
@@ -298,6 +298,28 @@ void object::test<19>()
"POLYGON ((20 90, 40 96, 56 95, 80 90, 90 70, 95 45, 90 20, 80 10, 60 15, 45 5, 20 10, 10 20, 5 40, 11 60, 20 70, 20 90), (40 80, 15 45, 21 30, 40 20, 70 20, 80 40, 80 60, 70 80, 40 80))" );
}
+//------------------------------------------------
+
+// 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
diff --git a/tests/unit/triangulate/DelaunayTest.cpp b/tests/unit/triangulate/DelaunayTest.cpp
index 5bae9331f..0b324c7df 100644
--- a/tests/unit/triangulate/DelaunayTest.cpp
+++ b/tests/unit/triangulate/DelaunayTest.cpp
@@ -5,6 +5,8 @@
// tut
#include <tut/tut.hpp>
// geos
+#include <geos/algorithm/ConvexHull.h>
+#include <geos/coverage/CoverageUnion.h>
#include <geos/triangulate/quadedge/QuadEdge.h>
#include <geos/triangulate/quadedge/QuadEdgeSubdivision.h>
#include <geos/triangulate/IncrementalDelaunayTriangulator.h>
@@ -39,6 +41,27 @@ 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> 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)
{
@@ -273,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
diff --git a/tests/unit/triangulate/quadedge/QuadEdgeSubdivisionTest.cpp b/tests/unit/triangulate/quadedge/QuadEdgeSubdivisionTest.cpp
index c6bd9ca1e..9521a8a7e 100644
--- a/tests/unit/triangulate/quadedge/QuadEdgeSubdivisionTest.cpp
+++ b/tests/unit/triangulate/quadedge/QuadEdgeSubdivisionTest.cpp
@@ -94,6 +94,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::
@@ -104,11 +109,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());
}
-----------------------------------------------------------------------
Summary of changes:
.../triangulate/IncrementalDelaunayTriangulator.h | 22 +++++
.../triangulate/quadedge/QuadEdgeSubdivision.h | 5 --
.../IncrementalDelaunayTriangulator.cpp | 82 +++++++++++++++---
src/triangulate/VoronoiDiagramBuilder.cpp | 5 ++
src/triangulate/quadedge/QuadEdgeSubdivision.cpp | 5 ++
tests/unit/algorithm/hull/ConcaveHullTest.cpp | 22 +++++
tests/unit/triangulate/DelaunayTest.cpp | 97 ++++++++++++++++++++++
.../quadedge/QuadEdgeSubdivisionTest.cpp | 9 +-
8 files changed, 229 insertions(+), 18 deletions(-)
hooks/post-receive
--
GEOS
More information about the geos-commits
mailing list