first attempt at surface reconstruction

docs/examples/reconstruction.py

@@ -0,0 +1,39 @@
+import math
+from pathlib import Path
+from compas.datastructures import Mesh
+from compas.geometry import Pointcloud
+from compas.geometry import Rotation, Scale
+from compas_view2.app import App
+from compas_cgal.reconstruction import poisson_surface_reconstruction
+
+HERE = Path(__file__).parent
+DATA = HERE / ".." / ".." / "data"
+FILE = DATA / "oni.xyz"
+
+points = []
+normals = []
+with open(FILE, "r") as f:
+    for line in f:
+        x, y, z, nx, ny, nz = line.strip().split()
+        points.append([float(x), float(y), float(z)])
+        normals.append([float(nx), float(ny), float(nz)])
+
+V, F = poisson_surface_reconstruction(points, normals)
+mesh = Mesh.from_vertices_and_faces(V, F)
+
+R = Rotation.from_axis_and_angle([1, 0, 0], math.radians(90))
+S = Scale.from_factors([5, 5, 5])
+T = R * S
+
+cloud = Pointcloud(V)
+cloud.transform(T)
+mesh.transform(T)
+
+viewer = App(width=1600, height=900)
+viewer.view.camera.position = [-5, -5, 1.5]
+viewer.view.camera.look_at([0, 0, 1.5])
+
+viewer.add(mesh)
+viewer.add(cloud)
+
+viewer.run()

docs/examples/reconstruction.rst

@@ -0,0 +1,11 @@
+********************************************************************************
+Poisson Surface Reconstruction
+********************************************************************************
+
+.. figure:: /_images/cgal_reconstruction.png
+    :figclass: figure
+    :class: figure-img img-fluid
+
+
+.. literalinclude:: reconstruction.py
+    :language: python

include/compas.h

@@ -8,31 +8,24 @@
 
 #include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
 #include <CGAL/Surface_mesh.h>
+#include <CGAL/Polyhedron_3.h>
+#include <CGAL/Polyhedron_incremental_builder_3.h>
+#include <CGAL/Polyhedron_items_with_id_3.h>
 
 using Kernel = CGAL::Exact_predicates_inexact_constructions_kernel;
-using Polyline = std::vector<Kernel::Point_3>;
+using Point = Kernel::Point_3;
+using Vector = Kernel::Vector_3;
+using Polyline = std::vector<Point>;
 using Polylines = std::list<Polyline>;
-using Mesh = CGAL::Surface_mesh<Kernel::Point_3>;
-
-// How to define and work with Polyhedra?
-// Define and register structured output.
+using Polyhedron = CGAL::Polyhedron_3<Kernel, CGAL::Polyhedron_items_with_id_3>;
+using Mesh = CGAL::Surface_mesh<Point>;
 
 namespace compas
 {
     using RowMatrixXd = Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
     using RowMatrixXi = Eigen::Matrix<int, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
 
-    struct Point
-    {
-        double x;
-        double y;
-        double z;
-    };
-
-    // Polyhedron
-    // polyhedron_from_vertices_and_faces(
-    //     const RowMatrixXd &V,
-    //     const RowMatrixXi &F);
+    Polyhedron polyhedron_from_vertices_and_faces(const RowMatrixXd &V, const RowMatrixXi &F);
 
     Mesh mesh_from_vertices_and_faces(const RowMatrixXd &V, const RowMatrixXi &F);
 
@@ -43,6 +36,8 @@ namespace compas
     std::tuple<compas::RowMatrixXd, compas::RowMatrixXi> quadmesh_to_vertices_and_faces(const Mesh &mesh);
 
     std::vector<compas::RowMatrixXd> polylines_to_lists_of_points(Polylines polylines);
+
+    std::tuple<compas::RowMatrixXd, compas::RowMatrixXi> polyhedron_to_vertices_and_faces(Polyhedron polyhedron);
 }
 
 #endif /* COMPAS_H */

setup.py

@@ -63,6 +63,7 @@ def get_library_dirs():
                 "src/subdivision.cpp",
                 "src/triangulations.cpp",
                 "src/skeletonization.cpp",
+                "src/reconstruction.cpp",
             ]
         ),
         include_dirs=["./include", get_eigen_include(), get_pybind_include()],

src/compas.cpp

@@ -1,6 +1,56 @@
 #include <compas.h>
 #include <pybind11/pybind11.h>
 
+template <class HDS>
+class Build_polyhedron : public CGAL::Modifier_base<HDS>
+{
+private:
+    const compas::RowMatrixXd &m_vertices;
+    const compas::RowMatrixXi &m_faces;
+
+public:
+    Build_polyhedron(const compas::RowMatrixXd &V, const compas::RowMatrixXi &F) : m_vertices(V), m_faces(F)
+    {
+    }
+
+    void operator()(HDS &hds)
+    {
+        typedef typename HDS::Vertex Vertex;
+        typedef typename Vertex::Point Point;
+
+        CGAL::Polyhedron_incremental_builder_3<HDS> Builder(hds, true);
+
+        Builder.begin_surface(m_vertices.rows(), m_faces.rows());
+
+        for (int i = 0; i < m_vertices.rows(); i++)
+        {
+            Builder.add_vertex(Point(m_vertices(i, 0), m_vertices(i, 1), m_vertices(i, 2)));
+        }
+
+        for (int i = 0; i < m_faces.rows(); i++)
+        {
+            Builder.begin_facet();
+            for (int j = 0; j < m_faces.cols(); j++)
+            {
+                Builder.add_vertex_to_facet(m_faces(i, j));
+            }
+            Builder.end_facet();
+        }
+
+        Builder.end_surface();
+    }
+};
+
+Polyhedron polyhedron_from_vertices_and_faces(
+    const compas::RowMatrixXd &V,
+    const compas::RowMatrixXi &F)
+{
+    Polyhedron polyhedron;
+    Build_polyhedron<Polyhedron::HalfedgeDS> build(V, F);
+    polyhedron.delegate(build);
+    return polyhedron;
+}
+
 Mesh compas::mesh_from_vertices_and_faces(
     const compas::RowMatrixXd &V,
     const compas::RowMatrixXi &F)
@@ -103,13 +153,13 @@ compas::quadmesh_to_vertices_and_faces(
     compas::RowMatrixXd V(v, 3);
     compas::RowMatrixXi F(f, 4);
 
-    Mesh::Property_map<Mesh::Vertex_index, Kernel::Point_3> location = mesh.points();
+    Mesh::Property_map<Mesh::Vertex_index, Kernel::Point_3> vertex_location = mesh.points();
 
     for (Mesh::Vertex_index vd : mesh.vertices())
     {
-        V(vd, 0) = (double)location[vd][0];
-        V(vd, 1) = (double)location[vd][1];
-        V(vd, 2) = (double)location[vd][2];
+        V(vd, 0) = (double)vertex_location[vd][0];
+        V(vd, 1) = (double)vertex_location[vd][1];
+        V(vd, 2) = (double)vertex_location[vd][2];
     }
 
     for (Mesh::Face_index fd : mesh.faces())
@@ -151,17 +201,58 @@ compas::polylines_to_lists_of_points(
     return pointsets;
 }
 
-// Polyhedron compas::polyhedron_from_vertices_and_faces(
-//     const RowMatrixXd &V,
-//     const RowMatrixXi &F);
-// {
-//     int v = V.rows();
-//     int f = F.rows();
+std::tuple<compas::RowMatrixXd, compas::RowMatrixXi>
+compas::polyhedron_to_vertices_and_faces(
+    Polyhedron polyhedron)
+{
+    int v = polyhedron.size_of_vertices();
+    int f = polyhedron.size_of_facets();
+
+    compas::RowMatrixXd V(v, 3);
+    compas::RowMatrixXi F(f, 3);
+
+    std::size_t i = 0;
+
+    // for (auto vi = polyhedron.vertices_begin(); vi != polyhedron.vertices_end(); ++vi)
+    // {
+    //     V(i, 0) = (double)vi->point().x();
+    //     V(i, 1) = (double)vi->point().y();
+    //     V(i, 2) = (double)vi->point().z();
+
+    //     vi->id() = i;
+
+    //     i++;
+    // }
+
+    for (Polyhedron::Vertex_handle vh : polyhedron.vertex_handles())
+    {
+        V(i, 0) = (double)vh->point().x();
+        V(i, 1) = (double)vh->point().y();
+        V(i, 2) = (double)vh->point().z();
+
+        vh->id() = i;
+
+        i++;
+    }
+
+    i = 0;
+
+    for (Polyhedron::Facet_handle fh : polyhedron.facet_handles())
+    {
+        std::size_t j = 0;
+
+        Polyhedron::Halfedge_handle start = fh->halfedge(), h = start;
+        do
+        {
+            F(i, j) = h->vertex()->id();
+            h = h->next();
+            j++;
 
-//     Polyhedron polyhedron;
+        } while (h != start);
 
-//     polyhedron.begin_surface();
-//     polyhedron.end_surface();
+        i++;
+    }
 
-//     return polyhedron;
-// }
+    std::tuple<compas::RowMatrixXd, compas::RowMatrixXi> result = std::make_tuple(V, F);
+    return result;
+}

src/compas_cgal.cpp

@@ -10,6 +10,7 @@ void init_measure(pybind11::module &);
 void init_triangulations(pybind11::module &);
 void init_subdivision(pybind11::module &);
 void init_skeletonization(pybind11::module &);
+void init_reconstruction(pybind11::module &);
 
 // the first parameter here ("_cgal") will be the name of the "so" or "pyd" file that will be produced by PyBind11
 // which is the entry point from where all other modules will be accessible.
@@ -26,4 +27,5 @@ PYBIND11_MODULE(_cgal, m)
     init_triangulations(m);
     init_subdivision(m);
     init_skeletonization(m);
+    init_reconstruction(m);
 }

src/compas_cgal/reconstruction.py

@@ -0,0 +1,9 @@
+from __future__ import annotations
+import numpy as np
+from compas_cgal._cgal import reconstruction
+
+
+def poisson_surface_reconstruction(points, normals):
+    P = np.asarray(points, dtype=np.float64)
+    N = np.asarray(normals, dtype=np.float64)
+    return reconstruction.poisson_surface_reconstruction(P, N)

src/reconstruction.cpp

@@ -1,9 +1,37 @@
 #include "reconstruction.h"
+#include <CGAL/poisson_surface_reconstruction.h>
+
+typedef std::pair<Point, Vector> PwN;
 
 std::tuple<compas::RowMatrixXd, compas::RowMatrixXi>
 poisson_surface_reconstruction(
-    Eigen::Ref<const compas::RowMatrixXd> &P)
+    Eigen::Ref<const compas::RowMatrixXd> &P,
+    Eigen::Ref<const compas::RowMatrixXd> &N)
 {
+    Polyhedron mesh;
+    std::vector<PwN> points;
+
+    for (int i = 0; i < P.rows(); i++)
+    {
+        points.push_back(PwN(
+            Point(P(i, 0), P(i, 1), P(i, 2)),
+            Vector(N(i, 0), N(i, 1), N(i, 2))));
+    }
+
+    double average_spacing = CGAL::compute_average_spacing<CGAL::Sequential_tag>(
+        points,
+        6,
+        CGAL::parameters::point_map(CGAL::First_of_pair_property_map<PwN>()));
+
+    CGAL::poisson_surface_reconstruction_delaunay(
+        points.begin(),
+        points.end(),
+        CGAL::First_of_pair_property_map<PwN>(),
+        CGAL::Second_of_pair_property_map<PwN>(),
+        mesh,
+        average_spacing);
+
+    return compas::polyhedron_to_vertices_and_faces(mesh);
 }
 
 // ===========================================================================
@@ -17,5 +45,6 @@ void init_reconstruction(pybind11::module &m)
     submodule.def(
         "poisson_surface_reconstruction",
         &poisson_surface_reconstruction,
-        pybind11::arg("P").noconvert());
+        pybind11::arg("P").noconvert(),
+        pybind11::arg("N").noconvert());
 };

src/reconstruction.h

@@ -5,7 +5,7 @@
 
 std::tuple<compas::RowMatrixXd, compas::RowMatrixXi>
 poisson_surface_reconstruction(
-    Eigen::Ref<const compas::RowMatrixXd> &V,
-    Eigen::Ref<const compas::RowMatrixXi> &F);
+    Eigen::Ref<const compas::RowMatrixXd> &P,
+    Eigen::Ref<const compas::RowMatrixXd> &N);
 
 #endif /* COMPAS_RECONSTRUCTION_H */