+ add Poisson surface reconstruction

src/Mod/ReverseEngineering/App/SurfaceTriangulation.cpp

@@ -38,6 +38,7 @@
 #include <pcl/point_traits.h>
 #include <pcl/surface/gp3.h>
 #include <pcl/surface/grid_projection.h>
+#include <pcl/surface/poisson.h>
 //#include <pcl/surface/convex_hull.h>
 //#include <pcl/surface/concave_hull.h>
 #include <pcl/surface/organized_fast_mesh.h>
@@ -108,11 +109,86 @@ void SurfaceTriangulation::perform(double searchRadius, double mu)
     PolygonMesh mesh;
     gp3.reconstruct (mesh);
 
+    MeshConversion::convert(mesh, myMesh);
+
+    // Additional vertex information
+    //std::vector<int> parts = gp3.getPartIDs();
+    //std::vector<int> states = gp3.getPointStates();
+}
+
+// ----------------------------------------------------------------------------
+
+// See
+// http://www.cs.jhu.edu/~misha/Code/PoissonRecon/Version8.0/
+PoissonReconstruction::PoissonReconstruction(const Points::PointKernel& pts, Mesh::MeshObject& mesh)
+  : myPoints(pts)
+  , myMesh(mesh)
+  , depth(-1)
+  , solverDivide(-1)
+  , samplesPerNode(-1.0f)
+{
+}
+
+void PoissonReconstruction::perform(int ksearch)
+{
+    PointCloud<PointXYZ>::Ptr cloud (new PointCloud<PointXYZ>);
+    PointCloud<PointNormal>::Ptr cloud_with_normals (new PointCloud<PointNormal>);
+    search::KdTree<PointXYZ>::Ptr tree;
+    search::KdTree<PointNormal>::Ptr tree2;
+
+    for (Points::PointKernel::const_iterator it = myPoints.begin(); it != myPoints.end(); ++it) {
+        cloud->push_back(PointXYZ(it->x, it->y, it->z));
+    }
+
+    // Create search tree
+    tree.reset (new search::KdTree<PointXYZ> (false));
+    tree->setInputCloud (cloud);
+
+    // Normal estimation
+    NormalEstimation<PointXYZ, Normal> n;
+    PointCloud<Normal>::Ptr normals (new PointCloud<Normal> ());
+    n.setInputCloud (cloud);
+    //n.setIndices (indices[B);
+    n.setSearchMethod (tree);
+    n.setKSearch (ksearch);
+    n.compute (*normals);
+
+    // Concatenate XYZ and normal information
+    pcl::concatenateFields (*cloud, *normals, *cloud_with_normals);
+      
+    // Create search tree
+    tree2.reset (new search::KdTree<PointNormal>);
+    tree2->setInputCloud (cloud_with_normals);
+
+    // Init objects
+    Poisson<PointNormal> poisson;
+
+    // Set parameters
+    poisson.setInputCloud (cloud_with_normals);
+    poisson.setSearchMethod (tree2);
+    if (depth >= 1)
+        poisson.setDepth(depth);
+    if (solverDivide >= 1)
+        poisson.setSolverDivide(solverDivide);
+    if (samplesPerNode >= 1.0f)
+        poisson.setSamplesPerNode(samplesPerNode);
+
+    // Reconstruct
+    PolygonMesh mesh;
+    poisson.reconstruct (mesh);
+
+    MeshConversion::convert(mesh, myMesh);
+}
+
+// ----------------------------------------------------------------------------
+
+void MeshConversion::convert(const pcl::PolygonMesh& pclMesh, Mesh::MeshObject& meshObject)
+{
     // number of points
-    size_t nr_points  = mesh.cloud.width * mesh.cloud.height;
-    size_t point_size = mesh.cloud.data.size () / nr_points;
+    size_t nr_points  = pclMesh.cloud.width * pclMesh.cloud.height;
+    size_t point_size = pclMesh.cloud.data.size () / nr_points;
     // number of faces for header
-    size_t nr_faces = mesh.polygons.size ();
+    size_t nr_faces = pclMesh.polygons.size ();
 
     MeshCore::MeshPointArray points;
     points.reserve(nr_points);
@@ -123,21 +199,21 @@ void SurfaceTriangulation::perform(double searchRadius, double mu)
     MeshCore::MeshPoint vertex;
     for (size_t i = 0; i < nr_points; ++i) {
         int xyz = 0;
-        for (size_t d = 0; d < mesh.cloud.fields.size(); ++d) {
+        for (size_t d = 0; d < pclMesh.cloud.fields.size(); ++d) {
             int c = 0;
             // adding vertex
-            if ((mesh.cloud.fields[d].datatype ==
+            if ((pclMesh.cloud.fields[d].datatype ==
 #if PCL_VERSION_COMPARE(>,1,6,0)
                  pcl::PCLPointField::FLOAT32) &&
 #else
                  sensor_msgs::PointField::FLOAT32) &&
 #endif
-                (mesh.cloud.fields[d].name == "x" ||
-                 mesh.cloud.fields[d].name == "y" ||
-                 mesh.cloud.fields[d].name == "z"))
+                (pclMesh.cloud.fields[d].name == "x" ||
+                 pclMesh.cloud.fields[d].name == "y" ||
+                 pclMesh.cloud.fields[d].name == "z"))
             {
                 float value;
-                memcpy (&value, &mesh.cloud.data[i * point_size + mesh.cloud.fields[d].offset + c * sizeof (float)], sizeof (float));
+                memcpy (&value, &pclMesh.cloud.data[i * point_size + pclMesh.cloud.fields[d].offset + c * sizeof (float)], sizeof (float));
                 vertex[xyz] = value;
                 if (++xyz == 3) {
                     points.push_back(vertex);
@@ -149,20 +225,16 @@ void SurfaceTriangulation::perform(double searchRadius, double mu)
     // get faces
     MeshCore::MeshFacet face;
     for (size_t i = 0; i < nr_faces; i++) {
-        face._aulPoints[0] = mesh.polygons[i].vertices[0];
-        face._aulPoints[1] = mesh.polygons[i].vertices[1];
-        face._aulPoints[2] = mesh.polygons[i].vertices[2];
+        face._aulPoints[0] = pclMesh.polygons[i].vertices[0];
+        face._aulPoints[1] = pclMesh.polygons[i].vertices[1];
+        face._aulPoints[2] = pclMesh.polygons[i].vertices[2];
         facets.push_back(face);
     }
 
     MeshCore::MeshKernel kernel;
     kernel.Adopt(points, facets, true);
-    myMesh.swap(kernel);
-    myMesh.harmonizeNormals();
-
-    // Additional vertex information
-    //std::vector<int> parts = gp3.getPartIDs();
-    //std::vector<int> states = gp3.getPointStates();
+    meshObject.swap(kernel);
+    meshObject.harmonizeNormals();
 }
 
 #endif // HAVE_PCL_SURFACE