diff --git a/CHANGELOG.md b/CHANGELOG.md index 80a79d8..672df69 100644 --- a/CHANGELOG.md +++ b/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 diff --git a/CMakeLists.txt b/CMakeLists.txt index 3ffac27..3a2d3d8 100644 --- a/CMakeLists.txt +++ b/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) diff --git a/docs/_images/example_geodesics.png b/docs/_images/example_geodesics.png new file mode 100644 index 0000000..82751cb Binary files /dev/null and b/docs/_images/example_geodesics.png differ diff --git a/docs/examples/example_geodesics.py b/docs/examples/example_geodesics.py new file mode 100644 index 0000000..f191632 --- /dev/null +++ b/docs/examples/example_geodesics.py @@ -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() diff --git a/docs/examples/example_geodesics.rst b/docs/examples/example_geodesics.rst new file mode 100644 index 0000000..adc578f --- /dev/null +++ b/docs/examples/example_geodesics.rst @@ -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 diff --git a/src/compas_cgal/geodesics.py b/src/compas_cgal/geodesics.py new file mode 100644 index 0000000..dcab545 --- /dev/null +++ b/src/compas_cgal/geodesics.py @@ -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)) diff --git a/src/geodesics.cpp b/src/geodesics.cpp new file mode 100644 index 0000000..d9149b4 --- /dev/null +++ b/src/geodesics.cpp @@ -0,0 +1,317 @@ +#include "geodesics.h" + +#include +#include +#include +#include +#include + +// Type definitions for heat method +using HeatMethod = CGAL::Heat_method_3::Surface_mesh_geodesic_distances_3; +using vertex_descriptor = boost::graph_traits::vertex_descriptor; + + +compas::RowMatrixXd +pmp_heat_geodesic_distances( + Eigen::Ref vertices, + Eigen::Ref faces, + const std::vector& 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( + "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 vertices, + Eigen::Ref faces) + : mesh_(compas::mesh_from_vertices_and_faces(vertices, faces)), + n_vertices_(mesh_.number_of_vertices()), + distance_pmap_(mesh_.add_property_map("v:distance", 0.0).first), + hm_(mesh_) + {} + + compas::RowMatrixXd solve(const std::vector& 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 distance_pmap_; + HeatMethod hm_; +}; + + +std::tuple, std::vector> +pmp_geodesic_isolines_split( + Eigen::Ref vertices, + Eigen::Ref faces, + const std::vector& sources, + const std::vector& isovalues) +{ + using vertex_descriptor = boost::graph_traits::vertex_descriptor; + using edge_descriptor = boost::graph_traits::edge_descriptor; + + // 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("v:distance", 0.0).first; + + // Create heat method solver and compute geodesic distances + HeatMethod hm(mesh); + for (int src : sources) { + if (src >= 0 && src < n_vertices) { + hm.add_source(vertex_descriptor(src)); + } + } + hm.estimate_geodesic_distances(distance_pmap); + + // Create edge constraint map to mark isoline edges + auto ecm = mesh.add_property_map("e:is_constrained", false).first; + + // Refine mesh along each isovalue + for (double isovalue : isovalues) { + CGAL::Polygon_mesh_processing::refine_mesh_at_isolevel( + mesh, distance_pmap, isovalue, + CGAL::parameters::edge_is_constrained_map(ecm)); + } + + // Triangulate any polygon faces created by refinement + // (faces with >2 vertices on isoline are not split by refine_mesh_at_isolevel) + CGAL::Polygon_mesh_processing::triangulate_faces(mesh); + + // Split the mesh into connected components bounded by isocurves + std::vector components; + CGAL::Polygon_mesh_processing::split_connected_components( + mesh, components, + CGAL::parameters::edge_is_constrained_map(ecm)); + + // Convert each component back to vertices/faces with proper index remapping + // CGAL mesh indices may have gaps, so we need to create contiguous indices + std::vector all_vertices; + std::vector all_faces; + all_vertices.reserve(components.size()); + all_faces.reserve(components.size()); + + for (const auto& comp : components) { + std::size_t nv = comp.number_of_vertices(); + std::size_t nf = comp.number_of_faces(); + + compas::RowMatrixXd V(nv, 3); + compas::RowMatrixXi F(nf, 3); + + // Build vertex index remapping (CGAL index -> contiguous index) + std::unordered_map vertex_remap; + int v_idx = 0; + for (auto vd : comp.vertices()) { + vertex_remap[vd.idx()] = v_idx; + auto pt = comp.point(vd); + V(v_idx, 0) = CGAL::to_double(pt.x()); + V(v_idx, 1) = CGAL::to_double(pt.y()); + V(v_idx, 2) = CGAL::to_double(pt.z()); + v_idx++; + } + + // Build faces with remapped vertex indices + int f_idx = 0; + for (auto fd : comp.faces()) { + int i = 0; + for (auto vd : vertices_around_face(comp.halfedge(fd), comp)) { + F(f_idx, i) = vertex_remap[vd.idx()]; + i++; + } + f_idx++; + } + + all_vertices.push_back(std::move(V)); + all_faces.push_back(std::move(F)); + } + + return std::make_tuple(std::move(all_vertices), std::move(all_faces)); +} + + +std::vector +pmp_geodesic_isolines( + Eigen::Ref vertices, + Eigen::Ref faces, + const std::vector& sources, + const std::vector& isovalues) +{ + using vertex_descriptor = boost::graph_traits::vertex_descriptor; + using edge_descriptor = boost::graph_traits::edge_descriptor; + + compas::Mesh mesh = compas::mesh_from_vertices_and_faces(vertices, faces); + int n_vertices = mesh.number_of_vertices(); + + auto distance_pmap = mesh.add_property_map("v:distance", 0.0).first; + + HeatMethod hm(mesh); + for (int src : sources) { + if (src >= 0 && src < n_vertices) { + hm.add_source(vertex_descriptor(src)); + } + } + hm.estimate_geodesic_distances(distance_pmap); + + auto ecm = mesh.add_property_map("e:is_constrained", false).first; + + for (double isovalue : isovalues) { + CGAL::Polygon_mesh_processing::refine_mesh_at_isolevel( + mesh, distance_pmap, isovalue, + CGAL::parameters::edge_is_constrained_map(ecm)); + } + + // Build a graph from constrained edges for polyline extraction + typedef boost::adjacency_list IsolineGraph; + typedef boost::graph_traits::vertex_descriptor ig_vertex; + + IsolineGraph iso_graph; + std::map v_map; + + for (auto e : mesh.edges()) { + if (get(ecm, e)) { + auto h = mesh.halfedge(e); + auto src = mesh.source(h); + auto tgt = mesh.target(h); + + if (v_map.find(src) == v_map.end()) { + v_map[src] = boost::add_vertex(src, iso_graph); + } + if (v_map.find(tgt) == v_map.end()) { + v_map[tgt] = boost::add_vertex(tgt, iso_graph); + } + boost::add_edge(v_map[src], v_map[tgt], iso_graph); + } + } + + // Visitor to collect polylines + struct PolylineVisitor { + const compas::Mesh& mesh; + std::vector>& polylines; + std::vector current; + + PolylineVisitor(const compas::Mesh& m, std::vector>& p) + : mesh(m), polylines(p) {} + + void start_new_polyline() { current.clear(); } + void add_node(ig_vertex v) { + vertex_descriptor vd = boost::get(boost::vertex_bundle, *graph_ptr)[v]; + current.push_back(mesh.point(vd)); + } + void end_polyline() { + if (current.size() >= 2) polylines.push_back(current); + } + const IsolineGraph* graph_ptr; + }; + + std::vector> polyline_points; + PolylineVisitor visitor(mesh, polyline_points); + visitor.graph_ptr = &iso_graph; + + CGAL::split_graph_into_polylines(iso_graph, visitor, CGAL::internal::IsTerminalDefault()); + + // Convert to output format + std::vector polylines; + for (const auto& pl : polyline_points) { + compas::RowMatrixXd pts(pl.size(), 3); + for (std::size_t i = 0; i < pl.size(); ++i) { + pts(i, 0) = CGAL::to_double(pl[i].x()); + pts(i, 1) = CGAL::to_double(pl[i].y()); + pts(i, 2) = CGAL::to_double(pl[i].z()); + } + polylines.push_back(std::move(pts)); + } + + return polylines; +} + + +NB_MODULE(_geodesics, m) { + m.def( + "heat_geodesic_distances", + &pmp_heat_geodesic_distances, + "Compute geodesic distances using CGAL heat method.", + "vertices"_a, "faces"_a, "sources"_a); + + nanobind::class_(m, "HeatGeodesicSolver", + "Precomputed heat method solver for repeated geodesic queries.") + .def(nanobind::init, Eigen::Ref>(), + "vertices"_a, "faces"_a) + .def("solve", &HeatGeodesicSolver::solve, "sources"_a) + .def_prop_ro("num_vertices", &HeatGeodesicSolver::num_vertices); + + m.def( + "geodesic_isolines_split", + &pmp_geodesic_isolines_split, + "Split mesh into components along geodesic isolines.", + "vertices"_a, "faces"_a, "sources"_a, "isovalues"_a); + + m.def( + "geodesic_isolines", + &pmp_geodesic_isolines, + "Extract isoline polylines from geodesic distance field.", + "vertices"_a, "faces"_a, "sources"_a, "isovalues"_a); +} diff --git a/src/geodesics.h b/src/geodesics.h new file mode 100644 index 0000000..3de7b3f --- /dev/null +++ b/src/geodesics.h @@ -0,0 +1,61 @@ +#pragma once + +#include "compas.h" + +// CGAL Heat Method for geodesic distances +#include +#include +#include + +/** + * @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 vertices, + Eigen::Ref faces, + const std::vector& sources); + +/** + * @brief Split mesh into components along geodesic isolines. + * + * Computes geodesic distances from sources, then refines the mesh along + * specified isolevel values and splits into connected components. + * + * @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 + * @param isovalues Vector of isovalue thresholds for splitting + * @return Tuple of (list of vertices arrays, list of faces arrays) for each mesh component + */ +std::tuple, std::vector> +pmp_geodesic_isolines_split( + Eigen::Ref vertices, + Eigen::Ref faces, + const std::vector& sources, + const std::vector& isovalues); + +/** + * @brief Extract isoline polylines from geodesic distance field. + * + * Computes geodesic distances and extracts polylines along specified isovalues. + * + * @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 + * @param isovalues Vector of isovalue thresholds + * @return Vector of polylines, each as Nx3 matrix of points + */ +std::vector +pmp_geodesic_isolines( + Eigen::Ref vertices, + Eigen::Ref faces, + const std::vector& sources, + const std::vector& isovalues); + +// HeatGeodesicSolver class is defined in geodesics.cpp and exposed via nanobind diff --git a/src/types_std.cpp b/src/types_std.cpp index f4c786b..97e5b78 100644 --- a/src/types_std.cpp +++ b/src/types_std.cpp @@ -6,6 +6,7 @@ NB_MODULE(_types_std, m) { nb::bind_vector>(m, "VectorInt"); nb::bind_vector>>(m, "VectorVectorInt"); nb::bind_vector>(m, "VectorRowMatrixXd"); + nb::bind_vector>(m, "VectorRowMatrixXi"); nb::bind_vector>>(m, "VectorVectorRowMatrixXd"); } \ No newline at end of file diff --git a/tests/test_geodesics.py b/tests/test_geodesics.py new file mode 100644 index 0000000..f65fa6e --- /dev/null +++ b/tests/test_geodesics.py @@ -0,0 +1,177 @@ +import pytest +import numpy as np + +from compas.geometry import Sphere +from compas.datastructures import Mesh +from compas_cgal.geodesics import geodesic_isolines +from compas_cgal.geodesics import geodesic_isolines_split +from compas_cgal.geodesics import heat_geodesic_distances +from compas_cgal.geodesics import HeatGeodesicSolver + + +@pytest.fixture +def sphere_mesh(): + """Create a test sphere mesh.""" + sphere = Sphere(1.0, point=[0, 0, 0]) + V, F = sphere.to_vertices_and_faces(u=16, v=16, triangulated=True) + return V, F + + +def test_heat_geodesic_distances_basic(sphere_mesh): + """Test basic geodesic distance computation.""" + V, F = sphere_mesh + sources = [0] + + distances = heat_geodesic_distances((V, F), sources) + + # Basic validation + assert isinstance(distances, np.ndarray) + assert distances.shape[0] == len(V) # One distance per vertex + assert distances.ndim == 1 # Should be 1D array + assert distances[0] == pytest.approx(0.0) # Distance from source to itself is 0 + assert np.all(distances >= 0) # All distances should be non-negative + assert np.all(np.isfinite(distances)) # No inf or nan values + + +def test_heat_geodesic_distances_multiple_sources(sphere_mesh): + """Test geodesic distance computation with multiple sources.""" + V, F = sphere_mesh + sources = [0, 10, 20] + + distances = heat_geodesic_distances((V, F), sources) + + # Validate shape and properties + assert isinstance(distances, np.ndarray) + assert distances.shape[0] == len(V) + assert np.all(distances >= 0) + + # All source vertices should have distance 0 + for source_idx in sources: + assert distances[source_idx] == pytest.approx(0.0, abs=1e-6) + + +def test_heat_geodesic_solver_basic(sphere_mesh): + """Test HeatGeodesicSolver for single query.""" + V, F = sphere_mesh + + solver = HeatGeodesicSolver((V, F)) + + # Check num_vertices property + assert solver.num_vertices == len(V) + + # Solve from single source + distances = solver.solve([0]) + + assert isinstance(distances, np.ndarray) + assert distances.shape[0] == len(V) + assert distances[0] == pytest.approx(0.0) + assert np.all(distances >= 0) + + +def test_heat_geodesic_solver_multiple_queries(sphere_mesh): + """Test HeatGeodesicSolver for multiple queries (efficiency use case).""" + V, F = sphere_mesh + + solver = HeatGeodesicSolver((V, F)) + + # Solve multiple times from different sources + d0 = solver.solve([0]) + d1 = solver.solve([5]) + d2 = solver.solve([10]) + + # Each should be valid + for distances in [d0, d1, d2]: + assert isinstance(distances, np.ndarray) + assert distances.shape[0] == len(V) + assert np.all(distances >= 0) + + # Different sources should yield different distance fields + assert not np.allclose(d0, d1) + assert not np.allclose(d1, d2) + + +def test_geodesic_distances_consistency(sphere_mesh): + """Test that function and solver produce same results.""" + V, F = sphere_mesh + sources = [0] + + # Compute with function + d_func = heat_geodesic_distances((V, F), sources) + + # Compute with solver + solver = HeatGeodesicSolver((V, F)) + d_solver = solver.solve(sources) + + # Should be approximately equal + assert np.allclose(d_func, d_solver, rtol=1e-5) + + +def test_geodesic_distances_mesh_conversion(sphere_mesh): + """Test compatibility with COMPAS Mesh objects.""" + V, F = sphere_mesh + + # Create a Mesh object and convert back + mesh = Mesh.from_vertices_and_faces(V, F) + V2, F2 = mesh.to_vertices_and_faces() + + # Compute distances + distances = heat_geodesic_distances((V2, F2), [0]) + + assert isinstance(distances, np.ndarray) + assert distances.shape[0] == len(V2) + assert distances[0] == pytest.approx(0.0) + + +def test_geodesic_isolines_split(sphere_mesh): + """Test splitting mesh along geodesic isolines.""" + V, F = sphere_mesh + + components = geodesic_isolines_split((V, F), [0], [1.0, 2.0]) + + # Should return list of (vertices, faces) tuples + assert isinstance(components, list) + assert len(components) > 1 # Should have multiple components + + # Each component should be valid mesh data + total_faces = 0 + for v, f in components: + assert isinstance(v, np.ndarray) + assert isinstance(f, np.ndarray) + assert v.shape[1] == 3 # 3D vertices + assert f.shape[1] == 3 # Triangle faces + assert f.min() >= 0 + assert f.max() < len(v) + total_faces += len(f) + + # Total faces should be >= original (refinement adds faces) + assert total_faces >= len(F) + + +def test_geodesic_isolines(sphere_mesh): + """Test extracting isoline polylines.""" + V, F = sphere_mesh + + isolines = geodesic_isolines((V, F), [0], [1.0, 2.0]) + + # Should return list of connected polylines + assert isinstance(isolines, list) + assert len(isolines) > 0 # Should have some isolines + + # Each polyline should have multiple points + for pts in isolines: + assert isinstance(pts, np.ndarray) + assert pts.ndim == 2 + assert pts.shape[1] == 3 # 3D points + assert pts.shape[0] >= 2 # At least 2 points + + +def test_geodesic_isolines_empty(sphere_mesh): + """Test isolines with no crossings.""" + V, F = sphere_mesh + + # Use isovalue outside the distance range + isolines = geodesic_isolines((V, F), [0], [100.0]) + + # Should return empty list + assert isinstance(isolines, list) + assert len(isolines) == 0