add limits to flip geodesic shortener

README.md

@@ -154,9 +154,11 @@ path_pts = path_solver.find_geodesic_path(v_start=14, v_end=22)
 - `EdgeFlipGeodesicSolver(V, F)` construct an instance of the solver 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).
-- `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.find_geodesic_path(v_start, v_end, max_iterations=None, max_relative_length_decrease=None)` 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, max_iterations=None, max_relative_length_decrease=None)` 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, 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.
 
 ### Point Cloud Distance & Vector Heat
 

src/cpp/mesh.cpp

@@ -186,7 +186,8 @@ class EdgeFlipGeodesicsManager {
   }
 
   // Generate a point-to-point geodesic by straightening a Dijkstra path
-  DenseMatrix<double> find_geodesic_path(int64_t startVert, int64_t endVert) {
+  DenseMatrix<double> find_geodesic_path(int64_t startVert, int64_t endVert, size_t maxIterations = INVALID_IND,
+                                         double maxRelativeLengthDecrease = 0.) {
 
     // Get an initial dijkstra path
     std::vector<Halfedge> dijkstraPath = shortestEdgePath(*geom, mesh->vertex(startVert), mesh->vertex(endVert));
@@ -202,7 +203,7 @@ class EdgeFlipGeodesicsManager {
     flipNetwork->reinitializePath({dijkstraPath});
 
     // Straighten the path to geodesic
-    flipNetwork->iterativeShorten();
+    flipNetwork->iterativeShorten(maxIterations, maxRelativeLengthDecrease);
 
     // Extract the path and store it in the vector
     std::vector<Vector3> path3D = flipNetwork->getPathPolyline3D().front();
@@ -220,7 +221,8 @@ class EdgeFlipGeodesicsManager {
   }
 
   // Generate a point-to-point geodesic by straightening a poly-geodesic path
-  DenseMatrix<double> find_geodesic_path_poly(std::vector<int64_t> verts) {
+  DenseMatrix<double> find_geodesic_path_poly(std::vector<int64_t> verts, size_t maxIterations = INVALID_IND,
+                                              double maxRelativeLengthDecrease = 0.) {
 
     // Convert to a list of vertices
     std::vector<Halfedge> halfedges;
@@ -245,7 +247,7 @@ class EdgeFlipGeodesicsManager {
     flipNetwork->reinitializePath({halfedges});
 
     // Straighten the path to geodesic
-    flipNetwork->iterativeShorten();
+    flipNetwork->iterativeShorten(maxIterations, maxRelativeLengthDecrease);
 
     // Extract the path and store it in the vector
     std::vector<Vector3> path3D = flipNetwork->getPathPolyline3D().front();
@@ -264,7 +266,8 @@ class EdgeFlipGeodesicsManager {
 
 
   // Generate a point-to-point geodesic loop by straightening a poly-geodesic path
-  DenseMatrix<double> find_geodesic_loop(std::vector<int64_t> verts) {
+  DenseMatrix<double> find_geodesic_loop(std::vector<int64_t> verts, size_t maxIterations = INVALID_IND,
+                                         double maxRelativeLengthDecrease = 0.) {
 
     // Convert to a list of vertices
     std::vector<Halfedge> halfedges;
@@ -289,7 +292,7 @@ class EdgeFlipGeodesicsManager {
     flipNetwork->reinitializePath({halfedges});
 
     // Straighten the path to geodesic
-    flipNetwork->iterativeShorten();
+    flipNetwork->iterativeShorten(maxIterations, maxRelativeLengthDecrease);
 
     // Extract the path and store it in the vector
     std::vector<Vector3> path3D = flipNetwork->getPathPolyline3D().front();
@@ -335,9 +338,9 @@ void bind_mesh(py::module& m) {
 
   py::class_<EdgeFlipGeodesicsManager>(m, "EdgeFlipGeodesicsManager")
         .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_path", &EdgeFlipGeodesicsManager::find_geodesic_path, py::arg("source_vert"), py::arg("target_vert"), py::arg("maxIterations"), py::arg("maxRelativeLengthDecrease"))
+        .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"));
 }

src/potpourri3d/mesh.py

@@ -63,14 +63,29 @@ def __init__(self, V, F, t_coef=1.):
         validate_mesh(V, F, force_triangular=True)
         self.bound_solver = pp3db.EdgeFlipGeodesicsManager(V, F)
 
-    def find_geodesic_path(self, v_start, v_end):
-        return self.bound_solver.find_geodesic_path(v_start, v_end)
+    def find_geodesic_path(self, v_start, v_end, max_iterations=None, max_relative_length_decrease=None):
+        if max_iterations is None:
+            max_iterations = 2**63-1
+        if max_relative_length_decrease is None:
+            max_relative_length_decrease = 0.
+
+        return self.bound_solver.find_geodesic_path(v_start, v_end, max_iterations, max_relative_length_decrease)
 
-    def find_geodesic_path_poly(self, v_list):
-        return self.bound_solver.find_geodesic_path_poly(v_list)
+    def find_geodesic_path_poly(self, v_list, max_iterations=None, max_relative_length_decrease=None):
+        if max_iterations is None:
+            max_iterations = 2**63-1
+        if max_relative_length_decrease is None:
+            max_relative_length_decrease = 0.
+
+        return self.bound_solver.find_geodesic_path_poly(v_list, max_iterations, max_relative_length_decrease)
 
-    def find_geodesic_loop(self, v_list):
-        return self.bound_solver.find_geodesic_loop(v_list)
+    def find_geodesic_loop(self, v_list, max_iterations=None, max_relative_length_decrease=None):
+        if max_iterations is None:
+            max_iterations = 2**63-1
+        if max_relative_length_decrease is None:
+            max_relative_length_decrease = 0.
+
+        return self.bound_solver.find_geodesic_loop(v_list, max_iterations, max_relative_length_decrease)
 
 
 def cotan_laplacian(V, F, denom_eps=0.):

test/potpourri3d_test.py

@@ -171,6 +171,7 @@ def test_mesh_flip_geodesic(self):
         path_pts = path_solver.find_geodesic_path(v_start=14, v_end=22)
         self.assertEqual(len(path_pts.shape), 2)
         self.assertEqual(path_pts.shape[1], 3)
+        path_pts = path_solver.find_geodesic_path(v_start=14, v_end=22, max_iterations=100, max_relative_length_decrease=0.5)
 
         # Do some more
         for i in range(5):
@@ -182,6 +183,7 @@ def test_mesh_flip_geodesic(self):
         path_pts = path_solver.find_geodesic_path_poly([1173, 148, 870, 898])
         self.assertEqual(len(path_pts.shape), 2)
         self.assertEqual(path_pts.shape[1], 3)
+        path_pts = path_solver.find_geodesic_path_poly([1173, 148, 870, 898], max_iterations=100, max_relative_length_decrease=0.5)
 
         # Do a loop
         loop_pts = path_solver.find_geodesic_loop([1173, 148, 870, 898])
@@ -193,6 +195,7 @@ def test_mesh_flip_geodesic(self):
         loop_pts = path_solver.find_geodesic_loop([307, 757, 190]) 
         self.assertEqual(len(loop_pts.shape), 2)
         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_point_cloud_distance(self):
 

test/sample.py

@@ -63,7 +63,7 @@
 loop_pts = path_solver.find_geodesic_loop([1173, 148, 870, 898])
 ps.register_curve_network("flip loop", loop_pts, edges='loop')
 
-loop_pts = path_solver.find_geodesic_loop([307, 757, 190]) # this one contracts to a point
+loop_pts = path_solver.find_geodesic_loop([307, 757, 190], max_relative_length_decrease=.9) # this one otherwise contracts to a point
 ps.register_curve_network("flip loop", loop_pts, edges='loop')