fix: make clip_with_plane idempotent and deduplicate the kernel round trip

src/mesh.cpp

@@ -8,6 +8,7 @@
 #include <algorithm>
 #include <array>
 #include <cmath>
+#include <limits>
 #include <CGAL/Polygon_mesh_processing/bbox.h>
 #include <CGAL/Polygon_mesh_processing/measure.h>
 #include <CGAL/Polygon_mesh_processing/intersection.h>
@@ -666,6 +667,40 @@ int TriMesh::clipWithPlane(double a, double b, double c, double d, bool use_exac
   if (_mesh.number_of_vertices() == 0 || _mesh.number_of_faces() == 0)
     throw std::invalid_argument("clipWithPlane: source mesh is empty");
 
+  // PMP::clip is not idempotent. Re-clipping a mesh that already lies entirely
+  // inside the halfspace removes no face, but still corefines the mesh against
+  // the plane and splits every edge along the existing seam. Repeating a plane
+  // therefore grows the mesh without bound -- and each exact round trip rounds
+  // the new seam points onto coincident doubles, which is what eventually takes
+  // PMP::clip down (Loop3D/loop-cgal#14). Nothing strictly outside the halfspace
+  // means nothing to remove, so return the documented no-op untouched.
+  const double normal_length = std::sqrt(a * a + b * b + c * c);
+  if (!(normal_length > 0.0))
+    throw std::invalid_argument("clipWithPlane: degenerate plane normal (a=b=c=0)");
+  double max_coord = 0.0;
+  double max_signed_distance = -std::numeric_limits<double>::infinity();
+  for (auto v : _mesh.vertices())
+  {
+    const auto &p = _mesh.point(v);
+    max_coord = std::max(max_coord, std::abs(p.x()));
+    max_coord = std::max(max_coord, std::abs(p.y()));
+    max_coord = std::max(max_coord, std::abs(p.z()));
+    const double signed_distance =
+        (a * p.x() + b * p.y() + c * p.z() + d) / normal_length;
+    max_signed_distance = std::max(max_signed_distance, signed_distance);
+  }
+  // Scale the tolerance to the mesh's own coordinate magnitude: the round trip
+  // displaces a seam vertex by a few ulps of that magnitude, never more.
+  const double on_plane_tol = 1e-12 * std::max(1.0, max_coord);
+  if (max_signed_distance <= on_plane_tol)
+  {
+    if (LoopCGAL::verbose)
+      std::cout << "clipWithPlane: mesh already inside the halfspace (max signed "
+                << "distance " << max_signed_distance << " <= " << on_plane_tol
+                << "), nothing to clip." << std::endl;
+    return 0;
+  }
+
   const int faces_before = static_cast<int>(_mesh.number_of_faces());
 
   // Capture constraints by geometry before the clip. The plane overload of

src/meshutils.cpp

@@ -269,63 +269,98 @@ NumpyMesh export_mesh(const TriangleMesh &tm, double area_threshold,
   result.triangles = triangles_array;
   return result;
 }
-Exact_Mesh convert_to_exact(const TriMesh& input) {
-  Exact_Mesh result;
-  // perf: Surface_mesh Vertex_index values are dense integers bounded by
-  // num_vertices() (vertex container size incl. removed slots). A flat vector
-  // indexed by the raw descriptor gives O(1) remap without per-vertex heap
-  // nodes; sizing by num_vertices() stays correct even if the mesh is not
-  // compact (removed vertices leave holes but never raise the index bound).
-  std::vector<Exact_Mesh::Vertex_index> vmap(input.get_mesh().num_vertices());
+namespace {
 
-  // perf: preallocate property arrays to the source's exact live counts so the
-  // 1:1 topology rebuild never incrementally reallocates as elements are added.
-  result.reserve(input.get_mesh().number_of_vertices(),
-                 input.get_mesh().number_of_edges(),
-                 input.get_mesh().number_of_faces());
-
-  for (auto v : vertices(input.get_mesh())) {
-    const auto& p = input.get_mesh().point(v);
-    Exact_K::Point_3 ep(p.x(), p.y(), p.z());
-    vmap[static_cast<std::size_t>(v)] = result.add_vertex(ep);
+// A mesh can end up with two vertices at the same point from either direction:
+// PMP::clip splits seam edges and rounding exact coordinates back to double
+// collapses points that differ below double precision, while write_to_file /
+// read_from_file bypass export_mesh's dedup so the duplicates survive a reload.
+// PMP::clip cannot cope with them -- the next clip along an edge joining such a
+// pair degenerates (Loop3D/loop-cgal#14). Rebuilding through this key merges on
+// exact double equality, which is lossless: the two vertices ARE one point.
+struct PointKey {
+  double x, y, z;
+  bool operator==(const PointKey &o) const noexcept {
+    return x == o.x && y == o.y && z == o.z;
   }
+};
 
-  for (auto f : faces(input.get_mesh())) {
-    std::vector<Exact_Mesh::Vertex_index> face_vertices;
-    for (auto v : vertices_around_face(input.get_mesh().halfedge(f), input.get_mesh())) {
-      face_vertices.push_back(vmap[static_cast<std::size_t>(v)]);
-    }
-    result.add_face(face_vertices);
+struct PointKeyHash {
+  std::size_t operator()(const PointKey &k) const noexcept {
+    std::size_t h = std::hash<double>{}(k.x);
+    h ^= std::hash<double>{}(k.y) + 0x9e3779b97f4a7c15ULL + (h << 6) + (h >> 2);
+    h ^= std::hash<double>{}(k.z) + 0x9e3779b97f4a7c15ULL + (h << 6) + (h >> 2);
+    return h;
   }
+};
+
+// The exact and double rebuilds differ only in how a point is converted, so
+// they share one 1:1 topology walk.
+template <typename OutMesh, typename InMesh, typename MakePoint>
+OutMesh rebuild_unique(const InMesh &in, MakePoint make_point) {
+  using OutVertex = typename OutMesh::Vertex_index;
+  OutMesh out;
 
-  return result;
-}
-TriangleMesh convert_to_double_mesh(const Exact_Mesh& input) {
-  TriangleMesh result;
   // perf: flat vector remap keyed by the dense Vertex_index integer, sized by
   // num_vertices() (safe upper bound incl. removed slots). Replaces std::map to
   // drop the O(log n) lookup and per-vertex heap allocation.
-  std::vector<TriangleMesh::Vertex_index> vmap(input.num_vertices());
+  std::vector<OutVertex> vmap(in.num_vertices());
 
   // perf: preallocate property arrays to the source's exact live counts so the
   // 1:1 topology rebuild never incrementally reallocates as elements are added.
-  result.reserve(input.number_of_vertices(),
-                 input.number_of_edges(),
-                 input.number_of_faces());
-
-  for (auto v : vertices(input)) {
-    const auto& p = input.point(v);
-    Point dp(CGAL::to_double(p.x()), CGAL::to_double(p.y()), CGAL::to_double(p.z()));
-    vmap[static_cast<std::size_t>(v)] = result.add_vertex(dp);
+  out.reserve(in.number_of_vertices(), in.number_of_edges(), in.number_of_faces());
+
+  std::unordered_map<PointKey, OutVertex, PointKeyHash> unique_points;
+  unique_points.reserve(in.number_of_vertices());
+
+  for (auto v : vertices(in)) {
+    const auto &p = in.point(v);
+    PointKey key{CGAL::to_double(p.x()), CGAL::to_double(p.y()),
+                 CGAL::to_double(p.z())};
+    auto it = unique_points.find(key);
+    if (it == unique_points.end()) {
+      it = unique_points.emplace(key, out.add_vertex(make_point(p))).first;
+    }
+    vmap[static_cast<std::size_t>(v)] = it->second;
   }
 
-  for (auto f : faces(input)) {
-    std::vector<TriangleMesh::Vertex_index> face_vertices;
-    for (auto v : vertices_around_face(input.halfedge(f), input)) {
+  for (auto f : faces(in)) {
+    std::vector<OutVertex> face_vertices;
+    for (auto v : vertices_around_face(in.halfedge(f), in)) {
       face_vertices.push_back(vmap[static_cast<std::size_t>(v)]);
     }
-    result.add_face(face_vertices);
+    // A merge can collapse two corners of a face onto one vertex. The face then
+    // has zero area and no orientation, so drop it instead of handing add_face
+    // a degenerate triangle.
+    bool degenerate = false;
+    for (std::size_t i = 0; i < face_vertices.size() && !degenerate; ++i) {
+      for (std::size_t j = i + 1; j < face_vertices.size(); ++j) {
+        if (face_vertices[i] == face_vertices[j]) {
+          degenerate = true;
+          break;
+        }
+      }
+    }
+    if (degenerate) {
+      continue;
+    }
+    out.add_face(face_vertices);
   }
 
-  return result;
-}
+  return out;
+}
+
+} // namespace
+
+Exact_Mesh convert_to_exact(const TriMesh& input) {
+  return rebuild_unique<Exact_Mesh>(
+      input.get_mesh(),
+      [](const Point &p) { return Exact_K::Point_3(p.x(), p.y(), p.z()); });
+}
+
+TriangleMesh convert_to_double_mesh(const Exact_Mesh& input) {
+  return rebuild_unique<TriangleMesh>(input, [](const Exact_K::Point_3 &p) {
+    return Point(CGAL::to_double(p.x()), CGAL::to_double(p.y()),
+                 CGAL::to_double(p.z()));
+  });
+}

tests/test_clip_idempotence.py

@@ -0,0 +1,118 @@
+"""Regression tests for Loop3D/loop-cgal#14.
+
+Two coupled defects, both exercised here:
+
+1. ``clip_with_plane`` was not idempotent. Re-clipping a mesh that already lay
+   entirely inside the halfspace removed no face but still corefined against the
+   plane, splitting every edge along the existing seam. Repeating a plane grew
+   the mesh without bound -- the "long runtime" half of the issue.
+2. The exact round trip rebuilt the mesh with no deduplication, so the surplus
+   seam points rounded onto coincident doubles. ``PMP::clip`` cannot cope with
+   those, which is the reported segfault (observed on Linux/x86-64 only).
+"""
+
+from __future__ import annotations
+
+import numpy as np
+import pytest
+from scipy.spatial import Delaunay
+
+import loop_cgal
+
+
+def raw_mesh(mesh):
+    """Vertices/triangles with the export dedup disabled, so this reports the
+    mesh's own vertex set rather than a cleaned copy of it."""
+    return mesh.get_vertices_and_triangles(area_threshold=0.0, duplicate_vertex_threshold=0.0)
+
+
+def n_coincident(vertices):
+    return len(vertices) - len(np.unique(vertices, axis=0))
+
+
+@pytest.fixture
+def bumpy_surface():
+    """Irregular, non-axis-aligned surface: clip intersections land on
+    high-bit-depth rationals, which is what makes the rounding collide."""
+    rng = np.random.default_rng(0)
+    pts = rng.uniform(-10, 10, size=(250, 2))
+    z = 1.5 * np.sin(pts[:, 0] * 0.7) + 1.1 * np.cos(pts[:, 1] * 0.5) + rng.normal(0, 0.05, 250)
+    return loop_cgal.TriMesh.from_vertices_and_triangles(
+        np.column_stack([pts, z]).astype(np.float64),
+        Delaunay(pts).simplices.astype(np.int32),
+    )
+
+
+# The plane that cuts the fixture in two.
+PLANE = (-0.5063653104002493, -0.45214174227658877, -0.7342765264628509, -0.3550784848736548)
+
+
+def test_first_clip_removes_faces(bumpy_surface):
+    """Guard must not suppress a clip that has something to remove."""
+    before = bumpy_surface.n_cells
+    assert bumpy_surface.clip_with_plane(*PLANE) > 0
+    assert bumpy_surface.n_cells < before
+
+
+@pytest.mark.parametrize("use_exact_kernel", [True, False])
+def test_repeated_clip_is_idempotent(bumpy_surface, use_exact_kernel):
+    """The same plane twice is a no-op: 0 faces removed, mesh untouched."""
+    bumpy_surface.clip_with_plane(*PLANE, use_exact_kernel)
+    settled_v, settled_t = raw_mesh(bumpy_surface)
+
+    for _ in range(20):
+        assert bumpy_surface.clip_with_plane(*PLANE, use_exact_kernel) == 0
+
+    v, t = raw_mesh(bumpy_surface)
+    np.testing.assert_array_equal(v, settled_v)
+    np.testing.assert_array_equal(t, settled_t)
+
+
+@pytest.mark.parametrize("use_exact_kernel", [True, False])
+def test_repeated_clip_breeds_no_coincident_vertices(bumpy_surface, use_exact_kernel):
+    for _ in range(20):
+        bumpy_surface.clip_with_plane(*PLANE, use_exact_kernel)
+        v, _ = raw_mesh(bumpy_surface)
+        assert n_coincident(v) == 0
+
+
+def test_cycled_planes_reach_a_fixed_point(bumpy_surface):
+    """The faultgen shape: a pool of near-tangent planes applied over and over.
+    Mesh size must settle rather than grow with every pass."""
+    rng = np.random.default_rng(3)
+    v0, _ = raw_mesh(bumpy_surface)
+    planes = []
+    for _ in range(8):
+        normal = rng.normal(size=3)
+        normal /= np.linalg.norm(normal)
+        offset = -float(np.quantile(v0 @ normal, 0.97))
+        planes.append(tuple(float(x) for x in normal) + (offset,))
+
+    for i in range(len(planes) * 2):
+        bumpy_surface.clip_with_plane(*planes[i % len(planes)])
+    settled = bumpy_surface.n_points
+
+    for i in range(len(planes) * 3):
+        bumpy_surface.clip_with_plane(*planes[i % len(planes)])
+        assert n_coincident(raw_mesh(bumpy_surface)[0]) == 0
+    assert bumpy_surface.n_points == settled
+
+
+def test_coincident_input_vertices_are_merged():
+    """A mesh that arrives already carrying duplicates -- read_from_file bypasses
+    export_mesh's dedup -- must be cleaned on the way into the exact kernel."""
+    verts = np.array(
+        [[0.0, 0, 0], [1, 0, 0], [0, 1, 0], [1, 1, 0], [1, 0, 0], [0, 1, 0]],
+        dtype=np.float64,
+    )  # rows 4 and 5 duplicate rows 1 and 2
+    tris = np.array([[0, 1, 2], [4, 3, 5]], dtype=np.int32)
+    mesh = loop_cgal.TriMesh.from_vertices_and_triangles(verts, tris)
+    assert n_coincident(raw_mesh(mesh)[0]) == 2
+
+    mesh.clip_with_plane(1.0, 0.0, 0.0, -0.75)  # keeps x <= 0.75
+    assert n_coincident(raw_mesh(mesh)[0]) == 0
+
+
+def test_degenerate_plane_normal_raises(bumpy_surface):
+    with pytest.raises(ValueError):
+        bumpy_surface.clip_with_plane(0.0, 0.0, 0.0, 1.0)