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 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()); +}