From 5ed79bfd9e5f7242b5170b1eceaa4280ab723ffa Mon Sep 17 00:00:00 2001 From: lachlangrose Date: Mon, 28 Sep 2026 20:48:54 +0930 Subject: [PATCH] feat: snap implicit cuts to existing vertices, add does_self_intersect Follow-on work to the mesh review, sitting between the review follow-ups and the clip fix. - cut_with_implicit_function now snaps a crossing to an existing vertex when it lands within snap_tol (as a fraction of the edge), instead of emitting a sliver an infinitesimal distance from it. A per-edge dedup guard stops the second visit to a shared edge pushing a duplicate crossing vertex, which was tearing the mesh along every seam and compounding over sequential cuts. The snap is recorded per slot (newtris_on_seam), not by writing the isovalue into the shared per-vertex array. The per-vertex form leaks a per-edge fact into every sub-triangle incident to that vertex: a legitimate keep-side sub-triangle then reads as entirely on-seam and an inclusive keep/discard test drops it from BOTH sides, silently deleting kept surface wherever the seam runs through a vertex. The filter now asks whether any original (non-crossing) vertex lies strictly on the keep side, and applies only to the one-sided modes so PRESERVE_INTERSECTION still keeps the whole surface. - does_self_intersect() exposes PMP::does_self_intersect. A self-intersecting surface is outside the domain of the boolean operations, so this is the check to run before handing a mesh to them. - Document the threading contract: most methods release the GIL, so concurrent access to one TriMesh is a data race in C++, not just a wrong answer. corefine additionally mutates its `other` argument, which was undocumented. - perf: convert_to_exact / convert_to_double_mesh use a flat vector remap keyed on the dense Vertex_index instead of std::map, and preallocate to the source's live counts; set_mesh gains an rvalue overload so the converters' temporaries move rather than deep-copy; snapshot_constraint_coords(interior_only) skips border edges the caller is about to re-derive anyway. Tests: the snap partition (both sides sum to the input area, and snapping does not move the cut) and does_self_intersect. 82 passing. Co-Authored-By: Claude Opus 5 (1M context) --- .gitignore | 7 +- loop_cgal/bindings.cpp | 40 ++++- src/mesh.cpp | 276 +++++++++++++++++++++----------- src/mesh.h | 15 +- src/meshutils.cpp | 32 +++- tests/test_implicit_cut.py | 15 +- tests/test_implicit_cut_snap.py | 50 ++++++ tests/test_mesh_operations.py | 29 ++++ 8 files changed, 358 insertions(+), 106 deletions(-) create mode 100644 tests/test_implicit_cut_snap.py diff --git a/.gitignore b/.gitignore index 225d4d8..42b1193 100644 --- a/.gitignore +++ b/.gitignore @@ -44,4 +44,9 @@ loop_cgal/__pycache__/__init__.cpython-311.pyc build/* *.vtk *.ply -build_win64/* \ No newline at end of file +build_win64/* + +# Local import shim: symlinks __init__.py / utils.py / build/_loop_cgal.so into +# one directory so `PYTHONPATH=.local-shim` imports the locally built extension. +# Needed because the source `loop_cgal/` package has no .so and shadows the build. +.local-shim/ diff --git a/loop_cgal/bindings.cpp b/loop_cgal/bindings.cpp index 4e6aa27..f5f98f3 100644 --- a/loop_cgal/bindings.cpp +++ b/loop_cgal/bindings.cpp @@ -19,7 +19,14 @@ PYBIND11_MODULE(_loop_cgal, m) .value("KEEP_POSITIVE_SIDE", ImplicitCutMode::KEEP_POSITIVE_SIDE) .value("KEEP_NEGATIVE_SIDE", ImplicitCutMode::KEEP_NEGATIVE_SIDE) .export_values(); - py::class_(m, "TriMesh") + py::class_(m, "TriMesh", + "A CGAL surface mesh of triangles.\n\n" + "NOT THREAD-SAFE. Most methods release the GIL for the duration of the " + "CGAL call, so two Python threads touching the same TriMesh concurrently " + "are a data race in C++ -- corrupt memory, not merely a wrong answer. That " + "includes passing one mesh as the `other`/`surface` argument of a call " + "running on another thread. Give each thread its own meshes, or serialise " + "access yourself. Separate processes (the usual pattern here) are unaffected.") .def(py::init &, const pybind11::array_t &>(), py::arg("vertices"), py::arg("triangles")) .def(py::init([](const TriMesh& other) { return other.clone(); }), @@ -28,13 +35,16 @@ PYBIND11_MODULE(_loop_cgal, m) .def("clip_with_plane", &TriMesh::clipWithPlane, py::arg("a"), py::arg("b"), py::arg("c"), py::arg("d"), py::arg("use_exact_kernel") = true, + py::call_guard(), "Clip the mesh with the halfspace ax+by+cz+d <= 0. " "Uses PMP::clip(mesh, Plane_3) directly — no corefinement, no skirt construction. " "Returns the number of faces removed (0 = no-op).") .def("corefine", &TriMesh::corefine, py::arg("other"), py::arg("use_exact_kernel") = true, + py::call_guard(), "Corefine this mesh with another, inserting the shared intersection " - "polyline into BOTH meshes in place. After the call the two meshes " + "polyline into BOTH meshes in place. NOTE this MUTATES `other` as well " + "as self, so `other` must not be shared with a concurrent call. After the call the two meshes " "carry coincident vertices/edges along the intersection curve, so " "patches taken from each side stitch together watertight. " "Returns the number of vertices added to this mesh. " @@ -44,23 +54,31 @@ PYBIND11_MODULE(_loop_cgal, m) py::arg("preserve_intersection") = false, py::arg("preserve_intersection_clipper") = false, py::arg("use_exact_kernel") = true, + py::arg("extend") = true, + py::call_guard(), "Cut mesh with a clipper surface. Returns the number of faces removed " - "(0 indicates the clipper did not intersect or did not extend beyond the mesh).") + "(0 indicates the clipper did not intersect or did not extend beyond the mesh). " + "extend=True grows a skirt so a small clipper cuts all the way through; pass " + "extend=False when the clipper already spans the domain. " + "Mutates self and reads `surface`; neither may be touched by a concurrent call.") .def("remesh", &TriMesh::remesh, py::arg("split_long_edges") = true, py::arg("target_edge_length") = 10.0, py::arg("number_of_iterations") = 3, py::arg("protect_constraints") = true, - py::arg("relax_constraints") = false) + py::arg("relax_constraints") = false, + py::call_guard()) .def("save", &TriMesh::save, py::arg("area_threshold") = 1e-6, py::arg("duplicate_vertex_threshold") = 1e-6) .def("reverse_face_orientation", &TriMesh::reverseFaceOrientation, "Reverse the face orientation of the mesh.") .def("add_fixed_edges", &TriMesh::add_fixed_edges, py::arg("pairs"), + py::call_guard(), "Vertex index pairs defining edges to be fixed in mesh when remeshing.") .def("cut_with_implicit_function", &TriMesh::cut_with_implicit_function, py::arg("property"), py::arg("value"),py::arg("cutmode") = ImplicitCutMode::KEEP_POSITIVE_SIDE, py::arg("snap_tol") = 1e-4, + py::call_guard(), "Cut the mesh with an implicit function defined by vertex properties. " "snap_tol reuses an existing vertex when a crossing lands within that " "fraction of an edge from it, preventing sliver triangles (0 = disable).") @@ -69,7 +87,19 @@ PYBIND11_MODULE(_loop_cgal, m) .def("_cgal_n_vertices", &TriMesh::n_vertices, "Number of vertices in the CGAL mesh.") .def("_cgal_points", &TriMesh::get_points, "Vertex coordinates as (n, 3) numpy array.") .def("overlaps", &TriMesh::overlaps, py::arg("other"), py::arg("bbox_tol") = 1e-6, - "Return True if this mesh intersects another TriMesh (AABB fast-reject + triangle-level check).") + py::call_guard(), + "Return True if this mesh intersects another TriMesh (AABB fast-reject + triangle-level check). " + "Reads both meshes; neither may be mutated by a concurrent call.") + .def("does_self_intersect", &TriMesh::does_self_intersect, + py::call_guard(), + "Return True if the mesh intersects itself.\n\n" + "Wraps CGAL PMP::does_self_intersect. A self-intersecting surface is " + "outside the domain of the boolean operations - corefine, cut_with_surface " + "and the clip paths may fail or return a corrupt result on one - so this is " + "the check to run before handing a mesh to them, and after a remesh or a fit " + "that could have folded it.\n\n" + "Cost is O(n log n) via an AABB tree over the faces: cheap enough to assert " + "on, not free in a tight loop.") .def("clone", &TriMesh::clone, "Return a deep copy of this TriMesh at the C++ level (no pyvista round-trip). " "Preserves the full fixed-edge set including any user-added edges.") diff --git a/src/mesh.cpp b/src/mesh.cpp index 3e25cf6..ae82afd 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -5,6 +5,9 @@ #include #include #include +#include +#include +#include #include #include #include @@ -12,6 +15,7 @@ #include #include #include +#include #include #include #include @@ -116,7 +120,8 @@ void TriMesh::init() _edge_is_constrained_map = CGAL::make_boolean_property_map(_fixedEdges); } -std::vector> TriMesh::snapshot_constraint_coords() const +std::vector> +TriMesh::snapshot_constraint_coords(bool interior_only) const { std::vector> coords; coords.reserve(_fixedEdges.size()); @@ -128,6 +133,9 @@ std::vector> TriMesh::snapshot_constraint_coords() const if (!_mesh.has_valid_index(e) || _mesh.is_removed(e)) continue; auto h = _mesh.halfedge(e); + if (interior_only && + (_mesh.is_border(h) || _mesh.is_border(_mesh.opposite(h)))) + continue; coords.emplace_back(_mesh.point(_mesh.source(h)), _mesh.point(_mesh.target(h))); } return coords; @@ -380,7 +388,8 @@ void TriMesh::reverseFaceOrientation() int TriMesh::cutWithSurface(TriMesh &clipper, bool preserve_intersection, bool preserve_intersection_clipper, - bool use_exact_kernel) + bool use_exact_kernel, + bool extend) { if (LoopCGAL::verbose) { @@ -419,82 +428,91 @@ int TriMesh::cutWithSurface(TriMesh &clipper, // Interior vertices and faces are completely untouched, so the cut location // is preserved exactly for both planar and curved (listric) clippers. // ----------------------------------------------------------------------- - CGAL::Bbox_3 target_bb = PMP::bbox(_mesh); - const double target_diag = std::sqrt( - CGAL::square(target_bb.xmax() - target_bb.xmin()) + - CGAL::square(target_bb.ymax() - target_bb.ymin()) + - CGAL::square(target_bb.zmax() - target_bb.zmin())); - // Copy the clipper and stitch any collocated border vertices first. // stitch_borders can remove/merge vertices, so all subsequent passes must // operate on extended_mesh (not clipper._mesh) to keep vertex descriptors valid. TriangleMesh extended_mesh = clipper._mesh; PMP::stitch_borders(extended_mesh); - // Pass 1: accumulate per-vertex outward directions from each border halfedge - // in extended_mesh (post-stitch). For a border halfedge h (source→target), - // the outward direction is fn × d, where fn is the adjacent face normal and - // d is the normalised edge direction. This lies in the face's tangent plane - // and points away from the interior. - std::map outward_sum; - for (auto he : extended_mesh.halfedges()) + // Grow the clipper by extruding a skirt of new triangles outward from each + // boundary edge, so a clipper smaller than the target still cuts all the way + // through. Skipped when extend=false: callers passing an already-domain- + // spanning clipper must NOT be + // re-extended — the skirt of a large surface folds on concave boundaries and + // corefinement then rejects the near-coplanar overlap. + if (extend) { - if (!extended_mesh.is_border(he)) continue; - - const Point &ps = extended_mesh.point(extended_mesh.source(he)); - const Point &pt = extended_mesh.point(extended_mesh.target(he)); - Vector d = pt - ps; - const double d_len = std::sqrt(d.squared_length()); - if (d_len < 1e-10) continue; - d = d / d_len; - - const auto adj_face = extended_mesh.face(extended_mesh.opposite(he)); - const Vector fn = PMP::compute_face_normal(adj_face, extended_mesh); - const Vector out = CGAL::cross_product(fn, d); - const double out_len = std::sqrt(out.squared_length()); - if (out_len < 1e-10) continue; - - auto vs = extended_mesh.source(he); - auto vt = extended_mesh.target(he); - outward_sum.try_emplace(vs, 0.0, 0.0, 0.0); - outward_sum.try_emplace(vt, 0.0, 0.0, 0.0); - outward_sum[vs] = outward_sum[vs] + out / out_len; - outward_sum[vt] = outward_sum[vt] + out / out_len; - } + CGAL::Bbox_3 target_bb = PMP::bbox(_mesh); + const double target_diag = std::sqrt( + CGAL::square(target_bb.xmax() - target_bb.xmin()) + + CGAL::square(target_bb.ymax() - target_bb.ymin()) + + CGAL::square(target_bb.zmax() - target_bb.zmin())); + + // Pass 1: accumulate per-vertex outward directions from each border halfedge + // in extended_mesh (post-stitch). For a border halfedge h (source→target), + // the outward direction is fn × d, where fn is the adjacent face normal and + // d is the normalised edge direction. This lies in the face's tangent plane + // and points away from the interior. + std::map outward_sum; + for (auto he : extended_mesh.halfedges()) + { + if (!extended_mesh.is_border(he)) continue; + + const Point &ps = extended_mesh.point(extended_mesh.source(he)); + const Point &pt = extended_mesh.point(extended_mesh.target(he)); + Vector d = pt - ps; + const double d_len = std::sqrt(d.squared_length()); + if (d_len < 1e-10) continue; + d = d / d_len; + + const auto adj_face = extended_mesh.face(extended_mesh.opposite(he)); + const Vector fn = PMP::compute_face_normal(adj_face, extended_mesh); + const Vector out = CGAL::cross_product(fn, d); + const double out_len = std::sqrt(out.squared_length()); + if (out_len < 1e-10) continue; + + auto vs = extended_mesh.source(he); + auto vt = extended_mesh.target(he); + outward_sum.try_emplace(vs, 0.0, 0.0, 0.0); + outward_sum.try_emplace(vt, 0.0, 0.0, 0.0); + outward_sum[vs] = outward_sum[vs] + out / out_len; + outward_sum[vt] = outward_sum[vt] + out / out_len; + } - // Pass 2: add one new vertex per boundary vertex pushed outward by target_diag, - // and stitch a skirt quad (two triangles) per boundary edge of extended_mesh. - // The winding order (source, target, target_new) produces normals consistent - // with the adjacent interior face. - std::map skirt_vertex; - for (auto &[v, dir_sum] : outward_sum) - { - const double len = std::sqrt(dir_sum.squared_length()); - if (len < 1e-10) continue; - const Vector dir = dir_sum / len; - const Point &p = extended_mesh.point(v); - skirt_vertex[v] = extended_mesh.add_vertex( - Point(p.x() + target_diag * dir.x(), - p.y() + target_diag * dir.y(), - p.z() + target_diag * dir.z())); - } + // Pass 2: add one new vertex per boundary vertex pushed outward by target_diag, + // and stitch a skirt quad (two triangles) per boundary edge of extended_mesh. + // The winding order (source, target, target_new) produces normals consistent + // with the adjacent interior face. + std::map skirt_vertex; + for (auto &[v, dir_sum] : outward_sum) + { + const double len = std::sqrt(dir_sum.squared_length()); + if (len < 1e-10) continue; + const Vector dir = dir_sum / len; + const Point &p = extended_mesh.point(v); + skirt_vertex[v] = extended_mesh.add_vertex( + Point(p.x() + target_diag * dir.x(), + p.y() + target_diag * dir.y(), + p.z() + target_diag * dir.z())); + } - for (auto he : extended_mesh.halfedges()) - { - if (!extended_mesh.is_border(he)) continue; - auto vs = extended_mesh.source(he); - auto vt = extended_mesh.target(he); - if (!skirt_vertex.count(vs) || !skirt_vertex.count(vt)) continue; - auto vs_new = skirt_vertex[vs]; - auto vt_new = skirt_vertex[vt]; - extended_mesh.add_face(vs, vt, vt_new); - extended_mesh.add_face(vs, vt_new, vs_new); - } + for (auto he : extended_mesh.halfedges()) + { + if (!extended_mesh.is_border(he)) continue; + auto vs = extended_mesh.source(he); + auto vt = extended_mesh.target(he); + if (!skirt_vertex.count(vs) || !skirt_vertex.count(vt)) continue; + auto vs_new = skirt_vertex[vs]; + auto vt_new = skirt_vertex[vt]; + extended_mesh.add_face(vs, vt, vt_new); + extended_mesh.add_face(vs, vt_new, vs_new); + } - if (LoopCGAL::verbose) - std::cout << " cutWithSurface: added skirt of " - << skirt_vertex.size() << " new vertices over target_diag=" - << target_diag << "\n"; + if (LoopCGAL::verbose) + std::cout << " cutWithSurface: added skirt of " + << skirt_vertex.size() << " new vertices over target_diag=" + << target_diag << "\n"; + } TriMesh scaled_clipper(std::move(extended_mesh)); @@ -727,6 +745,10 @@ void TriMesh::cut_with_implicit_function(const std::vector &property, do // renumbered 0..n-1 contiguously, so idx() is always in range. if (_mesh.has_garbage()) _mesh.collect_garbage(); + // The cut replaces _mesh wholesale, so the constraint set has to be re-derived + // afterwards (see the end of this function). Only user-added interior + // constraints need carrying across by coordinate; borders are re-derived. + auto saved_constraints = snapshot_constraint_coords(/*interior_only=*/true); if (property.size() != _mesh.number_of_vertices()) throw std::invalid_argument( "cut_with_implicit_function: property size (" + @@ -830,6 +852,12 @@ void TriMesh::cut_with_implicit_function(const std::vector &property, do std::map new_point_on_edge; std::vector> newtris(tri_array.begin(), tri_array.end()); + // Which slots of each triangle are seam crossings rather than original + // vertices. A crossing may reuse an EXISTING vertex index (see snapping + // below), so the index alone cannot tell the two apart, and the keep/discard + // filter needs to know: a crossing lies on the isosurface by construction and + // says nothing about which side the sub-triangle belongs to. + std::vector> newtris_on_seam(newtris.size(), {false, false, false}); if (LoopCGAL::verbose) { std::cout << "Starting main loop over " << tri_array.size() << " triangles." << std::endl; @@ -844,6 +872,15 @@ void TriMesh::cut_with_implicit_function(const std::vector &property, do // for each edge of tri, check if edge crosses for (auto eid : tri2edge[t]) { + // A crossed edge is shared by two straddling triangles, so it is visited + // twice. Without this guard the second visit pushes a second, coincident + // crossing vertex and overwrites new_point_on_edge[eid] with it, leaving + // the two triangles referencing different vertices at the same position: + // the mesh is torn along every seam, and the tear compounds over + // sequential cuts. The crossing is a function of the edge alone, so + // reusing the vertex recorded by the first visit is exact. + if (new_point_on_edge.count(eid)) + continue; auto ends = edge_array[eid]; double f0 = vertex_properties[ends.first]; double f1 = vertex_properties[ends.second]; @@ -871,6 +908,13 @@ void TriMesh::cut_with_implicit_function(const std::vector &property, do // dropped by the degeneracy check when the new mesh is assembled, so no // sliver is ever created. newverts starts as a copy of verts, so the // original vertex indices ends.first/ends.second are valid indices into it. + // Snapping declares that existing vertex to BE the crossing on THIS edge. + // That is a per-edge fact, not a per-vertex one: the same vertex is still + // strictly on its own side of the isosurface as far as every other + // triangle around it is concerned. Recording it in newtris_on_seam (per + // slot) rather than by overwriting newvals (per vertex) keeps the two + // apart — the latter leaks the assertion into every incident sub-triangle + // and deletes kept surface wherever the seam runs through a vertex. if (ratio <= snap_tol) { new_point_on_edge[eid] = ends.first; @@ -924,9 +968,21 @@ void TriMesh::cut_with_implicit_function(const std::vector &property, do else sub = {{{tri[0], tri[1], c12}, {tri[0], c12, c20}, {c20, c12, tri[2]}}}; + // Slot-for-slot with `sub` above: true where the entry is c01/c12/c20. + std::array, 3> sub_on_seam; + if (lone == 0) + sub_on_seam = {{{false, true, true}, {true, false, false}, {true, true, false}}}; + else if (lone == 1) + sub_on_seam = {{{false, true, false}, {true, true, false}, {true, false, true}}}; + else + sub_on_seam = {{{false, false, true}, {false, true, true}, {true, true, false}}}; + newtris[t] = sub[0]; + newtris_on_seam[t] = sub_on_seam[0]; newtris.push_back(sub[1]); + newtris_on_seam.push_back(sub_on_seam[1]); newtris.push_back(sub[2]); + newtris_on_seam.push_back(sub_on_seam[2]); } // Build new CGAL mesh from newverts and newtris @@ -935,42 +991,76 @@ void TriMesh::cut_with_implicit_function(const std::vector &property, do new_vhandles.reserve(newverts.size()); for (auto &p : newverts) new_vhandles.push_back(newmesh.add_vertex(p)); - for (auto &tri : newtris) + for (std::size_t ti = 0; ti < newtris.size(); ++ti) { + const auto &tri = newtris[ti]; + const auto &on_seam = newtris_on_seam[ti]; // skip degenerate if (tri[0] == tri[1] || tri[1] == tri[2] || tri[0] == tri[2]) continue; - // Drop sub-triangles on the unwanted side. Comparisons are inclusive of the - // seam: a wrong-side seam sub-triangle has two crossing vertices whose value - // is exactly `value` plus one vertex strictly on the wrong side, so a strict - // comparison would spare it and leave a one-triangle fringe past the - // isovalue. Using <=/>= drops it, so the cut stops exactly at the seam. - // (NaN values compare false either way, so off-extent triangles are kept.) - if (ImplicitCutMode::KEEP_NEGATIVE_SIDE == cutmode) + + // A sub-triangle belongs to the side its ORIGINAL (non-crossing) vertices + // lie on; by construction they all lie on the same side. Crossing vertices + // sit on the isosurface and carry no side information, so they are skipped + // rather than compared — which is what stops a wrong-side seam fringe (two + // crossings plus one strictly-wrong-side vertex) from surviving, and equally + // stops a legitimate keep-side sliver from being discarded. + // A NaN value means the implicit function was not evaluated there (off + // extent); such triangles are kept, as before. + // Only the one-sided modes discard anything; PRESERVE_INTERSECTION splits + // along the seam but keeps the whole surface. + if (ImplicitCutMode::KEEP_POSITIVE_SIDE == cutmode || + ImplicitCutMode::KEEP_NEGATIVE_SIDE == cutmode) { - double v0 = newvals[tri[0]]; - double v1 = newvals[tri[1]]; - double v2 = newvals[tri[2]]; - if (v0 >= value && v1 >= value && v2 >= value) + const bool keep_positive = (ImplicitCutMode::KEEP_POSITIVE_SIDE == cutmode); + bool on_keep_side = false; + bool undetermined = false; + for (int k = 0; k < 3; ++k) { - continue; + if (on_seam[k]) + continue; + const double v = newvals[tri[k]]; + if (std::isnan(v)) + undetermined = true; + else if (keep_positive ? (v > value) : (v < value)) + on_keep_side = true; } - } - if (ImplicitCutMode::KEEP_POSITIVE_SIDE == cutmode) - { - double v0 = newvals[tri[0]]; - double v1 = newvals[tri[1]]; - double v2 = newvals[tri[2]]; - if (v0 <= value && v1 <= value && v2 <= value) - { + if (!on_keep_side && !undetermined) continue; - } } + newmesh.add_face(new_vhandles[tri[0]], new_vhandles[tri[1]], new_vhandles[tri[2]]); } + // newverts carries every original vertex, including those used only by + // dropped triangles. Left in place they accumulate over sequential cuts, and + // the caller — which must size its property array to n_vertices — pays to + // evaluate its implicit function on all of them. + PMP::remove_isolated_vertices(newmesh); + if (newmesh.has_garbage()) + newmesh.collect_garbage(); + // Replace internal mesh _mesh = std::move(newmesh); + + // Re-derive the constraint set exactly as clipWithPlane does: border edges of + // the new mesh (including the fresh cut boundary) plus any user-added interior + // constraints that survived, re-resolved by coordinate. Without this, + // _fixedEdges keeps Edge_index values into the destroyed mesh, which alias + // unrelated edges of the new one and mark them protected during a later + // remesh(), and the cut boundary is never protected at all. + // + // Limitation (as for clipWithPlane): a constraint the seam passes *through* is + // not recovered, because the inserted crossing vertex breaks the adjacency of + // its original endpoints. + // + // With no interior constraints to carry across (the usual case) init() gives + // the same result and skips building a coordinate lookup over every vertex, + // which costs ~5 ms per cut on a 100k-vertex mesh. + if (saved_constraints.empty()) + init(); + else + rebuild_fixed_edges_from_coords(saved_constraints); } double TriMesh::area() const @@ -1122,3 +1212,9 @@ bool TriMesh::overlaps(const TriMesh& other, double bbox_tol) const } return PMP::do_intersect(_mesh, other._mesh); } + +// Does this mesh intersect itself? +bool TriMesh::does_self_intersect() const +{ + return PMP::does_self_intersect(_mesh); +} diff --git a/src/mesh.h b/src/mesh.h index f67362c..fc8c59f 100644 --- a/src/mesh.h +++ b/src/mesh.h @@ -40,7 +40,8 @@ class TriMesh int cutWithSurface(TriMesh &surface, bool preserve_intersection = false, bool preserve_intersection_clipper = false, - bool use_exact_kernel = true); + bool use_exact_kernel = true, + bool extend = true); // Clip the mesh with a halfspace defined by the plane ax+by+cz+d=0. // The negative side (ax+by+cz+d < 0) is kept. @@ -79,11 +80,16 @@ class TriMesh std::size_t n_vertices() const; pybind11::array_t get_points() const; bool overlaps(const TriMesh& other, double bbox_tol = 1e-6) const; + bool does_self_intersect() const; TriMesh clone() const; void write_to_file(const std::string& path) const; static TriMesh read_from_file(const std::string& path); const TriangleMesh& get_mesh() const { return _mesh; } void set_mesh(const TriangleMesh& mesh) { _mesh = mesh; } + // perf: rvalue overload avoids a deep Surface_mesh copy; all current + // callers pass a temporary (convert_to_double_mesh result). Matches the + // move pattern used by the private TriMesh(TriangleMesh) ctor. + void set_mesh(TriangleMesh&& mesh) { _mesh = std::move(mesh); } private: // Internal constructor used by cutWithSurface to wrap a scaled copy. explicit TriMesh(TriangleMesh m) : _mesh(std::move(m)) { init(); } @@ -92,7 +98,12 @@ class TriMesh // _fixedEdges that is still live in _mesh. Edge_index handles are not // stable across clip()/collect_garbage() or the exact-kernel round // trip, so constraints must be tracked by geometry, not by index. - std::vector> snapshot_constraint_coords() const; + // interior_only skips edges that are borders of the current _mesh: those + // are re-derived wholesale by the rebuild, so a caller that only needs + // to carry user-added interior constraints across can pass true and get + // an empty snapshot (and so skip the rebuild) in the common case. + std::vector> + snapshot_constraint_coords(bool interior_only = false) const; // Rebuild _fixedEdges (and re-bind the constrained-edge map) from a // coordinate snapshot: border edges of the current _mesh are always diff --git a/src/meshutils.cpp b/src/meshutils.cpp index 6b39e58..4d1491a 100644 --- a/src/meshutils.cpp +++ b/src/meshutils.cpp @@ -271,18 +271,29 @@ NumpyMesh export_mesh(const TriangleMesh &tm, double area_threshold, } Exact_Mesh convert_to_exact(const TriMesh& input) { Exact_Mesh result; - std::map vmap; + // 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()); + + // 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[v] = result.add_vertex(ep); + vmap[static_cast(v)] = result.add_vertex(ep); } 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[v]); + face_vertices.push_back(vmap[static_cast(v)]); } result.add_face(face_vertices); } @@ -291,18 +302,27 @@ Exact_Mesh convert_to_exact(const TriMesh& input) { } TriangleMesh convert_to_double_mesh(const Exact_Mesh& input) { TriangleMesh result; - std::map vmap; + // 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()); + + // 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[v] = result.add_vertex(dp); + vmap[static_cast(v)] = result.add_vertex(dp); } for (auto f : faces(input)) { std::vector face_vertices; for (auto v : vertices_around_face(input.halfedge(f), input)) { - face_vertices.push_back(vmap[v]); + face_vertices.push_back(vmap[static_cast(v)]); } result.add_face(face_vertices); } diff --git a/tests/test_implicit_cut.py b/tests/test_implicit_cut.py index 57d2bfe..e1b3244 100644 --- a/tests/test_implicit_cut.py +++ b/tests/test_implicit_cut.py @@ -117,13 +117,24 @@ def test_snap_tol_reuses_nearby_vertex(): tm = loop_cgal.TriMesh.from_vertices_and_triangles(verts, tris) n_before = tm.n_points # 4.0 is a grid line; 4.0 + 1e-6 is far inside snap_tol (1e-4 of a unit edge). + value = 4.0 + 1e-6 tm.cut_with_implicit_function( - verts[:, 0].tolist(), 4.0 + 1e-6, ImplicitCutMode.KEEP_POSITIVE_SIDE, 1e-4 + verts[:, 0].tolist(), value, ImplicitCutMode.KEEP_POSITIVE_SIDE, 1e-4 ) - assert tm.n_points == n_before, ( + # The count may only shrink: snapping reuses the grid-line vertices instead + # of inserting near-duplicates, and the cut then drops the vertices left + # isolated on the discarded side. + assert tm.n_points <= n_before, ( f"snap_tol should reuse existing vertices, but point count grew " f"{n_before} -> {tm.n_points}" ) + # The real invariant: every retained vertex still sits on a grid line — no + # near-duplicate was inserted a hair off x == 4.0 by the crossing itself. + xs = tm._cgal_points()[:, 0] + assert np.all(np.abs(xs - np.round(xs)) < 1e-12), ( + f"a near-duplicate crossing vertex was inserted off the grid line: " + f"{xs[np.abs(xs - np.round(xs)) >= 1e-12]}" + ) def test_keep_positive_side_is_tight(): diff --git a/tests/test_implicit_cut_snap.py b/tests/test_implicit_cut_snap.py new file mode 100644 index 0000000..0ecd456 --- /dev/null +++ b/tests/test_implicit_cut_snap.py @@ -0,0 +1,50 @@ +"""Regression test for the implicit-cut snap dropping kept surface. + +The snap recorded "this vertex is the crossing" by writing the isovalue into the +shared per-vertex array. That is a per-edge fact stored per-vertex, so every +sub-triangle incident to the vertex then read it as lying on the seam, and an +inclusive keep/discard test dropped legitimate keep-side sub-triangles from BOTH +sides -- silently deleting surface wherever the seam ran through a vertex. +""" + +from __future__ import annotations + +import numpy as np +import pytest + +import loop_cgal +from loop_cgal._loop_cgal import ImplicitCutMode + + +# --- 1. snapping must not delete kept surface ------------------------------ + +# Two triangles sharing edge 0-2. Vertex 0 sits a hair off the isosurface, and +# its neighbour 3 is far away, so the crossing on edge 0-3 lands within snap_tol +# of vertex 0 and snaps onto it. +SNAP_VERTS = np.array([[0.0, 0, 0], [1, 0, 0], [0, 1, 0], [-1, 0, 0]]) +SNAP_TRIS = np.array([[0, 1, 2], [0, 2, 3]], dtype=np.int32) +SNAP_PROP = [1e-3, 1e-3, -1e-3, -100.0] + + +def _cut_area(mode, snap_tol): + m = loop_cgal.TriMesh.from_vertices_and_triangles(SNAP_VERTS, SNAP_TRIS) + m.cut_with_implicit_function(SNAP_PROP, 0.0, mode, snap_tol) + return m.area + + +@pytest.mark.parametrize( + "mode", [ImplicitCutMode.KEEP_POSITIVE_SIDE, ImplicitCutMode.KEEP_NEGATIVE_SIDE] +) +def test_snapping_does_not_change_the_kept_area(mode): + """Snapping is a meshing detail; it must not move the cut.""" + assert _cut_area(mode, 1e-4) == pytest.approx(_cut_area(mode, 0.0), abs=1e-5) + + +def test_the_two_sides_partition_the_surface(): + """Every sub-triangle belongs to exactly one side — none dropped, none double + counted, even where the seam runs through an existing vertex.""" + total = 1.0 # the two input triangles + kept = _cut_area(ImplicitCutMode.KEEP_POSITIVE_SIDE, 1e-4) + _cut_area( + ImplicitCutMode.KEEP_NEGATIVE_SIDE, 1e-4 + ) + assert kept == pytest.approx(total, abs=1e-5) diff --git a/tests/test_mesh_operations.py b/tests/test_mesh_operations.py index bc485dd..998b496 100644 --- a/tests/test_mesh_operations.py +++ b/tests/test_mesh_operations.py @@ -409,3 +409,32 @@ def _add_collocated_seam(pv_mesh: pv.PolyData) -> pv.PolyData: # cut_with_implicit_function is covered thoroughly in tests/test_implicit_cut.py # (the previous smoke test here only asserted n_verts/n_faces >= 0, which is true # for any mesh and so could never fail). + + +# --------------------------------------------------------------------------- +# does_self_intersect +# --------------------------------------------------------------------------- + + +def test_does_self_intersect_false_for_a_clean_surface(unit_trimesh): + assert unit_trimesh.does_self_intersect() is False + + +def test_does_self_intersect_true_for_crossing_triangles(): + """Two triangles that pass through one another, sharing no vertex.""" + verts = np.array( + [[0.0, 0, 0], [2, 0, 0], [1, 1, 0], [1, -1, 1], [1, 2, -1], [1, 0.5, 0.0]] + ) + tris = np.array([[0, 1, 2], [3, 4, 5]], dtype=np.int32) + mesh = loop_cgal.TriMesh.from_vertices_and_triangles(verts, tris) + assert mesh.does_self_intersect() is True + + +def test_does_self_intersect_does_not_modify_the_mesh(unit_trimesh): + """The check is const — it must not compact, reorder or otherwise touch.""" + before_points = np.array(unit_trimesh.points, copy=True) + before = (unit_trimesh.n_points, unit_trimesh.n_cells, unit_trimesh.area) + unit_trimesh.does_self_intersect() + after = (unit_trimesh.n_points, unit_trimesh.n_cells, unit_trimesh.area) + assert before == after + np.testing.assert_array_equal(before_points, np.array(unit_trimesh.points))