FlipEdgeNetwork::iterativeShorten() hangs on a contractible closed loop

[designbynumbers]

First, thanks for geometry-central and the flip-geodesics implementation — it's
been a pleasure to build on.

This is probably a rare situation in practice, but we ran into an input on which
FlipEdgeNetwork::iterativeShorten() never returns (100% CPU, unbounded). It
seems worth reporting regardless of how you'd want to address it. Whatever the
preferred fix, the short version is just: iterativeShorten() does not
terminate on a contractible closed loop.

Reproducer

Self-contained — the mesh is generated in code and the build fetches
geometry-central, so there's no external asset. Sources are attached
(reproducer.cpp, CMakeLists.txt); reproduces on master at 019669d.

cmake -B build -DCMAKE_BUILD_TYPE=Release
cmake --build build
./build/reproducer        # prints "SPIN CONFIRMED" and never converges

It builds a triangulated disk (with one sharp cone apex raised at the center;
13 vertices), marks the innermost ring as a closed loop via
constructFromEdgeSet, and runs iterativeShorten() in fixed-size chunks. The
surface is a topological disk, so it is simply connected and every closed loop on
it is contractible — including this one, which encircles the cone apex. Output:

chunk  cumIters    length        dLength      queueSize
   1       2000    8.24621125   +2.246e+00           2
   2       4000    8.24621125   +0.000e+00           2
   ...
   8      16000    8.24621125   +0.000e+00           2

The signature: the length plateaus, the wedge queue never drains, and the
iteration counter climbs without bound. (nFlips is frozen during the
spin — the repeated work is degenerate path relabeling, not edge flips.)

What seems to be happening

A contractible closed curve has no nontrivial geodesic representative, so the
loop collapses all the way down to a single self-edge. From there,
locallyShortenAt() dispatches to processSingleEdgeLoop(), which re-routes the
self-edge onto the two opposite edges of its triangle. By the triangle
inequality that replacement is always longer, and unlike the generic
straightening path (which has if (newPathLength > initPathLength) return;),
processSingleEdgeLoop() has no such guard. The two then form a period-2 cycle:

single self-edge loop (len L)
   --processSingleEdgeLoop--> two-edge loop (len L' > L)     [length up, no guard]
   --generic shorten-------->  single self-edge loop (len L) [length back down]
   --> repeat forever

(The chunked length readout looks constant only because 2000 is even and samples
the period-2 oscillation at the same phase; an instrumented trace shows it
oscillating, e.g. 5.495 ↔ 8.246.)

We don't think open paths or non-contractible loops are affected: open paths
never reach processSingleEdgeLoop(), and non-contractible loops reach a
straight geodesic and terminate normally.

A conservative hotfix (linked PR)

We have provided a hotfix as a linked PR, #254. It's intentionally minimal: it gives
each loop one extra field recording the shortest length it has reached while
collapsed to a single self-edge, and only lets processSingleEdgeLoop() unfold
when the loop is strictly shorter than the last time it unfolded.

What makes it safe to drop in is that it is exactly as powerful as the current
code on any input that already terminates. It only declines to unfold when a
loop revisits a single-self-edge state without having strictly shortened since
the previous one — and under the current deterministic shortening, that is
already a non-terminating cycle (if it cycles once with no reduction, it never
reduces). So it can only convert a hang into a return; it never shortens,
lengthens, or otherwise changes a run that currently completes.

We'd stress that this is a hang→return hotfix, not a finished design: it leaves a
contractible loop as its minimal single-edge representative rather than reporting
it as trivial, so you may well prefer something more thorough. The PR has the
full details (including why the obvious one-line fix turned out not to work).
Either way, we mostly wanted to make sure the hang itself was on your radar.

---

reproducer.cpp

// ============================================================================
// Minimal standalone reproducer:
//   FlipEdgeNetwork::iterativeShorten() does NOT terminate when the network
//   contains a CONTRACTIBLE closed loop.
//
// A contractible closed curve on a surface has no nontrivial geodesic
// representative -- its length infimum is 0, achieved only in the limit as the
// loop shrinks to a point. iterativeShorten() keeps flipping to shorten the
// loop, eventually collapses it onto a single self-edge, and then
// processSingleEdgeLoop() oscillates forever: it replaces the single edge with
// the two opposite edges of the incident triangle and re-queues both wedges,
// with NO monotone-length-decrease guard (unlike the generic locallyShortenAt
// path, which bails via `if (newPathLength > initPathLength) return;`). The
// unfold is always longer than the single edge it replaces (triangle
// inequality) and the generic step legitimately shortens it back to the same
// self-edge, so the two alternate forever -- a period-2 cycle.
//
// This program builds a triangulated disk with a sharp cone apex at the center.
// Being a topological disk it is simply connected, so EVERY closed loop on it is
// contractible (this one encircles the apex). It marks the innermost interior
// ring as a closed loop, and runs iterativeShorten in fixed-size chunks. The
// signature of the bug: the path length plateaus while the shorten-iteration
// counter keeps climbing without bound -- iterations are consumed doing no
// work. A correct implementation would detect the trivialized loop and stop.
//
// Self-contained: the mesh is generated in code, no external mesh file needed.
// ============================================================================

#include "geometrycentral/surface/flip_geodesics.h"
#include "geometrycentral/surface/manifold_surface_mesh.h"
#include "geometrycentral/surface/surface_mesh_factories.h"
#include "geometrycentral/surface/vertex_position_geometry.h"

#include <cmath>
#include <cstdlib>
#include <iostream>
#include <vector>

using namespace geometrycentral;
using namespace geometrycentral::surface;

int main(int argc, char** argv) {
  // ---- Build a triangulated disk: 1 center vertex + R concentric rings of
  //      nSeg vertices each. It is a topological disk, hence simply connected,
  //      so any closed loop on it is contractible.
  const int R = (argc > 4) ? std::atoi(argv[4]) : 2;     // number of rings
  const int nSeg = (argc > 5) ? std::atoi(argv[5]) : 6;  // vertices per ring
  // Cone apex height at the center: >0 makes v0 a positive-curvature cone point
  // (angle defect), which is what keeps the collapsed single-edge loop's wedge
  // perpetually turned (the real meshes concentrate curvature; a perfectly flat
  // symmetric disk instead collapses to a straight wedge and converges). Set
  // via argv[2]; an irregularity seed (argv[3]) perturbs vertices to mimic the
  // non-symmetric triangle fans of a real intrinsic triangulation.
  const double coneH = (argc > 2) ? std::atof(argv[2]) : 4.0;
  const double jitter = (argc > 3) ? std::atof(argv[3]) : 0.0;
  auto pseudo = [](int i) {  // deterministic [-1,1] hash, no <random> needed
    double s = std::sin((double)i * 12.9898) * 43758.5453;
    return 2.0 * (s - std::floor(s)) - 1.0;
  };

  std::vector<Vector3> pos;
  pos.push_back(Vector3{0.0, 0.0, coneH});  // v0 = cone apex
  auto ringVert = [&](int r, int k) -> size_t {
    // r in [1..R], k in [0..nSeg)
    return 1 + (size_t)(r - 1) * nSeg + (size_t)(((k % nSeg) + nSeg) % nSeg);
  };
  for (int r = 1; r <= R; r++) {
    for (int k = 0; k < nSeg; k++) {
      double ang = 2.0 * M_PI * k / nSeg;
      double rad = (double)r;
      int id = r * 131 + k;
      double jr = 1.0 + jitter * pseudo(id);
      double ja = jitter * 0.5 * pseudo(id + 9991) * (2.0 * M_PI / nSeg);
      pos.push_back(Vector3{rad * jr * std::cos(ang + ja),
                            rad * jr * std::sin(ang + ja), 0.0});
    }
  }

  std::vector<std::vector<size_t>> faces;
  // Center fan: center -> ring 1
  for (int k = 0; k < nSeg; k++) {
    faces.push_back({0, ringVert(1, k), ringVert(1, k + 1)});
  }
  // Quad strips between ring r and ring r+1, split into 2 triangles each.
  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);
      size_t c = ringVert(r + 1, k), d = ringVert(r + 1, k + 1);
      faces.push_back({a, c, d});
      faces.push_back({a, d, b});
    }
  }

  std::unique_ptr<ManifoldSurfaceMesh> mesh;
  std::unique_ptr<VertexPositionGeometry> geom;
  std::tie(mesh, geom) = makeManifoldSurfaceMeshAndGeometry(faces, pos);
  std::cout << "mesh: " << mesh->nVertices() << " verts, " << mesh->nFaces()
            << " faces, " << mesh->nEdges() << " edges, Euler chi="
            << ((long)mesh->nVertices() - (long)mesh->nEdges() +
                (long)mesh->nFaces())
            << " (disk => 1)\n";

  // ---- Mark the innermost interior ring (ring 1) as the cut loop. These edges
  //      form a closed cycle around the center vertex; since the surface is a
  //      disk this loop is contractible. They are interior (each shared by an inner fan
  //      triangle and an outer strip triangle), hence flippable.
  geom->requireEdgeIndices();
  EdgeData<bool> inPath(*mesh, false);
  int nLoopEdges = 0;
  for (int k = 0; k < nSeg; k++) {
    Vertex va = mesh->vertex(ringVert(1, k));
    Vertex vb = mesh->vertex(ringVert(1, k + 1));
    // find the edge va--vb
    for (Halfedge he : va.outgoingHalfedges()) {
      if (he.tipVertex() == vb) {
        inPath[he.edge()] = true;
        nLoopEdges++;
        break;
      }
    }
  }
  VertexData<bool> extraMarked(*mesh, false);  // no marked endpoints => closed
  std::cout << "marked " << nLoopEdges << " ring-1 edges as a closed loop\n";

  // ---- Build the flip-edge network from the marked edge set (the same entry
  //      point pendel uses: constructFromEdgeSet, which infers loopiness).
  std::unique_ptr<FlipEdgeNetwork> net =
      FlipEdgeNetwork::constructFromEdgeSet(*mesh, *geom, inPath, extraMarked);
  if (!net || net->paths.empty()) {
    std::cerr << "FAILED to construct network\n";
    return 2;
  }
  net->addAllWedgesToAngleQueue();
  std::cout << "network: " << net->paths.size() << " path(s), path[0].isClosed="
            << (net->paths[0]->isClosed ? "true" : "false")
            << ", initial length=" << net->length() << "\n\n";

  // ---- Run iterativeShorten in chunks. If the loop were converging, the
  //      length would strictly decrease toward a geodesic and the wedge queue
  //      would drain (length stabilizes, iterations stop being consumed).
  //      Instead we will see: length plateaus, iterations keep climbing.
  const size_t chunk = 2000;
  const int nChunks = (argc > 1) ? std::atoi(argv[1]) : 8;
  double prevLen = net->length();
  long long prevIters = net->nShortenIters;
  int flatStalledChunks = 0;  // length unchanged AND iters fully consumed
  std::cout << "chunk  cumIters    length        dLength      queueSize\n";
  std::cout << "-----  ---------   -----------   ----------   ---------\n";
  for (int c = 0; c < nChunks; c++) {
    net->iterativeShorten(chunk);
    double len = net->length();
    long long iters = net->nShortenIters;
    std::printf("%4d   %8lld   %11.8f   %+10.3e   %9zu\n", c + 1, iters, len,
                len - prevLen, net->wedgeAngleQueue.size());
    bool itersFullyConsumed = (iters - prevIters) >= (long long)chunk;
    if (len == prevLen && itersFullyConsumed && !net->wedgeAngleQueue.empty())
      flatStalledChunks++;
    else
      flatStalledChunks = 0;
    prevLen = len;
    prevIters = iters;
  }

  std::cout << "final path[0].size()=" << net->paths[0]->size()
            << " isClosed=" << (net->paths[0]->isClosed ? "true" : "false")
            << " nFlips=" << net->nFlips << "\n\n";

  // Verdict: a spin = the last >=3 chunks each consumed the FULL iteration
  // budget while the length did not change and the wedge queue never drained.
  // A correct run drains the queue (queueSize -> 0, iters stop being consumed).
  bool spinning = flatStalledChunks >= 3;
  if (spinning) {
    std::cout << "VERDICT: SPIN CONFIRMED -- iterativeShorten consumed every "
                 "iteration with zero length change and a non-empty queue. This "
                 "is the non-terminating contractible-loop hang.\n";
    return 1;
  }
  std::cout << "VERDICT: converged (queue drained). Increase the cone height "
               "(argv[2]) to sharpen the curvature concentration.\n";
  return 0;
}

CMakeLists.txt

# Standalone build for the contractible-loop hang reproducer.
#
# Self-contained: fetches geometry-central (UPSTREAM by default -- the bug is in
# upstream flip_geodesics.cpp, not anything fork-specific). Override the repo/tag
# with -DGC_REPOSITORY=... -DGC_TAG=... to test a fix branch.
#
#   cmake -B build -DCMAKE_BUILD_TYPE=Release
#   cmake --build build
#   ./build/reproducer        # exits 1 and prints "SPIN CONFIRMED"
#
cmake_minimum_required(VERSION 3.16)
project(contractible_loop_hang CXX)

set(CMAKE_CXX_STANDARD 17)
set(CMAKE_CXX_STANDARD_REQUIRED ON)

set(GC_REPOSITORY "https://github.com/nmwsharp/geometry-central.git"
    CACHE STRING "geometry-central git repo")
set(GC_TAG "master" CACHE STRING "geometry-central git tag/branch/sha")

include(FetchContent)
FetchContent_Declare(
  geometrycentral
  GIT_REPOSITORY ${GC_REPOSITORY}
  GIT_TAG        ${GC_TAG}
)
FetchContent_MakeAvailable(geometrycentral)

add_executable(reproducer reproducer.cpp)
target_link_libraries(reproducer geometry-central)