From 22e2e866a2da70c91ef71584193ad94424403ad0 Mon Sep 17 00:00:00 2001 From: nzfeng Date: Tue, 17 Jun 2025 21:30:22 -0400 Subject: [PATCH 1/2] First commit of marching triangles --- .../surface/utilities/barycentric_vector.md | 6 +- .../surface/barycentric_vector.h | 1 + .../surface/barycentric_vector.ipp | 17 +++ .../surface/marching_triangles.h | 12 ++ .../surface/signed_heat_method.h | 1 - src/CMakeLists.txt | 1 + src/surface/fast_marching_method.cpp | 6 +- src/surface/marching_triangles.cpp | 129 ++++++++++++++++++ src/surface/signed_heat_method.cpp | 29 +--- 9 files changed, 174 insertions(+), 28 deletions(-) create mode 100644 include/geometrycentral/surface/marching_triangles.h create mode 100644 src/surface/marching_triangles.cpp diff --git a/docs/docs/surface/utilities/barycentric_vector.md b/docs/docs/surface/utilities/barycentric_vector.md index b83814ea..4db72a3d 100644 --- a/docs/docs/surface/utilities/barycentric_vector.md +++ b/docs/docs/surface/utilities/barycentric_vector.md @@ -25,7 +25,7 @@ which indicates what kind of vector it is. ??? func "`#!cpp BarycentricVector::BarycentricVector(Face f, Vector3 faceCoords)`" - Construct a barycentric vector that lies along the given face `f`, with the given barycentric coordinates `faceCoords`. + Construct a barycentric vector that lies in the given face `f`, with the given barycentric coordinates `faceCoords`. ??? func "`#!cpp BarycentricVector::BarycentricVector(Edge e, Vector2 edgeCoords)`" @@ -35,6 +35,10 @@ which indicates what kind of vector it is. Construct a (zero-length) barycentric vector that lies on the given vertex `v`. +??? func "`#!cpp BarycentricVector::BarycentricVector(Halfedge he, Face f)`" + + Construct a barycentric vector from the given halfedge, that lies in the given face `f`. + ??? func "`#!cpp BarycentricVector::BarycentricVector(SurfacePoint pA, SurfacePoint pB)`" Construct a barycentric vector from the given [surface points](../surface_point), i.e. barycentric points. The vector direction goes from `pA` to `pB`. diff --git a/include/geometrycentral/surface/barycentric_vector.h b/include/geometrycentral/surface/barycentric_vector.h index f55db8de..31781e19 100644 --- a/include/geometrycentral/surface/barycentric_vector.h +++ b/include/geometrycentral/surface/barycentric_vector.h @@ -20,6 +20,7 @@ struct BarycentricVector { BarycentricVector(Vertex v); BarycentricVector(Edge e, Vector2 edgeCoords); BarycentricVector(Face f, Vector3 faceCoords); + BarycentricVector(Halfedge he, Face f); BarycentricVectorType type = BarycentricVectorType::Face; diff --git a/include/geometrycentral/surface/barycentric_vector.ipp b/include/geometrycentral/surface/barycentric_vector.ipp index 637e3892..04f74d85 100644 --- a/include/geometrycentral/surface/barycentric_vector.ipp +++ b/include/geometrycentral/surface/barycentric_vector.ipp @@ -44,6 +44,23 @@ inline BarycentricVector::BarycentricVector(SurfacePoint pA, SurfacePoint pB) { faceCoords = pB_in_face.faceCoords - pA_in_face.faceCoords; } +inline BarycentricVector::BarycentricVector(Halfedge he_, Face f) : type(BarycentricVectorType::Face), face(f) { + + int eIdx = 0; + double sign = 0.; + for (Halfedge he : f.adjacentHalfedges()) { + if (he.edge() == he_.edge()) { + sign = (he.tailVertex() == he_.tailVertex() && he.tipVertex() == he_.tipVertex()) ? 1. : -1.; + break; + } + eIdx++; + } + faceCoords = {0, 0, 0}; + faceCoords[(eIdx + 1) % 3] = 1; + faceCoords[eIdx] = -1; + faceCoords *= sign; +} + // == Methods inline BarycentricVector BarycentricVector::inSomeFace() const { diff --git a/include/geometrycentral/surface/marching_triangles.h b/include/geometrycentral/surface/marching_triangles.h new file mode 100644 index 00000000..43075ddb --- /dev/null +++ b/include/geometrycentral/surface/marching_triangles.h @@ -0,0 +1,12 @@ +#pragma once + +#include "geometrycentral/surface/barycentric_vector.h" + +namespace geometrycentral { +namespace surface { + +std::vector> marchingTriangles(IntrinsicGeometryInterface& geom, const VertexData& u, + double isoval = 0.); + +} +} // namespace geometrycentral \ No newline at end of file diff --git a/include/geometrycentral/surface/signed_heat_method.h b/include/geometrycentral/surface/signed_heat_method.h index b4ede4d0..d29734bc 100644 --- a/include/geometrycentral/surface/signed_heat_method.h +++ b/include/geometrycentral/surface/signed_heat_method.h @@ -80,7 +80,6 @@ class SignedHeatSolver { double lengthOfSegment(const SurfacePoint& pA, const SurfacePoint& pB) const; SurfacePoint midSegmentSurfacePoint(const SurfacePoint& pA, const SurfacePoint& pB) const; std::complex projectedNormal(const SurfacePoint& pA, const SurfacePoint& pB, const Edge& e) const; - BarycentricVector barycentricVectorInFace(const Halfedge& he, const Face& f) const; FaceData sampleAtFaceBarycenters(const Vector>& Xt); Vector integrateWithZeroSetConstraint(const Vector& rhs, const std::vector& curves, const std::vector& points, diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 63a5be71..b74e243d 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -26,6 +26,7 @@ SET(SRCS surface/embedded_geometry_interface.cpp surface/edge_length_geometry.cpp surface/vertex_position_geometry.cpp + surface/marching_triangles.cpp surface/mutation_manager.cpp surface/mesh_graph_algorithms.cpp surface/direction_fields.cpp diff --git a/src/surface/fast_marching_method.cpp b/src/surface/fast_marching_method.cpp index af735907..6c8de707 100644 --- a/src/surface/fast_marching_method.cpp +++ b/src/surface/fast_marching_method.cpp @@ -92,7 +92,7 @@ VertexData FMMDistance(IntrinsicGeometryInterface& geometry, // 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.next().tipVertex()] = (he.vertex() == pA.vertex) ? -1 : 1; signs[he.twin().next().tipVertex()] = -signs[he.next().tipVertex()]; } else { Face commonFace = sharedFace(pA, pB); @@ -100,9 +100,9 @@ VertexData FMMDistance(IntrinsicGeometryInterface& geometry, 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); + BarycentricVector normal = -tangent.rotate90(geometry); for (Vertex v : commonFace.adjacentVertices()) { - BarycentricVector u(SurfacePoint(v), pA); + BarycentricVector u(pA, SurfacePoint(v)); signs[v] = (dot(geometry, normal, u) > 0) ? 1 : -1; } } diff --git a/src/surface/marching_triangles.cpp b/src/surface/marching_triangles.cpp new file mode 100644 index 00000000..5436d8fb --- /dev/null +++ b/src/surface/marching_triangles.cpp @@ -0,0 +1,129 @@ +#include "geometrycentral/surface/marching_triangles.h" + +namespace geometrycentral { +namespace surface { + +namespace { + +std::vector>> +getCurveComponents(SurfaceMesh& mesh, const std::vector& curveNodes, + const std::vector>& curveEdges) { + + // Note: This function assumes that `curveNodes` has been de-duplicated. + std::vector> edgesToAdd = curveEdges; + std::vector>> curves; + size_t nSegs = curveEdges.size(); + while (edgesToAdd.size() > 0) { + std::array startSeg = edgesToAdd.back(); + edgesToAdd.pop_back(); + curves.emplace_back(); + std::vector>& currCurve = curves.back(); + currCurve.push_back(startSeg); + + // Add segs to the front end until we can't. + std::array currSeg = startSeg; + while (true) { + const SurfacePoint& front = curveNodes[currSeg[1]]; + bool didWeFindOne = false; + for (size_t i = 0; i < edgesToAdd.size(); i++) { + std::array otherSeg = edgesToAdd[i]; + if (curveNodes[otherSeg[0]] == front) { + currSeg = otherSeg; + currCurve.push_back(otherSeg); + edgesToAdd.erase(edgesToAdd.begin() + i); + didWeFindOne = true; + break; + } + } + if (!didWeFindOne) break; + } + + // Add segs to the back end until we can't. + currSeg = startSeg; + while (true) { + const SurfacePoint& back = curveNodes[currSeg[0]]; + bool didWeFindOne = false; + for (size_t i = 0; i < edgesToAdd.size(); i++) { + std::array otherSeg = edgesToAdd[i]; + if (curveNodes[otherSeg[1]] == back) { + currSeg = otherSeg; + currCurve.insert(currCurve.begin(), otherSeg); + edgesToAdd.erase(edgesToAdd.begin() + i); + didWeFindOne = true; + break; + } + } + if (!didWeFindOne) break; + } + } + return curves; +} +} // namespace + +std::vector> marchingTriangles(IntrinsicGeometryInterface& geom, const VertexData& u, + double isoval) { + + SurfaceMesh& mesh = *u.getMesh(); + std::vector nodes; + std::vector> edges; + for (Face f : mesh.faces()) { + std::vector hits; + BarycentricVector gradient(f); + for (Halfedge he : f.adjacentHalfedges()) { + // Record edge crossings + Edge e = he.edge(); + Vertex v0 = e.firstVertex(); + Vertex v1 = e.secondVertex(); + double u0 = u[v0]; + double u1 = u[v1]; + double lB = std::min(u0, u1); + double uB = std::max(u0, u1); + if (lB == uB && lB == isoval) { + hits.clear(); + hits = {SurfacePoint(v0), SurfacePoint(v1)}; + break; + } + double t = (isoval - lB) / (uB - lB); + if (u0 > u1) t = 1. - t; + if (t <= 1. && t >= 0.) hits.emplace_back(e, t); + // Compute gradient of the scalar function. + BarycentricVector heVec(he.next(), f); + BarycentricVector ePerp = heVec.rotate90(geom); + gradient += ePerp * u[he.vertex()]; + } + if (hits.size() != 2) continue; + + // Orient segments so that smaller values are always on the "inside" of the curve. + std::array seg; + for (int i = 0; i < 2; i++) { + auto iter = std::find(nodes.begin(), nodes.end(), hits[i]); + if (iter != nodes.end()) { + seg[i] = iter - nodes.begin(); + } else { + nodes.push_back(hits[i]); + seg[i] = nodes.size() - 1; + } + } + BarycentricVector segTangent(hits[0], hits[1]); + BarycentricVector segNormal = -segTangent.inFace(f).rotate90(geom); + if (dot(geom, segNormal, gradient) > 0.) { + edges.push_back(seg); + } else { + edges.push_back({seg[1], seg[0]}); + } + } + std::vector>> components = getCurveComponents(mesh, nodes, edges); + size_t nCurves = components.size(); + std::vector> curves(nCurves); + for (size_t i = 0; i < nCurves; i++) { + size_t nEdges = components[i].size(); + for (size_t j = 0; j < nEdges; j++) { + curves[i].push_back(nodes[components[i][j][0]]); + } + curves[i].push_back(nodes[components[i][nEdges - 1][1]]); + } + return curves; +} + +} // namespace surface +} // namespace geometrycentral \ No newline at end of file diff --git a/src/surface/signed_heat_method.cpp b/src/surface/signed_heat_method.cpp index 6e005a41..4673618a 100644 --- a/src/surface/signed_heat_method.cpp +++ b/src/surface/signed_heat_method.cpp @@ -113,7 +113,7 @@ Vector> SignedHeatSolver::integrateVectorHeatFlow(const std for (Edge e : f.adjacentEdges()) { size_t eIdx = geom.edgeIndices[e]; double w = scalarCrouzeixRaviart(b, e); - BarycentricVector heVec = barycentricVectorInFace(e.halfedge(), f); + BarycentricVector heVec(e.halfedge(), f); heVec /= heVec.norm(geom); BarycentricVector heVecN = heVec.rotate90(geom); triplets.emplace_back(m, eIdx, w * dot(geom, heVec, segTangent)); @@ -162,8 +162,8 @@ VertexData SignedHeatSolver::integrateVectorField(const Vector SignedHeatSolver::projectedNormal(const SurfacePoint& pA, c return normal; } -BarycentricVector SignedHeatSolver::barycentricVectorInFace(const Halfedge& he_, const Face& f) const { - - int eIdx = 0; - double sign = 0.; - for (Halfedge he : f.adjacentHalfedges()) { - if (he.edge() == he_.edge()) { - sign = (he.tailVertex() == he_.tailVertex() && he.tipVertex() == he_.tipVertex()) ? 1. : -1.; - break; - } - eIdx++; - } - Vector3 faceCoords = {0, 0, 0}; - faceCoords[(eIdx + 1) % 3] = 1; - faceCoords[eIdx] = -1; - return BarycentricVector(f, sign * faceCoords); -} - FaceData SignedHeatSolver::sampleAtFaceBarycenters(const Vector>& Xt) { geom.requireEdgeIndices(); @@ -724,7 +707,7 @@ FaceData SignedHeatSolver::sampleAtFaceBarycenters(const Vect Vector3 faceCoords = {0, 0, 0}; for (Halfedge he : f.adjacentHalfedges()) { size_t eIdx = geom.edgeIndices[he.edge()]; - BarycentricVector e1 = barycentricVectorInFace(he, f); + BarycentricVector e1(he, f); if (!he.orientation()) e1 *= -1; BarycentricVector e2 = e1.rotate90(geom); e1 /= e1.norm(geom); From 1c25bb08e1fed11e08bc74240a0d6853333d9eef Mon Sep 17 00:00:00 2001 From: nzfeng Date: Tue, 17 Jun 2025 21:46:52 -0400 Subject: [PATCH 2/2] Fix bug in SurfacePoint's inEdge(). Yikes --- include/geometrycentral/surface/surface_point.ipp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/include/geometrycentral/surface/surface_point.ipp b/include/geometrycentral/surface/surface_point.ipp index 63a55c9c..3c4f0b0e 100644 --- a/include/geometrycentral/surface/surface_point.ipp +++ b/include/geometrycentral/surface/surface_point.ipp @@ -151,9 +151,9 @@ inline SurfacePoint SurfacePoint::inEdge(Edge targetEdge) const { switch (type) { case SurfacePointType::Vertex: if (vertex == targetEdge.halfedge().tailVertex()) { - return SurfacePoint(targetEdge, 0); - } else if (vertex == targetEdge.halfedge().tipVertex()) { return SurfacePoint(targetEdge, 1); + } else if (vertex == targetEdge.halfedge().tipVertex()) { + return SurfacePoint(targetEdge, 0); } break; case SurfacePointType::Edge: