fix: marching triangles isovalue degeneracies

src/surface/marching_triangles.cpp

@@ -1,10 +1,68 @@
 #include "geometrycentral/surface/marching_triangles.h"
 
+#include <algorithm>
+#include <cmath>
+
 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<SurfacePoint>& hits, const SurfacePoint& p) {
+  if (std::find(hits.begin(), hits.end(), p) == hits.end()) {
+    hits.push_back(p);
+  }
+}
+
+void addEdgeIsoHits(Halfedge he, const VertexData<double>& u, double isoval, std::vector<SurfacePoint>& 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<std::vector<std::array<size_t, 2>>>
 getCurveComponents(SurfaceMesh& mesh, const std::vector<SurfacePoint>& curveNodes,
                    const std::vector<std::array<size_t, 2>>& curveEdges) {
@@ -70,28 +128,21 @@ std::vector<std::vector<SurfacePoint>> marchingTriangles(IntrinsicGeometryInterf
     std::vector<SurfacePoint> 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<size_t, 2> seg;
@@ -126,4 +177,4 @@ std::vector<std::vector<SurfacePoint>> marchingTriangles(IntrinsicGeometryInterf
 }
 
 } // namespace surface
-} // namespace geometrycentral
+} // namespace geometrycentral

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

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 <memory>
+#include <tuple>
+#include <vector>
+
+using namespace geometrycentral;
+using namespace geometrycentral::surface;
+
+namespace {
+
+std::tuple<std::unique_ptr<ManifoldSurfaceMesh>, std::unique_ptr<VertexPositionGeometry>> buildSingleTriangleMesh() {
+  std::vector<Vector3> verts = {{0, 0, 0}, {1, 0, 0}, {0, 1, 0}};
+  std::vector<std::vector<size_t>> faces = {{0, 1, 2}};
+  return makeManifoldSurfaceMeshAndGeometry(faces, verts);
+}
+
+std::tuple<std::unique_ptr<ManifoldSurfaceMesh>, std::unique_ptr<VertexPositionGeometry>> buildTwoTriangleMesh() {
+  std::vector<Vector3> verts = {{0, 0, 0}, {1, 0, 0}, {1, 1, 0}, {0, 1, 0}};
+  std::vector<std::vector<size_t>> faces = {{0, 1, 2}, {0, 2, 3}};
+  return makeManifoldSurfaceMeshAndGeometry(faces, verts);
+}
+
+} // namespace
+
+TEST(MarchingTrianglesSuite, ExactIsoVertexCreatesSegment) {
+  std::unique_ptr<ManifoldSurfaceMesh> meshPtr;
+  std::unique_ptr<VertexPositionGeometry> geomPtr;
+  std::tie(meshPtr, geomPtr) = buildSingleTriangleMesh();
+
+  ManifoldSurfaceMesh& mesh = *meshPtr;
+  VertexData<double> values(mesh);
+  values[mesh.vertex(0)] = 0.;
+  values[mesh.vertex(1)] = 1.;
+  values[mesh.vertex(2)] = -1.;
+
+  std::vector<std::vector<SurfacePoint>> 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<ManifoldSurfaceMesh> meshPtr;
+  std::unique_ptr<VertexPositionGeometry> geomPtr;
+  std::tie(meshPtr, geomPtr) = buildTwoTriangleMesh();
+
+  ManifoldSurfaceMesh& mesh = *meshPtr;
+  VertexData<double> 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<std::vector<SurfacePoint>> curves = marchingTriangles(*geomPtr, values, -0.891336467690912);
+
+  ASSERT_EQ(curves.size(), 1);
+  EXPECT_EQ(curves[0].size(), 3);
+}
+
+TEST(MarchingTrianglesSuite, FullyIsoValuedFaceIsSkipped) {
+  std::unique_ptr<ManifoldSurfaceMesh> meshPtr;
+  std::unique_ptr<VertexPositionGeometry> geomPtr;
+  std::tie(meshPtr, geomPtr) = buildSingleTriangleMesh();
+
+  ManifoldSurfaceMesh& mesh = *meshPtr;
+  VertexData<double> values(mesh, 0.);
+
+  std::vector<std::vector<SurfacePoint>> curves = marchingTriangles(*geomPtr, values, 0.);
+
+  EXPECT_TRUE(curves.empty());
+}