Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

More robust Delaunay Triangulation frame size heuristic #728

Merged
merged 2 commits into from
Nov 8, 2022
Merged

More robust Delaunay Triangulation frame size heuristic #728

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
129552d
More robust Delaunay Triangulation frame size heuristic
dr-jts Nov 8, 2022
d215829
Add some doc
dr-jts Nov 8, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 3 additions & 1 deletion include/geos/triangulate/quadedge/QuadEdgeSubdivision.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -52,6 +52,9 @@ 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.
Expand Down Expand Up @@ -514,4 +517,3 @@ class GEOS_DLL QuadEdgeSubdivision {
} //namespace geos.triangulate.quadedge
} //namespace geos.triangulate
} //namespace goes

16 changes: 5 additions & 11 deletions src/triangulate/quadedge/QuadEdgeSubdivision.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -80,21 +80,15 @@ QuadEdgeSubdivision::createFrame(const geom::Envelope& env)

double deltaX = env.getWidth();
double deltaY = env.getHeight();
double offset = 0.0;
if(deltaX > deltaY) {
offset = deltaX * 10.0;
}
else {
offset = deltaY * 10.0;
}
double offset = std::max(deltaX, deltaY) * FRAME_SIZE_FACTOR;

frameVertex[0] = Vertex((env.getMaxX() + env.getMinX()) / 2.0, env
.getMaxY() + offset);
frameVertex[0] = Vertex((env.getMaxX() + env.getMinX()) / 2.0,
env.getMaxY() + offset);
frameVertex[1] = Vertex(env.getMinX() - offset, env.getMinY() - offset);
frameVertex[2] = Vertex(env.getMaxX() + offset, env.getMinY() - offset);

frameEnv = Envelope(frameVertex[0].getCoordinate(), frameVertex[1]
.getCoordinate());
frameEnv = Envelope(frameVertex[0].getCoordinate(),
frameVertex[1].getCoordinate());
frameEnv.expandToInclude(frameVertex[2].getCoordinate());
}
void
Expand Down
16 changes: 15 additions & 1 deletion tests/unit/triangulate/DelaunayTest.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -187,7 +187,7 @@ void object::test<8>
{
const char* wkt = "MULTIPOINT(-118.3964065 56.0557,-118.396406 56.0475,-118.396407 56.04,-118.3968 56)";
const char* expectedEdges =
"MULTILINESTRING ((-118.3964065 56.0557, -118.396406 56.0475), (-118.396407 56.04, -118.396406 56.0475), (-118.3968 56, -118.396407 56.04))";
"MULTILINESTRING ((-118.3968 56, -118.3964065 56.0557), (-118.3968 56, -118.396407 56.04), (-118.396407 56.04, -118.396406 56.0475), (-118.3964065 56.0557, -118.396406 56.0475), (-118.3968 56, -118.396406 56.0475))";

runDelaunay(wkt, false, expectedEdges, 0.001);
}
Expand Down Expand Up @@ -260,4 +260,18 @@ void object::test<12>
}
}

// failure case due to initial frame size too small
// see https://github.com/libgeos/geos/issues/719, https://github.com/locationtech/jts/pull/931
template<>
template<>
void object::test<13>
()
{
const char* wkt =
"MULTIPOINT ((0 194), (66 151), (203 80), (273 43), (340 0))";
const char* expected =
"GEOMETRYCOLLECTION (POLYGON ((0 194, 66 151, 203 80, 0 194)), POLYGON ((0 194, 203 80, 273 43, 0 194)), POLYGON ((273 43, 203 80, 340 0, 273 43)), POLYGON ((340 0, 203 80, 66 151, 340 0)))";
runDelaunay(wkt, true, expected);
}

} // namespace tut
3 changes: 2 additions & 1 deletion tests/unit/triangulate/VoronoiTest.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -114,7 +114,8 @@ void object::test<1>
std::unique_ptr<QuadEdgeSubdivision> subdiv = builder.getSubdivision();

ensure_equals(subdiv->getTolerance(), 0);
ensure_equals(subdiv->getEnvelope().toString(), "Env[-3540:4020,-3436:4050]");
//-- disable this test, since it depends on Subdiv frame size and so is brittle
//ensure_equals(subdiv->getEnvelope().toString(), "Env[-3540:4020,-3436:4050]");

}
// 1 - Case with a single point
Expand Down
20 changes: 9 additions & 11 deletions tests/unit/triangulate/quadedge/QuadEdgeSubdivisionTest.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -85,31 +85,31 @@ template<>
void object::test<2>
()
{
QuadEdgeSubdivision* subdiv;
auto sites = reader.read("MULTIPOINT ((150 200), (180 270), (275 163))");
auto sites = reader.read("MULTIPOINT ((100 100), (150 200), (200 100))");
auto siteCoords = DelaunayTriangulationBuilder::extractUniqueCoordinates(*sites);
Envelope Env = DelaunayTriangulationBuilder::envelope(*siteCoords);

double expandBy = std::max(Env.getWidth(), Env.getHeight());
Env.expandBy(expandBy);

IncrementalDelaunayTriangulator::VertexList vertices = DelaunayTriangulationBuilder::toVertices(*siteCoords);
subdiv = new quadedge::QuadEdgeSubdivision(Env, 0);
IncrementalDelaunayTriangulator triangulator(subdiv);
std::unique_ptr<QuadEdgeSubdivision> subdiv(new quadedge::QuadEdgeSubdivision(Env, 0));
IncrementalDelaunayTriangulator triangulator(subdiv.get());
triangulator.insertSites(vertices);

//Test for getVoronoiDiagram::
const GeometryFactory& geomFact(*GeometryFactory::getDefaultInstance());
std::unique_ptr<GeometryCollection> polys = subdiv->getVoronoiDiagram(geomFact);
const char* expected_str =
"GEOMETRYCOLLECTION (POLYGON ((-5849.974929324658 2268.0517257497568, -4529.9920486948895 2247.139449440667, 221.20588235294116 210.91176470588235, -684.4227119984187 -2848.644297291955, -5849.974929324658 2268.0517257497568)), POLYGON ((212.5 -3774.5, -684.4227119984187 -2848.644297291955, 221.20588235294116 210.91176470588235, 2448.7167655626645 2188.608343256571, 6235.0370264064295 2248.0370264064295, 212.5 -3774.5)), POLYGON ((-4529.9920486948895 2247.139449440667, 2448.7167655626645 2188.608343256571, 221.20588235294116 210.91176470588235, -4529.9920486948895 2247.139449440667)))";
// std::cout << polys->toString() << std::endl;
//std::cout << polys->toString() << std::endl;
ensure(polys->getNumGeometries()== 3);

auto expected = reader.read(expected_str);
// 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)))"
);
polys->normalize();
expected->normalize();
ensure(polys->equalsExact(expected.get(), 1e-7));
delete subdiv;
// ensure(polys->getCoordinateDimension() == expected->getCoordinateDimension());
}

Expand Down Expand Up @@ -158,5 +158,3 @@ template<> template<> void object::test<3>
}

} // namespace tut