add geodesic tracing

README.md

@@ -136,9 +136,12 @@ logmap = solver.compute_log_map(sourceV)
 - `MeshVectorHeatSolver.compute_log_map(v_ind)` compute the logarithmic map centered at the given source vertex
   - `v_ind` index of the source vertex
 
-### Mesh Geodesic Paths
 
-Use [edge flips to compute geodesic paths](https://nmwsharp.com/research/flip-geodesics/) on surfaces. These methods are especially useful when you want the path itself, rather than the distance. These routines use an iterative strategy which is quite fast, but note that it is not guaranteed to generate the globally-shortest geodesic (they sometimes find some other very short geodesic instead).
+### Mesh Geodesics Paths
+
+Use [edge flips to compute geodesic paths](https://nmwsharp.com/research/flip-geodesics/) on surfaces. These methods take an initial path, loop, or start & end points along the surface, and straighten the path out to be geodesic.
+
+This approach is mainly useful when you want the path itself, rather than the distance. These routines use an iterative strategy which is quite fast, but note that it is not guaranteed to generate a globally-shortest geodesic (they sometimes find some other very short geodesic instead if straightening falls into different local minimum).
 
 <img src="https://github.com/nmwsharp/potpourri3d/blob/master/media/elephant_geodesic.jpg" height="400">
 
@@ -159,6 +162,28 @@ path_pts = path_solver.find_geodesic_path(v_start=14, v_end=22)
 - `EdgeFlipGeodesicSolver.find_geodesic_loop(v_list, max_iterations=None, max_relative_length_decrease=None)` like `find_geodesic_path_poly()`, but connects the first to last point to find a closed geodesic loop.
 
 In the functions above, the optional argument `max_iterations` is an integer, giving the the maximum number of shortening iterations to perform (default: no limit). The optional argument `max_relative_length_decrease` is a float limiting the maximum decrease in length for the path, e.g. `0.5` would mean the resulting path is at least `0.5 * L` length, where `L` is the initial length.
+
+### Mesh Geodesic Tracing
+
+Given an initial point and direction/length, these routines trace out a geodesic path along the surface of the mesh and return it as a polyline.
+
+```python
+import potpourri3d as pp3d
+
+V, F = # your mesh
+tracer = pp3d.GeodesicTracer(V,F) # shares precomputation for repeated traces
+
+trace_pts = tracer.trace_geodesic_from_vertex(22, np.array((0.3, 0.5, 0.4)))
+# trace_pts is a Vx3 numpy array of points forming the path
+```
+
+- `GeodesicTracer(V, F)` construct an instance of the tracer class.
+  - `V` a Nx3 real numpy array of vertices
+  - `F` a Mx3 integer numpy array of faces, with 0-based vertex indices (must form a manifold, oriented triangle mesh).
+- `GeodesicTracer.trace_geodesic_from_vertex(start_vert, direction_xyz, max_iterations=None)` trace a geodesic from `start_vert`. `direction_xyz` is a length-3 vector giving the direction to walk trace in 3D xyz coordinates, it will be projected onto the tangent space of the vertex. The magnitude of `direction_xyz` determines the distance walked. Output is an `Nx3` numpy array of positions which define the path as a polyline along the surface.
+- `GeodesicTracer.trace_geodesic_from_face(start_face, bary_coords, direction_xyz, max_iterations=None)` similar to above, but from a point in a face. `bary_coords` is a length-3 vector of barycentric coordinates giving the location within the face to start from.
+
+Set `max_iterations` to terminate early after tracing the path through some number of faces/edges (default: no limit).
 
 ### Point Cloud Distance & Vector Heat
 

src/cpp/mesh.cpp

@@ -7,6 +7,7 @@
 #include "geometrycentral/surface/simple_polygon_mesh.h"
 #include "geometrycentral/surface/surface_mesh.h"
 #include "geometrycentral/surface/surface_mesh_factories.h"
+#include "geometrycentral/surface/trace_geodesic.h"
 #include "geometrycentral/surface/vector_heat_method.h"
 #include "geometrycentral/surface/vertex_position_geometry.h"
 #include "geometrycentral/utilities/eigen_interop_helpers.h"
@@ -315,6 +316,90 @@ class EdgeFlipGeodesicsManager {
   std::unique_ptr<FlipEdgeNetwork> flipNetwork;
 };
 
+class GeodesicTracer {
+
+public:
+  GeodesicTracer(DenseMatrix<double> verts, DenseMatrix<int64_t> faces) {
+
+    // Construct the internal mesh and geometry
+    mesh.reset(new ManifoldSurfaceMesh(faces));
+    geom.reset(new VertexPositionGeometry(*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);
+      }
+    }
+
+    geom->requireVertexTangentBasis();
+    geom->requireFaceTangentBasis();
+  }
+
+  // Generate a geodesic by tracing from a vertex along a tangent direction
+  DenseMatrix<double> trace_geodesic_worker(SurfacePoint start_point, Vector2 start_dir,
+                                            size_t max_iters = INVALID_IND) {
+
+    TraceOptions opts;
+    opts.includePath = true;
+    opts.errorOnProblem = false;
+    opts.barrierEdges = nullptr;
+    opts.maxIters = max_iters;
+
+    TraceGeodesicResult result = traceGeodesic(*geom, start_point, start_dir, opts);
+
+    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;
+  }
+
+  // Generate a geodesic by tracing from a vertex along a tangent direction
+  DenseMatrix<double> trace_geodesic_from_vertex(int64_t startVert, Eigen::Vector3d direction_xyz,
+                                                 size_t max_iters = INVALID_IND) {
+
+    // Project the input direction onto the tangent basis
+    Vertex vert = mesh->vertex(startVert);
+    Vector3 direction{direction_xyz(0), direction_xyz(1), direction_xyz(2)};
+    Vector3 bX = geom->vertexTangentBasis[vert][0];
+    Vector3 bY = geom->vertexTangentBasis[vert][1];
+
+    // magnitude matters! it determines the length
+    Vector2 tangent_dir{dot(direction, bX), dot(direction, bY)};
+
+    return trace_geodesic_worker(SurfacePoint(vert), tangent_dir, max_iters);
+  }
+
+  // Generate a geodesic by tracing from a face along a tangent direction
+  DenseMatrix<double> trace_geodesic_from_face(int64_t startFace, Eigen::Vector3d bary_coords,
+                                               Eigen::Vector3d direction_xyz, size_t max_iters = INVALID_IND) {
+
+    // Project the input direction onto the tangent basis
+    Face face = mesh->face(startFace);
+    Vector3 bary{bary_coords(0), bary_coords(1), bary_coords(2)};
+    Vector3 direction{direction_xyz(0), direction_xyz(1), direction_xyz(2)};
+    Vector3 bX = geom->faceTangentBasis[face][0];
+    Vector3 bY = geom->faceTangentBasis[face][1];
+
+    // magnitude matters! it determines the length
+    Vector2 tangent_dir{dot(direction, bX), dot(direction, bY)};
+
+    return trace_geodesic_worker(SurfacePoint(face, bary), tangent_dir, max_iters);
+  }
+
+private:
+  std::unique_ptr<ManifoldSurfaceMesh> mesh;
+  std::unique_ptr<VertexPositionGeometry> geom;
+};
+
 
 // Actual binding code
 // clang-format off
@@ -342,5 +427,10 @@ void bind_mesh(py::module& m) {
         .def("find_geodesic_path_poly", &EdgeFlipGeodesicsManager::find_geodesic_path_poly, py::arg("vert_list"), py::arg("maxIterations"), py::arg("maxRelativeLengthDecrease"))
         .def("find_geodesic_loop", &EdgeFlipGeodesicsManager::find_geodesic_loop, py::arg("vert_list"), py::arg("maxIterations"), py::arg("maxRelativeLengthDecrease"));
 
-  //m.def("read_mesh", &read_mesh, "Read a mesh from file.", py::arg("filename"));
+
+  py::class_<GeodesicTracer>(m, "GeodesicTracer")
+        .def(py::init<DenseMatrix<double>, DenseMatrix<int64_t>>())
+        .def("trace_geodesic_from_vertex", &GeodesicTracer::trace_geodesic_from_vertex, py::arg("start_vert"), py::arg("direction_xyz"), py::arg("max_iters"))
+        .def("trace_geodesic_from_face", &GeodesicTracer::trace_geodesic_from_face, py::arg("start_face"), py::arg("bary_coords"), py::arg("direction_xyz"), py::arg("max_iters"));
+
 }

src/potpourri3d/mesh.py

@@ -86,6 +86,29 @@ def find_geodesic_loop(self, v_list, max_iterations=None, max_relative_length_de
             max_relative_length_decrease = 0.
 
         return self.bound_solver.find_geodesic_loop(v_list, max_iterations, max_relative_length_decrease)
+
+class GeodesicTracer():
+
+    def __init__(self, V, F, t_coef=1.):
+        validate_mesh(V, F, force_triangular=True)
+        self.bound_tracer = pp3db.GeodesicTracer(V, F)
+
+    def trace_geodesic_from_vertex(self, start_vert, direction_xyz, max_iterations=None):
+        if max_iterations is None:
+            max_iterations = 2**63-1
+
+        direction_xyz = np.array(direction_xyz)
+
+        return self.bound_tracer.trace_geodesic_from_vertex(start_vert, direction_xyz, max_iterations)
+
+    def trace_geodesic_from_face(self, start_face, bary_coords, direction_xyz, max_iterations=None):
+        if max_iterations is None:
+            max_iterations = 2**63-1
+
+        bary_coords = np.array(bary_coords)
+        direction_xyz = np.array(direction_xyz)
+
+        return self.bound_tracer.trace_geodesic_from_face(start_face, bary_coords, direction_xyz, max_iterations)
 
 
 def cotan_laplacian(V, F, denom_eps=0.):

test/potpourri3d_test.py

@@ -197,6 +197,32 @@ def test_mesh_flip_geodesic(self):
         self.assertEqual(loop_pts.shape[1], 3)
         loop_pts = path_solver.find_geodesic_loop([307, 757, 190], max_iterations=100, max_relative_length_decrease=0.5)
 
+
+    def test_geodesic_trace(self):
+
+        V, F = pp3d.read_mesh(os.path.join(asset_path, "bunny_small.ply"))
+
+        # Test stateful version
+        tracer = pp3d.GeodesicTracer(V,F)
+
+        # Trace from a vertex
+        trace_pts = tracer.trace_geodesic_from_vertex(22, np.array((0.3, 0.5, 0.4)))
+        self.assertEqual(len(trace_pts.shape), 2)
+        self.assertEqual(trace_pts.shape[1], 3)
+
+        trace_pts = tracer.trace_geodesic_from_vertex(22, np.array((0.3, 0.5, 0.4)), max_iterations=10)
+        self.assertEqual(len(trace_pts.shape), 2)
+        self.assertEqual(trace_pts.shape[1], 3)
+
+        # Trace from a face
+        trace_pts = tracer.trace_geodesic_from_face(31, np.array((0.1, 0.4, 0.5)), np.array((0.3, 0.5, 0.4)))
+        self.assertEqual(len(trace_pts.shape), 2)
+        self.assertEqual(trace_pts.shape[1], 3)
+
+        trace_pts = tracer.trace_geodesic_from_face(31, np.array((0.1, 0.4, 0.5)), np.array((0.3, 0.5, 0.4)), max_iterations=10)
+        self.assertEqual(len(trace_pts.shape), 2)
+        self.assertEqual(trace_pts.shape[1], 3)
+
     def test_point_cloud_distance(self):
 
         P = generate_verts()

test/sample.py

@@ -67,6 +67,15 @@
 ps.register_curve_network("flip loop", loop_pts, edges='loop')
 
 
+# Trace geodesics
+tracer = pp3d.GeodesicTracer(V,F)
+
+trace_pts = tracer.trace_geodesic_from_vertex(22, np.array((0.3, 0.5, 0.4)))
+ps.register_curve_network("trace vertex geodesic", trace_pts, edges='line')
+
+trace_pts = tracer.trace_geodesic_from_face(31, np.array((0.1, 0.4, 0.5)), np.array((0.3, 0.5, 0.4)))
+ps.register_curve_network("trace face geodesic", trace_pts, edges='line')
+
 ## = Point cloud test
 P = V
 ps_cloud = ps.register_point_cloud("cloud", P)