From 1df721f579615733556fb1f22ec77e3b1bd65744 Mon Sep 17 00:00:00 2001 From: Liam White Date: Wed, 10 Jun 2026 16:04:02 -0400 Subject: [PATCH 1/2] fix(surface): handle exact isovalue vertices in marching triangles Classify scalar samples as below, equal, or above the isovalue before collecting contour hits. Deduplicate hits so an exact isovalue vertex shared by two incident edges does not produce a three-hit face and get skipped. Also skip ambiguous degenerate cases explicitly when a face has anything other than two unique contour hits. --- src/surface/marching_triangles.cpp | 87 +++++++++++++++++++++++------- 1 file changed, 69 insertions(+), 18 deletions(-) diff --git a/src/surface/marching_triangles.cpp b/src/surface/marching_triangles.cpp index 5436d8fb..826778ee 100644 --- a/src/surface/marching_triangles.cpp +++ b/src/surface/marching_triangles.cpp @@ -1,10 +1,68 @@ #include "geometrycentral/surface/marching_triangles.h" +#include +#include + namespace geometrycentral { namespace surface { namespace { +enum class IsoSign { Below, Equal, Above, Invalid }; + +IsoSign isoSign(double x, double iso) { + if (!std::isfinite(x) || !std::isfinite(iso)) return IsoSign::Invalid; + if (x < iso) return IsoSign::Below; + if (x > iso) return IsoSign::Above; + return IsoSign::Equal; +} + +void addUniqueHit(std::vector& hits, const SurfacePoint& p) { + if (std::find(hits.begin(), hits.end(), p) == hits.end()) { + hits.push_back(p); + } +} + +void addEdgeIsoHits(Halfedge he, const VertexData& u, double isoval, std::vector& hits) { + + Edge e = he.edge(); + Vertex v0 = e.firstVertex(); + Vertex v1 = e.secondVertex(); + + double u0 = u[v0]; + double u1 = u[v1]; + + IsoSign s0 = isoSign(u0, isoval); + IsoSign s1 = isoSign(u1, isoval); + + if (s0 == IsoSign::Invalid || s1 == IsoSign::Invalid) return; + + // Entire edge is on the contour. + if (s0 == IsoSign::Equal && s1 == IsoSign::Equal) { + addUniqueHit(hits, SurfacePoint(v0)); + addUniqueHit(hits, SurfacePoint(v1)); + return; + } + + // One endpoint lies exactly on the contour. + if (s0 == IsoSign::Equal) { + addUniqueHit(hits, SurfacePoint(v0)); + return; + } + + if (s1 == IsoSign::Equal) { + addUniqueHit(hits, SurfacePoint(v1)); + return; + } + + // Strict crossing. + if (s0 != s1) { + double t = (isoval - u0) / (u1 - u0); + t = std::max(0.0, std::min(1.0, t)); + addUniqueHit(hits, SurfacePoint(e, t)); + } +} + std::vector>> getCurveComponents(SurfaceMesh& mesh, const std::vector& curveNodes, const std::vector>& curveEdges) { @@ -70,28 +128,21 @@ std::vector> marchingTriangles(IntrinsicGeometryInterf 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); + + // Record edge crossings robustly. + addEdgeIsoHits(he, u, isoval, hits); + // 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; + if (hits.size() != 2) { + // hits.size() == 0: no contour. + // hits.size() == 1: contour only touches this face at one vertex. + // hits.size() == 3: whole face is iso-valued, so the 1D contour is ambiguous. + continue; + } // Orient segments so that smaller values are always on the "inside" of the curve. std::array seg; @@ -126,4 +177,4 @@ std::vector> marchingTriangles(IntrinsicGeometryInterf } } // namespace surface -} // namespace geometrycentral \ No newline at end of file +} // namespace geometrycentral From 214f588d7149501850f23b18cd5b86747a58440d Mon Sep 17 00:00:00 2001 From: Liam White Date: Wed, 10 Jun 2026 16:04:25 -0400 Subject: [PATCH 2/2] test(surface): cover marching triangles isovalue degeneracies Add regression coverage for a contour passing through an exact isovalue vertex, ensuring the face still emits a segment. Also cover fully iso-valued faces, which should be skipped because the one-dimensional contour is ambiguous. --- test/CMakeLists.txt | 1 + test/src/marching_triangles_test.cpp | 96 ++++++++++++++++++++++++++++ 2 files changed, 97 insertions(+) create mode 100644 test/src/marching_triangles_test.cpp diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index a55687fe..2e7aeae3 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -104,6 +104,7 @@ set(TEST_SRCS src/intrinsic_triangulation_test.cpp src/linear_algebra_test.cpp src/load_test_meshes.cpp + src/marching_triangles_test.cpp src/point_cloud_test.cpp src/poisson_disk_sampler_test.cpp src/polygon_operators_test.cpp diff --git a/test/src/marching_triangles_test.cpp b/test/src/marching_triangles_test.cpp new file mode 100644 index 00000000..fadcb9a9 --- /dev/null +++ b/test/src/marching_triangles_test.cpp @@ -0,0 +1,96 @@ +#include "geometrycentral/surface/marching_triangles.h" +#include "geometrycentral/surface/surface_mesh_factories.h" + +#include "gtest/gtest.h" + +#include +#include +#include + +using namespace geometrycentral; +using namespace geometrycentral::surface; + +namespace { + +std::tuple, std::unique_ptr> buildSingleTriangleMesh() { + std::vector verts = {{0, 0, 0}, {1, 0, 0}, {0, 1, 0}}; + std::vector> faces = {{0, 1, 2}}; + return makeManifoldSurfaceMeshAndGeometry(faces, verts); +} + +std::tuple, std::unique_ptr> buildTwoTriangleMesh() { + std::vector verts = {{0, 0, 0}, {1, 0, 0}, {1, 1, 0}, {0, 1, 0}}; + std::vector> faces = {{0, 1, 2}, {0, 2, 3}}; + return makeManifoldSurfaceMeshAndGeometry(faces, verts); +} + +} // namespace + +TEST(MarchingTrianglesSuite, ExactIsoVertexCreatesSegment) { + std::unique_ptr meshPtr; + std::unique_ptr geomPtr; + std::tie(meshPtr, geomPtr) = buildSingleTriangleMesh(); + + ManifoldSurfaceMesh& mesh = *meshPtr; + VertexData values(mesh); + values[mesh.vertex(0)] = 0.; + values[mesh.vertex(1)] = 1.; + values[mesh.vertex(2)] = -1.; + + std::vector> curves = marchingTriangles(*geomPtr, values, 0.); + + ASSERT_EQ(curves.size(), 1); + ASSERT_EQ(curves[0].size(), 2); + + bool hasIsoVertex = false; + bool hasCrossing = false; + for (const SurfacePoint& p : curves[0]) { + if (p.type == SurfacePointType::Vertex && p.vertex == mesh.vertex(0)) { + hasIsoVertex = true; + } + if (p.type == SurfacePointType::Edge) { + Vertex vA = p.edge.firstVertex(); + Vertex vB = p.edge.secondVertex(); + bool onOppositeEdge = (vA == mesh.vertex(1) && vB == mesh.vertex(2)) || + (vA == mesh.vertex(2) && vB == mesh.vertex(1)); + if (onOppositeEdge) { + hasCrossing = true; + EXPECT_NEAR(p.tEdge, 0.5, 1e-12); + } + } + } + + EXPECT_TRUE(hasIsoVertex); + EXPECT_TRUE(hasCrossing); +} + +TEST(MarchingTrianglesSuite, SharedEdgeCrossingsStitchIntoOneCurve) { + std::unique_ptr meshPtr; + std::unique_ptr geomPtr; + std::tie(meshPtr, geomPtr) = buildTwoTriangleMesh(); + + ManifoldSurfaceMesh& mesh = *meshPtr; + VertexData values(mesh); + values[mesh.vertex(0)] = -0.8929741240814866; + values[mesh.vertex(1)] = -0.8929741240814866; + values[mesh.vertex(2)] = 0.1403065708384703; + values[mesh.vertex(3)] = 0.1403065708384703; + + std::vector> curves = marchingTriangles(*geomPtr, values, -0.891336467690912); + + ASSERT_EQ(curves.size(), 1); + EXPECT_EQ(curves[0].size(), 3); +} + +TEST(MarchingTrianglesSuite, FullyIsoValuedFaceIsSkipped) { + std::unique_ptr meshPtr; + std::unique_ptr geomPtr; + std::tie(meshPtr, geomPtr) = buildSingleTriangleMesh(); + + ManifoldSurfaceMesh& mesh = *meshPtr; + VertexData values(mesh, 0.); + + std::vector> curves = marchingTriangles(*geomPtr, values, 0.); + + EXPECT_TRUE(curves.empty()); +}