Add PointCloudLocalTriangulation interface (#11)

* Add PointCloudLocalTriangulation interface

	Add class in src/point_cloud.cpp
	Expose class in point_cloud.py
	Add test based on cartwheel pointcloud

* Use std::max instead of MAX

Author: lkskstlr

src/cpp/point_cloud.cpp

@@ -3,6 +3,7 @@
 #include "geometrycentral/pointcloud/point_cloud_heat_solver.h"
 #include "geometrycentral/pointcloud/point_cloud_io.h"
 #include "geometrycentral/pointcloud/point_position_geometry.h"
+#include "geometrycentral/pointcloud/local_triangulation.h"
 #include "geometrycentral/utilities/eigen_interop_helpers.h"
 
 #include <pybind11/eigen.h>
@@ -11,6 +12,7 @@
 
 #include "Eigen/Dense"
 
+
 namespace py = pybind11;
 
 using namespace geometrycentral;
@@ -120,6 +122,59 @@ class PointCloudHeatSolverEigen {
   std::unique_ptr<PointCloudHeatSolver> solver;
 };
 
+// A class that exposes the local pointcloud trinagulation to python
+class PointCloudLocalTriangulation {
+public:
+  PointCloudLocalTriangulation(DenseMatrix<double> points, bool withDegeneracyHeuristic) : withDegeneracyHeuristic(withDegeneracyHeuristic) {
+
+    // Construct the internal cloud and geometry
+    cloud.reset(new PointCloud(points.rows()));
+    geom.reset(new PointPositionGeometry(*cloud));
+    for (size_t i = 0; i < cloud->nPoints(); i++) {
+      for (size_t j = 0; j < 3; j++) {
+        geom->positions[i][j] = points(i, j);
+      }
+    }
+  }
+
+
+  Eigen::Matrix<int, Eigen::Dynamic, Eigen::Dynamic> get_local_triangulation() {
+    PointData<std::vector<std::array<Point, 3>>> local_triangulation = buildLocalTriangulations(*cloud, *geom, withDegeneracyHeuristic);
+
+    int max_neigh = 0;
+
+    size_t idx = 0;
+    for (Point v : cloud->points()) {
+      max_neigh = std::max(max_neigh, static_cast<int>(local_triangulation[v].size()));
+      if (idx != v.getIndex()) {
+        py::print("Error. Index of points not consistent. (Idx, v.getIndex) = ", idx, v.getIndex());
+      }
+      idx++;
+    }
+
+    Eigen::Matrix<int, Eigen::Dynamic, Eigen::Dynamic> out(cloud->nPoints(), 3 * max_neigh);
+    out.setConstant(-1);
+
+    for (Point v : cloud->points()) {
+      int i = 0;
+      for (auto const &neighs : local_triangulation[v]) {
+        out(v.getIndex(), i + 0) = (int) neighs[0].getIndex();
+        out(v.getIndex(), i + 1) = (int) neighs[1].getIndex();
+        out(v.getIndex(), i + 2) = (int) neighs[2].getIndex();
+        i += 3;
+      }
+    }
+
+    return out;
+  }
+
+private:
+  const bool withDegeneracyHeuristic;
+  std::unique_ptr<PointCloud> cloud;
+  std::unique_ptr<PointPositionGeometry> geom;
+  std::unique_ptr<PointCloudHeatSolver> solver;
+};
+
 
 // Actual binding code
 // clang-format off
@@ -134,4 +189,8 @@ void bind_point_cloud(py::module& m) {
         .def("transport_tangent_vector", &PointCloudHeatSolverEigen::transport_tangent_vector, py::arg("source_point"), py::arg("vector"))
         .def("transport_tangent_vectors", &PointCloudHeatSolverEigen::transport_tangent_vectors, py::arg("source_points"), py::arg("vectors"))
         .def("compute_log_map", &PointCloudHeatSolverEigen::compute_log_map, py::arg("source_point"));
+
+  py::class_<PointCloudLocalTriangulation>(m, "PointCloudLocalTriangulation")
+    .def(py::init<DenseMatrix<double>, bool>())
+    .def("get_local_triangulation", &PointCloudLocalTriangulation::get_local_triangulation);
 }

src/potpourri3d/point_cloud.py

@@ -35,3 +35,24 @@ def transport_tangent_vectors(self, p_inds, vectors):
 
     def compute_log_map(self, p_ind):
         return self.bound_solver.compute_log_map(p_ind)
+
+
+class PointCloudLocalTriangulation():
+
+    def __init__(self, P, with_degeneracy_heuristic=True):
+        validate_points(P)
+        self.bound_triangulation = pp3db.PointCloudLocalTriangulation(P, with_degeneracy_heuristic)
+
+    def get_local_triangulation(self):
+        """Return the local point cloud triangulation
+        
+        The out matrix has the following convention:
+            size: num_points, max_neighs, 3. max_neighs is the maximum number of neighbors
+            out[point_idx, neigh_idx, :] are the indices of the 3 neighbors
+            -1 is used as the fill value for unused elements if num_neighs < max_neighs for a point
+        """
+        out = self.bound_triangulation.get_local_triangulation()
+        assert out.shape[-1] % 3 == 0
+        max_neighs = out.shape[-1] // 3
+        out = np.reshape(out, [-1, max_neighs, 3])
+        return out

test/potpourri3d_test.py

@@ -237,5 +237,43 @@ def test_point_cloud_vector_heat(self):
         self.assertEqual(logmap.shape[0], P.shape[0])
         self.assertEqual(logmap.shape[1], 2)
 
+    def test_point_cloud_local_triangulation(self):
+        # Test local triangulation for a "cartwheel" pointcloud
+
+        num = 31
+        t = np.linspace(0, 2*np.pi, num-1, endpoint=False)
+        points = np.concatenate([np.zeros([1, 3]), np.stack([np.cos(t), np.sin(t), 0*t], 1)], 0)
+
+        pcl_local_tri = pp3d.PointCloudLocalTriangulation(points)
+        idxs = pcl_local_tri.get_local_triangulation()
+
+        def next_id(i):
+            assert i != 0
+            if i == num-1:
+                return 1
+            return i+1
+
+        def prev_id(i):
+            assert i != 0
+            if i == 1:
+                return num-1
+            return i-1
+
+
+        # Explicitly check for cartwheel
+        # Use sets for order invariance
+
+        res0 = set(tuple(r) for r in idxs[0])
+        ref0 = set((0, j, next_id(j)) for j in range(1, num))
+        self.assertEqual(res0, ref0)
+
+        for i in range(1, num):
+            self.assertTrue(np.all(idxs[i, 2:] == -1))
+            
+            res = set(tuple(r) for r in idxs[i, :2])
+            ref = {(i, next_id(i), 0), (i, 0, prev_id(i))}
+            self.assertEqual(res, ref)
+
+
 if __name__ == '__main__':
     unittest.main()