From 51ee9a5c4bfa4bddc866223502078b01bbaf28c8 Mon Sep 17 00:00:00 2001 From: lachlangrose Date: Thu, 12 Mar 2026 08:35:13 +1030 Subject: [PATCH 1/2] Enhance cutWithSurface method to expand the clipper along the surface to ensure clippers are longer than surfaces being clipped. return number of faces post clip. --- loop_cgal/bindings.cpp | 4 +- src/mesh.cpp | 130 ++++++++++++++++++++++++++++++++++++----- src/mesh.h | 8 ++- 3 files changed, 125 insertions(+), 17 deletions(-) diff --git a/loop_cgal/bindings.cpp b/loop_cgal/bindings.cpp index 1546f4b..cd9fa21 100644 --- a/loop_cgal/bindings.cpp +++ b/loop_cgal/bindings.cpp @@ -25,7 +25,9 @@ PYBIND11_MODULE(_loop_cgal, m) .def("cut_with_surface", &TriMesh::cutWithSurface, py::arg("surface"), py::arg("preserve_intersection") = false, py::arg("preserve_intersection_clipper") = false, - py::arg("use_exact_kernel") = true) + py::arg("use_exact_kernel") = true, + "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).") .def("remesh", &TriMesh::remesh, py::arg("split_long_edges") = true, py::arg("target_edge_length") = 10.0, py::arg("number_of_iterations") = 3, diff --git a/src/mesh.cpp b/src/mesh.cpp index 157b2f3..74f83a4 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -5,9 +5,10 @@ #include #include #include -#include +#include #include #include +#include #include #include #include @@ -248,6 +249,7 @@ void TriMesh::remesh(bool split_long_edges, if (split_long_edges) if (LoopCGAL::verbose) std::cout << "Splitting long edges in iteration " << iter + 1 << ".\n"; + PMP::split_long_edges( edges(_mesh), target_edge_length, _mesh, CGAL::parameters::edge_is_constrained_map(_edge_is_constrained_map)); @@ -289,12 +291,11 @@ void TriMesh::reverseFaceOrientation() } } -void TriMesh::cutWithSurface(TriMesh &clipper, +int TriMesh::cutWithSurface(TriMesh &clipper, bool preserve_intersection, bool preserve_intersection_clipper, bool use_exact_kernel) { - if (LoopCGAL::verbose) { std::cout << "Cutting mesh with surface." << std::endl; @@ -304,29 +305,124 @@ void TriMesh::cutWithSurface(TriMesh &clipper, if (!CGAL::is_valid_polygon_mesh(_mesh, LoopCGAL::verbose)) { std::cerr << "Error: Source mesh is invalid!" << std::endl; - return; + return 0; } if (!CGAL::is_valid_polygon_mesh(clipper._mesh, LoopCGAL::verbose)) { std::cerr << "Error: Clipper mesh is invalid!" << std::endl; - return; + return 0; } if (_mesh.number_of_vertices() == 0 || _mesh.number_of_faces() == 0) { std::cerr << "Error: Source mesh is empty!" << std::endl; - return; + return 0; } if (clipper._mesh.number_of_vertices() == 0 || clipper._mesh.number_of_faces() == 0) { std::cerr << "Error: Clipper mesh is empty!" << std::endl; - return; + return 0; + } + + // Merge any collocated border vertices on the target before clipping. + // Collocated vertices produce zero-area faces whose degenerate bounding boxes + // can make PMP::do_intersect return false even when the meshes overlap. + PMP::stitch_borders(_mesh); + // Remove any isolated (orphan) vertices — they survive PMP::clip unchanged + // and would pollute the output point set with stale positions. + PMP::remove_isolated_vertices(_mesh); + + // ----------------------------------------------------------------------- + // Grow the clipper by extruding a skirt of new triangles outward from each + // boundary edge. Each skirt quad's normal matches its adjacent border face, + // so the extension is along the surface tangent — not a global scale. + // 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())); + + // Pass 1: accumulate per-ver tex outward directions from each adjacent border face. + // 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 its interior. + std::map outward_sum; + for (auto he : clipper._mesh.halfedges()) + { + if (!clipper._mesh.is_border(he)) continue; + + const Point &ps = clipper._mesh.point(clipper._mesh.source(he)); + const Point &pt = clipper._mesh.point(clipper._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 = clipper._mesh.face(clipper._mesh.opposite(he)); + const Vector fn = PMP::compute_face_normal(adj_face, clipper._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 = clipper._mesh.source(he); + auto vt = clipper._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: copy the clipper, merge any collocated border vertices (they + // produce degenerate faces whose normals are near-zero, which would cause + // skirt edges to be silently skipped and leave bridge gaps), then add one + // new vertex per boundary vertex pushed outward by target_diag, and stitch + // a skirt quad (two triangles) per boundary edge. The winding order + // (source, target, target_new) produces normals consistent with the + // adjacent interior face. + TriangleMesh extended_mesh = clipper._mesh; + PMP::stitch_borders(extended_mesh); + + 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 : clipper._mesh.halfedges()) + { + if (!clipper._mesh.is_border(he)) continue; + auto vs = clipper._mesh.source(he); + auto vt = clipper._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); } - bool intersection = PMP::do_intersect(_mesh, clipper._mesh); + 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)); + + const int faces_before = static_cast(_mesh.number_of_faces()); + + bool intersection = PMP::do_intersect(_mesh, scaled_clipper._mesh); if (intersection) { // Clip tm with clipper @@ -337,21 +433,18 @@ void TriMesh::cutWithSurface(TriMesh &clipper, try { - // bool flag = - // PMP::clip(_mesh, clipper._mesh, CGAL::parameters::clip_volume(false)); bool flag = false; try { if (use_exact_kernel){ - Exact_Mesh exact_clipper = convert_to_exact(clipper); + Exact_Mesh exact_clipper = convert_to_exact(scaled_clipper); Exact_Mesh exact_mesh = convert_to_exact(*this); flag = PMP::clip(exact_mesh, exact_clipper, CGAL::parameters::clip_volume(false)); - set_mesh(convert_to_double_mesh(exact_mesh)); + set_mesh(convert_to_double_mesh(exact_mesh)); } else{ - flag = PMP::clip(_mesh, clipper._mesh, CGAL::parameters::clip_volume(false)); + flag = PMP::clip(_mesh, scaled_clipper._mesh, CGAL::parameters::clip_volume(false)); } - } catch (const std::exception &e) { @@ -384,6 +477,15 @@ void TriMesh::cutWithSurface(TriMesh &clipper, << std::endl; } } + + const int faces_after = static_cast(_mesh.number_of_faces()); + if (faces_after >= faces_before && LoopCGAL::verbose) + { + std::cerr << "Warning: cutWithSurface removed no faces (before=" + << faces_before << ", after=" << faces_after + << "). Clipper may not extend beyond the mesh." << std::endl; + } + return faces_before - faces_after; } NumpyMesh TriMesh::save(double area_threshold, diff --git a/src/mesh.h b/src/mesh.h index c491b7e..a231bc7 100644 --- a/src/mesh.h +++ b/src/mesh.h @@ -28,8 +28,9 @@ class TriMesh TriMesh(const pybind11::array_t &vertices, const pybind11::array_t &triangles); - // Method to cut the mesh with another surface object - void cutWithSurface(TriMesh &surface, + // Method to cut the mesh with another surface object. + // Returns the number of faces removed (0 = no-op / bad cut). + int cutWithSurface(TriMesh &surface, bool preserve_intersection = false, bool preserve_intersection_clipper = false, bool use_exact_kernel = true); @@ -47,6 +48,9 @@ class TriMesh const TriangleMesh& get_mesh() const { return _mesh; } void set_mesh(const TriangleMesh& mesh) { _mesh = mesh; } private: + // Internal constructor used by cutWithSurface to wrap a scaled copy. + explicit TriMesh(TriangleMesh m) : _mesh(std::move(m)) { init(); } + std::set _fixedEdges; TriangleMesh _mesh; // The underlying CGAL surface mesh CGAL::Boolean_property_map> From 3a8648ae149dc225d3f3ba654869811823c4de71 Mon Sep 17 00:00:00 2001 From: lachlangrose Date: Thu, 12 Mar 2026 08:39:52 +1030 Subject: [PATCH 2/2] tests: add test for mesh clipping --- tests/test_mesh_operations.py | 167 ++++++++++++++++++++++++++++++++++ 1 file changed, 167 insertions(+) diff --git a/tests/test_mesh_operations.py b/tests/test_mesh_operations.py index 4eb9c5d..daf6fa6 100644 --- a/tests/test_mesh_operations.py +++ b/tests/test_mesh_operations.py @@ -60,6 +60,173 @@ def test_remesh_changes_vertices(square_surface, kwargs): # assert after_v >= 0.5 * before_v +def _make_curved_clipper(x_offset: float, y_half: float, z_top: float, z_bot: float, n: int = 6) -> pv.PolyData: + """Build a curved (listric) clipper surface that bows in the X direction. + + The surface lies roughly in the YZ plane at x=x_offset but curves so that + deeper vertices are offset further in X. This simulates a listric fault + being used as a clipper. + """ + ys = np.linspace(-y_half, y_half, n) + zs = np.linspace(z_top, z_bot, n) + verts = [] + for z in zs: + # curvature: bow increases with depth (small enough that cut stays within 50 m of X=0) + # The deepest target vertex is at z=-7000; max x-deviation = c * 7000 < tol=50 m → c < 0.007 + curve = x_offset + 0.005 * (z - z_top) + for y in ys: + verts.append([curve, y, z]) + verts = np.array(verts, dtype=float) + faces = [] + for iz in range(n - 1): + for iy in range(n - 1): + i0 = iz * n + iy + i1 = i0 + 1 + i2 = i0 + n + i3 = i2 + 1 + faces.extend([[3, i0, i1, i3], [3, i0, i3, i2]]) + faces = np.array(faces).flatten() + return pv.PolyData(verts, faces) + + +@pytest.mark.parametrize("clipper_factory,label", [ + ( + # Planar clipper: YZ plane at X=0, far smaller than the 4 km deep target + lambda: pv.Plane( + center=(0.0, 0.0, -250.0), + direction=(1.0, 0.0, 0.0), + i_size=500.0, + j_size=500.0, + i_resolution=4, + j_resolution=4, + ).triangulate(), + "planar", + ), + ( + # Curved (listric) clipper: also smaller than the target, bows in X with depth + lambda: _make_curved_clipper(x_offset=0.0, y_half=250.0, z_top=0.0, z_bot=-500.0), + "curved", + ), +]) +def test_no_bridge_when_clipper_shorter_than_target(clipper_factory, label): + """Clipper that doesn't extend fully through the target must not leave a bridge. + + Setup + ----- + Target : vertical fault surface in the XZ plane (Y=0), 10 km wide × 4 km deep. + Clipper : either a planar or curved (listric) surface at X≈0, only 500 m deep + — far smaller than the 4 km depth of the target. + + Without the fix, the lower portion (Z < -500 m) is left as a bridge. With + boundary extension, the clipper rim is pushed out to cover the full target so + the result must lie entirely on one side of X=0. The interior vertices of the + clipper are unchanged, so the cut location is preserved for both planar and + curved cases. + """ + # Large vertical fault in XZ plane (Y=0): 10 km wide, 4 km deep. + target_pv = pv.Plane( + center=(0.0, 0.0, -2000.0), + direction=(0.0, 1.0, 0.0), + i_size=10000.0, + j_size=4000.0, + i_resolution=20, + j_resolution=8, + ).triangulate() + + target = loop_cgal.TriMesh(target_pv) + clipper = loop_cgal.TriMesh(clipper_factory()) + + faces_removed = target.cut_with_surface(clipper) + result = target.to_pyvista() + + assert result.n_points > 0, f"[{label}] Clip removed the entire mesh" + assert faces_removed > 0, f"[{label}] Clip was a no-op — meshes may not intersect" + + # After a clean cut at X≈0, all remaining vertices must be on ONE side. + # A bridge would leave vertices significantly on both sides. + xs = result.points[:, 0] + tol = 50.0 # 50 m tolerance for vertices exactly on the cut boundary + on_positive = np.any(xs > tol) + on_negative = np.any(xs < -tol) + assert not (on_positive and on_negative), ( + f"[{label}] Bridge detected: vertices span both sides of the cut plane " + f"(X range [{xs.min():.1f}, {xs.max():.1f}] m)." + ) + + +@pytest.mark.parametrize("collocate_target,collocate_clipper,label", [ + (True, False, "collocated_target"), + (False, True, "collocated_clipper"), + (True, True, "collocated_both"), +]) +def test_no_artefact_with_collocated_vertices(collocate_target, collocate_clipper, label): + """Collocated (duplicate) vertices must not cause bridges or a no-op clip. + + Two failure modes are possible: + - Collocated vertices in the *target* create zero-area faces that can make + PMP::do_intersect return false, silently skipping the entire clip. + - Collocated vertices in the *clipper* create degenerate border faces whose + normals are near-zero, causing skirt edges to be skipped and leaving a + bridge gap at those positions. + + Both are fixed by calling PMP::stitch_borders before clipping. + """ + def _add_collocated_seam(pv_mesh: pv.PolyData) -> pv.PolyData: + """Duplicate a row of interior vertices along the mesh midline (Z=-2000). + + The duplicate vertices are collocated with the originals but topologically + distinct, creating zero-area triangles along the seam. + """ + pts = pv_mesh.points.copy() + # Find vertices near the midline (Z ≈ -2000) + seam = np.where(np.abs(pts[:, 2] + 2000.0) < 50.0)[0] + if len(seam) == 0: + return pv_mesh # nothing to duplicate + extra = pts[seam].copy() # exact same positions + new_pts = np.vstack([pts, extra]) + return pv.PolyData(new_pts, pv_mesh.faces) + + target_pv = pv.Plane( + center=(0.0, 0.0, -2000.0), + direction=(0.0, 1.0, 0.0), + i_size=10000.0, + j_size=4000.0, + i_resolution=20, + j_resolution=8, + ).triangulate() + + clipper_pv = pv.Plane( + center=(0.0, 0.0, -250.0), + direction=(1.0, 0.0, 0.0), + i_size=500.0, + j_size=500.0, + i_resolution=4, + j_resolution=4, + ).triangulate() + + if collocate_target: + target_pv = _add_collocated_seam(target_pv) + if collocate_clipper: + clipper_pv = _add_collocated_seam(clipper_pv) + + target = loop_cgal.TriMesh(target_pv) + clipper = loop_cgal.TriMesh(clipper_pv) + + faces_removed = target.cut_with_surface(clipper) + result = target.to_pyvista() + + assert result.n_points > 0, f"[{label}] Clip removed the entire mesh" + assert faces_removed > 0, f"[{label}] Clip was a no-op (do_intersect likely wrong due to degenerate faces)" + + xs = result.points[:, 0] + tol = 50.0 + on_positive = np.any(xs > tol) + on_negative = np.any(xs < -tol) + assert not (on_positive and on_negative), ( + f"[{label}] Bridge detected: X range [{xs.min():.1f}, {xs.max():.1f}] m" + ) + + def test_cut_with_implicit_function(square_surface): tm = loop_cgal.TriMesh(square_surface) # create a scalar property that varies across vertices