Add Geodesic Trace Function

README.md

@@ -157,6 +157,8 @@ 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.
+
+- `trace_geodesic(V, F, startVert, directionX, directionY)` computes a geodesic tracing from a starting vertex along a tangent-space direction.  Note that this is separate from the `EdgeFlipGeodesicSolver`.
 
 ### 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>
@@ -312,6 +314,40 @@ class EdgeFlipGeodesicsManager {
   std::unique_ptr<FlipEdgeNetwork> flipNetwork;
 };
 
+// Generate a geodesic by tracing from a vertex along a tangent direction
+DenseMatrix<double> trace_geodesic(DenseMatrix<double> verts, DenseMatrix<int64_t> faces, 
+                                     int64_t start_vert, double direction_x, double direction_y,
+                                     size_t max_iters = INVALID_IND) {
+
+  // Construct the internal mesh and geometry
+  ManifoldSurfaceMesh mesh(faces);
+  VertexPositionGeometry geom(mesh);
+  for (size_t i = 0; i < mesh.nVertices(); i++) {
+    for (size_t j = 0; j < 3; j++) {
+      geom.inputVertexPositions[i][j] = verts(i, j);
+    }
+  }
+
+  TraceGeodesicResult result =
+      traceGeodesic(geom, SurfacePoint(mesh.vertex(start_vert)), 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;
+}
+
 
 // Actual binding code
 // clang-format off
@@ -338,6 +374,9 @@ void bind_mesh(py::module& m) {
         .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"));
+
+  m.def("trace_geodesic", &trace_geodesic, py::arg("verts"), py::arg("faces"),
+    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

@@ -71,7 +71,11 @@ 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, V, F, start_vert, direction_x, direction_y, max_iters = 1000):
+    return pp3db.trace_geodesic(V, F, start_vert, direction_x, direction_y, max_iters)
+
 
 def cotan_laplacian(V, F, denom_eps=0.):
     validate_mesh(V, F, force_triangular=True)