add heat method geodesics module

CMakeLists.txt

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

src/compas_cgal/geodesics.py

@@ -0,0 +1,108 @@
+"""Geodesic distance computation using CGAL heat method."""
+
+from __future__ import annotations
+
+from typing import TYPE_CHECKING
+
+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 heat_geodesic_distances as _heat_geodesic_distances
+from compas_cgal._geodesics import HeatGeodesicSolver as _HeatGeodesicSolver
+
+if TYPE_CHECKING:
+    from compas.datastructures import Mesh
+
+__all__ = ["heat_geodesic_distances", "HeatGeodesicSolver"]
+
+
+def heat_geodesic_distances(mesh: Mesh, 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 : Mesh
+        A compas mesh (must be triangulated).
+    sources : list[int]
+        Source vertex indices.
+
+    Returns
+    -------
+    NDArray
+        Geodesic distances from the nearest source to each vertex.
+        Shape is (n_vertices,).
+
+    Examples
+    --------
+    >>> from compas.datastructures import Mesh
+    >>> from compas_cgal.geodesics import heat_geodesic_distances
+    >>> mesh = Mesh.from_obj('model.obj')
+    >>> distances = heat_geodesic_distances(mesh, [0])  # distances from vertex 0
+
+    """
+    V, F = mesh.to_vertices_and_faces()
+    vertices = np.asarray(V, dtype=np.float64)
+    faces = np.asarray(F, dtype=np.int32)
+
+    result = _heat_geodesic_distances(vertices, faces, sources)
+
+    # Return as 1D array
+    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 : Mesh
+        A compas mesh (must be triangulated).
+
+    Examples
+    --------
+    >>> from compas.datastructures import Mesh
+    >>> from compas_cgal.geodesics import HeatGeodesicSolver
+    >>> mesh = Mesh.from_obj('model.obj')
+    >>> 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: Mesh) -> None:
+        V, F = mesh.to_vertices_and_faces()
+        vertices = np.asarray(V, dtype=np.float64)
+        faces = np.asarray(F, dtype=np.int32)
+        self._solver = _HeatGeodesicSolver(vertices, faces)
+
+    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

src/geodesics.cpp

@@ -0,0 +1,147 @@
+#include "geodesics.h"
+
+#include <CGAL/Simple_cartesian.h>
+#include <CGAL/Surface_mesh.h>
+#include <CGAL/Heat_method_3/Surface_mesh_geodesic_distances_3.h>
+
+// Type definitions for heat method
+using HeatMethod = CGAL::Heat_method_3::Surface_mesh_geodesic_distances_3<compas::Mesh>;
+using vertex_descriptor = boost::graph_traits<compas::Mesh>::vertex_descriptor;
+
+
+compas::RowMatrixXd
+pmp_heat_geodesic_distances(
+    Eigen::Ref<compas::RowMatrixXd> vertices,
+    Eigen::Ref<compas::RowMatrixXi> faces,
+    const std::vector<int>& sources)
+{
+    // Convert input to CGAL mesh
+    compas::Mesh mesh = compas::mesh_from_vertices_and_faces(vertices, faces);
+
+    int n_vertices = mesh.number_of_vertices();
+
+    // Create property map for distances
+    auto distance_pmap = mesh.add_property_map<vertex_descriptor, double>(
+        "v:distance", 0.0).first;
+
+    // Create heat method solver (precomputation happens here)
+    HeatMethod hm(mesh);
+
+    // Add all sources
+    for (int src : sources) {
+        if (src >= 0 && src < n_vertices) {
+            vertex_descriptor vd = vertex_descriptor(src);
+            hm.add_source(vd);
+        }
+    }
+
+    // Compute geodesic distances
+    hm.estimate_geodesic_distances(distance_pmap);
+
+    // Copy results to output matrix
+    compas::RowMatrixXd result(n_vertices, 1);
+    for (vertex_descriptor vd : mesh.vertices()) {
+        result(vd.idx(), 0) = get(distance_pmap, vd);
+    }
+
+    return result;
+}
+
+
+/**
+ * @brief Heat method solver with precomputation for repeated queries.
+ *
+ * This class stores the mesh and precomputed data to allow efficient
+ * repeated geodesic distance computations from different source vertices.
+ */
+class HeatGeodesicSolver {
+public:
+    HeatGeodesicSolver(
+        Eigen::Ref<compas::RowMatrixXd> vertices,
+        Eigen::Ref<compas::RowMatrixXi> faces)
+        : mesh_(compas::mesh_from_vertices_and_faces(vertices, faces)),
+          n_vertices_(mesh_.number_of_vertices()),
+          distance_pmap_(mesh_.add_property_map<vertex_descriptor, double>("v:distance", 0.0).first),
+          hm_(mesh_)
+    {}
+
+    compas::RowMatrixXd solve(const std::vector<int>& sources) {
+        // Clear previous sources
+        hm_.clear_sources();
+
+        // Add all sources
+        for (int src : sources) {
+            if (src >= 0 && src < n_vertices_) {
+                vertex_descriptor vd = vertex_descriptor(src);
+                hm_.add_source(vd);
+            }
+        }
+
+        // Compute geodesic distances
+        hm_.estimate_geodesic_distances(distance_pmap_);
+
+        // Copy results to output matrix
+        compas::RowMatrixXd result(n_vertices_, 1);
+        for (vertex_descriptor vd : mesh_.vertices()) {
+            result(vd.idx(), 0) = get(distance_pmap_, vd);
+        }
+
+        return result;
+    }
+
+    int num_vertices() const { return n_vertices_; }
+
+private:
+    compas::Mesh mesh_;
+    int n_vertices_;
+    compas::Mesh::Property_map<vertex_descriptor, double> distance_pmap_;
+    HeatMethod hm_;
+};
+
+
+NB_MODULE(_geodesics, m) {
+    m.def(
+        "heat_geodesic_distances",
+        &pmp_heat_geodesic_distances,
+        "Compute geodesic distances using CGAL heat method.\n\n"
+        "Parameters\n"
+        "----------\n"
+        "vertices : array-like\n"
+        "    Vertex positions as Nx3 matrix (float64)\n"
+        "faces : array-like\n"
+        "    Face indices as Mx3 matrix (int32)\n"
+        "sources : list[int]\n"
+        "    Source vertex indices\n"
+        "\n"
+        "Returns\n"
+        "-------\n"
+        "array\n"
+        "    Geodesic distances from sources to all vertices (Nx1)",
+        "vertices"_a,
+        "faces"_a,
+        "sources"_a);
+
+    nanobind::class_<HeatGeodesicSolver>(m, "HeatGeodesicSolver",
+        "Precomputed heat method solver for repeated geodesic queries.\n\n"
+        "Use this class when computing geodesic distances from multiple\n"
+        "different sources on the same mesh. The expensive precomputation\n"
+        "is done once in the constructor, and solve() can be called many\n"
+        "times efficiently.")
+        .def(nanobind::init<Eigen::Ref<compas::RowMatrixXd>, Eigen::Ref<compas::RowMatrixXi>>(),
+             "vertices"_a, "faces"_a,
+             "Create solver with precomputation for the given mesh.")
+        .def("solve", &HeatGeodesicSolver::solve,
+             "sources"_a,
+             "Compute geodesic distances from source vertices.\n\n"
+             "Parameters\n"
+             "----------\n"
+             "sources : list[int]\n"
+             "    Source vertex indices\n"
+             "\n"
+             "Returns\n"
+             "-------\n"
+             "array\n"
+             "    Geodesic distances (Nx1)")
+        .def_prop_ro("num_vertices", &HeatGeodesicSolver::num_vertices,
+             "Number of vertices in the mesh.");
+}

src/geodesics.h

@@ -0,0 +1,22 @@
+#pragma once
+
+#include "compas.h"
+
+// CGAL Heat Method for geodesic distances
+#include <CGAL/Heat_method_3/Surface_mesh_geodesic_distances_3.h>
+
+/**
+ * @brief Compute geodesic distances from source vertices using CGAL heat method.
+ *
+ * @param vertices Matrix of vertex positions as Nx3 matrix (float64)
+ * @param faces Matrix of face indices as Mx3 matrix (int32)
+ * @param sources Vector of source vertex indices
+ * @return RowMatrixXd containing distances from sources to all vertices (Nx1)
+ */
+compas::RowMatrixXd
+pmp_heat_geodesic_distances(
+    Eigen::Ref<compas::RowMatrixXd> vertices,
+    Eigen::Ref<compas::RowMatrixXi> faces,
+    const std::vector<int>& sources);
+
+// HeatGeodesicSolver class is defined in geodesics.cpp and exposed via nanobind