feat(meshing): trimesh_remesh — protect_boundary + sharp-feature edge constraints

src/compas_cgal/meshing.py

@@ -14,6 +14,8 @@ def trimesh_remesh(
     target_edge_length: float,
     number_of_iterations: int = 10,
     do_project: bool = True,
+    protect_boundary: bool = False,
+    protect_sharp_edges_angle_deg: float = 0.0,
 ) -> VerticesFacesNumpy:
     """Remeshing of a triangle mesh.
 
@@ -27,15 +29,26 @@ def trimesh_remesh(
         Number of remeshing iterations.
     do_project : bool, optional
         If True, reproject vertices onto the input surface when they are created or displaced.
+    protect_boundary : bool, optional
+        If True, constrain all boundary edges so they are NOT split, collapsed,
+        or flipped during remeshing. Use this to preserve the input mesh's
+        boundary curve verbatim — including sharp corners that the default
+        smoothing pass would otherwise round.
+    protect_sharp_edges_angle_deg : float, optional
+        Dihedral threshold in degrees for interior feature detection. Edges
+        whose adjacent faces form a dihedral angle ≥ this value are marked
+        constrained and preserved. ``0.0`` disables interior-feature
+        detection (default).
 
     Returns
     -------
     VerticesFacesNumpy
 
     Notes
     -----
-    This remeshing function only constrains the edges on the boundary of the mesh.
-    Protecting specific features or edges is not implemented yet.
+    Without ``protect_boundary`` or ``protect_sharp_edges_angle_deg`` set,
+    boundary edges follow CGAL's default remeshing behaviour: re-sampled
+    to ``target_edge_length`` (visible corner rounding is the cost).
 
     Examples
     --------
@@ -52,7 +65,15 @@ def trimesh_remesh(
     V, F = mesh
     V = np.asarray(V, dtype=np.float64, order="C")
     F = np.asarray(F, dtype=np.int32, order="C")
-    return _meshing.pmp_trimesh_remesh(V, F, target_edge_length, number_of_iterations, do_project)
+    return _meshing.pmp_trimesh_remesh(
+        V,
+        F,
+        target_edge_length,
+        number_of_iterations,
+        do_project,
+        protect_boundary,
+        protect_sharp_edges_angle_deg,
+    )
 
 
 def trimesh_dual(

src/meshing.cpp

@@ -5,6 +5,7 @@
 #include <CGAL/boost/graph/Dual.h>
 #include <CGAL/boost/graph/helpers.h>
 #include <CGAL/Polygon_mesh_processing/compute_normal.h>
+#include <CGAL/Polygon_mesh_processing/detect_features.h>
 
 #include <iostream>
 #include <fstream>
@@ -43,21 +44,52 @@ pmp_trimesh_remesh(
     Eigen::Ref<const compas::RowMatrixXi> faces_a,
     double target_edge_length,
     unsigned int number_of_iterations,
-    bool do_project)
+    bool do_project,
+    bool protect_boundary,
+    double protect_sharp_edges_angle_deg)
 {
     // Convert input matrices to CGAL mesh and keep a copy for projection
     compas::Mesh original_mesh = compas::mesh_from_vertices_and_faces(vertices_a, faces_a);
     compas::Mesh mesh_a = compas::mesh_from_vertices_and_faces(vertices_a, faces_a);
 
+    // Build an edge-is-constrained property map. Edges marked True are
+    // not split / collapsed / flipped during isotropic_remeshing — this
+    // is how features are preserved.
+    auto ecm = mesh_a.add_property_map<
+        boost::graph_traits<compas::Mesh>::edge_descriptor, bool>(
+        "e:is_constrained", false).first;
+
+    // Constrain all boundary edges when requested. CGAL's
+    // isotropic_remeshing otherwise re-samples boundary edges per
+    // target_edge_length, which rounds visible corners.
+    if (protect_boundary) {
+        for (auto e : edges(mesh_a)) {
+            auto h = halfedge(e, mesh_a);
+            if (is_border(h, mesh_a) || is_border(opposite(h, mesh_a), mesh_a)) {
+                put(ecm, e, true);
+            }
+        }
+    }
+
+    // Detect sharp interior edges (dihedral > threshold) and constrain
+    // them too. 0.0 disables (default).
+    if (protect_sharp_edges_angle_deg > 0.0) {
+        CGAL::Polygon_mesh_processing::detect_sharp_edges(
+            mesh_a, protect_sharp_edges_angle_deg, ecm);
+    }
+
     // Perform isotropic remeshing
     CGAL::Polygon_mesh_processing::isotropic_remeshing(
         faces(mesh_a),
         target_edge_length,
         mesh_a,
         CGAL::Polygon_mesh_processing::parameters::number_of_iterations(number_of_iterations)
-                       .do_project(do_project));
+                       .do_project(do_project)
+                       .edge_is_constrained_map(ecm)
+                       .protect_constraints(protect_boundary
+                                            || protect_sharp_edges_angle_deg > 0.0));
+
 
-
     // Clean up the mesh
     mesh_a.collect_garbage();
 
@@ -888,7 +920,9 @@ NB_MODULE(_meshing, m) {
         "faces_a"_a,
         "target_edge_length"_a,
         "number_of_iterations"_a = 10,
-        "do_project"_a = true
+        "do_project"_a = true,
+        "protect_boundary"_a = false,
+        "protect_sharp_edges_angle_deg"_a = 0.0
     );
 
     m.def(

src/meshing.h

@@ -63,7 +63,9 @@ pmp_trimesh_remesh(
     Eigen::Ref<const compas::RowMatrixXi> faces_a,
     double target_edge_length,
     unsigned int number_of_iterations = 10,
-    bool do_project = true);
+    bool do_project = true,
+    bool protect_boundary = false,
+    double protect_sharp_edges_angle_deg = 0.0);
 
 
 /**

tests/test_meshing.py

@@ -59,6 +59,104 @@ def test_remesh(sample_mesh):
     assert remeshed_mesh.is_valid()
 
 
+@pytest.fixture
+def square_with_hole():
+    """Annular layer: square outer boundary, square hole in the middle.
+
+    Tests boundary preservation on a topology with TWO loops (outer + hole).
+    """
+    outer = [[0.0, 0.0, 0.0], [4.0, 0.0, 0.0], [4.0, 4.0, 0.0], [0.0, 4.0, 0.0]]
+    inner = [[1.0, 1.0, 0.0], [3.0, 1.0, 0.0], [3.0, 3.0, 0.0], [1.0, 3.0, 0.0]]
+    V = np.asarray(outer + inner, dtype=np.float64)
+    # Fan-triangulate the annulus: connect each outer corner to two inner
+    # corners. 8 triangles total.
+    F = np.asarray(
+        [
+            [0, 1, 5], [0, 5, 4],
+            [1, 2, 6], [1, 6, 5],
+            [2, 3, 7], [2, 7, 6],
+            [3, 0, 4], [3, 4, 7],
+        ],
+        dtype=np.int32,
+    )
+    return V, F
+
+
+def _corner_set(V, atol=1e-6):
+    return {tuple(np.round(p / atol).astype(int)) for p in V}
+
+
+def test_remesh_default_subdivides_boundary(square_with_hole):
+    """Without protect_boundary the boundary IS subdivided (default CGAL behaviour)."""
+    V, F = square_with_hole
+    V_new, F_new = trimesh_remesh((V, F), target_edge_length=0.5, number_of_iterations=5)
+    # Boundary subdivision adds vertices on outer + inner loops.
+    assert V_new.shape[0] > V.shape[0], (
+        f"default mode should add boundary verts; got {V_new.shape[0]} <= {V.shape[0]}"
+    )
+
+
+def test_remesh_protect_boundary_keeps_corners(square_with_hole):
+    """With protect_boundary=True every original corner survives verbatim."""
+    V, F = square_with_hole
+    V_new, F_new = trimesh_remesh(
+        (V, F),
+        target_edge_length=0.5,
+        number_of_iterations=5,
+        protect_boundary=True,
+    )
+    new_corners = _corner_set(V_new)
+    orig_corners = _corner_set(V)
+    missing = orig_corners - new_corners
+    assert not missing, f"protect_boundary lost original corners: {missing}"
+
+
+def test_remesh_protect_boundary_keeps_boundary_vertex_count(square_with_hole):
+    """protect_boundary=True must not insert NEW vertices along the boundary loops.
+
+    Bounds the boundary-vertex count to the original 8 (4 outer + 4 inner).
+    """
+    V, F = square_with_hole
+    V_new, F_new = trimesh_remesh(
+        (V, F),
+        target_edge_length=0.5,
+        number_of_iterations=5,
+        protect_boundary=True,
+    )
+    # Count vertices that lie on the original boundary segments. The 8
+    # original corners are all collinear with the outer/inner square edges.
+    # Any remeshed vertex on a boundary edge can be detected by checking
+    # if its half-edge has no twin in the new mesh.
+    he_set = set()
+    for face in F_new:
+        for k in range(3):
+            i, j = int(face[k]), int(face[(k + 1) % 3])
+            he_set.add((i, j))
+    boundary_verts = set()
+    for (i, j) in he_set:
+        if (j, i) not in he_set:
+            boundary_verts.add(i)
+            boundary_verts.add(j)
+    # Original boundary had 8 verts (4 outer + 4 inner); protect_boundary
+    # forbids splitting boundary edges so count must stay exactly 8.
+    assert len(boundary_verts) == 8, (
+        f"protect_boundary=True must keep 8 boundary verts; got {len(boundary_verts)}"
+    )
+
+
+def test_remesh_protect_sharp_edges_default_disabled():
+    """protect_sharp_edges_angle_deg=0.0 (default) must not affect output."""
+    # Two co-planar triangles sharing an interior edge — no sharp dihedral.
+    V = np.asarray([[0, 0, 0], [1, 0, 0], [1, 1, 0], [0, 1, 0]], dtype=np.float64)
+    F = np.asarray([[0, 1, 2], [0, 2, 3]], dtype=np.int32)
+    V_a, F_a = trimesh_remesh((V, F), 0.3, number_of_iterations=5)
+    V_b, F_b = trimesh_remesh(
+        (V, F), 0.3, number_of_iterations=5, protect_sharp_edges_angle_deg=0.0
+    )
+    np.testing.assert_array_equal(V_a.shape, V_b.shape)
+    np.testing.assert_array_equal(F_a.shape, F_b.shape)
+
+
 def test_dual(sample_mesh):
     """Test the dual functionality."""
     # Get mesh data