Begin adding Geodesic Trace Function

README.md

@@ -157,6 +157,7 @@ path_pts = path_solver.find_geodesic_path(v_start=14, v_end=22)
 - `EdgeFlipGeodesicSolver.find_geodesic_path(v_start, v_end)` compute a geodesic from `v_start` to `v_end`. Output is an `Nx3` numpy array of positions which define the path as a polyline along the surface.
 - `EdgeFlipGeodesicSolver.find_geodesic_path_poly(v_list)` like `find_geodesic_path()`, but takes as input a list of vertices `[v_start, v_a, v_b, ..., v_end]`, which is shorted to find a path from `v_start` to `v_end`. Useful for finding geodesics which are not shortest paths. The input vertices do not need to be connected; the routine internally constructs a piecwise-Dijkstra path between them. However, that path must not cross itself.
 - `EdgeFlipGeodesicSolver.find_geodesic_loop(v_list)` like `find_geodesic_path_poly()`, but connects the first to last point to find a closed geodesic loop.
+- `EdgeFlipGeodesicSolver.trace_geodesic(startVert, directionX, directionY)` computes a geodesic tracing from a starting vertex along a tangent-space direction.
 
 ### Point Cloud Distance & Vector Heat
 

src/cpp/mesh.cpp

@@ -9,6 +9,8 @@
 #include "geometrycentral/surface/surface_mesh_factories.h"
 #include "geometrycentral/surface/vector_heat_method.h"
 #include "geometrycentral/surface/vertex_position_geometry.h"
+#include "geometrycentral/surface/surface_point.h"
+#include "geometrycentral/surface/trace_geodesic.h"
 #include "geometrycentral/utilities/eigen_interop_helpers.h"
 
 #include <pybind11/eigen.h>
@@ -306,6 +308,29 @@ class EdgeFlipGeodesicsManager {
     return out;
   }
 
+  // Generate a geodesic by tracing from a vertex along a tangent direction
+  DenseMatrix<double> trace_geodesic(int64_t startVert, double direction_x, double direction_y,
+                                     size_t max_iters = INVALID_IND) {
+    TraceGeodesicResult result =
+        traceGeodesic(*geom, SurfacePoint(mesh->vertex(startVert)), Vector2{direction_x, direction_y},
+                      TraceOptions{true, false, NULL, max_iters});
+
+    if (!result.hasPath) {
+      throw std::runtime_error("geodesic trace encountered an error");
+    }
+
+    // Extract the path and store it in the vector
+    DenseMatrix<double> out(result.pathPoints.size(), 3);
+    for (size_t i = 0; i < result.pathPoints.size(); i++) {
+      Vector3 point = result.pathPoints[i].interpolate(geom->vertexPositions);
+      for (size_t j = 0; j < 3; j++) {
+        out(i, j) = point[j];
+      }
+    }
+
+    return out;
+  }
+
 private:
   std::unique_ptr<ManifoldSurfaceMesh> mesh;
   std::unique_ptr<VertexPositionGeometry> geom;
@@ -337,7 +362,8 @@ void bind_mesh(py::module& m) {
         .def(py::init<DenseMatrix<double>, DenseMatrix<int64_t>>())
         .def("find_geodesic_path", &EdgeFlipGeodesicsManager::find_geodesic_path, py::arg("source_vert"), py::arg("target_vert"))
         .def("find_geodesic_path_poly", &EdgeFlipGeodesicsManager::find_geodesic_path_poly, py::arg("vert_list"))
-        .def("find_geodesic_loop", &EdgeFlipGeodesicsManager::find_geodesic_loop, py::arg("vert_list"));
+        .def("find_geodesic_loop", &EdgeFlipGeodesicsManager::find_geodesic_loop, py::arg("vert_list"))
+        .def("trace_geodesic", &EdgeFlipGeodesicsManager::trace_geodesic, py::arg("start_vert"), py::arg("direction_x"), py::arg("direction_y"), py::arg("max_iters"));
 
   //m.def("read_mesh", &read_mesh, "Read a mesh from file.", py::arg("filename"));
 }

src/potpourri3d/mesh.py

@@ -72,6 +72,9 @@ def find_geodesic_path_poly(self, v_list):
     def find_geodesic_loop(self, v_list):
         return self.bound_solver.find_geodesic_loop(v_list)
 
+    def trace_geodesic(self, start_vert, direction_x, direction_y, max_iters = 1000):
+        return self.bound_solver.trace_geodesic(start_vert, direction_x, direction_y, max_iters)
+
 
 def cotan_laplacian(V, F, denom_eps=0.):
     validate_mesh(V, F, force_triangular=True)