From 7b2294018d9c86c3f6bdf1512265714c66074c56 Mon Sep 17 00:00:00 2001 From: Nicholas Sharp Date: Wed, 19 Jun 2024 12:39:34 -0700 Subject: [PATCH] add limits to flip geodesic shortener --- README.md | 8 +++++--- src/cpp/mesh.cpp | 21 ++++++++++++--------- src/potpourri3d/mesh.py | 27 +++++++++++++++++++++------ test/potpourri3d_test.py | 3 +++ test/sample.py | 2 +- 5 files changed, 42 insertions(+), 19 deletions(-) diff --git a/README.md b/README.md index d805466..167b927 100644 --- a/README.md +++ b/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 diff --git a/src/cpp/mesh.cpp b/src/cpp/mesh.cpp index e0d8f32..775f926 100644 --- a/src/cpp/mesh.cpp +++ b/src/cpp/mesh.cpp @@ -186,7 +186,8 @@ class EdgeFlipGeodesicsManager { } // Generate a point-to-point geodesic by straightening a Dijkstra path - DenseMatrix find_geodesic_path(int64_t startVert, int64_t endVert) { + DenseMatrix find_geodesic_path(int64_t startVert, int64_t endVert, size_t maxIterations = INVALID_IND, + double maxRelativeLengthDecrease = 0.) { // Get an initial dijkstra path std::vector 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 path3D = flipNetwork->getPathPolyline3D().front(); @@ -220,7 +221,8 @@ class EdgeFlipGeodesicsManager { } // Generate a point-to-point geodesic by straightening a poly-geodesic path - DenseMatrix find_geodesic_path_poly(std::vector verts) { + DenseMatrix find_geodesic_path_poly(std::vector verts, size_t maxIterations = INVALID_IND, + double maxRelativeLengthDecrease = 0.) { // Convert to a list of vertices std::vector 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 path3D = flipNetwork->getPathPolyline3D().front(); @@ -264,7 +266,8 @@ class EdgeFlipGeodesicsManager { // Generate a point-to-point geodesic loop by straightening a poly-geodesic path - DenseMatrix find_geodesic_loop(std::vector verts) { + DenseMatrix find_geodesic_loop(std::vector verts, size_t maxIterations = INVALID_IND, + double maxRelativeLengthDecrease = 0.) { // Convert to a list of vertices std::vector 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 path3D = flipNetwork->getPathPolyline3D().front(); @@ -335,9 +338,9 @@ void bind_mesh(py::module& m) { py::class_(m, "EdgeFlipGeodesicsManager") .def(py::init, DenseMatrix>()) - .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")); } diff --git a/src/potpourri3d/mesh.py b/src/potpourri3d/mesh.py index 2d1320c..53c6ed4 100644 --- a/src/potpourri3d/mesh.py +++ b/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.): diff --git a/test/potpourri3d_test.py b/test/potpourri3d_test.py index 32d4a67..1cfe605 100644 --- a/test/potpourri3d_test.py +++ b/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): diff --git a/test/sample.py b/test/sample.py index 00240ca..89c4d35 100644 --- a/test/sample.py +++ b/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')