From e7cf243df411c894be01696605dfb80616bdab85 Mon Sep 17 00:00:00 2001 From: nzfeng Date: Fri, 6 Jun 2025 23:45:55 -0400 Subject: [PATCH 1/6] All-vertex case working; fix bug in barycentric vector sharedFace --- .../surface/barycentric_vector.ipp | 2 +- .../surface/fast_marching_method.h | 6 + src/surface/fast_marching_method.cpp | 150 +++++++++++++++--- 3 files changed, 132 insertions(+), 26 deletions(-) diff --git a/include/geometrycentral/surface/barycentric_vector.ipp b/include/geometrycentral/surface/barycentric_vector.ipp index a0109767..637e3892 100644 --- a/include/geometrycentral/surface/barycentric_vector.ipp +++ b/include/geometrycentral/surface/barycentric_vector.ipp @@ -531,8 +531,8 @@ inline Face sharedFace(const BarycentricVector& u, const BarycentricVector& w) { case BarycentricVectorType::Vertex: for (Vertex v : u.face.adjacentVertices()) { if (v == w.vertex) return u.face; - break; } + break; } break; } diff --git a/include/geometrycentral/surface/fast_marching_method.h b/include/geometrycentral/surface/fast_marching_method.h index d56cdab6..0128aa05 100644 --- a/include/geometrycentral/surface/fast_marching_method.h +++ b/include/geometrycentral/surface/fast_marching_method.h @@ -1,6 +1,8 @@ #pragma once +#include "geometrycentral/surface/barycentric_vector.h" #include "geometrycentral/surface/intrinsic_geometry_interface.h" +#include "geometrycentral/surface/surface_point.h" #include "geometrycentral/utilities/utilities.h" #include @@ -11,6 +13,10 @@ namespace geometrycentral { namespace surface { +VertexData FMMDistance(IntrinsicGeometryInterface& geometry, + const std::vector>>& initialDistances, + bool sign = false); + VertexData FMMDistance(IntrinsicGeometryInterface& geometry, const std::vector>& initialDistances); diff --git a/src/surface/fast_marching_method.cpp b/src/surface/fast_marching_method.cpp index 65eed588..a2340bcc 100644 --- a/src/surface/fast_marching_method.cpp +++ b/src/surface/fast_marching_method.cpp @@ -11,7 +11,7 @@ namespace { // The super fun quadratic distance function in the Fast Marching Method on triangle meshes // TODO parameter c isn't actually defined in paper, so I guessed that it was an error -double eikonalDistanceSubroutine(double a, double b, double theta, double dA, double dB) { +double eikonalDistanceSubroutine(double a, double b, double theta, double dA, double dB, int sign) { if (theta <= PI / 2.0) { @@ -24,14 +24,14 @@ double eikonalDistanceSubroutine(double a, double b, double theta, double dA, do double quadB = 2 * b * u * (a * cTheta - b); double quadC = b * b * (u * u - a * a * sTheta2); double sqrtVal = std::sqrt(quadB * quadB - 4 * quadA * quadC); - // double tVals[] = {(-quadB + sqrtVal) / (2*quadA), // seems to always be the first one - // (-quadB - sqrtVal) / (2*quadA)}; - - double t = (-quadB + sqrtVal) / (2 * quadA); - if (u < t && a * cTheta < b * (t - u) / t && b * (t - u) / t < a / cTheta) { + double tVals[] = {(-quadB + sqrtVal) / (2 * quadA), (-quadB - sqrtVal) / (2 * quadA)}; + double t = sign > 0 ? tVals[0] : tVals[1]; + // double t = (-quadB + sqrtVal) / (2 * quadA); + double y = b * (t - u) / t; + if (u < sign * t && a * cTheta < y && y < a / cTheta) { return dA + t; } else { - return std::min(b + dA, a + dB); + return std::min(b * sign + dA, a * sign + dB); } } @@ -44,35 +44,122 @@ double eikonalDistanceSubroutine(double a, double b, double theta, double dA, do double altitude = 2 * area / c; // distance to base, must be inside triangle since obtuse double baseDist = maxDist + altitude; - return std::min({b + dA, a + dB, baseDist}); + return std::min({b * sign + dA, a * sign + dB, baseDist}); } } } // namespace - VertexData FMMDistance(IntrinsicGeometryInterface& geometry, - const std::vector>& initialDistances) { + const std::vector>>& initialDistances, + bool sign) { typedef std::pair Entry; SurfaceMesh& mesh = geometry.mesh; geometry.requireEdgeLengths(); geometry.requireCornerAngles(); - + // TODO this could handle nonmanifold geometry with a few small tweaks - if(!mesh.isManifold()) { + if (!mesh.isManifold()) { throw std::runtime_error("handling of nonmanifold mesh not yet implemented"); } // Initialize VertexData distances(mesh, std::numeric_limits::infinity()); + VertexData signs(mesh, 1); VertexData finalized(mesh, false); - std::priority_queue, std::greater> frontierPQ; - for (auto& x : initialDistances) { - frontierPQ.push(std::make_pair(x.second, x.first)); + auto cmp = [](Entry left, Entry right) { return std::abs(left.first) > std::abs(right.first); }; + std::priority_queue, decltype(cmp)> frontierPQ(cmp); + // Initialize signs + if (sign) { + for (Vertex v : mesh.vertices()) signs[v] = 0; + for (auto& curve : initialDistances) { + size_t nNodes = curve.size(); + for (size_t i = 0; i < nNodes - 1; i++) { + const SurfacePoint& pA = curve[i].first; + const SurfacePoint& pB = curve[i + 1].first; + Edge commonEdge = sharedEdge(pA, pB); + if (commonEdge != Edge()) { + Halfedge he = commonEdge.halfedge(); + signs[he.next().tipVertex()] = (he.vertex() == pA.vertex) ? 1 : -1; + signs[he.twin().next().tipVertex()] = (he.vertex() != pA.vertex) ? 1 : -1; + signs[pA.vertex] = 1; + signs[pB.vertex] = 1; + } else { + Face commonFace = sharedFace(pA, pB); + if (commonFace == Face()) { + throw std::logic_error("For signed fast marching distance, each curve segment must share a common face."); + } + BarycentricVector tangent(pA, pB); + BarycentricVector normal = tangent.rotate90(geometry); + for (Vertex v : commonFace.adjacentVertices()) { + BarycentricVector u(SurfacePoint(v), pA); + signs[v] = (dot(geometry, normal, u) > 0) ? 1 : -1; + } + } + } + // Fill in the signs of faces around the fan of any vertices. + for (size_t i = 0; i < nNodes; i++) { + const SurfacePoint& p = curve[i].first; + if (p.type != SurfacePointType::Vertex) continue; + Halfedge startHe = p.vertex.halfedge(); + Halfedge currHe = startHe.next().next().twin(); + while (currHe != startHe) { + if (signs[currHe.tipVertex()] == 0 || signs[currHe.tipVertex()] == signs[p.vertex] || + signs[currHe.next().tipVertex()] != 0) { + currHe = currHe.next().next().twin(); + continue; + } + signs[currHe.next().tipVertex()] = -1; + currHe = currHe.next().next().twin(); + } + } + } + } + // Initialize distances + for (auto& curve : initialDistances) { + for (auto& x : curve) { + const SurfacePoint& p = x.first; + switch (p.type) { + case (SurfacePointType::Vertex): { + frontierPQ.push(std::make_pair(x.second, p.vertex)); + break; + } + case (SurfacePointType::Edge): { + const Vertex& vA = p.edge.firstVertex(); + const Vertex& vB = p.edge.secondVertex(); + double l = geometry.edgeLengths[p.edge]; + frontierPQ.push(std::make_pair(x.second + signs[vA] * p.tEdge * l, vA)); + frontierPQ.push(std::make_pair(x.second + signs[vB] * (1. - p.tEdge) * l, vB)); + break; + } + case (SurfacePointType::Face): { + Halfedge he = p.face.halfedge(); + const Vertex& vA = he.vertex(); + const Vertex& vB = he.next().vertex(); + const Vertex& vC = he.next().next().vertex(); + double lAB = geometry.edgeLengths[he.edge()]; + double lBC = geometry.edgeLengths[he.next().edge()]; + double lCA = geometry.edgeLengths[he.next().next().edge()]; + double lAB2 = lAB * lAB; + double lBC2 = lBC * lBC; + double lCA2 = lCA * lCA; + double u = p.faceCoords[0]; + double v = p.faceCoords[1]; + double w = p.faceCoords[2]; + double dist2_A = lAB2 * (v * (1. - u)) + lCA2 * (w * (1. - u)) - lAB2 * v * w; // squared distance from p to vA + double dist2_B = lAB2 * (u * (1. - v)) + lBC2 * (w * (1. - v)) - lCA2 * u * w; // squared distance from p to vB + double dist2_C = lCA2 * (u * (1. - w)) + lBC2 * (v * (1. - w)) - lBC2 * u * v; // squared distance from p to vC + frontierPQ.push(std::make_pair(x.second + signs[vA] * std::sqrt(dist2_A), vA)); + frontierPQ.push(std::make_pair(x.second + signs[vB] * std::sqrt(dist2_B), vB)); + frontierPQ.push(std::make_pair(x.second + signs[vC] * std::sqrt(dist2_C), vC)); + break; + } + } + } } size_t nFound = 0; size_t nVert = mesh.nVertices(); @@ -95,16 +182,16 @@ VertexData FMMDistance(IntrinsicGeometryInterface& geometry, finalized[currV] = true; nFound++; - // Add any eligible neighbors for (Halfedge he : currV.incomingHalfedges()) { Vertex neighVert = he.vertex(); // Add with length if (!finalized[neighVert]) { - double newDist = currDist + geometry.edgeLengths[he.edge()]; - if (newDist < distances[neighVert]) { - frontierPQ.push(std::make_pair(currDist + geometry.edgeLengths[he.edge()], neighVert)); + if (signs[neighVert] == 0) signs[neighVert] = signs[currV]; + double newDist = currDist + signs[neighVert] * geometry.edgeLengths[he.edge()]; + if (std::abs(newDist) < std::abs(distances[neighVert])) { + frontierPQ.push(std::make_pair(newDist, neighVert)); distances[neighVert] = newDist; } continue; @@ -121,9 +208,9 @@ VertexData FMMDistance(IntrinsicGeometryInterface& geometry, double lenA = geometry.edgeLengths[he.next().edge()]; double distA = distances[neighVert]; double theta = geometry.cornerAngles[he.next().next().corner()]; - double newDist = eikonalDistanceSubroutine(lenA, lenB, theta, distA, distB); - - if (newDist < distances[newVert]) { + if (signs[newVert] == 0) signs[newVert] = (signs[currV] != 0) ? signs[currV] : signs[he.next().vertex()]; + double newDist = eikonalDistanceSubroutine(lenA, lenB, theta, distA, distB, signs[newVert]); + if (std::abs(newDist) < std::abs(distances[newVert])) { frontierPQ.push(std::make_pair(newDist, newVert)); distances[newVert] = newDist; } @@ -142,9 +229,9 @@ VertexData FMMDistance(IntrinsicGeometryInterface& geometry, double lenA = geometry.edgeLengths[heT.next().next().edge()]; double distA = distances[neighVert]; double theta = geometry.cornerAngles[heT.next().next().corner()]; - double newDist = eikonalDistanceSubroutine(lenA, lenB, theta, distA, distB); - - if (newDist < distances[newVert]) { + if (signs[newVert] == 0) signs[newVert] = (signs[currV] != 0) ? signs[currV] : signs[he.next().vertex()]; + double newDist = eikonalDistanceSubroutine(lenA, lenB, theta, distA, distB, signs[newVert]); + if (std::abs(newDist) < std::abs(distances[newVert])) { frontierPQ.push(std::make_pair(newDist, newVert)); distances[newVert] = newDist; } @@ -152,7 +239,20 @@ VertexData FMMDistance(IntrinsicGeometryInterface& geometry, } } } + // for (Vertex v : mesh.vertices()) distances[v] = signs[v]; + return distances; +} + +VertexData FMMDistance(IntrinsicGeometryInterface& geometry, + const std::vector>& initialDistances, bool sign) { + + std::vector>> initialConditions; + initialConditions.emplace_back(); + for (const auto& x : initialDistances) { + initialConditions.back().emplace_back(SurfacePoint(x.first), x.second); + } + VertexData distances = FMMDistance(geometry, initialConditions, sign); return distances; } From 383efa741926fb253b2e07b005ed12cccec43dc4 Mon Sep 17 00:00:00 2001 From: nzfeng Date: Sat, 7 Jun 2025 00:22:17 -0400 Subject: [PATCH 2/6] Now debugged for zero initial conditions --- src/surface/fast_marching_method.cpp | 33 ++++++++++++++++++---------- 1 file changed, 22 insertions(+), 11 deletions(-) diff --git a/src/surface/fast_marching_method.cpp b/src/surface/fast_marching_method.cpp index a2340bcc..5f3f94a8 100644 --- a/src/surface/fast_marching_method.cpp +++ b/src/surface/fast_marching_method.cpp @@ -83,11 +83,11 @@ VertexData FMMDistance(IntrinsicGeometryInterface& geometry, const SurfacePoint& pB = curve[i + 1].first; Edge commonEdge = sharedEdge(pA, pB); if (commonEdge != Edge()) { + // Assign +/- signs to the "third" vertices of each face straddling this edge. + // These vertices might themselves lie on the curve, in which case we overwrite them below. Halfedge he = commonEdge.halfedge(); signs[he.next().tipVertex()] = (he.vertex() == pA.vertex) ? 1 : -1; - signs[he.twin().next().tipVertex()] = (he.vertex() != pA.vertex) ? 1 : -1; - signs[pA.vertex] = 1; - signs[pB.vertex] = 1; + signs[he.twin().next().tipVertex()] = -signs[he.next().tipVertex()]; } else { Face commonFace = sharedFace(pA, pB); if (commonFace == Face()) { @@ -101,20 +101,31 @@ VertexData FMMDistance(IntrinsicGeometryInterface& geometry, } } } - // Fill in the signs of faces around the fan of any vertices. + } + // Vertices on the curve are always assumed to have positive sign. + for (auto& curve : initialDistances) { + size_t nNodes = curve.size(); + for (size_t i = 0; i < nNodes; i++) { + const SurfacePoint& p = curve[i].first; + if (p.type != SurfacePointType::Vertex) continue; + signs[p.vertex] = 1; + } + } + // Fill in the signs of faces around the fan of any vertices. + for (auto& curve : initialDistances) { + size_t nNodes = curve.size(); for (size_t i = 0; i < nNodes; i++) { const SurfacePoint& p = curve[i].first; if (p.type != SurfacePointType::Vertex) continue; Halfedge startHe = p.vertex.halfedge(); - Halfedge currHe = startHe.next().next().twin(); - while (currHe != startHe) { - if (signs[currHe.tipVertex()] == 0 || signs[currHe.tipVertex()] == signs[p.vertex] || - signs[currHe.next().tipVertex()] != 0) { - currHe = currHe.next().next().twin(); - continue; + Halfedge currHe = startHe; + while (true) { + if (signs[currHe.tipVertex()] != 0 && signs[currHe.tipVertex()] != signs[p.vertex] && + signs[currHe.next().tipVertex()] == 0) { + signs[currHe.next().tipVertex()] = -1; } - signs[currHe.next().tipVertex()] = -1; currHe = currHe.next().next().twin(); + if (currHe == startHe) break; } } } From 4b4e8f9f4d411a5a9d4aa7b2c33d7c12e12480d6 Mon Sep 17 00:00:00 2001 From: nzfeng Date: Sat, 7 Jun 2025 17:46:49 -0400 Subject: [PATCH 3/6] Thought I debugged everything; now debugging wavy isobands especially for negative regions --- .../surface/algorithms/geodesic_distance.md | 42 +++++++++++++++++++ .../surface/fast_marching_method.h | 2 +- src/surface/fast_marching_method.cpp | 31 +++++++++----- 3 files changed, 64 insertions(+), 11 deletions(-) diff --git a/docs/docs/surface/algorithms/geodesic_distance.md b/docs/docs/surface/algorithms/geodesic_distance.md index d42a41fc..f852b0c3 100644 --- a/docs/docs/surface/algorithms/geodesic_distance.md +++ b/docs/docs/surface/algorithms/geodesic_distance.md @@ -174,6 +174,48 @@ Computing exact geodesic paths also allows one to compute exact [log maps](/surf Returns the gradient of the distance function at `v` (i.e. the unit tangent vector at `v` which points away from the closest source) +## (Signed) Fast Marching Method for Distance + +These routines implement the _fast marching method_ (FMM) for geodesic distance on triangle meshes, as described by [Kimmel and Sethian 1998](https://www.pnas.org/doi/pdf/10.1073/pnas.95.15.8431). + +This implementation allows the user to specify the sign of the initial velocity, allowing one to obtain either signed or unsigned distance to a set of closed curves. + +Note that as a wave-based propagation method, if the source curves are not closed, then signed distance output is not guaranteed to be well-behaved; if you're looking for robust signed distance computation, consider the [Signed Heat Method routines](/surface/algorithms/signed_heat_method). If your curves are closed, then FMM will in general be more efficient. + +`#include "geometrycentral/surface/fast_marching_method.h"` + +??? func "`#!cpp VertexData FMMDistance(IntrinsicGeometryInterface& geometry, const std::vector>>& initialDistances, bool sign = false)`" + + Compute the distance from a source set of [Surface Points](/surface/utilities/surface_point/), initialized with the given distance values. + + If `sign` is false, computes unsigned distance. If `sign` is true, computes signed distance by interpreting the source points as a set of curves and propagating their orientation (gradient is continuous across the source). + +??? func "`#!cpp VertexData FMMDistance(IntrinsicGeometryInterface& geometry, const std::vector>& initialDistances, bool sign = false)`" + + Compute the distance from a set of source vertices, initialized with the given distance values. + + If `sign` is false, computes unsigned distance. If `sign` is true, computes signed distance by interpreting the source vertices as a set of curves and propagating their orientation (gradient is continuous across the source). + +Example +```cpp +#include "geometrycentral/surface/fast_marching_method.h" +#include "geometrycentral/surface/meshio.h" + +// Load a mesh +std::unique_ptr mesh; +std::unique_ptr geometry; +std::tie(mesh, geometry) = loadMesh(filename); + +// Pick some vertices and initial distances. +std::vector> initialDistances; +initialDistances.emplace_back(mesh->vertex(7), 0.); +initialDistances.emplace_back(mesh->vertex(8), 0.); + +// Compute distance +VertexData dist = FMMDistance(*geometry, initialDistances); +/* do something useful */ +``` + ## Heat Method for Distance These routines implement the [Heat Method for Geodesic Distance](http://www.cs.cmu.edu/~kmcrane/Projects/HeatMethod/paper.pdf). This algorithm uses short time heat flow to compute distance on surfaces. Because the main burden is simply solving linear systems of equations, it tends to be faster than polyhedral schemes, especially when computing distance multiple times on the same surface. In the computational geometry sense, this method is an approximation, as the result is not precisely equal to the polyhedral distance on the surface; nonetheless it is fast and well-suited for many applications. diff --git a/include/geometrycentral/surface/fast_marching_method.h b/include/geometrycentral/surface/fast_marching_method.h index 0128aa05..f48c678d 100644 --- a/include/geometrycentral/surface/fast_marching_method.h +++ b/include/geometrycentral/surface/fast_marching_method.h @@ -18,7 +18,7 @@ VertexData FMMDistance(IntrinsicGeometryInterface& geometry, bool sign = false); VertexData FMMDistance(IntrinsicGeometryInterface& geometry, - const std::vector>& initialDistances); + const std::vector>& initialDistances, bool sign = false); } // namespace surface diff --git a/src/surface/fast_marching_method.cpp b/src/surface/fast_marching_method.cpp index 5f3f94a8..19f5fbc9 100644 --- a/src/surface/fast_marching_method.cpp +++ b/src/surface/fast_marching_method.cpp @@ -31,20 +31,22 @@ double eikonalDistanceSubroutine(double a, double b, double theta, double dA, do if (u < sign * t && a * cTheta < y && y < a / cTheta) { return dA + t; } else { - return std::min(b * sign + dA, a * sign + dB); + if (sign) return std::min(b * sign + dA, a * sign + dB); + return std::max(b * sign + dA, a * sign + dB); } } // Custom by Nick to get acceptable results in obtuse triangles without fancy unfolding else { - - double maxDist = std::max(dA, dB); // all points on base are less than this far away, by convexity + // all points on base are less than this far away, by convexity + double maxDist = (sign > 0) ? std::min(dA, dB) : std::max(dA, dB); double c = std::sqrt(a * a + b * b - 2 * a * b * std::cos(theta)); double area = 0.5 * std::sin(theta) * a * b; double altitude = 2 * area / c; // distance to base, must be inside triangle since obtuse - double baseDist = maxDist + altitude; + double baseDist = maxDist + sign * altitude; - return std::min({b * sign + dA, a * sign + dB, baseDist}); + if (sign) return std::min({b * sign + dA, a * sign + dB, baseDist}); + return std::max({b * sign + dA, a * sign + dB, baseDist}); } } @@ -71,7 +73,12 @@ VertexData FMMDistance(IntrinsicGeometryInterface& geometry, VertexData signs(mesh, 1); VertexData finalized(mesh, false); - auto cmp = [](Entry left, Entry right) { return std::abs(left.first) > std::abs(right.first); }; + auto cmp = [&signs](Entry left, Entry right) { + if ((signs[left.second] == 1) && (signs[right.second] == 1)) { + return left.first > right.first; + } + return left.first < right.first; + }; std::priority_queue, decltype(cmp)> frontierPQ(cmp); // Initialize signs if (sign) { @@ -201,7 +208,10 @@ VertexData FMMDistance(IntrinsicGeometryInterface& geometry, if (!finalized[neighVert]) { if (signs[neighVert] == 0) signs[neighVert] = signs[currV]; double newDist = currDist + signs[neighVert] * geometry.edgeLengths[he.edge()]; - if (std::abs(newDist) < std::abs(distances[neighVert])) { + // Even though this is equivalent to (signs[neighVert * (newDist - distances[neighVert]) < 0), the latter seems + // more numerically unstable. + if ((signs[neighVert] && newDist < distances[neighVert]) || + (!signs[neighVert] && newDist > distances[neighVert])) { frontierPQ.push(std::make_pair(newDist, neighVert)); distances[neighVert] = newDist; } @@ -221,7 +231,8 @@ VertexData FMMDistance(IntrinsicGeometryInterface& geometry, double theta = geometry.cornerAngles[he.next().next().corner()]; if (signs[newVert] == 0) signs[newVert] = (signs[currV] != 0) ? signs[currV] : signs[he.next().vertex()]; double newDist = eikonalDistanceSubroutine(lenA, lenB, theta, distA, distB, signs[newVert]); - if (std::abs(newDist) < std::abs(distances[newVert])) { + // if (signs[newVert] * (newDist - distances[newVert]) < 0.) { + if ((signs[newVert] && newDist < distances[newVert]) || (!signs[newVert] && newDist > distances[newVert])) { frontierPQ.push(std::make_pair(newDist, newVert)); distances[newVert] = newDist; } @@ -242,7 +253,8 @@ VertexData FMMDistance(IntrinsicGeometryInterface& geometry, double theta = geometry.cornerAngles[heT.next().next().corner()]; if (signs[newVert] == 0) signs[newVert] = (signs[currV] != 0) ? signs[currV] : signs[he.next().vertex()]; double newDist = eikonalDistanceSubroutine(lenA, lenB, theta, distA, distB, signs[newVert]); - if (std::abs(newDist) < std::abs(distances[newVert])) { + // if (signs[newVert] * (newDist - distances[newVert]) < 0.) { + if ((signs[newVert] && newDist < distances[newVert]) || (!signs[newVert] && newDist > distances[newVert])) { frontierPQ.push(std::make_pair(newDist, newVert)); distances[newVert] = newDist; } @@ -250,7 +262,6 @@ VertexData FMMDistance(IntrinsicGeometryInterface& geometry, } } } - // for (Vertex v : mesh.vertices()) distances[v] = signs[v]; return distances; } From a4311231d8d41cdc00fc0d208306d3791a47e9c4 Mon Sep 17 00:00:00 2001 From: nzfeng Date: Sat, 7 Jun 2025 18:32:50 -0400 Subject: [PATCH 4/6] Specialize to zero set case --- src/surface/fast_marching_method.cpp | 62 +++++++++++----------------- 1 file changed, 25 insertions(+), 37 deletions(-) diff --git a/src/surface/fast_marching_method.cpp b/src/surface/fast_marching_method.cpp index 19f5fbc9..87035d55 100644 --- a/src/surface/fast_marching_method.cpp +++ b/src/surface/fast_marching_method.cpp @@ -11,7 +11,7 @@ namespace { // The super fun quadratic distance function in the Fast Marching Method on triangle meshes // TODO parameter c isn't actually defined in paper, so I guessed that it was an error -double eikonalDistanceSubroutine(double a, double b, double theta, double dA, double dB, int sign) { +double eikonalDistanceSubroutine(double a, double b, double theta, double dA, double dB) { if (theta <= PI / 2.0) { @@ -24,29 +24,24 @@ double eikonalDistanceSubroutine(double a, double b, double theta, double dA, do double quadB = 2 * b * u * (a * cTheta - b); double quadC = b * b * (u * u - a * a * sTheta2); double sqrtVal = std::sqrt(quadB * quadB - 4 * quadA * quadC); - double tVals[] = {(-quadB + sqrtVal) / (2 * quadA), (-quadB - sqrtVal) / (2 * quadA)}; - double t = sign > 0 ? tVals[0] : tVals[1]; - // double t = (-quadB + sqrtVal) / (2 * quadA); - double y = b * (t - u) / t; - if (u < sign * t && a * cTheta < y && y < a / cTheta) { + double t = (-quadB + sqrtVal) / (2 * quadA); // always the positive root (negative root corresponds to slope -1) + if (u < t && a * cTheta < b * (t - u) / t && b * (t - u) / t < a / cTheta) { return dA + t; } else { - if (sign) return std::min(b * sign + dA, a * sign + dB); - return std::max(b * sign + dA, a * sign + dB); + return std::min(b + dA, a + dB); } } // Custom by Nick to get acceptable results in obtuse triangles without fancy unfolding else { - // all points on base are less than this far away, by convexity - double maxDist = (sign > 0) ? std::min(dA, dB) : std::max(dA, dB); + + double maxDist = std::max(dA, dB); // all points on base are less than this far away, by convexity double c = std::sqrt(a * a + b * b - 2 * a * b * std::cos(theta)); double area = 0.5 * std::sin(theta) * a * b; double altitude = 2 * area / c; // distance to base, must be inside triangle since obtuse - double baseDist = maxDist + sign * altitude; + double baseDist = maxDist + altitude; - if (sign) return std::min({b * sign + dA, a * sign + dB, baseDist}); - return std::max({b * sign + dA, a * sign + dB, baseDist}); + return std::min({b + dA, a + dB, baseDist}); } } @@ -73,13 +68,7 @@ VertexData FMMDistance(IntrinsicGeometryInterface& geometry, VertexData signs(mesh, 1); VertexData finalized(mesh, false); - auto cmp = [&signs](Entry left, Entry right) { - if ((signs[left.second] == 1) && (signs[right.second] == 1)) { - return left.first > right.first; - } - return left.first < right.first; - }; - std::priority_queue, decltype(cmp)> frontierPQ(cmp); + std::priority_queue, std::greater> frontierPQ; // Initialize signs if (sign) { for (Vertex v : mesh.vertices()) signs[v] = 0; @@ -137,7 +126,9 @@ VertexData FMMDistance(IntrinsicGeometryInterface& geometry, } } } - // Initialize distances + // Initialize distances. It's also possible to propagate signed distances and define eikonalDistanceSubroutine() + // accordingly, but it gets tricky initializing distances of unvisited vertices as +/-inf. So just compute unsigned + // distances, then sign -- rather than try to propagate signed distance. for (auto& curve : initialDistances) { for (auto& x : curve) { const SurfacePoint& p = x.first; @@ -150,8 +141,8 @@ VertexData FMMDistance(IntrinsicGeometryInterface& geometry, const Vertex& vA = p.edge.firstVertex(); const Vertex& vB = p.edge.secondVertex(); double l = geometry.edgeLengths[p.edge]; - frontierPQ.push(std::make_pair(x.second + signs[vA] * p.tEdge * l, vA)); - frontierPQ.push(std::make_pair(x.second + signs[vB] * (1. - p.tEdge) * l, vB)); + frontierPQ.push(std::make_pair(x.second + p.tEdge * l, vA)); + frontierPQ.push(std::make_pair(x.second + (1. - p.tEdge) * l, vB)); break; } case (SurfacePointType::Face): { @@ -171,9 +162,9 @@ VertexData FMMDistance(IntrinsicGeometryInterface& geometry, double dist2_A = lAB2 * (v * (1. - u)) + lCA2 * (w * (1. - u)) - lAB2 * v * w; // squared distance from p to vA double dist2_B = lAB2 * (u * (1. - v)) + lBC2 * (w * (1. - v)) - lCA2 * u * w; // squared distance from p to vB double dist2_C = lCA2 * (u * (1. - w)) + lBC2 * (v * (1. - w)) - lBC2 * u * v; // squared distance from p to vC - frontierPQ.push(std::make_pair(x.second + signs[vA] * std::sqrt(dist2_A), vA)); - frontierPQ.push(std::make_pair(x.second + signs[vB] * std::sqrt(dist2_B), vB)); - frontierPQ.push(std::make_pair(x.second + signs[vC] * std::sqrt(dist2_C), vC)); + frontierPQ.push(std::make_pair(x.second + std::sqrt(dist2_A), vA)); + frontierPQ.push(std::make_pair(x.second + std::sqrt(dist2_B), vB)); + frontierPQ.push(std::make_pair(x.second + std::sqrt(dist2_C), vC)); break; } } @@ -207,11 +198,8 @@ VertexData FMMDistance(IntrinsicGeometryInterface& geometry, // Add with length if (!finalized[neighVert]) { if (signs[neighVert] == 0) signs[neighVert] = signs[currV]; - double newDist = currDist + signs[neighVert] * geometry.edgeLengths[he.edge()]; - // Even though this is equivalent to (signs[neighVert * (newDist - distances[neighVert]) < 0), the latter seems - // more numerically unstable. - if ((signs[neighVert] && newDist < distances[neighVert]) || - (!signs[neighVert] && newDist > distances[neighVert])) { + double newDist = currDist + geometry.edgeLengths[he.edge()]; + if (newDist < distances[neighVert]) { frontierPQ.push(std::make_pair(newDist, neighVert)); distances[neighVert] = newDist; } @@ -230,9 +218,8 @@ VertexData FMMDistance(IntrinsicGeometryInterface& geometry, double distA = distances[neighVert]; double theta = geometry.cornerAngles[he.next().next().corner()]; if (signs[newVert] == 0) signs[newVert] = (signs[currV] != 0) ? signs[currV] : signs[he.next().vertex()]; - double newDist = eikonalDistanceSubroutine(lenA, lenB, theta, distA, distB, signs[newVert]); - // if (signs[newVert] * (newDist - distances[newVert]) < 0.) { - if ((signs[newVert] && newDist < distances[newVert]) || (!signs[newVert] && newDist > distances[newVert])) { + double newDist = eikonalDistanceSubroutine(lenA, lenB, theta, distA, distB); + if (newDist < distances[newVert]) { frontierPQ.push(std::make_pair(newDist, newVert)); distances[newVert] = newDist; } @@ -252,9 +239,8 @@ VertexData FMMDistance(IntrinsicGeometryInterface& geometry, double distA = distances[neighVert]; double theta = geometry.cornerAngles[heT.next().next().corner()]; if (signs[newVert] == 0) signs[newVert] = (signs[currV] != 0) ? signs[currV] : signs[he.next().vertex()]; - double newDist = eikonalDistanceSubroutine(lenA, lenB, theta, distA, distB, signs[newVert]); - // if (signs[newVert] * (newDist - distances[newVert]) < 0.) { - if ((signs[newVert] && newDist < distances[newVert]) || (!signs[newVert] && newDist > distances[newVert])) { + double newDist = eikonalDistanceSubroutine(lenA, lenB, theta, distA, distB); + if (newDist < distances[newVert]) { frontierPQ.push(std::make_pair(newDist, newVert)); distances[newVert] = newDist; } @@ -262,6 +248,8 @@ VertexData FMMDistance(IntrinsicGeometryInterface& geometry, } } } + // Sign unsigned distances. This will probably only do something reasonable if the curves are closed. + for (Vertex v : mesh.vertices()) distances[v] *= signs[v]; return distances; } From bf64f314c4ab76e7afde2a0e1a39a76906e23411 Mon Sep 17 00:00:00 2001 From: nzfeng Date: Sat, 7 Jun 2025 18:54:46 -0400 Subject: [PATCH 5/6] Give up on the most general implementation of solving the eikonal equation, add to the docs --- docs/docs/surface/algorithms/geodesic_distance.md | 10 ++++++---- src/surface/fast_marching_method.cpp | 5 +++-- 2 files changed, 9 insertions(+), 6 deletions(-) diff --git a/docs/docs/surface/algorithms/geodesic_distance.md b/docs/docs/surface/algorithms/geodesic_distance.md index f852b0c3..da894bcd 100644 --- a/docs/docs/surface/algorithms/geodesic_distance.md +++ b/docs/docs/surface/algorithms/geodesic_distance.md @@ -174,11 +174,11 @@ Computing exact geodesic paths also allows one to compute exact [log maps](/surf Returns the gradient of the distance function at `v` (i.e. the unit tangent vector at `v` which points away from the closest source) -## (Signed) Fast Marching Method for Distance +## Fast Marching Method for Distance These routines implement the _fast marching method_ (FMM) for geodesic distance on triangle meshes, as described by [Kimmel and Sethian 1998](https://www.pnas.org/doi/pdf/10.1073/pnas.95.15.8431). -This implementation allows the user to specify the sign of the initial velocity, allowing one to obtain either signed or unsigned distance to a set of closed curves. +This implementation allows you to obtain unsigned distance to a set of points, or signed distance to a set of curves. Note that as a wave-based propagation method, if the source curves are not closed, then signed distance output is not guaranteed to be well-behaved; if you're looking for robust signed distance computation, consider the [Signed Heat Method routines](/surface/algorithms/signed_heat_method). If your curves are closed, then FMM will in general be more efficient. @@ -188,13 +188,13 @@ Note that as a wave-based propagation method, if the source curves are not close Compute the distance from a source set of [Surface Points](/surface/utilities/surface_point/), initialized with the given distance values. - If `sign` is false, computes unsigned distance. If `sign` is true, computes signed distance by interpreting the source points as a set of curves and propagating their orientation (gradient is continuous across the source). + If `sign` is false, computes unsigned distance. If `sign` is true, computes signed distance by interpreting the source points as a set of curves and propagating their orientation. ??? func "`#!cpp VertexData FMMDistance(IntrinsicGeometryInterface& geometry, const std::vector>& initialDistances, bool sign = false)`" Compute the distance from a set of source vertices, initialized with the given distance values. - If `sign` is false, computes unsigned distance. If `sign` is true, computes signed distance by interpreting the source vertices as a set of curves and propagating their orientation (gradient is continuous across the source). + If `sign` is false, computes unsigned distance. If `sign` is true, computes signed distance by interpreting the source vertices as a set of curves and propagating their orientation. Example ```cpp @@ -216,6 +216,8 @@ VertexData dist = FMMDistance(*geometry, initialDistances); /* do something useful */ ``` +If computing signed distance, distances should be initialized to 0 on the source points. + ## Heat Method for Distance These routines implement the [Heat Method for Geodesic Distance](http://www.cs.cmu.edu/~kmcrane/Projects/HeatMethod/paper.pdf). This algorithm uses short time heat flow to compute distance on surfaces. Because the main burden is simply solving linear systems of equations, it tends to be faster than polyhedral schemes, especially when computing distance multiple times on the same surface. In the computational geometry sense, this method is an approximation, as the result is not precisely equal to the polyhedral distance on the surface; nonetheless it is fast and well-suited for many applications. diff --git a/src/surface/fast_marching_method.cpp b/src/surface/fast_marching_method.cpp index 87035d55..69a7f42c 100644 --- a/src/surface/fast_marching_method.cpp +++ b/src/surface/fast_marching_method.cpp @@ -127,8 +127,9 @@ VertexData FMMDistance(IntrinsicGeometryInterface& geometry, } } // Initialize distances. It's also possible to propagate signed distances and define eikonalDistanceSubroutine() - // accordingly, but it gets tricky initializing distances of unvisited vertices as +/-inf. So just compute unsigned - // distances, then sign -- rather than try to propagate signed distance. + // accordingly, but it gets cumbersome managing directions of inequalities, and initializing distances of unvisited + // vertices as +/-inf. So just compute unsigned distances, then sign -- rather than try to actually solve the true + // BVP. for (auto& curve : initialDistances) { for (auto& x : curve) { const SurfacePoint& p = x.first; From 9df4338886a164cbb07e595a69dd3a37c0c66240 Mon Sep 17 00:00:00 2001 From: nzfeng Date: Wed, 11 Jun 2025 16:31:49 -0400 Subject: [PATCH 6/6] Finalize signed FMM --- .../surface/algorithms/geodesic_distance.md | 2 +- src/surface/fast_marching_method.cpp | 68 +++++++++++-------- 2 files changed, 41 insertions(+), 29 deletions(-) diff --git a/docs/docs/surface/algorithms/geodesic_distance.md b/docs/docs/surface/algorithms/geodesic_distance.md index da894bcd..20ba664a 100644 --- a/docs/docs/surface/algorithms/geodesic_distance.md +++ b/docs/docs/surface/algorithms/geodesic_distance.md @@ -178,7 +178,7 @@ Computing exact geodesic paths also allows one to compute exact [log maps](/surf These routines implement the _fast marching method_ (FMM) for geodesic distance on triangle meshes, as described by [Kimmel and Sethian 1998](https://www.pnas.org/doi/pdf/10.1073/pnas.95.15.8431). -This implementation allows you to obtain unsigned distance to a set of points, or signed distance to a set of curves. +This implementation allows you to obtain unsigned distance to a source set of points, or signed distance to a source set of curves, where the distance values on the source set may be initialized to any value. (When initial distances are not 0, "signed" means that the gradient of distance is continuous across the source curves.) Note that as a wave-based propagation method, if the source curves are not closed, then signed distance output is not guaranteed to be well-behaved; if you're looking for robust signed distance computation, consider the [Signed Heat Method routines](/surface/algorithms/signed_heat_method). If your curves are closed, then FMM will in general be more efficient. diff --git a/src/surface/fast_marching_method.cpp b/src/surface/fast_marching_method.cpp index 69a7f42c..af735907 100644 --- a/src/surface/fast_marching_method.cpp +++ b/src/surface/fast_marching_method.cpp @@ -11,7 +11,7 @@ namespace { // The super fun quadratic distance function in the Fast Marching Method on triangle meshes // TODO parameter c isn't actually defined in paper, so I guessed that it was an error -double eikonalDistanceSubroutine(double a, double b, double theta, double dA, double dB) { +double eikonalDistanceSubroutine(double a, double b, double theta, double dA, double dB, int sign) { if (theta <= PI / 2.0) { @@ -24,24 +24,25 @@ double eikonalDistanceSubroutine(double a, double b, double theta, double dA, do double quadB = 2 * b * u * (a * cTheta - b); double quadC = b * b * (u * u - a * a * sTheta2); double sqrtVal = std::sqrt(quadB * quadB - 4 * quadA * quadC); - double t = (-quadB + sqrtVal) / (2 * quadA); // always the positive root (negative root corresponds to slope -1) - if (u < t && a * cTheta < b * (t - u) / t && b * (t - u) / t < a / cTheta) { + double tVals[] = {(-quadB + sqrtVal) / (2 * quadA), (-quadB - sqrtVal) / (2 * quadA)}; + double t = sign > 0 ? tVals[0] : tVals[1]; + double y = b * (t - u) / t; + if (u < sign * t && a * cTheta < y && y < a / cTheta) { return dA + t; } else { - return std::min(b + dA, a + dB); + return std::min(b * sign + dA, a * sign + dB); } } // Custom by Nick to get acceptable results in obtuse triangles without fancy unfolding else { - - double maxDist = std::max(dA, dB); // all points on base are less than this far away, by convexity + double maxDist = std::max(sign * dA, sign * dB); // all points on base are less than this far away, by convexity double c = std::sqrt(a * a + b * b - 2 * a * b * std::cos(theta)); double area = 0.5 * std::sin(theta) * a * b; double altitude = 2 * area / c; // distance to base, must be inside triangle since obtuse double baseDist = maxDist + altitude; - return std::min({b + dA, a + dB, baseDist}); + return std::min({b * sign + dA, a * sign + dB, sign * baseDist}); } } @@ -67,8 +68,17 @@ VertexData FMMDistance(IntrinsicGeometryInterface& geometry, VertexData distances(mesh, std::numeric_limits::infinity()); VertexData signs(mesh, 1); VertexData finalized(mesh, false); + VertexData isSource(mesh, false); - std::priority_queue, std::greater> frontierPQ; + auto cmp = [&signs](Entry left, Entry right) { + if (signs[left.second] != signs[right.second]) { + // We're looking at an edge that intersects the source, which may not have initial distance 0. + // I *think* this can only happen if the positive vertex lies on the source, in which case it's the "closer" one. + return signs[left.second] != 1; + } + return signs[left.second] * left.first > signs[right.second] * right.first; + }; + std::priority_queue, decltype(cmp)> frontierPQ(cmp); // Initialize signs if (sign) { for (Vertex v : mesh.vertices()) signs[v] = 0; @@ -126,24 +136,24 @@ VertexData FMMDistance(IntrinsicGeometryInterface& geometry, } } } - // Initialize distances. It's also possible to propagate signed distances and define eikonalDistanceSubroutine() - // accordingly, but it gets cumbersome managing directions of inequalities, and initializing distances of unvisited - // vertices as +/-inf. So just compute unsigned distances, then sign -- rather than try to actually solve the true - // BVP. + // Initialize distances. for (auto& curve : initialDistances) { for (auto& x : curve) { const SurfacePoint& p = x.first; switch (p.type) { case (SurfacePointType::Vertex): { frontierPQ.push(std::make_pair(x.second, p.vertex)); + isSource[p.vertex] = true; break; } case (SurfacePointType::Edge): { const Vertex& vA = p.edge.firstVertex(); const Vertex& vB = p.edge.secondVertex(); double l = geometry.edgeLengths[p.edge]; - frontierPQ.push(std::make_pair(x.second + p.tEdge * l, vA)); - frontierPQ.push(std::make_pair(x.second + (1. - p.tEdge) * l, vB)); + frontierPQ.push(std::make_pair(x.second + signs[vA] * p.tEdge * l, vA)); + frontierPQ.push(std::make_pair(x.second + signs[vB] * (1. - p.tEdge) * l, vB)); + isSource[vA] = true; + isSource[vB] = true; break; } case (SurfacePointType::Face): { @@ -163,9 +173,12 @@ VertexData FMMDistance(IntrinsicGeometryInterface& geometry, double dist2_A = lAB2 * (v * (1. - u)) + lCA2 * (w * (1. - u)) - lAB2 * v * w; // squared distance from p to vA double dist2_B = lAB2 * (u * (1. - v)) + lBC2 * (w * (1. - v)) - lCA2 * u * w; // squared distance from p to vB double dist2_C = lCA2 * (u * (1. - w)) + lBC2 * (v * (1. - w)) - lBC2 * u * v; // squared distance from p to vC - frontierPQ.push(std::make_pair(x.second + std::sqrt(dist2_A), vA)); - frontierPQ.push(std::make_pair(x.second + std::sqrt(dist2_B), vB)); - frontierPQ.push(std::make_pair(x.second + std::sqrt(dist2_C), vC)); + frontierPQ.push(std::make_pair(x.second + signs[vA] * std::sqrt(dist2_A), vA)); + frontierPQ.push(std::make_pair(x.second + signs[vB] * std::sqrt(dist2_B), vB)); + frontierPQ.push(std::make_pair(x.second + signs[vC] * std::sqrt(dist2_C), vC)); + isSource[vA] = true; + isSource[vB] = true; + isSource[vC] = true; break; } } @@ -197,10 +210,10 @@ VertexData FMMDistance(IntrinsicGeometryInterface& geometry, Vertex neighVert = he.vertex(); // Add with length - if (!finalized[neighVert]) { + if (!finalized[neighVert] && !isSource[neighVert]) { if (signs[neighVert] == 0) signs[neighVert] = signs[currV]; - double newDist = currDist + geometry.edgeLengths[he.edge()]; - if (newDist < distances[neighVert]) { + double newDist = currDist + signs[neighVert] * geometry.edgeLengths[he.edge()]; + if (signs[neighVert] * newDist < signs[neighVert] * distances[neighVert] || std::isinf(distances[neighVert])) { frontierPQ.push(std::make_pair(newDist, neighVert)); distances[neighVert] = newDist; } @@ -210,7 +223,7 @@ VertexData FMMDistance(IntrinsicGeometryInterface& geometry, // Check the third point of the "left" triangle straddling this edge if (he.isInterior()) { Vertex newVert = he.next().next().vertex(); - if (!finalized[newVert]) { + if (!finalized[newVert] && !isSource[newVert]) { // Compute the distance double lenB = geometry.edgeLengths[he.next().next().edge()]; @@ -219,8 +232,8 @@ VertexData FMMDistance(IntrinsicGeometryInterface& geometry, double distA = distances[neighVert]; double theta = geometry.cornerAngles[he.next().next().corner()]; if (signs[newVert] == 0) signs[newVert] = (signs[currV] != 0) ? signs[currV] : signs[he.next().vertex()]; - double newDist = eikonalDistanceSubroutine(lenA, lenB, theta, distA, distB); - if (newDist < distances[newVert]) { + double newDist = eikonalDistanceSubroutine(lenA, lenB, theta, distA, distB, signs[newVert]); + if (signs[newVert] * newDist < signs[newVert] * distances[newVert] || std::isinf(distances[newVert])) { frontierPQ.push(std::make_pair(newDist, newVert)); distances[newVert] = newDist; } @@ -231,7 +244,7 @@ VertexData FMMDistance(IntrinsicGeometryInterface& geometry, Halfedge heT = he.twin(); if (heT.isInterior()) { Vertex newVert = heT.next().next().vertex(); - if (!finalized[newVert]) { + if (!finalized[newVert] && !isSource[newVert]) { // Compute the distance double lenB = geometry.edgeLengths[heT.next().edge()]; @@ -240,8 +253,8 @@ VertexData FMMDistance(IntrinsicGeometryInterface& geometry, double distA = distances[neighVert]; double theta = geometry.cornerAngles[heT.next().next().corner()]; if (signs[newVert] == 0) signs[newVert] = (signs[currV] != 0) ? signs[currV] : signs[he.next().vertex()]; - double newDist = eikonalDistanceSubroutine(lenA, lenB, theta, distA, distB); - if (newDist < distances[newVert]) { + double newDist = eikonalDistanceSubroutine(lenA, lenB, theta, distA, distB, signs[newVert]); + if (signs[newVert] * newDist < signs[newVert] * distances[newVert] || std::isinf(distances[newVert])) { frontierPQ.push(std::make_pair(newDist, newVert)); distances[newVert] = newDist; } @@ -249,8 +262,7 @@ VertexData FMMDistance(IntrinsicGeometryInterface& geometry, } } } - // Sign unsigned distances. This will probably only do something reasonable if the curves are closed. - for (Vertex v : mesh.vertices()) distances[v] *= signs[v]; + return distances; }