From dbc61b78ce24232603347a6ff52a77d3cf530c0c Mon Sep 17 00:00:00 2001 From: Mattias Fagerlund Date: Fri, 30 Oct 2026 06:37:04 +0100 Subject: [PATCH 1/3] Trace intrinsic edges between vertices robustly on degenerate meshes traceIntrinsicHalfedgeAlongInput traced an edge between two input vertices with traceGeodesic and kept the untrimmed trace when it missed the end vertex. On meshes with zero-area faces, or where the ray passes exactly through a vertex, the trace ends inside a face or turns off course, and the exported polyline is not a path between the two vertices (FlipEdgeNetwork::getPathPolyline on Thingi10K meshes). New traceGeodesicBetweenVertices: the trimmed trace is kept when it reaches the end vertex with the edge's length (1e-9 relative), so clean meshes give the same points bit for bit; otherwise the edge is traced again by unfolding its faces along the line (traceGeodesicByUnfolding: layout from positions where available, both sides tried at a vertex hit, re-aimed at the end vertex's unfolded position). surfacePathLength measures a surface point path. Tests: through a vertex of a zero-area face, from a corner of angle pi, FlipEdgeNetwork's polyline through a zero-area face, and bit-identity with the old trace on a clean mesh. --- .../geometrycentral/surface/trace_geodesic.h | 29 ++ .../signpost_intrinsic_triangulation.cpp | 6 +- src/surface/trace_geodesic.cpp | 353 ++++++++++++++++++ test/CMakeLists.txt | 1 + test/src/trace_geodesic_test.cpp | 131 +++++++ 5 files changed, 519 insertions(+), 1 deletion(-) create mode 100644 test/src/trace_geodesic_test.cpp diff --git a/include/geometrycentral/surface/trace_geodesic.h b/include/geometrycentral/surface/trace_geodesic.h index 33540194..97e06e77 100644 --- a/include/geometrycentral/surface/trace_geodesic.h +++ b/include/geometrycentral/surface/trace_geodesic.h @@ -52,5 +52,34 @@ TraceGeodesicResult traceGeodesic(IntrinsicGeometryInterface& geom, Face startFa // Return value indicates success. If true, the resulting ends in the 1-ring of targetVertex as expected. bool trimTraceResult(TraceGeodesicResult& traceResult, Vertex targetVertex); +// Trace a geodesic segment known to run from vertex startVert to vertex endVert (an intrinsic edge between two +// vertices of the input mesh): traceVec in startVert's tangent space, its norm the segment's length. The result +// runs from startVert through the edge crossings to endVert. The face-by-face trace above is used and trimmed +// with trimTraceResult; when it does not reach endVert's 1-ring, or the length of the trimmed path differs from +// norm(traceVec) by more than lengthTol relative (the trace turned off course at a degenerate face or at a vertex +// it passed through), the segment is traced again by unfolding its faces along the line, which has neither +// failure. If that fails too, the face-by-face trace is returned as before: trimmed if it trimmed, else untrimmed. +std::vector traceGeodesicBetweenVertices(IntrinsicGeometryInterface& geom, Vertex startVert, + Vector2 traceVec, Vertex endVert, + const TraceOptions& traceOptions = defaultTraceOptions, + double lengthTol = 1e-9); + +// Trace by unfolding: the faces the segment crosses are laid out one by one in the segment's frame (startVert at +// the origin, the segment along +x), each crossed edge chosen by the side of the line its new vertex lies on; +// from the vertex positions if geom is a VertexPositionGeometry (the edge lengths lose the width of a sliver). +// Where the line passes within rounding of a vertex both sides are tried. The line arrives where endVert is laid +// out at a point whose chord from the start has the segment's length (within arriveTol relative) and crosses +// every edge crossed so far; the crossings are then put on that chord, so the path ends exactly at endVert. +// traceVec only aims the line: where the line misses endVert (a direction carried through degenerate faces can +// be off while the length is right), it is turned toward endVert's unfolded position and unfolded again (up to +// 8 times). Returns an empty vector if it never arrives. +std::vector traceGeodesicByUnfolding(IntrinsicGeometryInterface& geom, Vertex startVert, + Vector2 traceVec, Vertex endVert, double arriveTol = 1e-6, + size_t maxIters = INVALID_IND); + +// Length of a path of surface points (from the vertex positions if geom is a VertexPositionGeometry), each step +// measured in a face the two points share; NaN if two consecutive points share no face. +double surfacePathLength(IntrinsicGeometryInterface& geom, const std::vector& path); + } // namespace surface } // namespace geometrycentral diff --git a/src/surface/signpost_intrinsic_triangulation.cpp b/src/surface/signpost_intrinsic_triangulation.cpp index 9ab4aee5..f46cef3e 100644 --- a/src/surface/signpost_intrinsic_triangulation.cpp +++ b/src/surface/signpost_intrinsic_triangulation.cpp @@ -75,10 +75,14 @@ std::vector SignpostIntrinsicTriangulation::traceIntrinsicHalfedge TraceOptions options; options.includePath = true; options.maxIters = mesh.nFaces() * 10; + Vertex endVert = he.twin().vertex(); + if (trimEnd && startP.type == SurfacePointType::Vertex && vertexLocations[endVert].type == SurfacePointType::Vertex) { + // between two input vertices: trimmed and checked, traced again by unfolding where the trace went off course + return traceGeodesicBetweenVertices(inputGeom, startP.vertex, traceVec, vertexLocations[endVert].vertex, options); + } TraceGeodesicResult result = traceGeodesic(inputGeom, startP, traceVec, options); // Trim off end crumbs if applicable - Vertex endVert = he.twin().vertex(); if (trimEnd && vertexLocations[endVert].type == SurfacePointType::Vertex) { bool success = trimTraceResult(result, endVert); if (success) { diff --git a/src/surface/trace_geodesic.cpp b/src/surface/trace_geodesic.cpp index 7a9f55ad..ad5d6921 100644 --- a/src/surface/trace_geodesic.cpp +++ b/src/surface/trace_geodesic.cpp @@ -6,6 +6,9 @@ #include #include +#include +#include +#include using std::cout; using std::endl; @@ -876,5 +879,355 @@ bool trimTraceResult(TraceGeodesicResult& traceResult, Vertex targetVertex) { } +namespace { +// Position of a surface point in face f's layout (f.halfedge()'s tail at the origin, as +// vertexCoordinatesInTriangle); false if the point is not on f. +bool pointInFace(IntrinsicGeometryInterface& geom, Face f, const SurfacePoint& p, Vector2& pos) { + std::array c = vertexCoordinatesInTriangle(geom, f); + Halfedge h = f.halfedge(); + switch (p.type) { + case SurfacePointType::Vertex: + for (int k = 0; k < 3; k++, h = h.next()) + if (h.tailVertex() == p.vertex) { + pos = c[k]; + return true; + } + return false; + case SurfacePointType::Edge: + for (int k = 0; k < 3; k++, h = h.next()) + if (h.edge() == p.edge) { + Vector2 tail = c[k], tip = c[(k + 1) % 3]; + if (h != p.edge.halfedge()) std::swap(tail, tip); + pos = tail + p.tEdge * (tip - tail); + return true; + } + return false; + case SurfacePointType::Face: + if (p.face != f) return false; + pos = c[0] * p.faceCoords.x + c[1] * p.faceCoords.y + c[2] * p.faceCoords.z; + return true; + } + return false; +} +} // namespace + +double surfacePathLength(IntrinsicGeometryInterface& geom, const std::vector& path) { + geom.requireHalfedgeVectorsInFace(); + VertexPositionGeometry* positions = dynamic_cast(&geom); + double L = 0.; + std::vector faces; + for (size_t i = 0; i + 1 < path.size(); i++) { + const SurfacePoint& p = path[i]; + faces.clear(); + switch (p.type) { + case SurfacePointType::Vertex: + for (Face f : p.vertex.adjacentFaces()) faces.push_back(f); + break; + case SurfacePointType::Edge: + if (p.edge.halfedge().isInterior()) faces.push_back(p.edge.halfedge().face()); + if (p.edge.halfedge().twin().isInterior()) faces.push_back(p.edge.halfedge().twin().face()); + break; + case SurfacePointType::Face: + faces.push_back(p.face); + break; + } + bool shared = false; + for (Face f : faces) { + Vector2 a, b; + if (pointInFace(geom, f, p, a) && pointInFace(geom, f, path[i + 1], b)) { + // from the positions if the geometry has them (the edge lengths lose the width of a sliver) + L += positions ? norm(path[i + 1].interpolate(positions->inputVertexPositions) - + p.interpolate(positions->inputVertexPositions)) + : norm(b - a); + shared = true; + break; + } + } + if (!shared) { + geom.unrequireHalfedgeVectorsInFace(); + return std::numeric_limits::quiet_NaN(); + } + } + geom.unrequireHalfedgeVectorsInFace(); + return L; +} + +namespace { +// One unfolding of the line from startVert at tangent-space angle theta (see traceGeodesicByUnfolding). +// Returns true with the path in out on arrival; otherwise aimOffset is the angle (in the start's corner, +// radians) of endVert's unfolded position as seen from the start, where it was laid out nearest the line +// (NaN if it never was). +bool unfoldLine(IntrinsicGeometryInterface& geom, VertexPositionGeometry* positions, Vertex startVert, + double theta, double length, Vertex endVert, double arriveTol, size_t maxIters, + std::vector& out, double& aimOffset) { + out.clear(); + aimOffset = std::numeric_limits::quiet_NaN(); + auto arg = [](Vector2 v) { + double t = std::atan2(v.y, v.x); + return t < 0. ? t + 2. * M_PI : t; + }; + + // The corner of startVert holding the direction, and the direction's angle from that corner's first side + const double sumScale = geom.vertexAngleSums[startVert] / (startVert.isBoundary() ? M_PI : 2. * M_PI); + Halfedge he; + double d0 = 0.; + { + double bestOut = std::numeric_limits::infinity(); + Halfedge curr = startVert.halfedge(); + do { + if (!curr.isInterior()) break; + const double width = geom.cornerAngles[curr.corner()] / sumScale; // in the tangent space's angles + double off = theta - arg(geom.halfedgeVectorsInVertex[curr]); + off -= 2. * M_PI * std::floor(off / (2. * M_PI)); + if (off <= width) { + he = curr; + d0 = off * sumScale; + break; + } + // closest corner, in case rounding leaves the direction in none + const double outBy = std::fmin(off - width, 2. * M_PI - off); + if (outBy < bestOut) { + bestOut = outBy; + he = curr; + d0 = (2. * M_PI - off < off - width) ? 0. : geom.cornerAngles[curr.corner()]; + } + curr = curr.next().next().twin(); + } while (curr != startVert.halfedge()); + } + if (he == Halfedge()) return false; + double c0 = geom.cornerAngles[he.corner()]; + if (positions) { // the corner from the positions: from the lengths it is off by about 1e-8 near 0 and pi + const Vector3 p0 = positions->inputVertexPositions[startVert]; + const Vector3 e1 = positions->inputVertexPositions[he.tipVertex()] - p0; + const Vector3 e2 = positions->inputVertexPositions[he.next().tipVertex()] - p0; + if (norm(e1) > 0. && norm(e2) > 0.) c0 = std::atan2(norm(cross(e1, e2)), dot(e1, e2)); + } + d0 = clamp(d0, 0., c0); + + // The current face: the line enters it at the corner T (first face: the start) or through an edge, and leaves + // it through heCross (u -> w, u right of the line, w left). Positions in the segment's frame. + struct State { + Halfedge heCross; + Vertex u, w, T; + Vector2 pu, pw, pT; + size_t nOut, steps; + }; + State st; + st.heCross = he.next(); + st.u = he.tipVertex(); + st.w = he.next().tipVertex(); + st.T = startVert; + const double lu = geom.edgeLengths[he.edge()], lw = geom.edgeLengths[he.next().next().edge()]; + st.pu = Vector2{lu * std::cos(d0), -lu * std::sin(d0)}; + st.pw = Vector2{lw * std::cos(c0 - d0), lw * std::sin(c0 - d0)}; + st.pT = Vector2{0., 0.}; + st.steps = 0; + out.push_back(SurfacePoint(startVert)); + st.nOut = out.size(); + + // Where the line passes within rounding of a vertex, the side taken; the other one is tried on failure + struct Decision { + State before; + Halfedge heUR, heRW; + Vertex r; + Vector2 pr; + }; + std::vector decisions; + const double hitTol = 1e-9 * length; + size_t attempts = 0; + const size_t maxAttempts = 64; + double nearest = std::numeric_limits::infinity(); + + // the edge each crossing lies on, in the segment's frame (out[k + 1] on sides[k]), and whether the crossed + // halfedge is its edge's halfedge (the crossing's t from u) + std::vector> sides; + std::vector sideForward; + + // endVert laid out at pB: arrived if the chord from the start to pB has the segment's length (within arriveTol) + // and crosses every edge crossed so far between its ends, i.e. is a straight line from startVert to endVert + // in these faces. The crossings are moved onto the chord, so the path ends exactly at endVert. Otherwise + // pB, if nearest the line so far, aims the next unfolding. + auto tryArrive = [&](Vector2 pB) { + const double lB = norm(pB); + if (std::fabs(lB - length) > arriveTol * length) return false; + std::vector ts(sides.size()); + for (size_t k = 0; k < sides.size(); k++) { + const Vector2 pu = sides[k].first, pw = sides[k].second, d = pw - pu; + const double den = cross(pB, d); + if (den == 0.) return false; + const double t = -cross(pB, pu) / den; // along u -> w + const double lam = cross(pu, d) / den; // along the chord + const double eps = 1e-9; + if (!(t >= -eps && t <= 1. + eps && lam >= -eps && lam <= 1. + eps)) return false; + ts[k] = clamp(t, 0., 1.); + } + for (size_t k = 0; k < sides.size(); k++) out[k + 1].tEdge = sideForward[k] ? ts[k] : 1. - ts[k]; + out.push_back(SurfacePoint(endVert)); + return true; + }; + auto noteEnd = [&](Vector2 pB) { + if (std::fabs(pB.y) < nearest) { + nearest = std::fabs(pB.y); + aimOffset = std::atan2(pB.y, pB.x); + } + }; + + // the line takes the side of r given by left (exit through u -> r) or right (through r -> w) + auto take = [&](State& s, Halfedge heUR, Halfedge heRW, Vertex r, Vector2 pr, bool left) { + if (left) { + s.T = s.w; + s.pT = s.pw; + s.heCross = heUR; + s.w = r; + s.pw = pr; + } else { + s.T = s.u; + s.pT = s.pu; + s.heCross = heRW; + s.u = r; + s.pu = pr; + } + }; + // on failure: back to the last open decision, other side; false if none is left + auto backtrack = [&]() { + if (decisions.empty() || ++attempts > maxAttempts) return false; + Decision d = decisions.back(); + decisions.pop_back(); + st = d.before; + out.resize(st.nOut); + sides.resize(st.nOut - 1); + sideForward.resize(st.nOut - 1); + take(st, d.heUR, d.heRW, d.r, d.pr, !(d.pr.y > 0.)); + return true; + }; + + while (true) { + // where the line leaves the current face + const double den = st.pu.y - st.pw.y; + const double s = clamp(den != 0. ? st.pu.y / den : 0.5, 0., 1.); + const double xc = st.pu.x + s * (st.pw.x - st.pu.x); + bool failed = false; + if (xc >= length * (1. - arriveTol)) { + // the line ends in this face (or on its exit edge): endVert should be a corner here + const Vertex vs[3] = {st.u, st.w, st.T}; + const Vector2 ps[3] = {st.pu, st.pw, st.pT}; + for (int k = 0; k < 3; k++) + if (vs[k] == endVert) { + if (tryArrive(ps[k])) return true; + noteEnd(ps[k]); + } + failed = xc > length * (1. + arriveTol); + } + Halfedge heE = st.heCross.twin(); + if (!failed && (!heE.isInterior() || ++st.steps > maxIters)) failed = true; + if (failed) { + if (!backtrack()) { + out.clear(); + return false; + } + continue; + } + const bool fwd = st.heCross == st.heCross.edge().halfedge(); + out.push_back(SurfacePoint(st.heCross.edge(), fwd ? s : 1. - s)); + sides.emplace_back(st.pu, st.pw); + sideForward.push_back(fwd); + + // the next face: r beyond u -> w, from the vertex positions if the geometry has them (the edge lengths + // lose the width of a sliver), else from the edge lengths + Halfedge heUR = heE.next(), heRW = heUR.next(); + Vertex r = heRW.tailVertex(); + Vector2 e = st.pw - st.pu; + const double L = norm(e); + const double a = geom.edgeLengths[heUR.edge()], b = geom.edgeLengths[heRW.edge()]; + Vector2 pr; + if (L > 0.) { + const Vector2 eh = e / L, nh{eh.y, -eh.x}; // nh: right of u -> w, away from the current face + double along = (a * a - b * b + L * L) / (2. * L); + double h = std::sqrt(std::fmax(0., a * a - along * along)); + if (positions) { + const Vector3 pU = positions->inputVertexPositions[st.u]; + const Vector3 uw = positions->inputVertexPositions[st.w] - pU, ur = positions->inputVertexPositions[r] - pU; + const double l3 = norm(uw); + if (l3 > 0.) { + along = dot(ur, uw) / l3; + h = norm(cross(uw, ur)) / l3; + } + } + pr = st.pu + along * eh + h * nh; + } else { + pr = st.pu + Vector2{a, 0.}; + } + if (r == endVert) { + if (tryArrive(pr)) return true; + noteEnd(pr); + } + const bool left = pr.y > 0.; + if (std::fabs(pr.y) <= hitTol) { + State before = st; + before.nOut = out.size(); // the crossing of heCross stays + decisions.push_back(Decision{before, heUR, heRW, r, pr}); + } + take(st, heUR, heRW, r, pr, left); + st.nOut = out.size(); + } +} +} // namespace + +std::vector traceGeodesicByUnfolding(IntrinsicGeometryInterface& geom, Vertex startVert, + Vector2 traceVec, Vertex endVert, double arriveTol, + size_t maxIters) { + std::vector out; + const double length = norm(traceVec); + if (!(length > 0) || !std::isfinite(length) || startVert == endVert) return out; + geom.requireEdgeLengths(); + geom.requireCornerAngles(); + geom.requireVertexAngleSums(); + geom.requireHalfedgeVectorsInVertex(); + if (maxIters == INVALID_IND) maxIters = 10 * startVert.getMesh()->nFaces(); + VertexPositionGeometry* positions = dynamic_cast(&geom); + + // The direction aims the line; where it misses endVert (a direction carried through degenerate + // faces can be off while the length is right), it is turned toward endVert's unfolded position and + // the line unfolded again + const double sumScale = geom.vertexAngleSums[startVert] / (startVert.isBoundary() ? M_PI : 2. * M_PI); + double theta = std::atan2(traceVec.y, traceVec.x); + const int maxShots = 8; + bool ok = false; + for (int shot = 0; shot < maxShots && !ok; shot++) { + double aim; + ok = unfoldLine(geom, positions, startVert, theta, length, endVert, arriveTol, maxIters, out, aim); + if (!ok) { + if (!std::isfinite(aim)) break; + theta += aim / sumScale; + } + } + if (!ok) out.clear(); + + geom.unrequireEdgeLengths(); + geom.unrequireCornerAngles(); + geom.unrequireVertexAngleSums(); + geom.unrequireHalfedgeVectorsInVertex(); + return out; +} + +std::vector traceGeodesicBetweenVertices(IntrinsicGeometryInterface& geom, Vertex startVert, + Vector2 traceVec, Vertex endVert, + const TraceOptions& traceOptions, double lengthTol) { + TraceOptions options = traceOptions; + options.includePath = true; + TraceGeodesicResult result = traceGeodesic(geom, SurfacePoint(startVert), traceVec, options); + const bool trimmed = trimTraceResult(result, endVert); + if (trimmed) { + result.pathPoints.push_back(SurfacePoint(endVert)); + const double len = norm(traceVec); + if (std::fabs(surfacePathLength(geom, result.pathPoints) - len) <= lengthTol * len) return result.pathPoints; + } + std::vector unfolded = traceGeodesicByUnfolding(geom, startVert, traceVec, endVert, 1e-6, options.maxIters); + if (!unfolded.empty()) return unfolded; + if (trimmed) return result.pathPoints; + return traceGeodesic(geom, SurfacePoint(startVert), traceVec, options).pathPoints; +} + + } // namespace surface } // namespace geometrycentral diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 2e7aeae3..f7734d62 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -109,6 +109,7 @@ set(TEST_SRCS src/poisson_disk_sampler_test.cpp src/polygon_operators_test.cpp src/flip_geodesic_test.cpp + src/trace_geodesic_test.cpp src/stl_reader_test.cpp src/surface_misc_test.cpp src/utilities_test.cpp diff --git a/test/src/trace_geodesic_test.cpp b/test/src/trace_geodesic_test.cpp new file mode 100644 index 00000000..5ab7a59a --- /dev/null +++ b/test/src/trace_geodesic_test.cpp @@ -0,0 +1,131 @@ +// traceGeodesicBetweenVertices: a segment between two vertices that passes exactly through a vertex of a +// zero-area face, or starts in such a face's corner of angle pi, ends at its end vertex with its length; +// on a clean mesh the result is the trimmed face-by-face trace, unchanged; FlipEdgeNetwork::getPathPolyline through +// a zero-area face runs from the path's first vertex to its last. + +#include "geometrycentral/surface/flip_geodesics.h" +#include "geometrycentral/surface/meshio.h" +#include "geometrycentral/surface/surface_mesh_factories.h" +#include "geometrycentral/surface/trace_geodesic.h" + +#include "gtest/gtest.h" + +#include +#include + +using namespace geometrycentral; +using namespace geometrycentral::surface; + +namespace { + +// A planar patch with a zero-area face: M = (1, 0) lies on the edge A-B, the cap (A, B, M) has a corner of +// angle pi at M; C above and D below. Vertices A 0, M 1, B 2, C 3, D 4. +std::tuple, std::unique_ptr> capPatch() { + std::vector verts = {{0, 0, 0}, {1, 0, 0}, {2, 0, 0}, {1, 1, 0}, {1, -1, 0}}; + std::vector> faces = {{0, 1, 3}, {1, 2, 3}, {0, 2, 1}, {2, 0, 4}}; + return makeManifoldSurfaceMeshAndGeometry(faces, verts); +} + +Halfedge outgoing(Vertex from, Vertex to) { + for (Halfedge he : from.outgoingHalfedges()) + if (he.tipVertex() == to) return he; + return Halfedge(); +} + +void expectSegment(VertexPositionGeometry& geom, const std::vector& path, Vertex a, Vertex b, + double length) { + ASSERT_GE(path.size(), 2u); + EXPECT_TRUE(path.front().type == SurfacePointType::Vertex && path.front().vertex == a); + EXPECT_TRUE(path.back().type == SurfacePointType::Vertex && path.back().vertex == b); + const double L = surfacePathLength(geom, path); // NaN if two consecutive points share no face + EXPECT_NEAR(L, length, 1e-12 * length); + // the polyline is straight: every point on the segment a -> b + const Vector3 pa = geom.inputVertexPositions[a], pb = geom.inputVertexPositions[b]; + for (const SurfacePoint& p : path) { + const Vector3 x = p.interpolate(geom.inputVertexPositions); + EXPECT_NEAR(norm(cross(x - pa, pb - pa)) / norm(pb - pa), 0., 1e-12); + } +} + +} // namespace + +TEST(TraceGeodesicBetweenVertices, ThroughVertexOfZeroAreaFace) { + std::unique_ptr mesh; + std::unique_ptr geom; + std::tie(mesh, geom) = capPatch(); + geom->requireHalfedgeVectorsInVertex(); + Vertex C = mesh->vertex(3), D = mesh->vertex(4); + // C -> D runs straight down through M + Vector2 dir = geom->halfedgeVectorsInVertex[outgoing(C, mesh->vertex(1))].normalize(); + std::vector path = traceGeodesicBetweenVertices(*geom, C, 2. * dir, D); + expectSegment(*geom, path, C, D, 2.); +} + +TEST(TraceGeodesicBetweenVertices, FromCornerOfAnglePi) { + std::unique_ptr mesh; + std::unique_ptr geom; + std::tie(mesh, geom) = capPatch(); + geom->requireHalfedgeVectorsInVertex(); + geom->requireCornerAngles(); + Vertex M = mesh->vertex(1), D = mesh->vertex(4); + Halfedge ma = outgoing(M, mesh->vertex(0)); // the cap's corner at M spans from M -> A to M -> B + ASSERT_NEAR(geom->cornerAngles[ma.corner()], M_PI, 1e-12); + // M -> D leaves through the middle of that corner; M's angle sum is 2 pi, so the tangent space is unscaled + Vector2 dir = geom->halfedgeVectorsInVertex[ma].normalize() * Vector2::fromAngle(M_PI / 2.); + std::vector path = traceGeodesicBetweenVertices(*geom, M, dir, D); + expectSegment(*geom, path, M, D, 1.); +} + +TEST(TraceGeodesicBetweenVertices, CleanMeshUnchanged) { + // FlipOut paths on the bunny: each intrinsic edge traced as before (traceGeodesic, trimTraceResult, end + // vertex appended), bit for bit + std::unique_ptr mesh; + std::unique_ptr geom; + std::tie(mesh, geom) = readManifoldSurfaceMesh(std::string(GC_TEST_ASSETS_ABS_PATH) + "/bob_small.ply"); + std::mt19937 rng(7); + std::uniform_int_distribution pick(0, mesh->nVertices() - 1); + size_t traced = 0; + for (int q = 0; q < 10; q++) { + Vertex a = mesh->vertex(pick(rng)), b = mesh->vertex(pick(rng)); + if (a == b) continue; + auto net = FlipEdgeNetwork::constructFromDijkstraPath(*mesh, *geom, a, b); + net->iterativeShorten(); + for (Halfedge he : net->paths[0]->getHalfedgeList()) { + if (net->tri->edgeIsOriginal[he.edge()]) continue; + Vertex va = net->tri->vertexLocations[he.vertex()].vertex, vb = net->tri->vertexLocations[he.twin().vertex()].vertex; + // SignpostIntrinsicTriangulation::halfedgeVector + net->tri->requireVertexAngleSums(); + const double scale = net->tri->vertexAngleSums[he.vertex()] / (he.vertex().isBoundary() ? M_PI : 2. * M_PI); + Vector2 vec = Vector2::fromAngle(net->tri->signpostAngle[he] / scale) * net->tri->edgeLengths[he.edge()]; + TraceOptions options; + options.includePath = true; + options.maxIters = mesh->nFaces() * 10; + TraceGeodesicResult old = traceGeodesic(*geom, SurfacePoint(va), vec, options); + ASSERT_TRUE(trimTraceResult(old, vb)); + old.pathPoints.push_back(SurfacePoint(vb)); + std::vector now = traceGeodesicBetweenVertices(*geom, va, vec, vb, options); + ASSERT_EQ(now.size(), old.pathPoints.size()); + for (size_t k = 0; k < now.size(); k++) { + EXPECT_TRUE(now[k] == old.pathPoints[k]); + if (now[k].type == SurfacePointType::Edge) EXPECT_EQ(now[k].tEdge, old.pathPoints[k].tEdge); + } + traced++; + } + } + EXPECT_GT(traced, 10u); +} + +TEST(TraceGeodesicBetweenVertices, FlipOutPolylineThroughZeroAreaFace) { + // FlipOut's shortest path C -> D is the straight segment through M; its polyline must run from C to D on the + // surface, every step inside one face, with the path's length + std::unique_ptr mesh; + std::unique_ptr geom; + std::tie(mesh, geom) = capPatch(); + Vertex C = mesh->vertex(3), D = mesh->vertex(4); + auto net = FlipEdgeNetwork::constructFromDijkstraPath(*mesh, *geom, C, D); + net->iterativeShorten(); + std::vector> polyline = net->getPathPolyline(); + ASSERT_EQ(polyline.size(), 1u); + EXPECT_NEAR(net->length(), 2., 1e-6); // FlipOut stops within its angle tolerance of straight + expectSegment(*geom, polyline[0], C, D, 2.); +} From 9e96207bf0df8aca76c0b5b0457e60ff919bc937 Mon Sep 17 00:00:00 2001 From: Mattias Fagerlund Date: Fri, 30 Oct 2026 07:12:15 +0100 Subject: [PATCH 2/3] surfacePathLength: no face layout when the geometry has positions With vertex positions the length is summed from the positions; the faces are only needed to check that consecutive points share one, which is topological. Same result bit for bit, without requiring halfedgeVectorsInFace per call: the per-edge check in traceGeodesicBetweenVertices is now cheap on clean meshes. --- src/surface/trace_geodesic.cpp | 54 ++++++++++++++++++++++++++++++---- 1 file changed, 49 insertions(+), 5 deletions(-) diff --git a/src/surface/trace_geodesic.cpp b/src/surface/trace_geodesic.cpp index ad5d6921..a6303fdb 100644 --- a/src/surface/trace_geodesic.cpp +++ b/src/surface/trace_geodesic.cpp @@ -911,9 +911,56 @@ bool pointInFace(IntrinsicGeometryInterface& geom, Face f, const SurfacePoint& p } } // namespace +namespace { +// Whether surface point p lies on face f (vertex of f, edge of f, or in f): pointInFace without the layout. +bool onFace(const SurfacePoint& p, Face f) { + Halfedge h = f.halfedge(); + switch (p.type) { + case SurfacePointType::Vertex: + for (int k = 0; k < 3; k++, h = h.next()) + if (h.tailVertex() == p.vertex) return true; + return false; + case SurfacePointType::Edge: + for (int k = 0; k < 3; k++, h = h.next()) + if (h.edge() == p.edge) return true; + return false; + case SurfacePointType::Face: + return p.face == f; + } + return false; +} + +// Whether two consecutive path points share a face, over the faces of a in the order surfacePathLength tries them. +bool shareFace(const SurfacePoint& a, const SurfacePoint& b) { + switch (a.type) { + case SurfacePointType::Vertex: + for (Face f : a.vertex.adjacentFaces()) + if (onFace(b, f)) return true; + return false; + case SurfacePointType::Edge: + if (a.edge.halfedge().isInterior() && onFace(b, a.edge.halfedge().face())) return true; + if (a.edge.halfedge().twin().isInterior() && onFace(b, a.edge.halfedge().twin().face())) return true; + return false; + case SurfacePointType::Face: + return onFace(b, a.face); + } + return false; +} +} // namespace + double surfacePathLength(IntrinsicGeometryInterface& geom, const std::vector& path) { - geom.requireHalfedgeVectorsInFace(); VertexPositionGeometry* positions = dynamic_cast(&geom); + if (positions) { + // from the positions (the edge lengths lose the width of a sliver); a step's face is only checked, so no layout + const VertexData& pos = positions->inputVertexPositions; + double L = 0.; + for (size_t i = 0; i + 1 < path.size(); i++) { + if (!shareFace(path[i], path[i + 1])) return std::numeric_limits::quiet_NaN(); + L += norm(path[i + 1].interpolate(pos) - path[i].interpolate(pos)); + } + return L; + } + geom.requireHalfedgeVectorsInFace(); double L = 0.; std::vector faces; for (size_t i = 0; i + 1 < path.size(); i++) { @@ -935,10 +982,7 @@ double surfacePathLength(IntrinsicGeometryInterface& geom, const std::vectorinputVertexPositions) - - p.interpolate(positions->inputVertexPositions)) - : norm(b - a); + L += norm(b - a); shared = true; break; } From c800b93925f47dae5a4558be841497c4b558c7c0 Mon Sep 17 00:00:00 2001 From: Mattias Fagerlund Date: Fri, 30 Oct 2026 07:48:04 +0100 Subject: [PATCH 3/3] traceGeodesicBetweenVertices: keep a trace that ends at the end vertex without measuring it A trace whose untrimmed end lies within lengthTol of endVert (bounded from the edge lengths, O(1)) and that crossed only edge interiors on the way is kept as it is; only the others are measured with surfacePathLength and, if off, traced again by unfolding. At a vertex the trace picks a new direction and can still end at endVert by another route (near-coincident vertices), so such traces are still measured. The result is moved out instead of copied. Same polylines as before on every pair tested; on clean meshes the export costs the same as the old one. --- .../geometrycentral/surface/trace_geodesic.h | 10 ++-- src/surface/trace_geodesic.cpp | 56 +++++++++++++++++-- 2 files changed, 58 insertions(+), 8 deletions(-) diff --git a/include/geometrycentral/surface/trace_geodesic.h b/include/geometrycentral/surface/trace_geodesic.h index 97e06e77..1af4299d 100644 --- a/include/geometrycentral/surface/trace_geodesic.h +++ b/include/geometrycentral/surface/trace_geodesic.h @@ -55,10 +55,12 @@ bool trimTraceResult(TraceGeodesicResult& traceResult, Vertex targetVertex); // Trace a geodesic segment known to run from vertex startVert to vertex endVert (an intrinsic edge between two // vertices of the input mesh): traceVec in startVert's tangent space, its norm the segment's length. The result // runs from startVert through the edge crossings to endVert. The face-by-face trace above is used and trimmed -// with trimTraceResult; when it does not reach endVert's 1-ring, or the length of the trimmed path differs from -// norm(traceVec) by more than lengthTol relative (the trace turned off course at a degenerate face or at a vertex -// it passed through), the segment is traced again by unfolding its faces along the line, which has neither -// failure. If that fails too, the face-by-face trace is returned as before: trimmed if it trimmed, else untrimmed. +// with trimTraceResult. It is kept at once when it ends within lengthTol * norm(traceVec) of endVert (bounded from +// the edge lengths) and crossed only edge interiors on the way; otherwise it is kept if the trimmed path's length +// (surfacePathLength) is within lengthTol relative of norm(traceVec). When neither holds (the trace did not reach +// endVert's 1-ring, or turned off course at a degenerate face or at a vertex it passed through), the segment is +// traced again by unfolding its faces along the line, which has neither failure. If that fails too, the +// face-by-face trace is returned as before: trimmed if it trimmed, else untrimmed. std::vector traceGeodesicBetweenVertices(IntrinsicGeometryInterface& geom, Vertex startVert, Vector2 traceVec, Vertex endVert, const TraceOptions& traceOptions = defaultTraceOptions, diff --git a/src/surface/trace_geodesic.cpp b/src/surface/trace_geodesic.cpp index a6303fdb..f3518f9a 100644 --- a/src/surface/trace_geodesic.cpp +++ b/src/surface/trace_geodesic.cpp @@ -954,9 +954,12 @@ double surfacePathLength(IntrinsicGeometryInterface& geom, const std::vector& pos = positions->inputVertexPositions; double L = 0.; + Vector3 prev = path.empty() ? Vector3::zero() : path[0].interpolate(pos); for (size_t i = 0; i + 1 < path.size(); i++) { if (!shareFace(path[i], path[i + 1])) return std::numeric_limits::quiet_NaN(); - L += norm(path[i + 1].interpolate(pos) - path[i].interpolate(pos)); + const Vector3 next = path[i + 1].interpolate(pos); + L += norm(next - prev); + prev = next; } return L; } @@ -1254,21 +1257,66 @@ std::vector traceGeodesicByUnfolding(IntrinsicGeometryInterface& g return out; } +namespace { +// An upper bound on the distance from surface point p to vertex v, from the edge lengths: exact on an edge of v, +// by the triangle inequality over the corners in a face of v; infinite if p is not on an edge or face of v. +double endDistanceBound(IntrinsicGeometryInterface& geom, const SurfacePoint& p, Vertex v) { + const double inf = std::numeric_limits::infinity(); + switch (p.type) { + case SurfacePointType::Vertex: + return p.vertex == v ? 0. : inf; + case SurfacePointType::Edge: { + Halfedge h = p.edge.halfedge(); + geom.requireEdgeLengths(); + const double l = geom.edgeLengths[p.edge]; + geom.unrequireEdgeLengths(); + if (h.tailVertex() == v) return std::fabs(p.tEdge) * l; + if (h.tipVertex() == v) return std::fabs(1. - p.tEdge) * l; + return inf; + } + case SurfacePointType::Face: { + Halfedge h = p.face.halfedge(); + const double b[3] = {p.faceCoords.x, p.faceCoords.y, p.faceCoords.z}; // at h's tail, next's tail, next's tip + geom.requireEdgeLengths(); + double d = inf; + for (int k = 0; k < 3; k++, h = h.next()) { + if (h.tailVertex() != v) continue; + // v at corner k: the other two corners along h (corner k + 1) and along the previous halfedge (k + 2) + d = std::fabs(b[(k + 1) % 3]) * geom.edgeLengths[h.edge()] + + std::fabs(b[(k + 2) % 3]) * geom.edgeLengths[h.next().next().edge()]; + } + geom.unrequireEdgeLengths(); + return d; + } + } + return inf; +} +} // namespace + std::vector traceGeodesicBetweenVertices(IntrinsicGeometryInterface& geom, Vertex startVert, Vector2 traceVec, Vertex endVert, const TraceOptions& traceOptions, double lengthTol) { TraceOptions options = traceOptions; options.includePath = true; TraceGeodesicResult result = traceGeodesic(geom, SurfacePoint(startVert), traceVec, options); + const double len = norm(traceVec); + // A trace that ends at endVert and crossed only edge interiors on the way went right: then the trimmed trace is + // kept without measuring it. At a vertex (or an edge's end) the trace picks a new direction, and can still end at + // endVert by another route (near-coincident vertices), so such a trace is measured. + bool atEnd = !result.pathPoints.empty() && endDistanceBound(geom, result.pathPoints.back(), endVert) <= lengthTol * len; + for (size_t i = 1; atEnd && i + 1 < result.pathPoints.size(); i++) { + const SurfacePoint& p = result.pathPoints[i]; + if (p.type != SurfacePointType::Edge || !(p.tEdge > 0. && p.tEdge < 1.)) atEnd = false; + } const bool trimmed = trimTraceResult(result, endVert); if (trimmed) { result.pathPoints.push_back(SurfacePoint(endVert)); - const double len = norm(traceVec); - if (std::fabs(surfacePathLength(geom, result.pathPoints) - len) <= lengthTol * len) return result.pathPoints; + if (atEnd) return std::move(result.pathPoints); + if (std::fabs(surfacePathLength(geom, result.pathPoints) - len) <= lengthTol * len) return std::move(result.pathPoints); } std::vector unfolded = traceGeodesicByUnfolding(geom, startVert, traceVec, endVert, 1e-6, options.maxIters); if (!unfolded.empty()) return unfolded; - if (trimmed) return result.pathPoints; + if (trimmed) return std::move(result.pathPoints); return traceGeodesic(geom, SurfacePoint(startVert), traceVec, options).pathPoints; }