diff --git a/CHANGELOG.md b/CHANGELOG.md index 0fb82834..324c4033 100644 --- a/CHANGELOG.md +++ b/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`. diff --git a/docs/assets/images/example_meshing.png b/docs/assets/images/example_meshing.png index 78235a56..b775bfe5 100644 Binary files a/docs/assets/images/example_meshing.png and b/docs/assets/images/example_meshing.png differ diff --git a/docs/examples/example_meshing.md b/docs/examples/example_meshing.md index 54943087..6ad3a084 100644 --- a/docs/examples/example_meshing.md +++ b/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" diff --git a/docs/examples/example_meshing.py b/docs/examples/example_meshing.py index 54081172..ddc7fc12 100644 --- a/docs/examples/example_meshing.py +++ b/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() diff --git a/src/compas_cgal/meshing.py b/src/compas_cgal/meshing.py index 3c30fa17..1dd02f78 100644 --- a/src/compas_cgal/meshing.py +++ b/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,6 +30,23 @@ 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 ------- @@ -34,8 +54,9 @@ def trimesh_remesh( 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( diff --git a/src/meshing.cpp b/src/meshing.cpp index 3ac44f48..d973c095 100644 --- a/src/meshing.cpp +++ b/src/meshing.cpp @@ -5,6 +5,7 @@ #include #include #include +#include #include #include @@ -43,21 +44,90 @@ pmp_trimesh_remesh( Eigen::Ref 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::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::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( diff --git a/src/meshing.h b/src/meshing.h index 1dafebd5..b59a5c9d 100644 --- a/src/meshing.h +++ b/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 containing: * - New vertices as Rx3 matrix (float64) * - New faces as Sx3 matrix (int32) @@ -63,7 +68,10 @@ pmp_trimesh_remesh( Eigen::Ref 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)); /** diff --git a/tests/test_meshing.py b/tests/test_meshing.py index 60a55775..f1104588 100644 --- a/tests/test_meshing.py +++ b/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