[Feature] Expose geodesic trace as barycentric coordinates

src/cpp/mesh.cpp

@@ -601,7 +601,7 @@ class GeodesicTracer {
 
   // 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) {
+                                            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<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];
+    if (!return_bary){
+      // 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;
     }
 
+    // Extract the path as points and barycentric coordinates
+    // [x,y,z,type,id,t1,t2]
+    DenseMatrix<double> 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<double> 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<double> 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_<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"));
+        .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"));
 
 }

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)

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()