[mfagerlund]

What goes wrong

On meshes with zero-area faces or near-coincident vertices, FlipEdgeNetwork::getPathPolyline() can return a
polyline that is not a path on the surface between the two endpoints. It stops short of the target vertex, two
consecutive points share no face, or its length differs from FlipEdgeNetwork::length(). The path FlipOut
computed is fine; only the export is wrong.

Why

SignpostIntrinsicTriangulation::traceIntrinsicHalfedgeAlongInput traces each intrinsic edge of the path from
its start vertex along its signpost direction for its length with traceGeodesic, then trims the trace to the
end vertex's 1-ring with trimTraceResult. When the trim fails, it keeps the untrimmed trace, which does not end
at the end vertex. When the trim succeeds, the result is not checked.

The trace goes wrong in four ways on degenerate input:

• Zero-area faces. The barycentric step in a face of zero area solves a singular system, takes a near-zero
  step and ends inside the face.
• A ray through a vertex. The trace enters a face exactly at a corner (tEdge 0 or 1), the direction is
  clamped and the trace turns off course. On the meshes below this happens at vertices of zero-area faces.
• Near-coincident vertices. The trace reaches the end vertex's position, but the 1-ring test fails or the
  trimmed path takes a detour through the duplicate.
• Drifted signposts. On near-degenerate intrinsic triangles the signpost angles drift (1e-8 to
  1e-5 rad on the pairs below) while the edge lengths stay right. On a long edge the trace then misses the end vertex.

The fix

New traceGeodesicBetweenVertices(geom, startVert, traceVec, endVert, options, lengthTol = 1e-9) in
trace_geodesic.h, used by traceIntrinsicHalfedgeAlongInput when both ends of the intrinsic edge are input
vertices:

1. Trace with traceGeodesic and trim with trimTraceResult, as before. If the trim succeeds and the trace is
   right, return it: this is the old result, bit for bit. The trace is right if

   • its untrimmed end lies within 1e-9 of the edge length from the end vertex (an O(1) bound from the edge
     lengths of the end point's edge or face) and it crossed only edge interiors on the way, or else
   • the trimmed path's length (surfacePathLength, from the vertex positions) is within 1e-9 relative of the
     edge length.

   The first test is what makes the check free on clean meshes. The second catches a trace that hit a vertex,
   turned, and still ended at the end vertex by another route (Thingi10K 106830 below, near-coincident
   vertices).

2. Otherwise, trace the edge again with traceGeodesicByUnfolding. It lays the faces out one by one along the
   line in the segment's frame, from the vertex positions when the geometry has them (edge lengths lose the
   width of a sliver), and picks each exit edge by the side of the line the new vertex lies on. Where the line
   passes within 1e-9 of a vertex, both sides are tried. It arrives where the end vertex is laid out at a point
   whose chord from the start has the edge's length (1e-6 relative) and crosses every edge crossed so far; the
   crossings are then put on that chord, so the path ends exactly at the end vertex. The direction only aims
   the line: where it misses, it is turned toward the end vertex's unfolded position and unfolded again, up to
   8 times.

3. If that fails too, return the old result (trimmed if the trim succeeded, else the untrimmed trace).

surfacePathLength(geom, path) measures a path of surface points and returns NaN if two consecutive points
share no face. With vertex positions it needs no face layout (the shared face is a topological test).

Clean meshes are unchanged. Step 1 returns the old trace whenever it was already a path to the end vertex
with the edge's length, which is every edge on a clean mesh. The test CleanMeshUnchanged compares the old
trace and the new function point by point on FlipOut paths on bob_small.ply. On 14 standard meshes (300
random pairs each, 4200 paths) the polylines are bit-identical before and after, and so are FlipOut's lengths.

Cost. On the same 14 meshes (300 pairs each, single thread, old and new export alternated in one process,
min of 3 per pair): getPathPolyline takes 0.999x the old export's time (geometric mean over meshes; 0.979x to
1.025x per mesh), and iterativeShorten plus getPathPolyline 1.0000x (0.9998x to 1.0003x). The result is moved
out instead of copied, which pays for the check. Measuring every trace with surfacePathLength instead (no
first test) costs 1.20x on getPathPolyline: about 45 ns per path point against about 700 ns per point for
the trace itself.

Tests

test/src/trace_geodesic_test.cpp, 4 tests in TraceGeodesicBetweenVertices:

• ThroughVertexOfZeroAreaFace: a planar patch where vertex M = (1, 0, 0) lies on edge A-B, so face (A, B, M)
  has zero area; the segment C -> D runs straight down through M.
• FromCornerOfAnglePi: the segment M -> D, leaving M through the zero-area face's corner of angle pi.
• FlipOutPolylineThroughZeroAreaFace: the public path. FlipEdgeNetwork::constructFromDijkstraPath C -> D on
  the same patch, iterativeShorten, getPathPolyline; the polyline must run from C to D, step inside faces
  and have length 2.
• CleanMeshUnchanged: bit-identity with the old trace on a clean mesh.

Upstream's suite with the change: 205 of 205 pass (201 existing, 4 new).

Before the fix (the test file built against master; master lacks the two new functions, so the test file was
compiled with a shim that defines traceGeodesicBetweenVertices as master's export, i.e. traceGeodesic +
trimTraceResult + end vertex appended, and surfacePathLength as in the fix), 3 of the 4 fail:

[ RUN      ] TraceGeodesicBetweenVertices.ThroughVertexOfZeroAreaFace
trace_geodesic_test.cpp(40): error: Value of: path.back().type == SurfacePointType::Vertex && path.back().vertex == b
  Actual: false
Expected: true
trace_geodesic_test.cpp(42): error: The difference between L and length is 0.99999999899999992, which exceeds 1e-12 * length, where
L evaluates to 1.0000000010000001,
length evaluates to 2
[  FAILED  ] TraceGeodesicBetweenVertices.ThroughVertexOfZeroAreaFace
[ RUN      ] TraceGeodesicBetweenVertices.FromCornerOfAnglePi
trace_geodesic_test.cpp(40): error: Value of: path.back().type == SurfacePointType::Vertex && path.back().vertex == b
  Actual: false
Expected: true
trace_geodesic_test.cpp(42): error: The difference between L and length is 1, which exceeds 1e-12 * length, where
L evaluates to 0,
length evaluates to 1
[  FAILED  ] TraceGeodesicBetweenVertices.FromCornerOfAnglePi
[ RUN      ] TraceGeodesicBetweenVertices.CleanMeshUnchanged
[       OK ] TraceGeodesicBetweenVertices.CleanMeshUnchanged
[ RUN      ] TraceGeodesicBetweenVertices.FlipOutPolylineThroughZeroAreaFace
trace_geodesic_test.cpp(40): error: Value of: path.back().type == SurfacePointType::Vertex && path.back().vertex == b
  Actual: false
Expected: true
trace_geodesic_test.cpp(42): error: The difference between L and length is 0.99999998946328805, which exceeds 1e-12 * length, where
L evaluates to 1.0000000105367119,
length evaluates to 2
trace_geodesic_test.cpp(47): error: The difference between norm(cross(x - pa, pb - pa)) / norm(pb - pa) and 0. is 1.0536712058417663e-08, which exceeds 1e-12
[  FAILED  ] TraceGeodesicBetweenVertices.FlipOutPolylineThroughZeroAreaFace

On master, FlipOut's polyline from C to D stops at length 1 (at M, the vertex of the zero-area face) instead
of reaching D at length 2. With the fix all 4 pass.

Thingi10K

419 vertex pairs on 188 Thingi10K models where getPathPolyline was invalid on master. The meshes are the
Thingi10K models with vertices of identical coordinates merged, vertices sorted lexicographically and unused
vertices removed (that fixes the vertex indices below, 0-based); no other repair. Per pair: the shipped
FlipOut as a user runs it (constructFromDijkstraPath, iterativeShorten, getPathPolyline, default
settings), on master and on the branch. A polyline is valid if it starts at v0, ends at v1, every two
consecutive points share a face, and its length from the vertex positions equals length() to 1e-9
relative.

polyline	master	branch
valid	0	300
does not end at v1	96	8
two consecutive points share no face	211	25
length off by 1e-9 .. 1e-8	27	68
length off by 1e-8 .. 1e-6	13	13
length off by 1e-6 or more	72	5

FlipOut's reported length is bit-identical before and after on all 419 pairs; only the export changed.

Every number here is the same with and without the O(1) test in step 1: the 419 polylines are identical point
for point. The 419 split in two sets:

• 214 pairs where FlipOut converged (its length agrees with an independent geodesic solver to 1e-9).
  Invalid 214 -> 7.
  The 7 left (Thingi10K 109089, 270395, 285376, 331335, 521601, 331795, 331796) are intrinsic edges with no
  straight trace on the input surface: the end vertex is never laid out near the line, or only at a chord
  1.9e-4 longer than the intrinsic length. 6 of the 7 become valid when vertices closer than 1e-10 of the
  bounding box diagonal are merged first.
• 205 pairs where FlipOut stops on a longer path (it refuses a flip at a degenerate triangle). Invalid
  205 -> 112. Of the 112, 81 are length mismatches below 1e-6: there the polyline is a straight chord through
  the input faces and FlipOut's summed intrinsic edge lengths have drifted from the surface distance by
  1e-9 .. 1e-6 relative on near-degenerate intrinsic triangles. The other 31 are untraceable edges as above,
  27 of them valid once near-coincident vertices are merged.

Representative pairs:

Thingi10K file	v0	v1	cause	FlipOut length	master	branch
225354	60	480	zero-area face	50.277717	points share no face	valid (1e-15)
81131	61	483	zero-area face	13.152384	points share no face	valid (6e-11)
225354	740	680	ray through a vertex of a zero-area face	21.873890	stops short of v1 (polyline 18.81)	valid (5e-14)
225354	994	191	ray through a vertex of a zero-area face	100.244701	stops short of v1 (polyline 13.19)	valid (1e-15)
72228	805	718	near-coincident vertices (edge 3e-18 of the bbox diagonal)	12.773048	length off by 9.1e-5	valid (4e-15)
65394	427	909	near-coincident vertices (edge 9e-19 of the diagonal)	73.568874	length off by 9.2e-4	valid (2e-14)
106830	5209	4789	near-coincident vertices (edge 3e-19 of the diagonal); the trace hits one, turns by 0.04 rad and still ends at v1	16.419142	length off by 8.4e-3	valid (9e-14)
47732	2233	2378	drifted signpost (7e-8 rad)	14.585373	misses v1 by 3.5e-5 relative	valid (3e-15)
62989	444	946	drifted signpost (8e-6 rad)	66.973759	length off by 7.5e-9	valid (8e-12)
521601	1941	1935	untraceable intrinsic edge	33.925854	points share no face	points share no face (unchanged)

In brackets for "valid": the relative difference between the polyline's length and length(). Causes are
taken from the first intrinsic edge that the old export got wrong: the faces its trace touched and the two end
fans (a face with 2 area / longest edge^2 below 1e-12 is zero-area; an edge below 1e-10 of the bounding box
diagonal is near-coincident; a trace through a corner or an edge end is a vertex hit; otherwise the angle
between the old trace's first step and the new one).

How to reproduce

#include "geometrycentral/surface/flip_geodesics.h"
#include "geometrycentral/surface/meshio.h"
using namespace geometrycentral;
using namespace geometrycentral::surface;

std::unique_ptr<ManifoldSurfaceMesh> mesh;
std::unique_ptr<VertexPositionGeometry> geom;
std::tie(mesh, geom) = readManifoldSurfaceMesh("225354.ply"); // Thingi10K 225354, prepared as above
auto net = FlipEdgeNetwork::constructFromDijkstraPath(*mesh, *geom, mesh->vertex(994), mesh->vertex(191));
net->iterativeShorten();
std::vector<SurfacePoint> pts = net->getPathPolyline()[0];
// master: pts.back() is not vertex 191; the polyline is 13.19 long, net->length() is 100.24

Or, without Thingi10K, the patch from the tests: vertices (0,0,0), (1,0,0), (2,0,0), (1,1,0), (1,-1,0), faces
{0,1,3}, {1,2,3}, {0,2,1}, {2,0,4}, FlipOut from vertex 3 to vertex 4.

Note

The FlipOutContext export in #259 traces intrinsic edges the same way and will call
traceGeodesicBetweenVertices once this lands.