diff --git a/README.md b/README.md index 167b927..ac78c1b 100644 --- a/README.md +++ b/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). @@ -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 diff --git a/src/cpp/mesh.cpp b/src/cpp/mesh.cpp index 775f926..1a98eb6 100644 --- a/src/cpp/mesh.cpp +++ b/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 flipNetwork; }; +class GeodesicTracer { + +public: + GeodesicTracer(DenseMatrix verts, DenseMatrix 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 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 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 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 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 mesh; + std::unique_ptr 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_(m, "GeodesicTracer") + .def(py::init, DenseMatrix>()) + .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")); + } diff --git a/src/potpourri3d/mesh.py b/src/potpourri3d/mesh.py index 53c6ed4..8972386 100644 --- a/src/potpourri3d/mesh.py +++ b/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.): diff --git a/test/potpourri3d_test.py b/test/potpourri3d_test.py index 1cfe605..cd10da3 100644 --- a/test/potpourri3d_test.py +++ b/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() diff --git a/test/sample.py b/test/sample.py index 89c4d35..6077dcb 100644 --- a/test/sample.py +++ b/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)