diff --git a/src/cpp/mesh.cpp b/src/cpp/mesh.cpp index b466efe..aee7186 100644 --- a/src/cpp/mesh.cpp +++ b/src/cpp/mesh.cpp @@ -601,7 +601,7 @@ class GeodesicTracer { // 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) { + size_t max_iters = INVALID_IND, bool return_bary = false) { TraceOptions opts; opts.includePath = true; @@ -615,21 +615,61 @@ class GeodesicTracer { 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]; + if (!return_bary){ + // 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; } + // Extract the path as points and barycentric coordinates + // [x,y,z,type,id,t1,t2] + DenseMatrix out(result.pathPoints.size(), 7); // 7 columns + for (size_t i = 0; i < result.pathPoints.size(); i++) { + SurfacePoint spoint = result.pathPoints[i]; + Vector3 point = spoint.interpolate(geom->vertexPositions); + + if (spoint.type == SurfacePointType::Vertex) { + out(i, 0) = point.x; + out(i, 1) = point.y; + out(i, 2) = point.z; + out(i, 3) = 0; // type + out(i, 4) = spoint.vertex.getIndex(); + out(i, 5) = -1; // unused + out(i, 6) = -1; // unused + } + else if (spoint.type == SurfacePointType::Edge) { + out(i, 0) = point.x; + out(i, 1) = point.y; + out(i, 2) = point.z; + out(i, 3) = 1; + out(i, 4) = spoint.edge.getIndex(); + out(i, 5) = spoint.tEdge; + out(i, 6) = -1; + } + else if (spoint.type == SurfacePointType::Face) { + out(i, 0) = point.x; + out(i, 1) = point.y; + out(i, 2) = point.z; + out(i, 3) = 2; + out(i, 4) = spoint.face.getIndex(); + out(i, 5) = spoint.faceCoords.x; + out(i, 6) = spoint.faceCoords.y; + } + } + + 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) { + size_t max_iters = INVALID_IND, bool return_bary = false) { // Project the input direction onto the tangent basis Vertex vert = mesh->vertex(startVert); @@ -640,12 +680,12 @@ class GeodesicTracer { // 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); + return trace_geodesic_worker(SurfacePoint(vert), tangent_dir, max_iters, return_bary); } // 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) { + Eigen::Vector3d direction_xyz, size_t max_iters = INVALID_IND, bool return_bary = false) { // Project the input direction onto the tangent basis Face face = mesh->face(startFace); @@ -657,7 +697,7 @@ class GeodesicTracer { // 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); + return trace_geodesic_worker(SurfacePoint(face, bary), tangent_dir, max_iters, return_bary); } private: @@ -724,7 +764,7 @@ void bind_mesh(py::module& m) { 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")); + .def("trace_geodesic_from_vertex", &GeodesicTracer::trace_geodesic_from_vertex, py::arg("start_vert"), py::arg("direction_xyz"), py::arg("max_iters"), py::arg("return_bary")) + .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"), py::arg("return_bary")); } diff --git a/src/potpourri3d/mesh.py b/src/potpourri3d/mesh.py index e03e9e3..22824d3 100644 --- a/src/potpourri3d/mesh.py +++ b/src/potpourri3d/mesh.py @@ -162,23 +162,43 @@ 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): + def trace_geodesic_from_vertex(self, start_vert, direction_xyz, max_iterations=None, return_bary=False): 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) + trace_result = self.bound_tracer.trace_geodesic_from_vertex(start_vert, direction_xyz, max_iterations, return_bary) - def trace_geodesic_from_face(self, start_face, bary_coords, direction_xyz, max_iterations=None): + if not return_bary: return trace_result + return self.extract_xyz_and_bary(trace_result) + + def trace_geodesic_from_face(self, start_face, bary_coords, direction_xyz, max_iterations=None, return_bary=False): 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) - + trace_result = self.bound_tracer.trace_geodesic_from_face(start_face, bary_coords, direction_xyz, max_iterations, return_bary) + + if not return_bary: return trace_result + return self.extract_xyz_and_bary(trace_result) + + def extract_xyz_and_bary(self,arr): + xyz = arr[:, :3].copy() # Nx3 + + types = arr[:, 3].astype(int) + ids = arr[:, 4].astype(int) + params1 = arr[:, 5] + params2 = arr[:, 6] + + bary = [ + (id_, [] if t == 0 else [params1[i]] if t == 1 else [params1[i], params2[i]]) + for i, (t, id_) in enumerate(zip(types, ids)) + ] + + return xyz, bary def cotan_laplacian(V, F, denom_eps=0.): validate_mesh(V, F, force_triangular=True) diff --git a/test/potpourri3d_test.py b/test/potpourri3d_test.py index 33d9450..f506b20 100644 --- a/test/potpourri3d_test.py +++ b/test/potpourri3d_test.py @@ -238,6 +238,75 @@ def test_geodesic_trace(self): self.assertEqual(len(trace_pts.shape), 2) self.assertEqual(trace_pts.shape[1], 3) + + def test_geodesic_trace_bary(self): + V, F = pp3d.read_mesh(os.path.join(asset_path, "bunny_small.ply")) + tracer = pp3d.GeodesicTracer(V, F) + + # --- Trace from a vertex --- + trace_pts, barys = tracer.trace_geodesic_from_vertex( + 22, + np.array((0.3, 0.5, 0.4)), + return_bary=True + ) + self.assertEqual(len(trace_pts.shape), 2) + self.assertEqual(trace_pts.shape[1], 3) + self._check_bary_structure(trace_pts, barys) + + trace_pts, barys = tracer.trace_geodesic_from_vertex( + 22, + np.array((0.3, 0.5, 0.4)), + max_iterations=10, + return_bary=True + ) + self.assertEqual(len(trace_pts.shape), 2) + self.assertEqual(trace_pts.shape[1], 3) + self._check_bary_structure(trace_pts, barys) + + # --- Trace from a face --- + trace_pts, barys = tracer.trace_geodesic_from_face( + 31, + np.array((0.1, 0.4, 0.5)), + np.array((0.3, 0.5, 0.4)), + return_bary=True + ) + self.assertEqual(len(trace_pts.shape), 2) + self.assertEqual(trace_pts.shape[1], 3) + self._check_bary_structure(trace_pts, barys) + + trace_pts, barys = 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, + return_bary=True + ) + self.assertEqual(len(trace_pts.shape), 2) + self.assertEqual(trace_pts.shape[1], 3) + self._check_bary_structure(trace_pts, barys) + + + def _check_bary_structure(self, trace_pts, barys): + """Helper to verify that barycentric return has the expected structure.""" + # barys should exist and match trace length + self.assertTrue(hasattr(barys, '__len__')) + self.assertEqual(len(barys), len(trace_pts)) + + for entry in barys: + # each entry must be a 2-tuple: (primitive_index, coords) + self.assertIsInstance(entry, (list, tuple)) + self.assertEqual(len(entry), 2) + + primitive_id, coords = entry + + # primitive_id: integer index + self.assertIsInstance(primitive_id, (int, np.integer)) + + # coords: list or ndarray of numbers (possibly empty) + self.assertTrue(isinstance(coords, (list, np.ndarray))) + for c in coords: + self.assertTrue(isinstance(c, (float, np.floating, int, np.integer))) + def test_point_cloud_distance(self): P = generate_verts()