From 84b6f154f66dae0189eb669f02842b5340c48120 Mon Sep 17 00:00:00 2001 From: nzfeng Date: Sun, 12 Oct 2025 17:47:34 -0400 Subject: [PATCH 1/2] First commit for preprocessing curves --- src/surface/signed_heat_method.cpp | 175 +++++++++++++++++++++++++++++ 1 file changed, 175 insertions(+) diff --git a/src/surface/signed_heat_method.cpp b/src/surface/signed_heat_method.cpp index 4673618a..5a23a7a4 100644 --- a/src/surface/signed_heat_method.cpp +++ b/src/surface/signed_heat_method.cpp @@ -32,6 +32,7 @@ VertexData SignedHeatSolver::computeDistance(const std::vector& c const std::vector& points, const SignedHeatOptions& options) { + std::vector processedCurves = preprocessCurves(curves); Vector> Xt = integrateVectorHeatFlow(curves, points, options); return integrateVectorField(Xt, curves, points, options); } @@ -56,6 +57,180 @@ void SignedHeatSolver::setDiffusionTimeCoefficient(double tCoef_) { doubleVectorOp = doubleMassMat + shortTime * doubleConnectionLaplacian; } +std::vector SignedHeatSolver::preprocessCurves(const std::vector& curves) { + + // If two SurfacePoints are in a common or adjacent faces, return the faces. + auto adjacentFaces = [](const SurfacePoint& pA, const SurfacePoint& pB, Edge& commonEdge) -> std::pair { + switch (pA.type) { + case SurfacePointType::Vertex: + switch (pB.type) { + case SurfacePointType::Vertex: + for (Face f : pA.vertex.adjacentFaces()) { + for (Edge e : f.adjacentEdges()) { + for (Face g : e.adjacentFaces()) { + for (Vertex v : g.adjacentVertices()) { + if (v == pB.vertex) { + commonEdge = e; + return std::make_pair(f, g); + } + } + } + } + } + break; + case SurfacePointType::Edge: + for (Face f : pA.vertex.adjacentFaces()) { + for (Edge e : f.adjacentEdges()) { + for (Face g : e.adjacentFaces()) { + for (Edge h : g.adjacentEdges()) { + if (h == pB.edge) { + commonEdge = e; + return std::make_pair(f, g); + } + } + } + } + } + break; + case SurfacePointType::Face: + for (Face f : pA.vertex.adjacentFaces()) { + for (Edge e : f.adjacentEdges()) { + for (Face g : e.adjacentFaces()) { + if (g == pB.face) { + commonEdge = e; + return std::make_pair(f, g); + } + } + } + } + break; + } + break; + case SurfacePointType::Edge: + switch (pB.type) { + case SurfacePointType::Vertex: + for (Face f : pA.edge.adjacentFaces()) { + for (Edge e : f.adjacentEdges()) { + for (Face g : e.adjacentFaces()) { + for (Vertex v : g.adjacentVertices()) { + if (v == pB.vertex) { + commonEdge = e; + return std::make_pair(f, g); + } + } + } + } + } + break; + case SurfacePointType::Edge: + for (Face f : pA.edge.adjacentFaces()) { + for (Edge e : f.adjacentEdges()) { + for (Face g : e.adjacentFaces()) { + for (Edge h : g.adjacentEdges()) { + if (h == pB.edge) { + commonEdge = e; + return std::make_pair(f, g); + } + } + } + } + } + break; + case SurfacePointType::Face: + for (Face f : pA.edge.adjacentFaces()) { + for (Edge e : f.adjacentEdges()) { + for (Face g : e.adjacentFaces()) { + if (g == pB.face) { + commonEdge = e; + return std::make_pair(f, g); + } + } + } + } + break; + } + break; + case SurfacePointType::Face: + switch (pB.type) { + case SurfacePointType::Vertex: + for (Edge e : pA.face.adjacentFaces()) { + for (Face g : e.adjacentFaces()) { + for (Vertex v : g.adjacentVertices()) { + if (v == pB.vertex) { + commonEdge = e; + return std::make_pair(pA.face, g); + } + } + } + } + break; + case SurfacePointType::Edge: + for (Edge e : pA.face.adjacentFaces()) { + for (Face g : e.adjacentFaces()) { + for (Edge h : g.adjacentEdges()) { + if (h == pB.edge) { + commonEdge = e; + return std::make_pair(pA.face, g); + } + } + } + } + break; + case SurfacePointType::Face: + for (Edge e : pA.face.adjacentFaces()) { + for (Face g : e.adjacentFaces()) { + if (g == pB.face) { + commonEdge = e; + return std::make_pair(pA.face, g); + } + } + } + break; + } + break; + } + commonEdge = Edge(); + return std::make_pair(Face(), Face()); + }; + + // If there are gaps in the curves (i.e. curves are not sampled densely along the mesh), then there is an ambiguity of + // how to connect up the points into curves. Rather than, for example, automatically connecting up curves using a + // geodesic path -- which would impose curves that the user might not have meant -- we simply break up the input + // curves into components that do reflect how they're sampled. + std::vector newCurves; + for (const Curve& curve : curves) { + newCurves.emplace_back(); + newCurves.back().isSigned = curve.isSigned; + size_t nNodes = curve.nodes.size(); + for (size_t i = 0; i < nNodes - 1; i++) { + const SurfacePoint& pA = curve.nodes[i]; + const SurfacePoint& pB = curve.nodes[i + 1]; + newCurves.back().nodes.push_back(pA); + Face commonFace = sharedFace(pA, pB); + if (commonFace == Face()) { + if (curve.isSigned) { + // If the curve betwen pA and pB is unclear, but pA and pB are still in adjacent faces, then we can at + // least figure out the edge crossing. But don't do this for optimization for unsigned curves, which are + // still supposed to be restricted to mesh edges. + Edge commonEdge; + std::pair facePair = adjacentFaces(pA, pB, commonEdge); + if (facePair.first != Face() && facePair.second != Face()) { + // We want to solve for the edge crossing p s.t. = -. + BarycentricVector eVec(commonEdge, Vector2(1, -1)); // pick arbitrary direction for edge + } + } else { + // Create new curve + newCurves.emplace_back(); + newCurves.back().isSigned = curve.isSigned; + } + } + } + // Don't forget the last point + newCurves.back().nodes.push_back(curve.nodes[nNodes - 1]); + } + return newCurves; +} + Vector> SignedHeatSolver::integrateVectorHeatFlow(const std::vector& curves, const std::vector& points, const SignedHeatOptions& options) { From e2a737e339de7f179c69593b4de19250674e2a5b Mon Sep 17 00:00:00 2001 From: nzfeng Date: Mon, 13 Oct 2025 12:07:22 -0400 Subject: [PATCH 2/2] Preprocess curves by partitioning them into well-sampled components --- .../surface/algorithms/signed_heat_method.md | 2 + .../surface/signed_heat_method.h | 1 + src/surface/signed_heat_method.cpp | 166 ++---------------- 3 files changed, 13 insertions(+), 156 deletions(-) diff --git a/docs/docs/surface/algorithms/signed_heat_method.md b/docs/docs/surface/algorithms/signed_heat_method.md index afb8a366..c645beb8 100644 --- a/docs/docs/surface/algorithms/signed_heat_method.md +++ b/docs/docs/surface/algorithms/signed_heat_method.md @@ -64,6 +64,8 @@ The diffusion time can also be changed using the following function. Re-computes the time used for short time heat flow `tCoef * h^2`, where `h` is the mean distance between nodes of the mesh. +Note: On triangle meshes, each mesh edge crossed by an input curve should be explicitly represented by a SurfacePoint (otherwise the curve's geometry is not fully specified). If this is not the case, the code internally partitions input curves so that this condition is satisfied; the robustness of the heat method should fill in small gaps. Depending on how sparsely your curve is sampled, this may produce curve components with fewer than two nodes, which will be ignored; if many of your curves have fewer than two nodes, consider [FlipOut geodesics](../flip_geodesics) for shortest-path completion of curves. + ## Helper Types ### Curves diff --git a/include/geometrycentral/surface/signed_heat_method.h b/include/geometrycentral/surface/signed_heat_method.h index d29734bc..a6b44570 100644 --- a/include/geometrycentral/surface/signed_heat_method.h +++ b/include/geometrycentral/surface/signed_heat_method.h @@ -61,6 +61,7 @@ class SignedHeatSolver { SparseMatrix massMat, doubleMassMat, doubleConnectionLaplacian, doubleVectorOp; // Helpers + std::vector preprocessCurves(const std::vector& curves) const; Vector> integrateVectorHeatFlow(const std::vector& curves, const std::vector& points, const SignedHeatOptions& options); diff --git a/src/surface/signed_heat_method.cpp b/src/surface/signed_heat_method.cpp index 5a23a7a4..33cb5a71 100644 --- a/src/surface/signed_heat_method.cpp +++ b/src/surface/signed_heat_method.cpp @@ -33,8 +33,8 @@ VertexData SignedHeatSolver::computeDistance(const std::vector& c const SignedHeatOptions& options) { std::vector processedCurves = preprocessCurves(curves); - Vector> Xt = integrateVectorHeatFlow(curves, points, options); - return integrateVectorField(Xt, curves, points, options); + Vector> Xt = integrateVectorHeatFlow(processedCurves, points, options); + return integrateVectorField(Xt, processedCurves, points, options); } VertexData SignedHeatSolver::computeDistance(const std::vector& curves, @@ -57,146 +57,12 @@ void SignedHeatSolver::setDiffusionTimeCoefficient(double tCoef_) { doubleVectorOp = doubleMassMat + shortTime * doubleConnectionLaplacian; } -std::vector SignedHeatSolver::preprocessCurves(const std::vector& curves) { - - // If two SurfacePoints are in a common or adjacent faces, return the faces. - auto adjacentFaces = [](const SurfacePoint& pA, const SurfacePoint& pB, Edge& commonEdge) -> std::pair { - switch (pA.type) { - case SurfacePointType::Vertex: - switch (pB.type) { - case SurfacePointType::Vertex: - for (Face f : pA.vertex.adjacentFaces()) { - for (Edge e : f.adjacentEdges()) { - for (Face g : e.adjacentFaces()) { - for (Vertex v : g.adjacentVertices()) { - if (v == pB.vertex) { - commonEdge = e; - return std::make_pair(f, g); - } - } - } - } - } - break; - case SurfacePointType::Edge: - for (Face f : pA.vertex.adjacentFaces()) { - for (Edge e : f.adjacentEdges()) { - for (Face g : e.adjacentFaces()) { - for (Edge h : g.adjacentEdges()) { - if (h == pB.edge) { - commonEdge = e; - return std::make_pair(f, g); - } - } - } - } - } - break; - case SurfacePointType::Face: - for (Face f : pA.vertex.adjacentFaces()) { - for (Edge e : f.adjacentEdges()) { - for (Face g : e.adjacentFaces()) { - if (g == pB.face) { - commonEdge = e; - return std::make_pair(f, g); - } - } - } - } - break; - } - break; - case SurfacePointType::Edge: - switch (pB.type) { - case SurfacePointType::Vertex: - for (Face f : pA.edge.adjacentFaces()) { - for (Edge e : f.adjacentEdges()) { - for (Face g : e.adjacentFaces()) { - for (Vertex v : g.adjacentVertices()) { - if (v == pB.vertex) { - commonEdge = e; - return std::make_pair(f, g); - } - } - } - } - } - break; - case SurfacePointType::Edge: - for (Face f : pA.edge.adjacentFaces()) { - for (Edge e : f.adjacentEdges()) { - for (Face g : e.adjacentFaces()) { - for (Edge h : g.adjacentEdges()) { - if (h == pB.edge) { - commonEdge = e; - return std::make_pair(f, g); - } - } - } - } - } - break; - case SurfacePointType::Face: - for (Face f : pA.edge.adjacentFaces()) { - for (Edge e : f.adjacentEdges()) { - for (Face g : e.adjacentFaces()) { - if (g == pB.face) { - commonEdge = e; - return std::make_pair(f, g); - } - } - } - } - break; - } - break; - case SurfacePointType::Face: - switch (pB.type) { - case SurfacePointType::Vertex: - for (Edge e : pA.face.adjacentFaces()) { - for (Face g : e.adjacentFaces()) { - for (Vertex v : g.adjacentVertices()) { - if (v == pB.vertex) { - commonEdge = e; - return std::make_pair(pA.face, g); - } - } - } - } - break; - case SurfacePointType::Edge: - for (Edge e : pA.face.adjacentFaces()) { - for (Face g : e.adjacentFaces()) { - for (Edge h : g.adjacentEdges()) { - if (h == pB.edge) { - commonEdge = e; - return std::make_pair(pA.face, g); - } - } - } - } - break; - case SurfacePointType::Face: - for (Edge e : pA.face.adjacentFaces()) { - for (Face g : e.adjacentFaces()) { - if (g == pB.face) { - commonEdge = e; - return std::make_pair(pA.face, g); - } - } - } - break; - } - break; - } - commonEdge = Edge(); - return std::make_pair(Face(), Face()); - }; - +std::vector SignedHeatSolver::preprocessCurves(const std::vector& curves) const { // If there are gaps in the curves (i.e. curves are not sampled densely along the mesh), then there is an ambiguity of - // how to connect up the points into curves. Rather than, for example, automatically connecting up curves using a - // geodesic path -- which would impose curves that the user might not have meant -- we simply break up the input - // curves into components that do reflect how they're sampled. + // how to connect up the points into curves. + // Rather than, for example, automatically connecting up curves using a geodesic path -- which would impose curves + // that the user might not have meant -- we simply break up the input curves into components that do reflect how + // they're sampled. (This may, however, create curves with fewer than 2 nodes, which would get ignored.) std::vector newCurves; for (const Curve& curve : curves) { newCurves.emplace_back(); @@ -208,21 +74,9 @@ std::vector SignedHeatSolver::preprocessCurves(const std::vector& newCurves.back().nodes.push_back(pA); Face commonFace = sharedFace(pA, pB); if (commonFace == Face()) { - if (curve.isSigned) { - // If the curve betwen pA and pB is unclear, but pA and pB are still in adjacent faces, then we can at - // least figure out the edge crossing. But don't do this for optimization for unsigned curves, which are - // still supposed to be restricted to mesh edges. - Edge commonEdge; - std::pair facePair = adjacentFaces(pA, pB, commonEdge); - if (facePair.first != Face() && facePair.second != Face()) { - // We want to solve for the edge crossing p s.t. = -. - BarycentricVector eVec(commonEdge, Vector2(1, -1)); // pick arbitrary direction for edge - } - } else { - // Create new curve - newCurves.emplace_back(); - newCurves.back().isSigned = curve.isSigned; - } + // Create new curve + newCurves.emplace_back(); + newCurves.back().isSigned = curve.isSigned; } } // Don't forget the last point