ADD skeletonization and slicer.

Author: petrasvestartas

CMakeLists.txt

@@ -230,3 +230,6 @@ add_nanobind_module(meshing_ext src/meshing.cpp)
 add_nanobind_module(intersections_ext src/intersections.cpp)
 add_nanobind_module(measure_ext src/measure.cpp)
 add_nanobind_module(reconstruction_ext src/reconstruction.cpp)
+add_nanobind_module(skeletonization_ext src/skeletonization.cpp)
+add_nanobind_module(slicer_ext src/slicer.cpp)
+

data/elephant.off

@@ -0,0 +1,56 @@
+import math
+from pathlib import Path
+
+from compas.datastructures import Mesh
+from compas.geometry import Polyline
+from compas.geometry import Rotation
+from compas.geometry import Scale
+from compas.geometry import Translation
+from compas_viewer import Viewer
+from line_profiler import profile
+
+from compas_cgal.skeletonization import mesh_skeleton
+
+
+@profile
+def main():
+    """Skeletonize a mesh."""
+
+    input_file = Path(__file__).parent.parent.parent / "data" / "elephant.off"
+
+    rotation_x = Rotation.from_axis_and_angle([1, 0, 0], math.radians(60))
+    rotation_y = Rotation.from_axis_and_angle([0, 1, 0], math.radians(5))
+    rotation = rotation_y * rotation_x
+    scale = Scale.from_factors([5, 5, 5])
+    translation = Translation.from_vector([0, 0, 2])
+
+    mesh = Mesh.from_off(input_file).transformed(translation * rotation * scale)
+    # mesh = mesh.subdivided("loop")
+    v, f = mesh.to_vertices_and_faces(triangulated=True)
+
+    skeleton_edges = mesh_skeleton((v, f))
+
+    polylines = []
+    for start_point, end_point in skeleton_edges:
+        polyline = Polyline([start_point, end_point])
+        polylines.append(polyline)
+
+    return mesh, polylines
+
+
+mesh, polylines = main()
+
+# ==============================================================================
+# Visualize
+# ==============================================================================
+
+viewer = Viewer()
+viewer.renderer.camera.target = [0, 0, 1.5]
+viewer.renderer.camera.position = [-5, -5, 1.5]
+
+viewer.scene.add(mesh, opacity=0.5, show_points=False)
+
+for polyline in polylines:
+    viewer.scene.add(polyline, linewidth=5, show_points=True)
+
+viewer.show()

docs/examples/example_skeletonization.py

@@ -0,0 +1,18 @@
+Mesh Skeletonization
+====================
+
+This example demonstrates how to compute the geometric skeleton of a triangle mesh using COMPAS CGAL.
+
+Key Features:
+
+* Loading and transforming OFF mesh files
+* Mesh subdivision using Loop scheme
+* Skeleton computation using mean curvature flow
+* Visualization of mesh and skeleton polylines
+
+.. figure:: /_images/example_skeletonization.png
+    :figclass: figure
+    :class: figure-img img-fluid
+
+.. literalinclude:: example_skeletonization.py
+   :language: python

docs/examples/example_skeletonization.rst

@@ -0,0 +1,64 @@
+import math
+from pathlib import Path
+
+import numpy as np
+from compas.datastructures import Mesh
+from compas.geometry import Plane
+from compas.geometry import Point
+from compas.geometry import Polyline
+from compas.geometry import Vector
+from compas_viewer import Viewer
+from compas_viewer.config import Config
+from line_profiler import profile
+
+from compas_cgal.slicer import slice_mesh
+
+
+@profile
+def main():
+    # Get Mesh from STL
+    FILE = Path(__file__).parent.parent.parent / "data" / "3DBenchy.stl"
+    benchy = Mesh.from_stl(FILE)
+
+    V, F = benchy.to_vertices_and_faces()
+
+    # Get Slice planes from the bounding box
+    bbox = benchy.aabb()
+    normal = Vector(0, 0, 1)
+    planes = []
+    for i in np.linspace(bbox.zmin, bbox.zmax, 50):
+        plane = Plane(Point(0, 0, i), normal)
+        planes.append(plane)
+
+    # Slice
+    slicer_polylines = slice_mesh((V, F), planes)
+
+    # Convert edges to polylines
+    polylines = []
+    for polyline in slicer_polylines:
+        points = []
+        for point in polyline:
+            points.append(Point(*point))
+        polylines.append(Polyline(points))
+
+    return benchy, polylines
+
+
+mesh, polylines = main()
+
+# ==============================================================================
+# Visualize
+# ==============================================================================
+
+config = Config()
+config.camera.target = [0, 100, 0]
+config.camera.position = [0, -75, 50]
+config.camera.scale = 10
+
+viewer = Viewer(config=config)
+
+viewer.scene.add(mesh, opacity=0.5, show_lines=False, show_points=False)
+for polyline in polylines:
+    viewer.scene.add(polyline, linewidth=2, show_points=False)
+
+viewer.show()

docs/examples/example_slicer.py

@@ -0,0 +1,19 @@
+Mesh Slicing
+============
+
+This example demonstrates how to slice a mesh with multiple planes using COMPAS CGAL.
+
+Key Features:
+
+* Loading STL mesh files
+* Creating multiple slice planes from bounding box
+* Slicing mesh with multiple planes
+* Converting slice results to polylines
+* Visualization with semi-transparent mesh and slice curves
+
+.. figure:: /_images/example_slicer.png
+    :figclass: figure
+    :class: figure-img img-fluid
+
+.. literalinclude:: example_slicer.py
+   :language: python

docs/examples/example_slicer.rst

@@ -0,0 +1,58 @@
+import numpy as np
+from compas.plugins import plugin
+
+from compas_cgal import skeletonization_ext
+
+from .types import PolylinesNumpySkeleton
+from .types import VerticesFaces
+
+
+@plugin(category="mesh")
+def mesh_skeleton(mesh: VerticesFaces) -> PolylinesNumpySkeleton:
+    """Compute the geometric skeleton of a triangle mesh using mean curvature flow.
+
+    Parameters
+    ----------
+    mesh : VerticesFaces
+        A tuple containing:
+        * vertices: Nx3 array of vertex coordinates
+        * faces: Mx3 array of vertex indices
+
+    Returns
+    -------
+    PolylinesNumpySkeleton
+        List of polylines representing the skeleton edges.
+        Each polyline is a tuple of start and end point coordinates.
+
+    Raises
+    ------
+    TypeError
+        If the input mesh is not a tuple of vertices and faces.
+    ValueError
+        If the vertices array is not Nx3.
+        If the faces array is not Mx3.
+        If the face indices are out of range.
+        If the mesh is not manifold and closed.
+    RuntimeError
+        If the mesh contraction fails to converge.
+
+    Notes
+    -----
+    The input mesh must be manifold and closed.
+    The skeleton is computed using mean curvature flow.
+    """
+    V, F = mesh
+    V_numpy = np.asarray(V, dtype=np.float64, order="C")  # Ensure C-contiguous
+    F_numpy = np.asarray(F, dtype=np.int32, order="C")  # Ensure C-contiguous
+
+    # Get start and end points as flattened vectorS    
+    start_points, end_points = skeletonization_ext.mesh_skeleton(V_numpy, F_numpy)
+
+    # Convert flattened vectors to list of point coordinates
+    edges = []
+    for i in range(0, len(start_points), 3):
+        start = [start_points[i], start_points[i + 1], start_points[i + 2]]
+        end = [end_points[i], end_points[i + 1], end_points[i + 2]]
+        edges.append((start, end))
+
+    return edges

src/compas_cgal/skeletonization.py

@@ -0,0 +1,53 @@
+import numpy as np
+from compas.geometry import Plane
+from compas.plugins import plugin
+
+from compas_cgal import slicer_ext
+
+from .types import PolylinesNumpy
+from .types import VerticesFaces
+
+
+@plugin(category="trimesh", pluggable_name="trimesh_slice")
+def slice_mesh_planes(mesh: VerticesFaces, planes: list[Plane]) -> PolylinesNumpy:
+    """Slice a mesh by a list of planes.
+
+    Parameters
+    ----------
+    mesh : :attr:`compas_cgal.types.VerticesFaces`
+        The mesh to slice.
+    planes : :attr:`compas_cgal.types.Planes`
+        The slicing planes.
+
+    Returns
+    -------
+    :attr:`compas_cgal.types.PolylinesNumpy`
+        A list of slice polylines, with each polyline an array of points.
+
+    Examples
+    --------
+    >>> from compas.geometry import Sphere, Plane, Polyline
+    >>> from compas.itertools import linspace
+    >>> from compas_cgal.slicer import slice_mesh
+
+    >>> sphere = Sphere(1.0, point=[0, 0, 1.0])
+    >>> mesh = sphere.to_vertices_and_faces(u=32, v=32, triangulated=True)
+
+    >>> planes = [Plane([0, 0, z], [0, 0, 1.0]) for z in linspace(0.1, 1.9, 19)]
+
+    >>> result = slice_mesh(mesh, planes)
+    >>> polylines = [Polyline(points) for points in result]
+
+    """
+    vertices, faces = mesh
+    points, normals = zip(*planes)
+    V = np.asarray(vertices, dtype=np.float64, order="C")
+    F = np.asarray(faces, dtype=np.int32, order="C")
+    P = np.array(points, dtype=np.float64, order="C")
+    N = np.array(normals, dtype=np.float64, order="C")
+
+    pointsets = slicer_ext.slice_mesh(V, F, P, N)
+    return pointsets
+
+
+slice_mesh = slice_mesh_planes

src/compas_cgal/slicer.py

@@ -0,0 +1,60 @@
+#include "skeletonization.h"
+
+typedef CGAL::Mean_curvature_flow_skeletonization<compas::Mesh> Skeletonization;
+typedef Skeletonization::Skeleton Skeleton;
+typedef boost::graph_traits<Skeleton>::vertex_descriptor SkeletonVertex;
+typedef boost::graph_traits<Skeleton>::edge_descriptor SkeletonEdge;
+
+std::tuple<std::vector<double>, std::vector<double>>
+pmp_mesh_skeleton(
+    Eigen::Ref<const compas::RowMatrixXd> vertices,
+    Eigen::Ref<const compas::RowMatrixXi> faces)
+{
+    compas::Mesh mesh = compas::mesh_from_vertices_and_faces(vertices, faces);
+
+    Skeleton skeleton;
+    Skeletonization mcs(mesh);
+
+    mcs.contract_until_convergence();
+    mcs.convert_to_skeleton(skeleton);
+
+    // Initialize vectors to store start and end points
+    std::vector<double> start_points;
+    std::vector<double> end_points;
+
+    // Reserve space for efficiency
+    size_t num_edges = boost::num_edges(skeleton);
+    start_points.reserve(num_edges * 3);  // Each point has 3 coordinates
+    end_points.reserve(num_edges * 3);
+
+    // Extract skeleton edges
+    for (SkeletonEdge edge : CGAL::make_range(edges(skeleton))) {
+        const compas::Kernel::Point_3& start = skeleton[source(edge, skeleton)].point;
+        const compas::Kernel::Point_3& end = skeleton[target(edge, skeleton)].point;
+
+        // Add start point coordinates
+        start_points.push_back(start.x());
+        start_points.push_back(start.y());
+        start_points.push_back(start.z());
+
+        // Add end point coordinates
+        end_points.push_back(end.x());
+        end_points.push_back(end.y());
+        end_points.push_back(end.z());
+    }
+
+    return std::make_tuple(start_points, end_points);
+};
+
+NB_MODULE(skeletonization_ext, m) {
+
+    nb::bind_vector<std::vector<double>>(m, "VectorDouble");
+
+    m.def(
+        "mesh_skeleton",
+        &pmp_mesh_skeleton,
+        "Create a geometric skeleton from a mesh using mean curvature flow",
+        "vertices"_a,
+        "faces"_a
+    );
+}

src/skeletonization.cpp

@@ -0,0 +1,21 @@
+#pragma once
+
+#include "compas.h"
+
+// CGAL skeletonization
+#include <CGAL/Mean_curvature_flow_skeletonization.h>
+#include <CGAL/boost/graph/split_graph_into_polylines.h>
+
+/**
+ * @brief Compute the geometric skeleton of a triangle mesh using mean curvature flow.
+ * 
+ * @param vertices Matrix of vertex positions as Nx3 matrix in row-major order (float64)
+ * @param faces Matrix of face indices as Mx3 matrix in row-major order (int32)
+ * @return std::tuple<std::vector<double>, std::vector<double>> containing:
+ *         - Start points of skeleton edges as vector of 3D coordinates (float64)
+ *         - End points of skeleton edges as vector of 3D coordinates (float64)
+ */
+std::tuple<std::vector<double>, std::vector<double>>
+pmp_mesh_skeleton(
+    Eigen::Ref<const compas::RowMatrixXd> vertices,
+    Eigen::Ref<const compas::RowMatrixXi> faces);