ADD geodistance.

Author: petrasvestartas

docs/examples/geodistance.py

@@ -20,7 +20,7 @@
 # ==============================================================================
 
 source = trimesh.vertex_sample(size=1)[0]
-distance = igl.trimesh_geodistance(trimesh.to_vertices_and_faces(), source, method="heat")
+distance = igl.trimesh_geodistance(trimesh.to_vertices_and_faces(), source, method="exact")
 
 # ==============================================================================
 # Visualize

src/compas.hpp

@@ -23,12 +23,15 @@
 #include <Eigen/Geometry>
 
 // libigl includes
+#include <igl/boundary_loop.h>
 #include <igl/cotmatrix.h>
 #include <igl/massmatrix.h>
 #include <igl/invert_diag.h>
 #include <igl/gaussian_curvature.h>
 #include <igl/principal_curvature.h>
-#include <igl/boundary_loop.h>
+#include <igl/avg_edge_length.h>
+#include <igl/exact_geodesic.h>
+#include <igl/heat_geodesics.h>
 #include <igl/map_vertices_to_circle.h>
 #include <igl/harmonic.h>
 #include <igl/doublearea.h>

src/compas_libigl/__init__.py

@@ -3,6 +3,7 @@
 from ._nanobind import add, __doc__
 from .boundaries import trimesh_boundaries
 from .curvature import trimesh_gaussian_curvature, trimesh_principal_curvature
+from .geodistance import trimesh_geodistance
 
 # from .boundaries import trimesh_boundaries
 # from .curvature import trimesh_gaussian_curvature, trimesh_principal_curvature
@@ -106,6 +107,7 @@ def get_armadillo():
     trimesh_boundaries,
     trimesh_gaussian_curvature,
     trimesh_principal_curvature,
+    trimesh_geodistance,
     # "get",
     # "get_beetle",
     # "get_armadillo",

src/compas_libigl/geodistance.py

@@ -0,0 +1,43 @@
+import numpy as np
+from compas.plugins import plugin
+
+from compas_libigl import _geodistance
+
+
+@plugin(category="trimesh")
+def trimesh_geodistance(M, source, method="exact"):
+    """Compute the geodesic distance from every vertex of the mesh to a source vertex.
+
+    Parameters
+    ----------
+    M : (list, list)
+        A mesh represented by a list of vertices and a list of faces.
+    source : int
+        The index of the vertex from where the geodesic distances should be calculated.
+    method : {'exact', 'heat'}
+        The method for calculating the distances.
+        Default is `'exact'`.
+
+    Returns
+    -------
+    list of float
+        A list of geodesic distances from the source vertex.
+
+    Raises
+    ------
+    NotImplementedError
+        If ``method`` is not one of ``{'exact', 'heat'}``.
+
+    Examples
+    --------
+    >>>
+
+    """
+    V, F = M
+    V = np.asarray(V, dtype=np.float64)
+    F = np.asarray(F, dtype=np.int32)
+    if method == "exact":
+        return _geodistance.trimesh_geodistance_exact(V, F, source)
+    if method == "heat":
+        return _geodistance.trimesh_geodistance_heat(V, F, source)
+    raise NotImplementedError

src/geodistance.cpp

@@ -1,13 +1,64 @@
 #include "geodistance.hpp"
 
-void function(const Eigen::MatrixXd& V, const Eigen::MatrixXi& F) {}
+Eigen::VectorXd
+trimesh_geodistance_exact(
+    const compas::RowMatrixXd& V,
+    const compas::RowMatrixXi& F,
+    int vid)
+{
+	Eigen::VectorXd D;
+	Eigen::VectorXi VS, FS, VT, FT;
+
+	VS.resize(1);
+	VS << vid;
+
+	VT.setLinSpaced(V.rows(), 0, V.rows() - 1);
+
+	igl::exact_geodesic(V, F, VS, FS, VT, FT, D);
+
+    return D;
+}
+
+Eigen::VectorXd
+trimesh_geodistance_heat(
+    const compas::RowMatrixXd& V,
+    const compas::RowMatrixXi& F,
+    int vid)
+{
+    Eigen::VectorXi gamma;
+    gamma.resize(1);
+    gamma << vid;
+
+    igl::HeatGeodesicsData<double> data;
+    double t = std::pow(igl::avg_edge_length(V, F), 2);
+    igl::heat_geodesics_precompute(V, F, t, data);
+
+    Eigen::VectorXd D = Eigen::VectorXd::Zero(data.Grad.cols());
+    D(vid) = 1;
+
+    igl::heat_geodesics_solve(data, gamma, D);
+
+    return D;
+}
+
+
+
 
 NB_MODULE(_geodistance, m) {
 
     m.def(
-        "function_name",
-        &function,
+        "trimesh_geodistance_exact",
+        &trimesh_geodistance_exact,
+        "Description.",
+        "V"_a,
+        "F"_a,
+        "vid"_a);
+
+    m.def(
+        "trimesh_geodistance_heat",
+        &trimesh_geodistance_heat,
         "Description.",
         "V"_a,
-        "F"_a);
+        "F"_a,
+        "vid"_a);
 }