Add geodesic isolines and mesh splitting functionality

Author: pv

docs/_images/example_geodesics.png

@@ -0,0 +1,102 @@
+from pathlib import Path
+
+import numpy as np
+from compas.colors import Color, ColorMap
+from compas.datastructures import Mesh
+from compas.geometry import Box, Point, Polyline, Translation
+from compas_viewer import Viewer
+from compas_viewer.config import Config
+
+from compas_cgal.geodesics import (
+    HeatGeodesicSolver,
+    geodesic_isolines,
+    geodesic_isolines_split,
+    heat_geodesic_distances,
+)
+
+# Load mesh
+FILE = Path(__file__).parent.parent.parent / "data" / "elephant.off"
+mesh = Mesh.from_off(FILE)
+mesh.quads_to_triangles()
+V, F = mesh.to_vertices_and_faces()
+V_np = np.array(V)
+
+# Config
+X_OFF, Y_OFF = 0.75, 1.0
+ISOVALUES = [i/50 for i in range(50)]
+SPLIT_ISOVALUES = [i/20 for i in range(20)]
+COLORS = [Color.red(), Color.orange(), Color.yellow(), Color.green(),
+          Color.cyan(), Color.blue(), Color.purple(), Color.magenta()]
+cmap = ColorMap.from_two_colors(Color.blue(), Color.red())
+
+
+def make_mesh(V, F, offset):
+    m = Mesh.from_vertices_and_faces(V, F)
+    m.transform(Translation.from_vector([offset[0], offset[1], 0]))
+    return m
+
+
+def make_vertex_colors(distances):
+    return {i: cmap(d, minval=distances.min(), maxval=distances.max()) for i, d in enumerate(distances)}
+
+
+def make_isolines(sources, offset):
+    polylines = []
+    for pts in geodesic_isolines((V, F), sources, ISOVALUES):
+        points = [[pts[i, 0] + offset[0], pts[i, 1] + offset[1], pts[i, 2]] for i in range(len(pts))]
+        polylines.append(Polyline(points))
+    return polylines
+
+
+def make_split_meshes(sources, offset):
+    meshes = []
+    for i, (v, f) in enumerate(geodesic_isolines_split((V, F), sources, SPLIT_ISOVALUES)):
+        m = Mesh.from_vertices_and_faces(v.tolist(), f.tolist())
+        m.transform(Translation.from_vector([offset[0], offset[1], 0]))
+        meshes.append((m, COLORS[i % len(COLORS)]))
+    return meshes
+
+
+def make_source_points(sources, offset):
+    return [Point(V[s][0] + offset[0], V[s][1] + offset[1], V[s][2]) for s in sources]
+
+
+# Row 1: Single source
+src1 = [0]
+dist1 = heat_geodesic_distances((V, F), src1)
+
+# Row 2: Multiple sources (bbox corners)
+solver = HeatGeodesicSolver((V, F))
+bbox = Box.from_points(V_np)
+src2 = list(dict.fromkeys([int(np.argmin(np.linalg.norm(V_np - c, axis=1))) for c in bbox.vertices]))
+dist2 = solver.solve(src2)
+
+# Viewer
+config = Config()
+config.camera.target = [X_OFF, -Y_OFF / 2, 0]
+config.camera.position = [X_OFF, -2.0, 0.8]
+viewer = Viewer(config=config)
+
+# Row 1: Single Source
+g1 = viewer.scene.add_group("Single Source")
+g1.add(make_mesh(V, F, (0, 0)), use_vertexcolors=True, vertexcolor=make_vertex_colors(dist1), show_lines=False)
+for pt in make_source_points(src1, (0, 0)):
+    g1.add(pt, pointcolor=Color.black(), pointsize=20)
+# g1.add(make_mesh(V, F, (X_OFF, 0)), facecolor=Color.grey(), show_lines=True)
+for pl in make_isolines(src1, (X_OFF, 0)):
+    g1.add(pl, linecolor=Color.red(), lineswidth=5)
+for m, c in make_split_meshes(src1, (2 * X_OFF, 0)):
+    g1.add(m, facecolor=c, show_lines=False)
+
+# Row 2: Multiple Sources
+g2 = viewer.scene.add_group("Multiple Sources")
+g2.add(make_mesh(V, F, (0, -Y_OFF)), use_vertexcolors=True, vertexcolor=make_vertex_colors(dist2), show_lines=False)
+for pt in make_source_points(src2, (0, -Y_OFF)):
+    g2.add(pt, pointcolor=Color.black(), pointsize=20)
+# g2.add(make_mesh(V, F, (X_OFF, -Y_OFF)), facecolor=Color.grey(), show_lines=True)
+for pl in make_isolines(src2, (X_OFF, -Y_OFF)):
+    g2.add(pl, linecolor=Color.red(), lineswidth=5)
+for m, c in make_split_meshes(src2, (2 * X_OFF, -Y_OFF)):
+    g2.add(m, facecolor=c, show_lines=False)
+
+viewer.show()

docs/examples/example_geodesics.py

@@ -0,0 +1,29 @@
+Geodesic Distances and Isolines
+===============================
+
+This example demonstrates geodesic distance computation, isoline extraction, and mesh splitting using COMPAS CGAL.
+
+The visualization shows two rows:
+
+* **Row 1 (Single Source)**: Geodesic distances from a single vertex
+* **Row 2 (Multiple Sources)**: Geodesic distances from all bounding box corners
+
+Each row displays three columns:
+
+* **Heat Field**: Mesh colored by geodesic distance (blue → red gradient) with source points marked in black
+* **Isolines**: Extracted isoline polylines at regular distance intervals
+* **Split Mesh**: Mesh split into components along the isolines, each colored differently
+
+Key Features:
+
+* Heat method geodesic distance computation via ``heat_geodesic_distances``
+* Reusable solver via ``HeatGeodesicSolver`` for multiple queries
+* Isoline extraction as polylines via ``geodesic_isolines``
+* Mesh splitting along isolines via ``geodesic_isolines_split``
+
+.. figure:: /_images/example_geodesics.png
+    :figclass: figure
+    :class: figure-img img-fluid
+
+.. literalinclude:: example_geodesics.py
+   :language: python

docs/examples/example_geodesics.rst

@@ -6,11 +6,15 @@
 from numpy.typing import NDArray
 
 from compas_cgal import _types_std  # noqa: F401  # Load vector type bindings
+from compas_cgal._geodesics import geodesic_isolines as _geodesic_isolines
+from compas_cgal._geodesics import geodesic_isolines_split as _geodesic_isolines_split
 from compas_cgal._geodesics import heat_geodesic_distances as _heat_geodesic_distances
 from compas_cgal._geodesics import HeatGeodesicSolver as _HeatGeodesicSolver
+from compas_cgal.types import PolylinesNumpy
 from compas_cgal.types import VerticesFaces
+from compas_cgal.types import VerticesFacesNumpy
 
-__all__ = ["heat_geodesic_distances", "HeatGeodesicSolver"]
+__all__ = ["heat_geodesic_distances", "HeatGeodesicSolver", "geodesic_isolines_split", "geodesic_isolines"]
 
 
 def heat_geodesic_distances(mesh: VerticesFaces, sources: List[int]) -> NDArray:
@@ -102,3 +106,78 @@ def solve(self, sources: List[int]) -> NDArray:
     def num_vertices(self) -> int:
         """Number of vertices in the mesh."""
         return self._solver.num_vertices
+
+
+def geodesic_isolines_split(
+    mesh: VerticesFaces,
+    sources: List[int],
+    isovalues: List[float],
+) -> List[VerticesFacesNumpy]:
+    """Split mesh into components along geodesic isolines.
+
+    Computes geodesic distances from sources, refines the mesh along
+    specified isovalue thresholds, and splits into connected components.
+
+    Parameters
+    ----------
+    mesh : :attr:`compas_cgal.types.VerticesFaces`
+        A triangulated mesh as a tuple of vertices and faces.
+    sources : List[int]
+        Source vertex indices for geodesic distance computation.
+    isovalues : List[float]
+        Isovalue thresholds for splitting. The mesh will be refined
+        along curves where the geodesic distance equals each isovalue,
+        then split into connected components.
+
+    Returns
+    -------
+    List[:attr:`compas_cgal.types.VerticesFacesNumpy`]
+        List of mesh components as (vertices, faces) tuples.
+
+    Examples
+    --------
+    >>> from compas.geometry import Sphere
+    >>> from compas_cgal.geodesics import geodesic_isolines_split
+    >>> sphere = Sphere(1.0)
+    >>> mesh = sphere.to_vertices_and_faces(u=32, v=32, triangulated=True)
+    >>> components = geodesic_isolines_split(mesh, [0], [0.5, 1.0, 1.5])
+    >>> len(components)  # Number of mesh strips
+
+    """
+    V, F = mesh
+    V = np.asarray(V, dtype=np.float64, order="C")
+    F = np.asarray(F, dtype=np.int32, order="C")
+
+    vertices_list, faces_list = _geodesic_isolines_split(V, F, sources, isovalues)
+    return list(zip(vertices_list, faces_list))
+
+
+def geodesic_isolines(
+    mesh: VerticesFaces,
+    sources: List[int],
+    isovalues: List[float],
+) -> PolylinesNumpy:
+    """Extract isoline polylines from geodesic distance field.
+
+    Computes geodesic distances and extracts polylines along specified isovalues.
+
+    Parameters
+    ----------
+    mesh : :attr:`compas_cgal.types.VerticesFaces`
+        A triangulated mesh as a tuple of vertices and faces.
+    sources : List[int]
+        Source vertex indices for geodesic distance computation.
+    isovalues : List[float]
+        Isovalue thresholds for isoline extraction.
+
+    Returns
+    -------
+    :attr:`compas_cgal.types.PolylinesNumpy`
+        List of polyline segments as Nx3 arrays of points.
+
+    """
+    V, F = mesh
+    V = np.asarray(V, dtype=np.float64, order="C")
+    F = np.asarray(F, dtype=np.int32, order="C")
+
+    return list(_geodesic_isolines(V, F, sources, isovalues))