Skip to content

Commit

Permalink
Fix bug in handling hanging edges when triangulating pseudo-polygons
Browse files Browse the repository at this point in the history 
Add a regression test
  • Loading branch information
@artem-ogre
artem-ogre committed Nov 13, 2023
1 parent 46d12c7 commit ffb69c6
Show file tree
Hide file tree
Showing 4 changed files with 67 additions and 63 deletions.
20 changes: 10 additions & 10 deletions CDT/include/Triangulation.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -520,7 +520,7 @@ class CDT_EXPORT Triangulation

/// State for an iteration of triangulate pseudo-polygon
typedef tuple<IndexSizeType, IndexSizeType, TriInd, TriInd, Index>
TriangulatePseudopolygonTask;
TriangulatePseudoPolygonTask;

/**
* Insert an edge into constraint Delaunay triangulation
Expand All @@ -538,7 +538,7 @@ class CDT_EXPORT Triangulation
Edge edge,
Edge originalEdge,
EdgeVec& remaining,
std::vector<TriangulatePseudopolygonTask>& tppIterations);
std::vector<TriangulatePseudoPolygonTask>& tppIterations);

/**
* Insert an edge or its part into constraint Delaunay triangulation
Expand All @@ -556,7 +556,7 @@ class CDT_EXPORT Triangulation
Edge edge,
Edge originalEdge,
EdgeVec& remaining,
std::vector<TriangulatePseudopolygonTask>& tppIterations);
std::vector<TriangulatePseudoPolygonTask>& tppIterations);

/// State for iteration of conforming to edge
typedef tuple<Edge, EdgeVec, BoundaryOverlapCount> ConformToEdgeTask;
Expand Down Expand Up @@ -629,16 +629,16 @@ class CDT_EXPORT Triangulation
VertInd iVedge1,
VertInd iVedge2,
TriInd newNeighbor);
void triangulatePseudopolygon(
void triangulatePseudoPolygon(
const std::vector<VertInd>& poly,
const std::vector<TriInd>& outerTris,
unordered_map<Edge, TriInd>& outerTris,
TriInd iT,
TriInd iN,
std::vector<TriangulatePseudopolygonTask>& iterations);
void triangulatePseudopolygonIteration(
std::vector<TriangulatePseudoPolygonTask>& iterations);
void triangulatePseudoPolygonIteration(
const std::vector<VertInd>& poly,
const std::vector<TriInd>& outerTris,
std::vector<TriangulatePseudopolygonTask>& iterations);
unordered_map<Edge, TriInd>& outerTris,
std::vector<TriangulatePseudoPolygonTask>& iterations);
IndexSizeType findDelaunayPoint(
const std::vector<VertInd>& poly,
IndexSizeType iA,
Expand Down Expand Up @@ -859,7 +859,7 @@ void Triangulation<T, TNearPointLocator>::insertEdges(
if(isFinalized())
throw FinalizedError(CDT_SOURCE_LOCATION);

std::vector<TriangulatePseudopolygonTask> tppIterations;
std::vector<TriangulatePseudoPolygonTask> tppIterations;
EdgeVec remaining;
for(; first != last; ++first)
{
Expand Down
92 changes: 39 additions & 53 deletions CDT/include/Triangulation.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -521,7 +521,7 @@ void Triangulation<T, TNearPointLocator>::insertEdgeIteration(
const Edge edge,
const Edge originalEdge,
EdgeVec& remaining,
std::vector<TriangulatePseudopolygonTask>& tppIterations)
std::vector<TriangulatePseudoPolygonTask>& tppIterations)
{
const VertInd iA = edge.v1();
VertInd iB = edge.v2();
Expand Down Expand Up @@ -556,15 +556,15 @@ void Triangulation<T, TNearPointLocator>::insertEdgeIteration(
Triangle t = triangles[iT];
std::vector<TriInd> intersected(1, iT);
std::vector<VertInd> polyL, polyR;
std::vector<TriInd> outerTrisL, outerTrisR;
polyL.reserve(2);
polyL.push_back(iA);
polyL.push_back(iVL);
outerTrisL.push_back(edgeNeighbor(t, iA, iVL));
polyR.reserve(2);
polyR.push_back(iA);
polyR.push_back(iVR);
outerTrisR.push_back(edgeNeighbor(t, iA, iVR));
unordered_map<Edge, TriInd> outerTris;
outerTris[Edge(iA, iVL)] = edgeNeighbor(t, iA, iVL);
outerTris[Edge(iA, iVR)] = edgeNeighbor(t, iA, iVR);
VertInd iV = iA;

while(!t.containsVertex(iB))
Expand Down Expand Up @@ -616,14 +616,20 @@ void Triangulation<T, TNearPointLocator>::insertEdgeIteration(
locatePointLine(vertices[iVopo], a, b, distanceTolerance);
if(loc == PtLineLocation::Left)
{
outerTrisL.push_back(edgeNeighbor(tOpo, polyL.back(), iVopo));
const Edge e(polyL.back(), iVopo);
const TriInd outer = edgeNeighbor(tOpo, e.v1(), e.v2());
if(!outerTris.insert(std::make_pair(e, outer)).second)
outerTris[e] = noNeighbor; // hanging edge detected
polyL.push_back(iVopo);
iV = iVL;
iVL = iVopo;
}
else if(loc == PtLineLocation::Right)
{
outerTrisR.push_back(edgeNeighbor(tOpo, polyR.back(), iVopo));
const Edge e(polyR.back(), iVopo);
const TriInd outer = edgeNeighbor(tOpo, e.v1(), e.v2());
if(!outerTris.insert(std::make_pair(e, outer)).second)
outerTris[e] = noNeighbor; // hanging edge detected
polyR.push_back(iVopo);
iV = iVR;
iVR = iVopo;
Expand All @@ -635,8 +641,8 @@ void Triangulation<T, TNearPointLocator>::insertEdgeIteration(
iT = iTopo;
t = triangles[iT];
}
outerTrisL.push_back(edgeNeighbor(t, polyL.back(), iB));
outerTrisR.push_back(edgeNeighbor(t, polyR.back(), iB));
outerTris[Edge(polyL.back(), iB)] = edgeNeighbor(t, polyL.back(), iB);
outerTris[Edge(polyR.back(), iB)] = edgeNeighbor(t, polyR.back(), iB);
polyL.push_back(iB);
polyR.push_back(iB);

Expand All @@ -647,26 +653,16 @@ void Triangulation<T, TNearPointLocator>::insertEdgeIteration(
pivotVertexTriangleCW(iA);
if(m_vertTris[iB] == intersected.back())
pivotVertexTriangleCW(iB);
// Handle outer triangles for the cases of hanging edges
typedef std::vector<TriInd>::iterator TriIndIt;
std::sort(intersected.begin(), intersected.end());
for(TriIndIt it = outerTrisR.begin(); it != outerTrisR.end(); ++it)
if(std::binary_search(intersected.begin(), intersected.end(), *it))
*it = noNeighbor;
for(TriIndIt it = outerTrisL.begin(); it != outerTrisL.end(); ++it)
if(std::binary_search(intersected.begin(), intersected.end(), *it))
*it = noNeighbor;
// Remove intersected triangles
typedef std::vector<TriInd>::const_iterator TriIndCit;
for(TriIndCit it = intersected.begin(); it != intersected.end(); ++it)
makeDummy(*it);
{ // Triangulate pseudo-polygons on both sides
std::reverse(polyR.begin(), polyR.end());
std::reverse(outerTrisR.begin(), outerTrisR.end());
const TriInd iTL = addTriangle();
const TriInd iTR = addTriangle();
triangulatePseudopolygon(polyL, outerTrisL, iTL, iTR, tppIterations);
triangulatePseudopolygon(polyR, outerTrisR, iTR, iTL, tppIterations);
triangulatePseudoPolygon(polyL, outerTris, iTL, iTR, tppIterations);
triangulatePseudoPolygon(polyR, outerTris, iTR, iTL, tppIterations);
}

if(iB != edge.v2()) // encountered point on the edge
Expand All @@ -688,7 +684,7 @@ void Triangulation<T, TNearPointLocator>::insertEdge(
Edge edge,
const Edge originalEdge,
EdgeVec& remaining,
std::vector<TriangulatePseudopolygonTask>& tppIterations)
std::vector<TriangulatePseudoPolygonTask>& tppIterations)
{
// use iteration over recursion to avoid stack overflows
remaining.clear();
Expand Down Expand Up @@ -776,8 +772,7 @@ void Triangulation<T, TNearPointLocator>::conformToEdgeIteration(
// don't count super-triangle vertices
e1 = Edge(e1.v1() - m_nTargetVerts, e1.v2() - m_nTargetVerts);
e2 = Edge(e2.v1() - m_nTargetVerts, e2.v2() - m_nTargetVerts);
throw IntersectingConstraintsError(
e1, e2, CDT_SOURCE_LOCATION);
throw IntersectingConstraintsError(e1, e2, CDT_SOURCE_LOCATION);
}
break;
case IntersectingConstraintEdges::TryResolve: {
Expand Down Expand Up @@ -1597,23 +1592,15 @@ void Triangulation<T, TNearPointLocator>::changeNeighbor(
}

template <typename T, typename TNearPointLocator>
void Triangulation<T, TNearPointLocator>::triangulatePseudopolygon(
void Triangulation<T, TNearPointLocator>::triangulatePseudoPolygon(
const std::vector<VertInd>& poly,
const std::vector<TriInd>& outerTris,
const TriInd iT,
const TriInd iN,
std::vector<TriangulatePseudopolygonTask>& iterations)
unordered_map<Edge, TriInd>& outerTris,
TriInd iT,
TriInd iN,
std::vector<TriangulatePseudoPolygonTask>& iterations)
{
assert(poly.size() > 2);

// note: needed for proper linking with outer triangles
// during pseudo-polygon triangulation, vertex triangle
// will be set back, see asserts at the end
for(std::size_t i = 1; i < outerTris.size(); ++i)
if(outerTris[i] == noNeighbor)
m_vertTris[poly[i]] = noNeighbor;

// note: uses interation instead of recursion to avoid stack overflows
// note: uses iteration instead of recursion to avoid stack overflows
iterations.clear();
iterations.push_back(make_tuple(
IndexSizeType(0),
Expand All @@ -1623,19 +1610,15 @@ void Triangulation<T, TNearPointLocator>::triangulatePseudopolygon(
Index(0)));
while(!iterations.empty())
{
triangulatePseudopolygonIteration(poly, outerTris, iterations);
triangulatePseudoPolygonIteration(poly, outerTris, iterations);
}

// make sure adjacent triangles were restored
for(std::size_t i = 0; i < poly.size(); ++i)
assert(m_vertTris[poly[i]] != noNeighbor);
}

template <typename T, typename TNearPointLocator>
void Triangulation<T, TNearPointLocator>::triangulatePseudopolygonIteration(
void Triangulation<T, TNearPointLocator>::triangulatePseudoPolygonIteration(
const std::vector<VertInd>& poly,
const std::vector<TriInd>& outerTris,
std::vector<TriangulatePseudopolygonTask>& iterations)
unordered_map<Edge, TriInd>& outerTris,
std::vector<TriangulatePseudoPolygonTask>& iterations)
{
IndexSizeType iA, iB;
TriInd iT, iParent;
Expand All @@ -1651,11 +1634,10 @@ void Triangulation<T, TNearPointLocator>::triangulatePseudopolygonIteration(
const VertInd a = poly[iA];
const VertInd b = poly[iB];
const VertInd c = poly[iC];

// split pseudo-polygon in two parts and triangulate them
//
// note: first part needs to be pushed on stack last
// in order to be processed first
//
// note: second part needs to be pushed on stack first to be processed first

// second part: points after the Delaunay point
if(iB - iC > 1)
{
Expand All @@ -1664,13 +1646,16 @@ void Triangulation<T, TNearPointLocator>::triangulatePseudopolygonIteration(
}
else // pseudo-poly is reduced to a single outer edge
{
const TriInd outerTri = outerTris[iC];
const Edge outerEdge(b, c);
const TriInd outerTri = outerTris.at(outerEdge);
if(outerTri != noNeighbor)
{
assert(outerTri != iT);
t.neighbors[1] = outerTri;
changeNeighbor(outerTri, c, b, iT);
}
else
outerTris.at(outerEdge) = iT;
}
// first part: points before the Delaunay point
if(iC - iA > 1)
Expand All @@ -1680,22 +1665,23 @@ void Triangulation<T, TNearPointLocator>::triangulatePseudopolygonIteration(
}
else
{ // pseudo-poly is reduced to a single outer edge
const TriInd outerTri =
outerTris[iA] != noNeighbor ? outerTris[iA] : m_vertTris[c];
const Edge outerEdge(c, a);
const TriInd outerTri = outerTris.at(outerEdge);
if(outerTri != noNeighbor)
{
assert(outerTri != iT);
t.neighbors[2] = outerTri;
changeNeighbor(outerTri, c, a, iT);
}
else
outerTris.at(outerEdge) = iT;
}
// Finalize triangle
// note: only when triangle is finalized to we add it as a neighbor to
// parent to maintain triangulation topology consistency
triangles[iParent].neighbors[iInParent] = iT;
t.neighbors[0] = iParent;
t.vertices = detail::arr3(a, b, c);
setAdjacentTriangle(a, iT);
setAdjacentTriangle(c, iT);
}

Expand Down
9 changes: 9 additions & 0 deletions CDT/tests/cdt.test.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -936,3 +936,12 @@ TEST_CASE("Regression test issue #154 (3)", "")
REQUIRE(topologyString(cdt) == topologyString(outFile));
}
}

TEST_CASE("Regression test: hanging edge in pseudo-poly", "")
{
auto [vv, ee] = readInputFromFile<double>("inputs/hanging3.txt");
auto cdt = Triangulation<double>{};
cdt.insertVertices(vv);
cdt.insertEdges(ee);
REQUIRE(CDT::verifyTopology(cdt));
}
9 changes: 9 additions & 0 deletions CDT/tests/inputs/hanging3.txt
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,9 @@
7 1
1534.79 789.063
-785.078 788.629
789.094 533.067
1034.16 789.067
785.067 513.067
784 664.004
513.064 789.067
0 1

0 comments on commit ffb69c6

Please sign in to comment.