ADD skeletonization with vertex mapping

Author: zygi

CHANGELOG.md

@@ -8,6 +8,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 ## Unreleased
 
 ### Added
+* Added  `compas_cgal.skeletonization.mesh_skeleton_with_mapping`.
 
 ### Changed
 

docs/_images/example_skeletonization_with_mapping.png

@@ -0,0 +1,54 @@
+import math
+from pathlib import Path
+
+from compas.datastructures import Mesh
+from compas.geometry import Polyline, Point
+from compas.geometry import Rotation, Scale, Translation
+from compas_viewer import Viewer
+
+from compas_cgal.skeletonization import mesh_skeleton_with_mapping
+
+# Load and transform mesh
+input_file = Path(__file__).parent.parent.parent / "data" / "elephant.off"
+
+rotation = Rotation.from_axis_and_angle([0, 1, 0], math.radians(5)) * Rotation.from_axis_and_angle([1, 0, 0], math.radians(60))
+transform = Translation.from_vector([0, 0, 2]) * rotation * Scale.from_factors([5, 5, 5])
+
+mesh = Mesh.from_off(input_file).transformed(transform)
+v, f = mesh.to_vertices_and_faces(triangulated=True)
+
+# Compute skeleton with vertex mapping
+skeleton_edges, vertex_indices = mesh_skeleton_with_mapping((v, f))
+
+# Create polylines for skeleton edges
+polylines = [Polyline([start, end]) for start, end in skeleton_edges]
+
+# Select edge to highlight (100th or last if fewer edges)
+edge_idx = min(100, len(vertex_indices) - 1)
+start_indices, end_indices = vertex_indices[edge_idx]
+
+print(f"Mesh: {len(v)} vertices, Skeleton: {len(skeleton_edges)} edges")
+print(f"Edge {edge_idx}: {len(start_indices)} vertices → start, {len(end_indices)} vertices → end")
+
+# Visualize
+viewer = Viewer()
+viewer.renderer.camera.target = [0, 0, 1.5]
+viewer.renderer.camera.position = [-5, -5, 1.5]
+
+# Show mesh as backdrop
+viewer.scene.add(mesh, opacity=0.2, show_points=False, facecolor=[0.9, 0.9, 0.9])
+
+# Show vertices that map to selected skeleton edge
+for idx in start_indices:
+    viewer.scene.add(Point(*v[idx]), pointcolor=[1.0, 0.0, 0.0], pointsize=15)  # red = start
+for idx in end_indices:
+    viewer.scene.add(Point(*v[idx]), pointcolor=[0.0, 0.0, 1.0], pointsize=15)  # blue = end
+
+# Show skeleton edges
+for i, polyline in enumerate(polylines):
+    if i == edge_idx:
+        viewer.scene.add(polyline, linewidth=10, linecolor=[1.0, 1.0, 0.0], show_points=True, pointsize=20)  # yellow
+    else:
+        viewer.scene.add(polyline, linewidth=3, linecolor=[0.5, 0.5, 0.5], show_points=False)  # gray
+
+viewer.show()

docs/examples/example_skeletonization_with_mapping.py

@@ -0,0 +1,18 @@
+Mesh Skeletonization with Vertex Mapping
+=========================================
+
+This example demonstrates how to compute the geometric skeleton of a triangle mesh with vertex correspondence mapping using COMPAS CGAL.
+
+Key Features:
+
+* Computing skeleton with vertex correspondence information
+* Understanding which original mesh vertices map to each skeleton vertex
+* Colored visualization showing vertex mapping
+* Printing mapping statistics to understand the correspondence
+
+.. figure:: /_images/example_skeletonization_with_mapping.png
+    :figclass: figure
+    :class: figure-img img-fluid
+
+.. literalinclude:: example_skeletonization_with_mapping.py
+   :language: python

docs/examples/example_skeletonization_with_mapping.rst

@@ -5,6 +5,7 @@
 from compas_cgal import _types_std  # noqa: F401
 
 from .types import PolylinesNumpySkeleton
+from .types import SkeletonVertexMapping
 from .types import VerticesFaces
 
 
@@ -46,8 +47,8 @@ def mesh_skeleton(mesh: VerticesFaces) -> PolylinesNumpySkeleton:
     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.mesh_skeleton(V_numpy, F_numpy)
+    # Get start and end points as flattened vectors
+    start_points, end_points, _, _ = _skeletonization.mesh_skeleton(V_numpy, F_numpy)
 
     # Convert flattened vectors to list of point coordinates
     edges = []
@@ -57,3 +58,69 @@ def mesh_skeleton(mesh: VerticesFaces) -> PolylinesNumpySkeleton:
         edges.append((start, end))
 
     return edges
+
+
+@plugin(category="mesh")
+def mesh_skeleton_with_mapping(mesh: VerticesFaces) -> tuple[PolylinesNumpySkeleton, SkeletonVertexMapping]:
+    """Compute the geometric skeleton of a triangle mesh with vertex correspondence mapping.
+
+    Parameters
+    ----------
+    mesh : VerticesFaces
+        A tuple containing:
+        * vertices: Nx3 array of vertex coordinates
+        * faces: Mx3 array of vertex indices
+
+    Returns
+    -------
+    PolylinesSkeletonWithMapping
+        A tuple containing:
+        * edges: List of polylines representing the skeleton edges.
+          Each polyline is a tuple of start and end point coordinates.
+        * vertex_indices: List of tuples, each containing two lists of
+          vertex indices corresponding to the start and end vertices of
+          each skeleton edge. These are the original mesh vertices that
+          contracted to form each skeleton vertex.
+
+    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.
+    Each skeleton vertex corresponds to a set of original mesh vertices
+    that were contracted to that point during the skeletonization process.
+    (The set might be empty for some skeleton vertices that don't correspond
+    to any original vertex.)
+    """
+    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 and vertex indices
+    start_points, end_points, start_vertex_indices, end_vertex_indices = _skeletonization.mesh_skeleton(V_numpy, F_numpy)
+
+    # Convert flattened vectors to list of point coordinates
+    edges = []
+    vertex_indices = []
+
+    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))
+
+    # Process vertex indices - convert VectorInt to Python lists
+    for start_indices, end_indices in zip(start_vertex_indices, end_vertex_indices):
+        vertex_indices.append((list(start_indices), list(end_indices)))
+
+    return edges, vertex_indices

src/compas_cgal/skeletonization.py

@@ -42,6 +42,11 @@
 PolylinesNumpySkeleton = List[Tuple[List[float], List[float]]]
 """A list of polylines, where each polyline is represented by a tuple of start and end point coordinates."""
 
+SkeletonVertexMapping = List[Tuple[List[int], List[int]]]
+"""Vertex correspondence mapping for skeleton edges. Each tuple contains two lists:
+the first list has indices of original mesh vertices that contracted to the skeleton edge's start vertex,
+the second list has indices that contracted to the skeleton edge's end vertex."""
+
 Planes = Union[
     Sequence[compas.geometry.Plane],
     Sequence[Tuple[compas.geometry.Point, compas.geometry.Vector]],

src/compas_cgal/types.py

@@ -5,7 +5,7 @@ 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>>
+std::tuple<std::vector<double>, std::vector<double>, std::vector<std::vector<int>>, std::vector<std::vector<int>>>
 pmp_mesh_skeleton(
     Eigen::Ref<const compas::RowMatrixXd> vertices,
     Eigen::Ref<const compas::RowMatrixXi> faces)
@@ -15,22 +15,35 @@ pmp_mesh_skeleton(
     Skeleton skeleton;
     Skeletonization mcs(mesh);
 
+    // Create a vertex index map for the mesh
+    std::map<compas::Mesh::Vertex_index, int> vertex_index_map;
+    for (auto v : mesh.vertices()) {
+        vertex_index_map[v] = v.idx();
+    }
+
     mcs.contract_until_convergence();
     mcs.convert_to_skeleton(skeleton);
 
-    // Initialize vectors to store start and end points
+    // Initialize vectors to store start and end points and vertex indices
     std::vector<double> start_points;
     std::vector<double> end_points;
+    std::vector<std::vector<int>> start_vertex_indices;
+    std::vector<std::vector<int>> end_vertex_indices;
 
     // 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);
+    start_vertex_indices.reserve(num_edges);
+    end_vertex_indices.reserve(num_edges);
 
     // 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;
+        SkeletonVertex source_vertex = source(edge, skeleton);
+        SkeletonVertex target_vertex = target(edge, skeleton);
+
+        const compas::Kernel::Point_3& start = skeleton[source_vertex].point;
+        const compas::Kernel::Point_3& end = skeleton[target_vertex].point;
 
         // Add start point coordinates
         start_points.push_back(start.x());
@@ -41,9 +54,23 @@ pmp_mesh_skeleton(
         end_points.push_back(end.x());
         end_points.push_back(end.y());
         end_points.push_back(end.z());
+
+        // Extract vertex indices for start vertex
+        std::vector<int> start_indices;
+        for (auto v : skeleton[source_vertex].vertices) {
+            start_indices.push_back(vertex_index_map[v]);
+        }
+        start_vertex_indices.push_back(start_indices);
+
+        // Extract vertex indices for end vertex
+        std::vector<int> end_indices;
+        for (auto v : skeleton[target_vertex].vertices) {
+            end_indices.push_back(vertex_index_map[v]);
+        }
+        end_vertex_indices.push_back(end_indices);
     }
 
-    return std::make_tuple(start_points, end_points);
+    return std::make_tuple(start_points, end_points, start_vertex_indices, end_vertex_indices);
 };
 
 NB_MODULE(_skeletonization, m) {

src/skeletonization.cpp

@@ -8,14 +8,16 @@
 
 /**
  * @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:
+ * @return std::tuple<std::vector<double>, std::vector<double>, std::vector<std::vector<int>>, std::vector<std::vector<int>>> containing:
  *         - Start points of skeleton edges as vector of 3D coordinates (float64)
  *         - End points of skeleton edges as vector of 3D coordinates (float64)
+ *         - Start vertex indices: for each skeleton edge, indices of original mesh vertices that contracted to the start vertex
+ *         - End vertex indices: for each skeleton edge, indices of original mesh vertices that contracted to the end vertex
  */
-std::tuple<std::vector<double>, std::vector<double>>
+std::tuple<std::vector<double>, std::vector<double>, std::vector<std::vector<int>>, std::vector<std::vector<int>>>
 pmp_mesh_skeleton(
     Eigen::Ref<const compas::RowMatrixXd> vertices,
     Eigen::Ref<const compas::RowMatrixXi> faces);

src/skeletonization.h

@@ -1,5 +1,5 @@
 from compas.geometry import Box
-from compas_cgal.skeletonization import mesh_skeleton
+from compas_cgal.skeletonization import mesh_skeleton, mesh_skeleton_with_mapping
 
 
 def test_mesh_skeleton():
@@ -17,3 +17,42 @@ def test_mesh_skeleton():
     assert isinstance(edges[0], tuple)
     assert len(edges[0]) == 2
     assert len(edges[0][0]) == 3  # 3D points
+
+
+def test_mesh_skeleton_with_mapping():
+    """Test mesh skeletonization with vertex mapping."""
+    # Create test box mesh
+    box = Box.from_width_height_depth(2.0, 2.0, 2.0)
+    mesh = box.to_vertices_and_faces(triangulated=True)
+    vertices, faces = mesh
+
+    # Get skeleton with mapping
+    edges, vertex_indices = mesh_skeleton_with_mapping(mesh)
+
+    # Basic validation of edges
+    assert isinstance(edges, list)
+    assert len(edges) > 0
+    assert isinstance(edges[0], tuple)
+    assert len(edges[0]) == 2
+    assert len(edges[0][0]) == 3  # 3D points
+
+    # Validate vertex indices
+    assert isinstance(vertex_indices, list)
+    assert len(vertex_indices) == len(edges)
+
+    # Each edge should have vertex indices for start and end
+    for indices in vertex_indices:
+        assert isinstance(indices, tuple)
+        assert len(indices) == 2
+        start_indices, end_indices = indices
+
+        # Each skeleton vertex maps to one or more original vertices
+        assert isinstance(start_indices, list)
+        assert isinstance(end_indices, list)
+
+        # Indices should be valid
+        num_vertices = len(vertices)
+        for idx in start_indices:
+            assert 0 <= idx < num_vertices
+        for idx in end_indices:
+            assert 0 <= idx < num_vertices