diff --git a/src/mesh.cpp b/src/mesh.cpp index ae82afd..c5c01c9 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -8,6 +8,7 @@ #include #include #include +#include #include #include #include @@ -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::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(_mesh.number_of_faces()); // Capture constraints by geometry before the clip. The plane overload of diff --git a/src/meshutils.cpp b/src/meshutils.cpp index 4d1491a..a589ad4 100644 --- a/src/meshutils.cpp +++ b/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 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(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 face_vertices; - for (auto v : vertices_around_face(input.get_mesh().halfedge(f), input.get_mesh())) { - face_vertices.push_back(vmap[static_cast(v)]); - } - result.add_face(face_vertices); +struct PointKeyHash { + std::size_t operator()(const PointKey &k) const noexcept { + std::size_t h = std::hash{}(k.x); + h ^= std::hash{}(k.y) + 0x9e3779b97f4a7c15ULL + (h << 6) + (h >> 2); + h ^= std::hash{}(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 +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 vmap(input.num_vertices()); + std::vector 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(v)] = result.add_vertex(dp); + out.reserve(in.number_of_vertices(), in.number_of_edges(), in.number_of_faces()); + + std::unordered_map 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(v)] = it->second; } - for (auto f : faces(input)) { - std::vector face_vertices; - for (auto v : vertices_around_face(input.halfedge(f), input)) { + for (auto f : faces(in)) { + std::vector face_vertices; + for (auto v : vertices_around_face(in.halfedge(f), in)) { face_vertices.push_back(vmap[static_cast(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; -} \ No newline at end of file + return out; +} + +} // namespace + +Exact_Mesh convert_to_exact(const TriMesh& input) { + return rebuild_unique( + 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(input, [](const Exact_K::Point_3 &p) { + return Point(CGAL::to_double(p.x()), CGAL::to_double(p.y()), + CGAL::to_double(p.z())); + }); +} diff --git a/tests/test_clip_idempotence.py b/tests/test_clip_idempotence.py new file mode 100644 index 0000000..e325b82 --- /dev/null +++ b/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)