[geos-commits] [SCM] GEOS branch 3.11 updated. 16938a636c628f3e5ec81b2719bfd0929cf3609a
git at osgeo.org
git at osgeo.org
Wed Jan 31 12:51:10 PST 2024
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.11 has been updated
via 16938a636c628f3e5ec81b2719bfd0929cf3609a (commit)
from 2d2febc50e3e61af5d0ebe56d689613417d7a36a (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 16938a636c628f3e5ec81b2719bfd0929cf3609a
Author: Mike Taves <mwtoews at gmail.com>
Date: Thu Feb 1 09:50:36 2024 +1300
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 at gmail.com>
diff --git a/NEWS.md b/NEWS.md
index c368a0ad0..bd44ea16f 100644
--- a/NEWS.md
+++ b/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)
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 c9cfc6084..01a711bb2 100644
--- a/src/triangulate/VoronoiDiagramBuilder.cpp
+++ b/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);
}
diff --git a/src/triangulate/quadedge/QuadEdgeSubdivision.cpp b/src/triangulate/quadedge/QuadEdgeSubdivision.cpp
index 9cbd953e5..1c7bc8b96 100644
--- a/src/triangulate/quadedge/QuadEdgeSubdivision.cpp
+++ b/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])
diff --git a/tests/unit/algorithm/hull/ConcaveHullTest.cpp b/tests/unit/algorithm/hull/ConcaveHullTest.cpp
index 170762b03..830b1d833 100644
--- a/tests/unit/algorithm/hull/ConcaveHullTest.cpp
+++ b/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
diff --git a/tests/unit/triangulate/DelaunayTest.cpp b/tests/unit/triangulate/DelaunayTest.cpp
index e0ddaddb8..1faf4322a 100644
--- a/tests/unit/triangulate/DelaunayTest.cpp
+++ b/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
diff --git a/tests/unit/triangulate/quadedge/QuadEdgeSubdivisionTest.cpp b/tests/unit/triangulate/quadedge/QuadEdgeSubdivisionTest.cpp
index 66049d9b9..865a703f7 100644
--- a/tests/unit/triangulate/quadedge/QuadEdgeSubdivisionTest.cpp
+++ b/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());
}
-----------------------------------------------------------------------
Summary of changes:
NEWS.md | 1 +
.../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 | 99 ++++++++++++++++++++++
.../quadedge/QuadEdgeSubdivisionTest.cpp | 9 +-
9 files changed, 232 insertions(+), 18 deletions(-)
hooks/post-receive
--
GEOS
More information about the geos-commits
mailing list