Merge pull request #58 from jf---/jf/geodesics-module

add heat method geodesics module

Author: tomvanmele

CHANGELOG.md

@@ -8,7 +8,11 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 ## Unreleased
 
 ### Added
-* Added  `compas_cgal.skeletonization.mesh_skeleton_with_mapping`.
+* Added  `compas_cgal.skeletonization.mesh_skeleton_with_mapping`
+- `HeatGeodesicSolver` class with precomputation for repeated queries
+- `heat_geodesic_distances` function for single-shot usage
+    - uses CGAL Heat_method_3 with intrinsic Delaunay triangulation
+    - ~30% faster than `libigl` heat in `compas_slicer` workflow
 
 - Added `simplify_polyline` and `simplify_polylines` functions for polyline simplification using Douglas-Peucker algorithm
 - Added `closest_points_on_polyline` function for batch closest point queries on polylines

CMakeLists.txt

@@ -238,3 +238,4 @@ add_nanobind_module(_triangulation src/triangulation.cpp)
 add_nanobind_module(_subdivision src/subdivision.cpp)
 add_nanobind_module(_polylines src/polylines.cpp)
 
+add_nanobind_module(_geodesics src/geodesics.cpp)

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

@@ -0,0 +1,183 @@
+"""Geodesic distance computation using CGAL heat method."""
+
+from typing import List
+
+import numpy as np
+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", "geodesic_isolines_split", "geodesic_isolines"]
+
+
+def heat_geodesic_distances(mesh: VerticesFaces, sources: List[int]) -> NDArray:
+    """Compute geodesic distances from source vertices using CGAL heat method.
+
+    Uses CGAL's Heat_method_3 with intrinsic Delaunay triangulation for
+    accurate geodesic distance computation.
+
+    Parameters
+    ----------
+    mesh : :attr:`compas_cgal.types.VerticesFaces`
+        A triangulated mesh as a tuple of vertices and faces.
+    sources : List[int]
+        Source vertex indices.
+
+    Returns
+    -------
+    NDArray
+        Geodesic distances from the nearest source to each vertex.
+        Shape is (n_vertices,).
+
+    Examples
+    --------
+    >>> from compas.geometry import Box
+    >>> from compas_cgal.geodesics import heat_geodesic_distances
+    >>> box = Box(1)
+    >>> mesh = box.to_vertices_and_faces(triangulated=True)
+    >>> distances = heat_geodesic_distances(mesh, [0])  # distances from vertex 0
+
+    """
+    V, F = mesh
+    V = np.asarray(V, dtype=np.float64, order="C")
+    F = np.asarray(F, dtype=np.int32, order="C")
+
+    result = _heat_geodesic_distances(V, F, sources)
+    return result.flatten()
+
+
+class HeatGeodesicSolver:
+    """Precomputed heat method solver for repeated geodesic queries.
+
+    Use this class when computing geodesic distances from multiple
+    different sources on the same mesh. The expensive precomputation
+    is done once in the constructor, and solve() can be called many
+    times efficiently.
+
+    Parameters
+    ----------
+    mesh : :attr:`compas_cgal.types.VerticesFaces`
+        A triangulated mesh as a tuple of vertices and faces.
+
+    Examples
+    --------
+    >>> from compas.geometry import Sphere
+    >>> from compas_cgal.geodesics import HeatGeodesicSolver
+    >>> sphere = Sphere(1.0)
+    >>> mesh = sphere.to_vertices_and_faces(u=32, v=32, triangulated=True)
+    >>> solver = HeatGeodesicSolver(mesh)  # precomputation happens here
+    >>> d0 = solver.solve([0])  # distances from vertex 0
+    >>> d1 = solver.solve([1])  # distances from vertex 1 (fast, reuses precomputation)
+
+    """
+
+    def __init__(self, mesh: VerticesFaces) -> None:
+        V, F = mesh
+        V = np.asarray(V, dtype=np.float64, order="C")
+        F = np.asarray(F, dtype=np.int32, order="C")
+        self._solver = _HeatGeodesicSolver(V, F)
+
+    def solve(self, sources: List[int]) -> NDArray:
+        """Compute geodesic distances from source vertices.
+
+        Parameters
+        ----------
+        sources : List[int]
+            Source vertex indices.
+
+        Returns
+        -------
+        NDArray
+            Geodesic distances from the nearest source to each vertex.
+            Shape is (n_vertices,).
+
+        """
+        result = self._solver.solve(sources)
+        return result.flatten()
+
+    @property
+    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))