diff --git a/NEWS.md b/NEWS.md index c368a0ad09..bd44ea16f4 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 95f14d243a..942d62d33c 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 @@ class GEOS_DLL IncrementalDelaunayTriangulator { typedef std::vector 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 MUST 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 diff --git a/include/geos/triangulate/quadedge/QuadEdgeSubdivision.h b/include/geos/triangulate/quadedge/QuadEdgeSubdivision.h index edd2aa3ca5..24492d51c0 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 d4cf79d1a3..a742b7ee12 100644 --- a/src/triangulate/IncrementalDelaunayTriangulator.cpp +++ b/src/triangulate/IncrementalDelaunayTriangulator.cpp @@ -21,18 +21,27 @@ #include #include #include +#include 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 diff --git a/src/triangulate/VoronoiDiagramBuilder.cpp b/src/triangulate/VoronoiDiagramBuilder.cpp index c9cfc60849..01a711bb2b 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 9cbd953e56..1c7bc8b968 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 170762b03c..830b1d8334 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 e0ddaddb84..1faf4322aa 100644 --- a/tests/unit/triangulate/DelaunayTest.cpp +++ b/tests/unit/triangulate/DelaunayTest.cpp @@ -5,6 +5,9 @@ // tut #include // geos +#include +// #include +#include #include #include #include @@ -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 tris = builder.getTriangles(geomFact); + + // std::unique_ptr hullTris = geos::coverage::CoverageUnion::Union(tris.get()); + std::unique_ptr hullTris = geos::operation::overlayng::CoverageUnion::geomunion(tris.get()); + std::unique_ptr 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 66049d9b93..865a703f79 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 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()); }