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

CHANGELOG.md

@@ -9,6 +9,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ### Added
 
+* Added `protect_boundary`, `protect_sharp_edges_angle_deg`, and `keep_points` parameters to `compas_cgal.meshing.trimesh_remesh`. `protect_boundary` constrains all boundary edges (no split/collapse/flip) so an open mesh's border curve — including sharp corners — is preserved verbatim instead of being re-sampled to `target_edge_length`; `protect_sharp_edges_angle_deg` marks interior edges whose dihedral angle ≥ the threshold as constrained (`0.0` disables interior-feature detection); `keep_points` (an Nx3 array) pins only the mesh vertices coincident with the given points (matched by coordinate within tolerance, via `vertex_is_constrained_map`), so specific vertices — e.g. the four corners of an open patch — survive while the rest of the boundary is still re-sampled.
+* Reworked `docs/examples/example_meshing.py` / `.md` to demonstrate all three boundary modes on the RhinoVault shell side by side (default, `protect_boundary=True`, and `keep_points` with the four corners), replacing the example image.
 * Added `compas_cgal.straight_skeleton_2.extrude` wrapping CGAL's `extrude_skeleton`, turning a 2D polygon (optionally with holes) into a closed 3D roof mesh, with control over roof pitch via taper `angles` or straight skeleton `weights` and an optional `maximum_height`. Internally coplanar extrusion faces are merged into polygons (via `remesh_planar_patches`) to recover the true roof planes, and `extrude` returns a `(mesh, lines)` tuple: a triangulated, ready-to-display mesh (so non-convex roof faces render correctly) and the roof outline (eaves, hips and ridges) as a list of `compas.geometry.Line`.
 * Added the `extrude_straight_skeleton` binding in `src/straight_skeleton_2.cpp` / `.h`.
 * Added `docs/examples/example_straight_skeleton_2_extrude.py` and `.md` (with image `docs/assets/images/example_straight_skeleton_2_extrude.png`) demonstrating roof generation for a few corner cases — a footprint without holes, one with a single hole, and one with multiple holes — and wired the example into the Examples nav in `mkdocs.yml`.

docs/examples/example_meshing.md

@@ -2,14 +2,23 @@
 
 ![Mesh Remeshing](../assets/images/example_meshing.png)
 
-This example demonstrates how to remesh a triangle mesh using COMPAS CGAL.
+This example demonstrates how to remesh a triangle mesh with COMPAS CGAL, and the three ways `trimesh_remesh` can treat the boundary of an open mesh.
+
+A RhinoVault funicular shell is triangulated and remeshed to the same coarse target edge length three ways, shown side by side (each drawn over the faint original for reference):
+
+* **Left — `protect_boundary=False`** (default): boundary edges are re-sampled to the target length, so the shell's open perimeter is coarsened and its corners rounded.
+* **Middle — `protect_boundary=True`**: every boundary edge is constrained, so the whole perimeter is preserved verbatim while the interior is coarsened.
+* **Right — `keep_points=<4 corners>`**: only the four corner vertices are pinned (matched by coordinate). The rest of the boundary is still re-sampled, but those four points (shown in red) survive exactly.
+
+`protect_boundary` and `keep_points` are complementary: the first keeps the entire boundary curve, the second keeps only the specific vertices you name. The companion parameter `protect_sharp_edges_angle_deg` similarly constrains interior feature edges whose dihedral angle exceeds a threshold.
 
 Key Features:
 
 * Loading PLY mesh files
-* Mesh transformation and centering
-* Remeshing with target edge length
-* Side-by-side visualization of original and remeshed models
+* Remeshing with a target edge length via `trimesh_remesh`
+* Preserving the whole boundary with `protect_boundary`
+* Pinning only specific vertices with `keep_points`
+* Side-by-side visualization of the three results over the original
 
 ```python
 ---8<--- "docs/examples/example_meshing.py"

docs/examples/example_meshing.py

@@ -1,7 +1,10 @@
 from pathlib import Path
 
 from compas.datastructures import Mesh
+from compas.geometry import Pointcloud
+from compas.geometry import Translation
 from compas_viewer import Viewer
+from compas_viewer.scene.tagobject import Tag
 
 from compas_cgal.meshing import trimesh_remesh
 
@@ -17,20 +20,70 @@
 VF = mesh.to_vertices_and_faces()
 
 # =============================================================================
-# Remesh
+# Footprint bounds + the four base corners
 # =============================================================================
 
-V1, F1 = trimesh_remesh(VF, target_edge_length=1, number_of_iterations=10)
+coords = mesh.vertices_attributes("xyz")
+xmin = min(x for x, _, _ in coords)
+xmax = max(x for x, _, _ in coords)
+ymin = min(y for _, y, _ in coords)
+ymax = max(y for _, y, _ in coords)
+zmax = max(z for _, _, z in coords)
 
-remeshed = Mesh.from_vertices_and_faces(V1, F1)
+# The four base corners: the boundary vertices nearest the xy-extremes of the
+# footprint. Their coordinates are what we hand to ``keep_points``.
+boundary_xyz = [mesh.vertex_coordinates(v) for v in mesh.vertices_on_boundary()]
 
-# ==============================================================================
-# Visualize
-# ==============================================================================
+
+def nearest_corner(tx, ty):
+    return min(boundary_xyz, key=lambda p: (p[0] - tx) ** 2 + (p[1] - ty) ** 2)
+
+
+corners = [
+    nearest_corner(xmin, ymin),
+    nearest_corner(xmax, ymin),
+    nearest_corner(xmax, ymax),
+    nearest_corner(xmin, ymax),
+]
+
+# =============================================================================
+# Remesh three ways to a coarse target (larger than the input boundary spacing)
+# =============================================================================
+
+TARGET_EDGE_LENGTH = 3
+ITERATIONS = 10
+
+# 1) Default: the open perimeter is coarsened and its corners rounded.
+V0, F0 = trimesh_remesh(VF, TARGET_EDGE_LENGTH, ITERATIONS, protect_boundary=False)
+
+# 2) protect_boundary=True: every boundary edge is constrained, so the whole
+#    perimeter is preserved verbatim while the interior is coarsened.
+V1, F1 = trimesh_remesh(VF, TARGET_EDGE_LENGTH, ITERATIONS, protect_boundary=True)
+
+# 3) keep_points: only the four corners are pinned — the rest of the boundary
+#    is still re-sampled, but those four vertices survive exactly.
+V2, F2 = trimesh_remesh(VF, TARGET_EDGE_LENGTH, ITERATIONS, keep_points=corners)
+
+# =============================================================================
+# Visualize — three panels side by side, each over the faint original
+# =============================================================================
+
+dx = 1.3 * (xmax - xmin)  # offset between panels
+xmid = 0.5 * (xmin + xmax)
 
 viewer = Viewer(width=1600, height=900)
 
-viewer.scene.add(mesh, show_points=False, opacity=0.25)
-viewer.scene.add(remeshed, show_points=True)
+labels = ["protect_boundary = False", "protect_boundary = True", "keep_points = 4 corners"]
+for i, (V, F) in enumerate([(V0, F0), (V1, F1), (V2, F2)]):
+    shift = Translation.from_vector([i * dx, 0.0, 0.0])
+    remeshed = Mesh.from_vertices_and_faces(V.tolist(), F.tolist()).transformed(shift)
+    original = mesh.transformed(shift)
+    viewer.scene.add(original, show_points=False, opacity=0.25)
+    viewer.scene.add(remeshed, show_points=True)
+    viewer.scene.add(Tag(labels[i], (xmid + i * dx, ymin - 2.0, zmax + 1.0), height=30))
+
+# Emphasize the four kept corners on the third panel.
+kept = Pointcloud(corners).transformed(Translation.from_vector([2 * dx, 0.0, 0.0]))
+viewer.scene.add(kept, pointcolor=(1.0, 0.0, 0.0), pointsize=20, show_points=True)
 
 viewer.show()

src/compas_cgal/meshing.py

@@ -14,6 +14,9 @@ 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,
+    keep_points=None,
 ) -> VerticesFacesNumpy:
     """Remeshing of a triangle mesh.
 
@@ -27,15 +30,33 @@ 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).
+    keep_points : array-like, optional
+        An Nx3 array of points. The mesh vertex coincident with each point
+        (matched by coordinate, within tolerance) is pinned: it is neither
+        moved nor removed during remeshing, while the edges between such
+        points are still re-sampled. Use this to keep only specific vertices
+        — e.g. the four corners of an open patch — rather than the whole
+        boundary. ``None`` disables (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 +73,19 @@ 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)
+    if keep_points is None:
+        keep_points = np.empty((0, 3), dtype=np.float64)
+    keep_points = np.asarray(keep_points, dtype=np.float64, order="C").reshape(-1, 3)
+    return _meshing.pmp_trimesh_remesh(
+        V,
+        F,
+        target_edge_length,
+        number_of_iterations,
+        do_project,
+        protect_boundary,
+        protect_sharp_edges_angle_deg,
+        keep_points,
+    )
 
 
 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,90 @@ 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,
+    const compas::RowMatrixXd& keep_points)
 {
     // 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);
 
-    // 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));
+    // 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);
+    }
+
+    // Build a vertex-is-constrained property map from user-supplied points.
+    // A mesh vertex whose coordinates match a provided point (within
+    // tolerance) is pinned: isotropic_remeshing will neither move nor remove
+    // it, while the edges between such points are still re-sampled. This is
+    // the coordinate-snap behaviour used by pmp_trimesh_remesh_dual.
+    auto vcm = mesh_a.add_property_map<
+        boost::graph_traits<compas::Mesh>::vertex_descriptor, bool>(
+        "v:is_constrained", false).first;
+
+    bool has_kept_vertices = false;
+    if (keep_points.rows() > 0) {
+        for (Eigen::Index i = 0; i < keep_points.rows(); ++i) {
+            compas::Kernel::Point_3 target(keep_points(i, 0), keep_points(i, 1), keep_points(i, 2));
+            for (auto v : vertices(mesh_a)) {
+                if (CGAL::squared_distance(mesh_a.point(v), target) < 1e-6) {
+                    put(vcm, v, true);
+                    has_kept_vertices = true;
+                    break;
+                }
+            }
+        }
+    }
+
+    // Perform isotropic remeshing. The vertex-is-constrained map is only
+    // attached when points were supplied — otherwise CGAL's default corner
+    // handling (derived from the constrained edges) must be left untouched.
+    const bool protect = protect_boundary || protect_sharp_edges_angle_deg > 0.0;
+    if (has_kept_vertices) {
+        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)
+                           .edge_is_constrained_map(ecm)
+                           .vertex_is_constrained_map(vcm)
+                           .protect_constraints(protect));
+    } else {
+        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)
+                           .edge_is_constrained_map(ecm)
+                           .protect_constraints(protect));
+    }
+
 
-
     // Clean up the mesh
     mesh_a.collect_garbage();
 
@@ -888,7 +958,10 @@ 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,
+        "keep_points"_a = compas::RowMatrixXd(0, 3)
     );
 
     m.def(

src/meshing.h

@@ -53,6 +53,11 @@ namespace compas {
  * @param target_edge_length Desired length for mesh edges
  * @param number_of_iterations Number of remeshing iterations
  * @param do_project Whether to project vertices onto the input surface
+ * @param protect_boundary Constrain all boundary edges so they are preserved
+ * @param protect_sharp_edges_angle_deg Dihedral threshold for interior feature edges (0 disables)
+ * @param keep_points Nx3 matrix of points; the mesh vertex coincident with each
+ *                    (within tolerance) is pinned and survives remeshing verbatim.
+ *                    Empty (0x3) disables (default).
  * @return std::tuple<RowMatrixXd, RowMatrixXi> containing:
  *         - New vertices as Rx3 matrix (float64)
  *         - New faces as Sx3 matrix (int32)
@@ -63,7 +68,10 @@ 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,
+    const compas::RowMatrixXd& keep_points = compas::RowMatrixXd(0, 3));
 
 
 /**