From 3034e55317b496bda5da9393cdffb52d3a5b1744 Mon Sep 17 00:00:00 2001 From: Jason Cantarella Date: Sun, 21 Jun 2026 13:24:33 -0600 Subject: [PATCH] flip_geodesics: terminate iterativeShorten() on contractible closed loops A contractible closed loop has no geodesic representative, so iterativeShorten() would collapse it to a single self-edge and spin forever. processSingleEdgeLoop() re-routes the self-edge onto the two opposite edges of its triangle -- always longer by the triangle inequality -- with no progress guard, forming a period-2 limit cycle with the generic shortening step (unfold -> reshorten -> same self-edge -> ...). Add a per-loop guard: FlipEdgePath::lastSingleEdgeLoopLen records the shortest length a loop reaches while collapsed to a single self-edge, and processSingleEdgeLoop() only unfolds when strictly shorter than last time. A contractible loop returns to the same self-edge each round (no progress) and stops at its minimal single-edge form; a non-contractible loop strictly shortens toward its positive-length geodesic each round and is unaffected. Note: a blanket "never unfold if it lengthens" guard is NOT sufficient -- non-contractible loops can also pass through a single self-edge en route to their geodesic and need the (always-lengthening) unfold to straighten the final corner; blanket guarding strands them at a kink (verified: minAngle 2.45 vs geodesic 3.50 on cone-concentrated tori). Adds three deterministic regression tests (contractible terminates; non-contractible reaches geodesic; non-contractible-through-self-edge reaches geodesic). Full suite passes (226/226). Prepared with Claude Code (Opus 4.8) --- .../geometrycentral/surface/flip_geodesics.h | 6 + src/surface/flip_geodesics.cpp | 20 +++ test/src/flip_geodesic_test.cpp | 165 ++++++++++++++++++ 3 files changed, 191 insertions(+) diff --git a/include/geometrycentral/surface/flip_geodesics.h b/include/geometrycentral/surface/flip_geodesics.h index 2f820d7c..5ee4b2dc 100644 --- a/include/geometrycentral/surface/flip_geodesics.h +++ b/include/geometrycentral/surface/flip_geodesics.h @@ -4,6 +4,7 @@ #include "geometrycentral/surface/signpost_intrinsic_triangulation.h" #include "geometrycentral/surface/vertex_position_geometry.h" +#include #include #include #include @@ -47,6 +48,11 @@ class FlipEdgePath { FlipEdgeNetwork& network; // the parent network of which this path is a part bool isClosed; + // Shortest length this loop has reached while collapsed to a single self-edge (see + // FlipEdgeNetwork::processSingleEdgeLoop). Used to guarantee termination on contractible + // closed loops, which have no geodesic representative. Infinity until first reached. + double lastSingleEdgeLoopLen = std::numeric_limits::infinity(); + // Data about the segments in path. For a path segment id, contains the corresponding halfedge, as well as the ID of // the previous/next segment (or INVALID_IND if they are path endpoints) std::unordered_map> diff --git a/src/surface/flip_geodesics.cpp b/src/surface/flip_geodesics.cpp index 43c11655..e414de35 100644 --- a/src/surface/flip_geodesics.cpp +++ b/src/surface/flip_geodesics.cpp @@ -880,6 +880,26 @@ void FlipEdgeNetwork::processSingleEdgeLoop(FlipPathSegment& pathSegment, Segmen Halfedge he; std::tie(he, UNUSED1, UNUSED2) = edgePath.pathHeInfo[id]; + // Termination guard for contractible closed loops. + // + // A single self-edge loop is the maximally-collapsed form of a closed loop. The only move + // available here -- replacing the self-edge with the two opposite edges of its triangle -- + // ALWAYS lengthens the loop (triangle inequality), so it is only worthwhile if the subsequent + // straightening recovers that length and then some (net progress toward a geodesic). + // + // For a CONTRACTIBLE loop that never happens: the loop has no geodesic representative (its + // length infimum is zero, reached only as it shrinks to a point), so the unfold + reshorten + // returns to the very same self-edge and the process loops forever. We detect this by + // remembering the shortest length this loop has reached while collapsed to a single edge, and + // only unfolding when we are strictly shorter than last time. A non-contractible loop strictly + // shortens toward its (positive-length) geodesic on each such round and always passes this + // test; a contractible loop stalls and we stop, leaving it as its minimal single-edge form. + double curLoopLen = tri->edgeLengths[he.edge()]; + if (curLoopLen >= edgePath.lastSingleEdgeLoopLen * (1.0 - 1e-12)) { + return; // no net progress since the last unfold of this loop; stop to guarantee termination + } + edgePath.lastSingleEdgeLoopLen = curLoopLen; + // == Replace the old segment with the two opposite edges of the triangle diff --git a/test/src/flip_geodesic_test.cpp b/test/src/flip_geodesic_test.cpp index 93f19f70..da174da1 100644 --- a/test/src/flip_geodesic_test.cpp +++ b/test/src/flip_geodesic_test.cpp @@ -168,3 +168,168 @@ TEST_F(FlipGeodesicSuite, PiecewiseDijkstraMultipleOverlapAtEnd) { ASSERT_EQ(network->paths.size(), 1); EXPECT_EQ(network->paths[0]->getHalfedgeList().size(), 1); } + +// ============================================================ +// =============== Closed-loop shortening (termination) +// ============================================================ +// +// A contractible closed loop has no nontrivial geodesic representative, so flip-shortening can +// collapse it to a single self-edge and then spin forever unfolding/reshortening it. The fix in +// processSingleEdgeLoop must (a) make such loops TERMINATE, while (b) NOT regressing +// non-contractible loops, which can legitimately pass through a single self-edge on the way to a +// real geodesic (and need the unfold to straighten that last corner). + +namespace { + +// Deterministic [-1,1] hash (no needed), shared by the mesh builders below. +double flipLoopPseudo(int i) { + double s = std::sin((double)i * 12.9898) * 43758.5453; + return 2.0 * (s - std::floor(s)) - 1.0; +} + +// Flat triangulated disk with a sharp cone apex at the center. Intrinsically Euclidean away from +// the apex, so EVERY closed loop is contractible. Returns the mesh plus the ring-1 loop edges. +std::tuple, std::unique_ptr> +buildConeDisk(int R, int nSeg, double coneH) { + auto ringVert = [&](int r, int k) -> size_t { + return 1 + (size_t)(r - 1) * nSeg + (size_t)(((k % nSeg) + nSeg) % nSeg); + }; + std::vector pos; + pos.push_back(Vector3{0.0, 0.0, coneH}); // apex + for (int r = 1; r <= R; r++) + for (int k = 0; k < nSeg; k++) { + double ang = 2.0 * M_PI * k / nSeg; + pos.push_back(Vector3{(double)r * std::cos(ang), (double)r * std::sin(ang), 0.0}); + } + std::vector> faces; + for (int k = 0; k < nSeg; k++) faces.push_back({0, ringVert(1, k), ringVert(1, k + 1)}); + for (int r = 1; r < R; r++) + for (int k = 0; k < nSeg; k++) { + size_t a = ringVert(r, k), b = ringVert(r, k + 1), c = ringVert(r + 1, k), d = ringVert(r + 1, k + 1); + faces.push_back({a, c, d}); + faces.push_back({a, d, b}); + } + return makeManifoldSurfaceMeshAndGeometry(faces, pos); +} + +// Torus, optionally with one vertex spiked radially (a cone point) and the tube pinched. On +// cone-concentrated/pinched tori, a non-contractible ring collapses to a sub-pi self-edge whose +// unfold is required to reach the geodesic -- the case that distinguishes a correct fix from a +// blanket no-op. +std::tuple, std::unique_ptr> +buildTorus(int nU, int nV, double Rmaj, double rmin, double jitter, int coneVert, double coneAmt, double pinch) { + auto vid = [&](int i, int j) { return (size_t)((i % nU) * nV + (j % nV)); }; + std::vector pos(nU * nV); + for (int i = 0; i < nU; i++) + for (int j = 0; j < nV; j++) { + double u = 2 * M_PI * i / nU, v = 2 * M_PI * j / nV; + double rm = rmin * (1.0 + pinch * std::cos(u)); + double jit = 1.0 + jitter * flipLoopPseudo(i * 131 + j * 17); + double Rr = (Rmaj + rm * std::cos(v)) * jit; + Vector3 P{Rr * std::cos(u), Rr * std::sin(u), rm * std::sin(v) * jit}; + size_t id = vid(i, j); + if ((int)id == coneVert) P *= (1.0 + coneAmt); + pos[id] = P; + } + std::vector> faces; + for (int i = 0; i < nU; i++) + for (int j = 0; j < nV; j++) + faces.push_back({vid(i, j), vid(i + 1, j), vid(i + 1, j + 1)}), + faces.push_back({vid(i, j), vid(i + 1, j + 1), vid(i, j + 1)}); + return makeManifoldSurfaceMeshAndGeometry(faces, pos); +} + +// Mark the edge va--vb in the path set; returns false if not found. +bool markLoopEdge(ManifoldSurfaceMesh& mesh, EdgeData& inPath, size_t va, size_t vb) { + for (Halfedge he : mesh.vertex(va).outgoingHalfedges()) + if (he.tipVertex() == mesh.vertex(vb)) { + inPath[he.edge()] = true; + return true; + } + return false; +} + +} // namespace + +// A contractible closed loop must TERMINATE (not spin). We cap iterations generously so that a +// regression cannot hang the test suite; a correct run drains the queue in a handful of steps. +TEST_F(FlipGeodesicSuite, ContractibleClosedLoopTerminates) { + const int nSeg = 6; + std::unique_ptr meshPtr; + std::unique_ptr geomPtr; + std::tie(meshPtr, geomPtr) = buildConeDisk(/*R=*/2, nSeg, /*coneH=*/4.0); + ManifoldSurfaceMesh& mesh = *meshPtr; + VertexPositionGeometry& geom = *geomPtr; + + EdgeData inPath(mesh, false); + for (int k = 0; k < nSeg; k++) ASSERT_TRUE(markLoopEdge(mesh, inPath, 1 + k, 1 + (k + 1) % nSeg)); + VertexData extraMarked(mesh, false); + + std::unique_ptr network = + FlipEdgeNetwork::constructFromEdgeSet(mesh, geom, inPath, extraMarked); + ASSERT_NE(network, nullptr); + ASSERT_EQ(network->paths.size(), 1); + ASSERT_TRUE(network->paths[0]->isClosed); + network->addAllWedgesToAngleQueue(); // queue the closed-loop seam wedge (ctor skips the first segment) + + const size_t cap = 100000; + network->iterativeShorten(cap); + + // Termination signature: the queue drained and we did NOT exhaust the cap (a spin would do both + // the opposite). This assertion cannot hang even if the guard is removed -- the cap bounds it. + EXPECT_TRUE(network->wedgeAngleQueue.empty()); + EXPECT_LT(network->nShortenIters, (long long)cap); +} + +// A non-contractible loop on a plain torus must converge to a closed geodesic. +TEST_F(FlipGeodesicSuite, NonContractibleClosedLoopReachesGeodesic) { + const int nU = 12, nV = 8; + std::unique_ptr meshPtr; + std::unique_ptr geomPtr; + std::tie(meshPtr, geomPtr) = buildTorus(nU, nV, 3.0, 1.0, 0.0, -1, 0.0, 0.0); + ManifoldSurfaceMesh& mesh = *meshPtr; + VertexPositionGeometry& geom = *geomPtr; + + EdgeData inPath(mesh, false); + auto vid = [&](int i, int j) { return (size_t)((i % nU) * nV + (j % nV)); }; + for (int j = 0; j < nV; j++) ASSERT_TRUE(markLoopEdge(mesh, inPath, vid(0, j), vid(0, j + 1))); + VertexData extraMarked(mesh, false); + + std::unique_ptr network = + FlipEdgeNetwork::constructFromEdgeSet(mesh, geom, inPath, extraMarked); + ASSERT_NE(network, nullptr); + ASSERT_TRUE(network->paths[0]->isClosed); + network->addAllWedgesToAngleQueue(); // queue the closed-loop seam wedge (ctor skips the first segment) + + network->iterativeShorten(100000); + EXPECT_TRUE(network->wedgeAngleQueue.empty()); + EXPECT_GE(network->minAngle(), M_PI - network->EPS_ANGLE); // straight => geodesic +} + +// The discriminating case: a non-contractible ring on a cone-concentrated, pinched torus whose +// final corner collapses to a sub-pi single self-edge. Reaching the geodesic REQUIRES unfolding +// that self-edge -- so a blanket "never unfold" guard leaves it stuck at a kink (minAngle < pi), +// while the round-trip-progress guard straightens it. This test fails on the naive fix. +TEST_F(FlipGeodesicSuite, NonContractibleLoopThroughSelfEdgeReachesGeodesic) { + const int nU = 4, nV = 4; + std::unique_ptr meshPtr; + std::unique_ptr geomPtr; + std::tie(meshPtr, geomPtr) = buildTorus(nU, nV, 2.0, 1.2, 0.8, /*coneVert=*/6, /*coneAmt=*/3.5, /*pinch=*/0.7); + ManifoldSurfaceMesh& mesh = *meshPtr; + VertexPositionGeometry& geom = *geomPtr; + + EdgeData inPath(mesh, false); + auto vid = [&](int i, int j) { return (size_t)((i % nU) * nV + (j % nV)); }; + for (int j = 0; j < nV; j++) ASSERT_TRUE(markLoopEdge(mesh, inPath, vid(0, j), vid(0, j + 1))); + VertexData extraMarked(mesh, false); + + std::unique_ptr network = + FlipEdgeNetwork::constructFromEdgeSet(mesh, geom, inPath, extraMarked); + ASSERT_NE(network, nullptr); + ASSERT_TRUE(network->paths[0]->isClosed); + network->addAllWedgesToAngleQueue(); // queue the closed-loop seam wedge (ctor skips the first segment) + + network->iterativeShorten(100000); + EXPECT_TRUE(network->wedgeAngleQueue.empty()); + EXPECT_GE(network->minAngle(), M_PI - network->EPS_ANGLE); +}