diff --git a/include/geometrycentral/surface/trace_geodesic.h b/include/geometrycentral/surface/trace_geodesic.h index 33540194..1af4299d 100644 --- a/include/geometrycentral/surface/trace_geodesic.h +++ b/include/geometrycentral/surface/trace_geodesic.h @@ -52,5 +52,36 @@ 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. 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, + 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..f3518f9a 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,447 @@ 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 + +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) { + 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.; + 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(); + const Vector3 next = path[i + 1].interpolate(pos); + L += norm(next - prev); + prev = next; + } + return L; + } + geom.requireHalfedgeVectorsInFace(); + 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)) { + L += 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; +} + +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)); + 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 std::move(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.); +}