+ add several surface reconstruction methods from pcl to Reen module

src/Mod/ReverseEngineering/App/AppReverseEngineering.cpp

@@ -59,12 +59,24 @@ public:
             "Iterations=5,Correction=True,PatchFactor=1.0"
         );
 #if defined(HAVE_PCL_SURFACE)
-        add_varargs_method("triangulate",&Module::triangulate,
+        add_keyword_method("triangulate",&Module::triangulate,
             "triangulate(PointKernel,searchRadius[,mu=2.5])."
         );
         add_keyword_method("poissonReconstruction",&Module::poissonReconstruction,
             "poissonReconstruction(PointKernel)."
         );
+        add_keyword_method("viewTriangulation",&Module::viewTriangulation,
+            "viewTriangulation(PointKernel, width, height)."
+        );
+        add_keyword_method("gridProjection",&Module::gridProjection,
+            "gridProjection(PointKernel)."
+        );
+        add_keyword_method("marchingCubesRBF",&Module::marchingCubesRBF,
+            "marchingCubesRBF(PointKernel)."
+        );
+        add_keyword_method("marchingCubesHoppe",&Module::marchingCubesHoppe,
+            "marchingCubesHoppe(PointKernel)."
+        );
 #endif
 #if defined(HAVE_PCL_OPENNURBS)
         add_keyword_method("fitBSpline",&Module::fitBSpline,
@@ -199,47 +211,270 @@ private:
         }
     }
 #if defined(HAVE_PCL_SURFACE)
-    Py::Object triangulate(const Py::Tuple& args)
+    /*
+import ReverseEngineering as Reen
+import Points
+import Mesh
+import random
+
+r=random.Random()
+
+p=Points.Points()
+pts=[]
+for i in range(21):
+  for j in range(21):
+    pts.append(App.Vector(i,j,r.gauss(5,0.05)))
+
+p.addPoints(pts)
+m=Reen.triangulate(Points=p,SearchRadius=2.2)
+Mesh.show(m)
+    */
+    Py::Object triangulate(const Py::Tuple& args, const Py::Dict& kwds)
     {
-        PyObject *pcObj;
+        PyObject *pts;
         double searchRadius;
+        PyObject *vec = 0;
+        int ksearch=5;
         double mu=2.5;
-        if (!PyArg_ParseTuple(args.ptr(), "O!d|d", &(Points::PointsPy::Type), &pcObj, &searchRadius, &mu))
+
+        static char* kwds_greedy[] = {"Points", "SearchRadius", "Mu", "KSearch",
+                                      "Normals", NULL};
+        if (!PyArg_ParseTupleAndKeywords(args.ptr(), kwds.ptr(), "O!d|diO", kwds_greedy,
+                                        &(Points::PointsPy::Type), &pts,
+                                        &searchRadius, &mu, &ksearch, &vec))
             throw Py::Exception();
 
-        Points::PointsPy* pPoints = static_cast<Points::PointsPy*>(pcObj);
+        Points::PointKernel* points = static_cast<Points::PointsPy*>(pts)->getPointKernelPtr();
-        Points::PointKernel* points = pPoints->getPointKernelPtr();
+
-        
         Mesh::MeshObject* mesh = new Mesh::MeshObject();
         SurfaceTriangulation tria(*points, *mesh);
-        tria.perform(searchRadius, mu);
+        tria.setMu(mu);
+        tria.setSearchRadius(searchRadius);
+        if (vec) {
+            Py::Sequence list(vec);
+            std::vector<Base::Vector3f> normals;
+            normals.reserve(list.size());
+            for (Py::Sequence::iterator it = list.begin(); it != list.end(); ++it) {
+                Base::Vector3d v = Py::Vector(*it).toVector();
+                normals.push_back(Base::convertTo<Base::Vector3f>(v));
+            }
+            tria.perform(normals);
+        }
+        else {
+            tria.perform(ksearch);
+        }
 
         return Py::asObject(new Mesh::MeshPy(mesh));
     }
     Py::Object poissonReconstruction(const Py::Tuple& args, const Py::Dict& kwds)
     {
-        PyObject *pcObj;
+        PyObject *pts;
+        PyObject *vec = 0;
         int ksearch=5;
         int octreeDepth=-1;
         int solverDivide=-1;
         double samplesPerNode=-1.0;
 
         static char* kwds_poisson[] = {"Points", "KSearch", "OctreeDepth", "SolverDivide",
-                                      "SamplesPerNode", NULL};
+                                      "SamplesPerNode", "Normals", NULL};
-        if (!PyArg_ParseTupleAndKeywords(args.ptr(), kwds.ptr(), "O!|iiid", kwds_poisson,
+        if (!PyArg_ParseTupleAndKeywords(args.ptr(), kwds.ptr(), "O!|iiidO", kwds_poisson,
-                                        &(Points::PointsPy::Type), &pcObj,
+                                        &(Points::PointsPy::Type), &pts,
-                                        &ksearch, &octreeDepth, &solverDivide, &samplesPerNode))
+                                        &ksearch, &octreeDepth, &solverDivide, &samplesPerNode, &vec))
             throw Py::Exception();
 
-        Points::PointsPy* pPoints = static_cast<Points::PointsPy*>(pcObj);
+        Points::PointKernel* points = static_cast<Points::PointsPy*>(pts)->getPointKernelPtr();
-        Points::PointKernel* points = pPoints->getPointKernelPtr();
 
         Mesh::MeshObject* mesh = new Mesh::MeshObject();
         Reen::PoissonReconstruction poisson(*points, *mesh);
         poisson.setDepth(octreeDepth);
         poisson.setSolverDivide(solverDivide);
         poisson.setSamplesPerNode(samplesPerNode);
-        poisson.perform(ksearch);
+        if (vec) {
+            Py::Sequence list(vec);
+            std::vector<Base::Vector3f> normals;
+            normals.reserve(list.size());
+            for (Py::Sequence::iterator it = list.begin(); it != list.end(); ++it) {
+                Base::Vector3d v = Py::Vector(*it).toVector();
+                normals.push_back(Base::convertTo<Base::Vector3f>(v));
+            }
+            poisson.perform(normals);
+        }
+        else {
+            poisson.perform(ksearch);
+        }
+
+        return Py::asObject(new Mesh::MeshPy(mesh));
+    }
+    /*
+import ReverseEngineering as Reen
+import Points
+import Mesh
+import random
+import math
+
+r=random.Random()
+
+p=Points.Points()
+pts=[]
+for i in range(21):
+  for j in range(21):
+    pts.append(App.Vector(i,j,r.random()))
+
+p.addPoints(pts)
+m=Reen.viewTriangulation(p,21,21)
+Mesh.show(m)
+
+def boxmueller():
+  r1,r2=random.random(),random.random()
+  return math.sqrt(-2*math.log(r1))*math.cos(2*math.pi*r2)
+
+p=Points.Points()
+pts=[]
+for i in range(21):
+  for j in range(21):
+    pts.append(App.Vector(i,j,r.gauss(5,0.05)))
+
+p.addPoints(pts)
+m=Reen.viewTriangulation(p,21,21)
+Mesh.show(m)
+    */
+    Py::Object viewTriangulation(const Py::Tuple& args, const Py::Dict& kwds)
+    {
+        PyObject *pts;
+        PyObject *vec = 0;
+        int width;
+        int height;
+
+        static char* kwds_greedy[] = {"Points", "Width", "Height", NULL};
+        if (!PyArg_ParseTupleAndKeywords(args.ptr(), kwds.ptr(), "O!|ii", kwds_greedy,
+                                        &(Points::PointsPy::Type), &pts,
+                                        &width, &height))
+            throw Py::Exception();
+
+        Points::PointKernel* points = static_cast<Points::PointsPy*>(pts)->getPointKernelPtr();
+
+        try {
+            Mesh::MeshObject* mesh = new Mesh::MeshObject();
+            ImageTriangulation view(width, height, *points, *mesh);
+            view.perform();
+
+            return Py::asObject(new Mesh::MeshPy(mesh));
+        }
+        catch (const Base::Exception& e) {
+            throw Py::RuntimeError(e.what());
+        }
+    }
+    Py::Object gridProjection(const Py::Tuple& args, const Py::Dict& kwds)
+    {
+        PyObject *pts;
+        PyObject *vec = 0;
+        int ksearch=5;
+
+        static char* kwds_greedy[] = {"Points", "KSearch", "Normals", NULL};
+        if (!PyArg_ParseTupleAndKeywords(args.ptr(), kwds.ptr(), "O!|iO", kwds_greedy,
+                                        &(Points::PointsPy::Type), &pts,
+                                        &ksearch, &vec))
+            throw Py::Exception();
+
+        Points::PointKernel* points = static_cast<Points::PointsPy*>(pts)->getPointKernelPtr();
+
+        Mesh::MeshObject* mesh = new Mesh::MeshObject();
+        GridReconstruction tria(*points, *mesh);
+        if (vec) {
+            Py::Sequence list(vec);
+            std::vector<Base::Vector3f> normals;
+            normals.reserve(list.size());
+            for (Py::Sequence::iterator it = list.begin(); it != list.end(); ++it) {
+                Base::Vector3d v = Py::Vector(*it).toVector();
+                normals.push_back(Base::convertTo<Base::Vector3f>(v));
+            }
+            tria.perform(normals);
+        }
+        else {
+            tria.perform(ksearch);
+        }
+
+        return Py::asObject(new Mesh::MeshPy(mesh));
+    }
+    Py::Object marchingCubesRBF(const Py::Tuple& args, const Py::Dict& kwds)
+    {
+        PyObject *pts;
+        PyObject *vec = 0;
+        int ksearch=5;
+
+        static char* kwds_greedy[] = {"Points", "KSearch", "Normals", NULL};
+        if (!PyArg_ParseTupleAndKeywords(args.ptr(), kwds.ptr(), "O!|iO", kwds_greedy,
+                                        &(Points::PointsPy::Type), &pts,
+                                        &ksearch, &vec))
+            throw Py::Exception();
+
+        Points::PointKernel* points = static_cast<Points::PointsPy*>(pts)->getPointKernelPtr();
+
+        Mesh::MeshObject* mesh = new Mesh::MeshObject();
+        MarchingCubesRBF tria(*points, *mesh);
+        if (vec) {
+            Py::Sequence list(vec);
+            std::vector<Base::Vector3f> normals;
+            normals.reserve(list.size());
+            for (Py::Sequence::iterator it = list.begin(); it != list.end(); ++it) {
+                Base::Vector3d v = Py::Vector(*it).toVector();
+                normals.push_back(Base::convertTo<Base::Vector3f>(v));
+            }
+            tria.perform(normals);
+        }
+        else {
+            tria.perform(ksearch);
+        }
+
+        return Py::asObject(new Mesh::MeshPy(mesh));
+    }
+    /*
+import ReverseEngineering as Reen
+import Points
+import Mesh
+import random
+
+r=random.Random()
+
+p=Points.Points()
+pts=[]
+for i in range(21):
+  for j in range(21):
+    pts.append(App.Vector(i,j,r.gauss(5,0.05)))
+
+p.addPoints(pts)
+m=Reen.marchingCubesHoppe(Points=p)
+Mesh.show(m)
+    */
+    Py::Object marchingCubesHoppe(const Py::Tuple& args, const Py::Dict& kwds)
+    {
+        PyObject *pts;
+        PyObject *vec = 0;
+        int ksearch=5;
+
+        static char* kwds_greedy[] = {"Points", "KSearch", "Normals", NULL};
+        if (!PyArg_ParseTupleAndKeywords(args.ptr(), kwds.ptr(), "O!|iO", kwds_greedy,
+                                        &(Points::PointsPy::Type), &pts,
+                                        &ksearch, &vec))
+            throw Py::Exception();
+
+        Points::PointKernel* points = static_cast<Points::PointsPy*>(pts)->getPointKernelPtr();
+
+        Mesh::MeshObject* mesh = new Mesh::MeshObject();
+        MarchingCubesHoppe tria(*points, *mesh);
+        if (vec) {
+            Py::Sequence list(vec);
+            std::vector<Base::Vector3f> normals;
+            normals.reserve(list.size());
+            for (Py::Sequence::iterator it = list.begin(); it != list.end(); ++it) {
+                Base::Vector3d v = Py::Vector(*it).toVector();
+                normals.push_back(Base::convertTo<Base::Vector3f>(v));
+            }
+            tria.perform(normals);
+        }
+        else {
+            tria.perform(ksearch);
+        }
 
         return Py::asObject(new Mesh::MeshPy(mesh));
     }
@@ -247,7 +482,7 @@ private:
 #if defined(HAVE_PCL_OPENNURBS)
     Py::Object fitBSpline(const Py::Tuple& args, const Py::Dict& kwds)
     {
-        PyObject *pcObj;
+        PyObject *pts;
         int degree = 2;
         int refinement = 4;
         int iterations = 10;
@@ -259,14 +494,13 @@ private:
         static char* kwds_approx[] = {"Points", "Degree", "Refinement", "Iterations",
                                       "InteriorSmoothness", "InteriorWeight", "BoundarySmoothness", "BoundaryWeight", NULL};
         if (!PyArg_ParseTupleAndKeywords(args.ptr(), kwds.ptr(), "O!|iiidddd", kwds_approx,
-                                        &(Points::PointsPy::Type), &pcObj,
+                                        &(Points::PointsPy::Type), &pts,
                                         &degree, &refinement, &iterations,
                                         &interiorSmoothness, &interiorWeight,
                                         &boundarySmoothness, &boundaryWeight))
             throw Py::Exception();
 
-        Points::PointsPy* pPoints = static_cast<Points::PointsPy*>(pcObj);
+        Points::PointKernel* points = static_cast<Points::PointsPy*>(pts)->getPointKernelPtr();
-        Points::PointKernel* points = pPoints->getPointKernelPtr();
 
         BSplineFitting fit(points->getBasicPoints());
         fit.setOrder(degree+1);

src/Mod/ReverseEngineering/App/SurfaceTriangulation.cpp

@@ -26,8 +26,10 @@
 #include "SurfaceTriangulation.h"
 #include <Mod/Points/App/Points.h>
 #include <Mod/Mesh/App/Mesh.h>
+#include <Mod/Mesh/App/Core/Algorithm.h>
 #include <Mod/Mesh/App/Core/Elements.h>
 #include <Mod/Mesh/App/Core/MeshKernel.h>
+#include <Base/Exception.h>
 
 // http://svn.pointclouds.org/pcl/tags/pcl-1.5.1/test/
 #if defined(HAVE_PCL_SURFACE)
@@ -42,6 +44,8 @@
 //#include <pcl/surface/convex_hull.h>
 //#include <pcl/surface/concave_hull.h>
 #include <pcl/surface/organized_fast_mesh.h>
+#include <pcl/surface/marching_cubes_rbf.h>
+#include <pcl/surface/marching_cubes_hoppe.h>
 #include <pcl/surface/ear_clipping.h>
 #include <pcl/common/common.h>
 #include <boost/random.hpp>
@@ -55,18 +59,25 @@ using namespace pcl::io;
 using namespace std;
 using namespace Reen;
 
+// See
+// http://www.ics.uci.edu/~gopi/PAPERS/Euro00.pdf
+// http://www.ics.uci.edu/~gopi/PAPERS/CGMV.pdf
 SurfaceTriangulation::SurfaceTriangulation(const Points::PointKernel& pts, Mesh::MeshObject& mesh)
-  : myPoints(pts), myMesh(mesh)
+  : myPoints(pts)
+  , myMesh(mesh)
+  , mu(0)
+  , searchRadius(0)
 {
 }
 
-void SurfaceTriangulation::perform(double searchRadius, double mu)
+void SurfaceTriangulation::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;
-    
+
+    cloud->reserve(myPoints.size());
     for (Points::PointKernel::const_iterator it = myPoints.begin(); it != myPoints.end(); ++it) {
         cloud->push_back(PointXYZ(it->x, it->y, it->z));
     }
@@ -81,12 +92,12 @@ void SurfaceTriangulation::perform(double searchRadius, double mu)
     n.setInputCloud (cloud);
     //n.setIndices (indices[B);
     n.setSearchMethod (tree);
-    n.setKSearch (20);
+    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);
@@ -104,6 +115,61 @@ void SurfaceTriangulation::perform(double searchRadius, double mu)
     gp3.setMinimumAngle(M_PI/18); // 10 degrees
     gp3.setMaximumAngle(2*M_PI/3); // 120 degrees
     gp3.setNormalConsistency(false);
+    gp3.setConsistentVertexOrdering(true);
+
+    // Reconstruct
+    PolygonMesh mesh;
+    gp3.reconstruct (mesh);
+
+    MeshConversion::convert(mesh, myMesh);
+
+    // Additional vertex information
+    //std::vector<int> parts = gp3.getPartIDs();
+    //std::vector<int> states = gp3.getPointStates();
+}
+
+void SurfaceTriangulation::perform(const std::vector<Base::Vector3f>& normals)
+{
+    if (myPoints.size() != normals.size())
+        throw Base::RuntimeError("Number of points doesn't match with number of normals");
+
+    PointCloud<PointNormal>::Ptr cloud_with_normals (new PointCloud<PointNormal>);
+    search::KdTree<PointNormal>::Ptr tree;
+
+    cloud_with_normals->reserve(myPoints.size());
+    std::size_t num_points = myPoints.size();
+    const std::vector<Base::Vector3f>& points = myPoints.getBasicPoints();
+    for (std::size_t index=0; index<num_points; index++) {
+        const Base::Vector3f& p = points[index];
+        const Base::Vector3f& n = normals[index];
+        PointNormal pn;
+        pn.x = p.x;
+        pn.y = p.y;
+        pn.z = p.z;
+        pn.normal_x = n.x;
+        pn.normal_y = n.y;
+        pn.normal_z = n.z;
+        cloud_with_normals->push_back(pn);
+    }
+
+    // Create search tree
+    tree.reset (new search::KdTree<PointNormal>);
+    tree->setInputCloud (cloud_with_normals);
+
+    // Init objects
+    GreedyProjectionTriangulation<PointNormal> gp3;
+
+    // Set parameters
+    gp3.setInputCloud (cloud_with_normals);
+    gp3.setSearchMethod (tree);
+    gp3.setSearchRadius (searchRadius);
+    gp3.setMu (mu);
+    gp3.setMaximumNearestNeighbors (100);
+    gp3.setMaximumSurfaceAngle(M_PI/4); // 45 degrees
+    gp3.setMinimumAngle(M_PI/18); // 10 degrees
+    gp3.setMaximumAngle(2*M_PI/3); // 120 degrees
+    gp3.setNormalConsistency(true);
+    gp3.setConsistentVertexOrdering(true);
 
     // Reconstruct
     PolygonMesh mesh;
@@ -136,6 +202,7 @@ void PoissonReconstruction::perform(int ksearch)
     search::KdTree<PointXYZ>::Ptr tree;
     search::KdTree<PointNormal>::Ptr tree2;
 
+    cloud->reserve(myPoints.size());
     for (Points::PointKernel::const_iterator it = myPoints.begin(); it != myPoints.end(); ++it) {
         cloud->push_back(PointXYZ(it->x, it->y, it->z));
     }
@@ -155,7 +222,7 @@ void PoissonReconstruction::perform(int ksearch)
 
     // 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);
@@ -180,6 +247,449 @@ void PoissonReconstruction::perform(int ksearch)
     MeshConversion::convert(mesh, myMesh);
 }
 
+void PoissonReconstruction::perform(const std::vector<Base::Vector3f>& normals)
+{
+    if (myPoints.size() != normals.size())
+        throw Base::RuntimeError("Number of points doesn't match with number of normals");
+
+    PointCloud<PointNormal>::Ptr cloud_with_normals (new PointCloud<PointNormal>);
+    search::KdTree<PointNormal>::Ptr tree;
+
+    cloud_with_normals->reserve(myPoints.size());
+    std::size_t num_points = myPoints.size();
+    const std::vector<Base::Vector3f>& points = myPoints.getBasicPoints();
+    for (std::size_t index=0; index<num_points; index++) {
+        const Base::Vector3f& p = points[index];
+        const Base::Vector3f& n = normals[index];
+        PointNormal pn;
+        pn.x = p.x;
+        pn.y = p.y;
+        pn.z = p.z;
+        pn.normal_x = n.x;
+        pn.normal_y = n.y;
+        pn.normal_z = n.z;
+        cloud_with_normals->push_back(pn);
+    }
+
+    // Create search tree
+    tree.reset (new search::KdTree<PointNormal>);
+    tree->setInputCloud (cloud_with_normals);
+
+    // Init objects
+    Poisson<PointNormal> poisson;
+
+    // Set parameters
+    poisson.setInputCloud (cloud_with_normals);
+    poisson.setSearchMethod (tree);
+    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);
+}
+
+// ----------------------------------------------------------------------------
+
+GridReconstruction::GridReconstruction(const Points::PointKernel& pts, Mesh::MeshObject& mesh)
+  : myPoints(pts)
+  , myMesh(mesh)
+{
+}
+
+void GridReconstruction::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;
+
+    cloud->reserve(myPoints.size());
+    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
+    GridProjection<PointNormal> grid;
+
+    // Set parameters
+    grid.setResolution(0.005); 
+    grid.setPaddingSize(3); 
+    grid.setNearestNeighborNum(100); 
+    grid.setMaxBinarySearchLevel(10);
+    grid.setInputCloud (cloud_with_normals);
+    grid.setSearchMethod (tree2);
+
+    // Reconstruct
+    PolygonMesh mesh;
+    grid.reconstruct (mesh);
+
+    MeshConversion::convert(mesh, myMesh);
+}
+
+void GridReconstruction::perform(const std::vector<Base::Vector3f>& normals)
+{
+    if (myPoints.size() != normals.size())
+        throw Base::RuntimeError("Number of points doesn't match with number of normals");
+
+    PointCloud<PointNormal>::Ptr cloud_with_normals (new PointCloud<PointNormal>);
+    search::KdTree<PointNormal>::Ptr tree;
+
+    cloud_with_normals->reserve(myPoints.size());
+    std::size_t num_points = myPoints.size();
+    const std::vector<Base::Vector3f>& points = myPoints.getBasicPoints();
+    for (std::size_t index=0; index<num_points; index++) {
+        const Base::Vector3f& p = points[index];
+        const Base::Vector3f& n = normals[index];
+        PointNormal pn;
+        pn.x = p.x;
+        pn.y = p.y;
+        pn.z = p.z;
+        pn.normal_x = n.x;
+        pn.normal_y = n.y;
+        pn.normal_z = n.z;
+        cloud_with_normals->push_back(pn);
+    }
+
+    // Create search tree
+    tree.reset (new search::KdTree<PointNormal>);
+    tree->setInputCloud (cloud_with_normals);
+
+    // Init objects
+    GridProjection<PointNormal> grid;
+
+    // Set parameters
+    grid.setResolution(0.005); 
+    grid.setPaddingSize(3); 
+    grid.setNearestNeighborNum(100); 
+    grid.setMaxBinarySearchLevel(10);
+    grid.setInputCloud (cloud_with_normals);
+    grid.setSearchMethod (tree);
+
+    // Reconstruct
+    PolygonMesh mesh;
+    grid.reconstruct (mesh);
+
+    MeshConversion::convert(mesh, myMesh);
+}
+
+// ----------------------------------------------------------------------------
+
+ImageTriangulation::ImageTriangulation(int width, int height, const Points::PointKernel& pts, Mesh::MeshObject& mesh)
+  : width(width)
+  , height(height)
+  , myPoints(pts)
+  , myMesh(mesh)
+{
+}
+
+void ImageTriangulation::perform()
+{
+    if (myPoints.size() != width * height)
+        throw Base::RuntimeError("Number of points doesn't match with given width and height");
+
+    //construct dataset
+    pcl::PointCloud<pcl::PointXYZ>::Ptr cloud_organized (new pcl::PointCloud<pcl::PointXYZ> ());
+    cloud_organized->width = width;
+    cloud_organized->height = height;
+    cloud_organized->points.resize (cloud_organized->width * cloud_organized->height);
+
+    const std::vector<Base::Vector3f>& points = myPoints.getBasicPoints();
+
+    int npoints = 0;
+    for (size_t i = 0; i < cloud_organized->height; i++) {
+        for (size_t j = 0; j < cloud_organized->width; j++) {
+            const Base::Vector3f& p = points[npoints];
+            cloud_organized->points[npoints].x = p.x;
+            cloud_organized->points[npoints].y = p.y;
+            cloud_organized->points[npoints].z = p.z;
+            npoints++;
+        }
+    }
+
+    OrganizedFastMesh<PointXYZ> ofm;
+
+    // Set parameters
+    ofm.setInputCloud (cloud_organized);
+    // This parameter is not yet implmented (pcl 1.7)
+    ofm.setMaxEdgeLength (1.5);
+    ofm.setTrianglePixelSize (1);
+    ofm.setTriangulationType (OrganizedFastMesh<PointXYZ>::TRIANGLE_ADAPTIVE_CUT);
+    ofm.storeShadowedFaces(true);
+
+    // Reconstruct
+    PolygonMesh mesh;
+    ofm.reconstruct (mesh);
+
+    MeshConversion::convert(mesh, myMesh);
+
+    // remove invalid points
+    //
+    MeshCore::MeshKernel& kernel = myMesh.getKernel();
+    const MeshCore::MeshFacetArray& face = kernel.GetFacets();
+    MeshCore::MeshAlgorithm meshAlg(kernel);
+    meshAlg.SetPointFlag(MeshCore::MeshPoint::INVALID);
+    std::vector<unsigned long> validPoints;
+    validPoints.reserve(face.size()*3);
+    for (MeshCore::MeshFacetArray::_TConstIterator it = face.begin(); it != face.end(); ++it) {
+        validPoints.push_back(it->_aulPoints[0]);
+        validPoints.push_back(it->_aulPoints[1]);
+        validPoints.push_back(it->_aulPoints[2]);
+    }
+
+    // remove duplicates
+    std::sort(validPoints.begin(), validPoints.end());
+    validPoints.erase(std::unique(validPoints.begin(), validPoints.end()), validPoints.end());
+    meshAlg.ResetPointsFlag(validPoints, MeshCore::MeshPoint::INVALID);
+
+    unsigned long countInvalid = meshAlg.CountPointFlag(MeshCore::MeshPoint::INVALID);
+    if (countInvalid > 0) {
+        std::vector<unsigned long> invalidPoints;
+        invalidPoints.reserve(countInvalid);
+        meshAlg.GetPointsFlag(invalidPoints, MeshCore::MeshPoint::INVALID);
+
+        kernel.DeletePoints(invalidPoints);
+    }
+}
+
+// ----------------------------------------------------------------------------
+
+Reen::MarchingCubesRBF::MarchingCubesRBF(const Points::PointKernel& pts, Mesh::MeshObject& mesh)
+  : myPoints(pts)
+  , myMesh(mesh)
+{
+}
+
+void Reen::MarchingCubesRBF::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;
+
+    cloud->reserve(myPoints.size());
+    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
+    pcl::MarchingCubesRBF<PointNormal> rbf;
+
+    // Set parameters
+    rbf.setIsoLevel (0);
+    rbf.setGridResolution (60, 60, 60);
+    rbf.setPercentageExtendGrid (0.1f);
+    rbf.setOffSurfaceDisplacement (0.02f);
+
+    rbf.setInputCloud (cloud_with_normals);
+    rbf.setSearchMethod (tree2);
+
+    // Reconstruct
+    PolygonMesh mesh;
+    rbf.reconstruct (mesh);
+
+    MeshConversion::convert(mesh, myMesh);
+}
+
+void Reen::MarchingCubesRBF::perform(const std::vector<Base::Vector3f>& normals)
+{
+    if (myPoints.size() != normals.size())
+        throw Base::RuntimeError("Number of points doesn't match with number of normals");
+
+    PointCloud<PointNormal>::Ptr cloud_with_normals (new PointCloud<PointNormal>);
+    search::KdTree<PointNormal>::Ptr tree;
+
+    cloud_with_normals->reserve(myPoints.size());
+    std::size_t num_points = myPoints.size();
+    const std::vector<Base::Vector3f>& points = myPoints.getBasicPoints();
+    for (std::size_t index=0; index<num_points; index++) {
+        const Base::Vector3f& p = points[index];
+        const Base::Vector3f& n = normals[index];
+        PointNormal pn;
+        pn.x = p.x;
+        pn.y = p.y;
+        pn.z = p.z;
+        pn.normal_x = n.x;
+        pn.normal_y = n.y;
+        pn.normal_z = n.z;
+        cloud_with_normals->push_back(pn);
+    }
+
+    // Create search tree
+    tree.reset (new search::KdTree<PointNormal>);
+    tree->setInputCloud (cloud_with_normals);
+
+
+    // Init objects
+    pcl::MarchingCubesRBF<PointNormal> rbf;
+
+    // Set parameters
+    rbf.setIsoLevel (0);
+    rbf.setGridResolution (60, 60, 60);
+    rbf.setPercentageExtendGrid (0.1f);
+    rbf.setOffSurfaceDisplacement (0.02f);
+
+    rbf.setInputCloud (cloud_with_normals);
+    rbf.setSearchMethod (tree);
+
+    // Reconstruct
+    PolygonMesh mesh;
+    rbf.reconstruct (mesh);
+
+    MeshConversion::convert(mesh, myMesh);
+}
+
+// ----------------------------------------------------------------------------
+
+Reen::MarchingCubesHoppe::MarchingCubesHoppe(const Points::PointKernel& pts, Mesh::MeshObject& mesh)
+  : myPoints(pts)
+  , myMesh(mesh)
+{
+}
+
+void Reen::MarchingCubesHoppe::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;
+
+    cloud->reserve(myPoints.size());
+    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
+    pcl::MarchingCubesHoppe<PointNormal> hoppe;
+
+    // Set parameters
+    hoppe.setIsoLevel (0);
+    hoppe.setGridResolution (60, 60, 60);
+    hoppe.setPercentageExtendGrid (0.1f);
+
+    hoppe.setInputCloud (cloud_with_normals);
+    hoppe.setSearchMethod (tree2);
+
+    // Reconstruct
+    PolygonMesh mesh;
+    hoppe.reconstruct (mesh);
+
+    MeshConversion::convert(mesh, myMesh);
+}
+
+void Reen::MarchingCubesHoppe::perform(const std::vector<Base::Vector3f>& normals)
+{
+    if (myPoints.size() != normals.size())
+        throw Base::RuntimeError("Number of points doesn't match with number of normals");
+
+    PointCloud<PointNormal>::Ptr cloud_with_normals (new PointCloud<PointNormal>);
+    search::KdTree<PointNormal>::Ptr tree;
+
+    cloud_with_normals->reserve(myPoints.size());
+    std::size_t num_points = myPoints.size();
+    const std::vector<Base::Vector3f>& points = myPoints.getBasicPoints();
+    for (std::size_t index=0; index<num_points; index++) {
+        const Base::Vector3f& p = points[index];
+        const Base::Vector3f& n = normals[index];
+        PointNormal pn;
+        pn.x = p.x;
+        pn.y = p.y;
+        pn.z = p.z;
+        pn.normal_x = n.x;
+        pn.normal_y = n.y;
+        pn.normal_z = n.z;
+        cloud_with_normals->push_back(pn);
+    }
+
+    // Create search tree
+    tree.reset (new search::KdTree<PointNormal>);
+    tree->setInputCloud (cloud_with_normals);
+
+
+    // Init objects
+    pcl::MarchingCubesHoppe<PointNormal> hoppe;
+
+    // Set parameters
+    hoppe.setIsoLevel (0);
+    hoppe.setGridResolution (60, 60, 60);
+    hoppe.setPercentageExtendGrid (0.1f);
+
+    hoppe.setInputCloud (cloud_with_normals);
+    hoppe.setSearchMethod (tree);
+
+    // Reconstruct
+    PolygonMesh mesh;
+    hoppe.reconstruct (mesh);
+
+    MeshConversion::convert(mesh, myMesh);
+}
+
 // ----------------------------------------------------------------------------
 
 void MeshConversion::convert(const pcl::PolygonMesh& pclMesh, Mesh::MeshObject& meshObject)

src/Mod/ReverseEngineering/App/SurfaceTriangulation.h

@@ -24,6 +24,9 @@
 #ifndef REEN_SURFACETRIANGULATION_H
 #define REEN_SURFACETRIANGULATION_H
 
+#include <Base/Vector3D.h>
+#include <vector>
+
 namespace Points {class PointKernel;}
 namespace Mesh {class MeshObject;}
 namespace pcl {struct PolygonMesh;}
@@ -40,18 +43,48 @@ class SurfaceTriangulation
 {
 public:
     SurfaceTriangulation(const Points::PointKernel&, Mesh::MeshObject&);
-    void perform(double searchRadius, double mu);
+    /** \brief Set the number of k nearest neighbors to use for the normal estimation.
-    
+      * \param[in] k the number of k-nearest neighbors
+      */
+    void perform(int ksearch);
+    /** \brief Pass the normals to the points given in the constructor.
+      * \param[in] normals the normals to the given points.
+      */
+    void perform(const std::vector<Base::Vector3f>& normals);
+
+    /** \brief Set the multiplier of the nearest neighbor distance to obtain the final search radius for each point
+      * (this will make the algorithm adapt to different point densities in the cloud).
+      * \param[in] mu the multiplier
+      */
+    inline void 
+    setMu (double mu) { this->mu = mu; }
+
+    /** \brief Set the sphere radius that is to be used for determining the k-nearest neighbors used for triangulating.
+      * \param[in] radius the sphere radius that is to contain all k-nearest neighbors
+      * \note This distance limits the maximum edge length!
+      */
+    inline void 
+    setSearchRadius (double radius) { this->searchRadius = radius; }
+
 private:
     const Points::PointKernel& myPoints;
     Mesh::MeshObject& myMesh;
+    double mu;
+    double searchRadius;
 };
 
 class PoissonReconstruction
 {
 public:
     PoissonReconstruction(const Points::PointKernel&, Mesh::MeshObject&);
+    /** \brief Set the number of k nearest neighbors to use for the normal estimation.
+      * \param[in] k the number of k-nearest neighbors
+      */
     void perform(int ksearch=5);
+    /** \brief Pass the normals to the points given in the constructor.
+      * \param[in] normals the normals to the given points.
+      */
+    void perform(const std::vector<Base::Vector3f>& normals);
 
     /** \brief Set the maximum depth of the tree that will be used for surface reconstruction.
       * \note Running at depth d corresponds to solving on a voxel grid whose resolution is no larger than
@@ -88,6 +121,72 @@ private:
     float samplesPerNode;
 };
 
+class GridReconstruction
+{
+public:
+    GridReconstruction(const Points::PointKernel&, Mesh::MeshObject&);
+    /** \brief Set the number of k nearest neighbors to use for the normal estimation.
+      * \param[in] k the number of k-nearest neighbors
+      */
+    void perform(int ksearch=5);
+    /** \brief Pass the normals to the points given in the constructor.
+      * \param[in] normals the normals to the given points.
+      */
+    void perform(const std::vector<Base::Vector3f>& normals);
+
+private:
+    const Points::PointKernel& myPoints;
+    Mesh::MeshObject& myMesh;
+};
+
+class ImageTriangulation
+{
+public:
+    ImageTriangulation(int width, int height, const Points::PointKernel&, Mesh::MeshObject&);
+    void perform();
+
+private:
+    int width, height;
+    const Points::PointKernel& myPoints;
+    Mesh::MeshObject& myMesh;
+};
+
+class MarchingCubesRBF
+{
+public:
+    MarchingCubesRBF(const Points::PointKernel&, Mesh::MeshObject&);
+    /** \brief Set the number of k nearest neighbors to use for the normal estimation.
+      * \param[in] k the number of k-nearest neighbors
+      */
+    void perform(int ksearch=5);
+    /** \brief Pass the normals to the points given in the constructor.
+      * \param[in] normals the normals to the given points.
+      */
+    void perform(const std::vector<Base::Vector3f>& normals);
+
+private:
+    const Points::PointKernel& myPoints;
+    Mesh::MeshObject& myMesh;
+};
+
+class MarchingCubesHoppe
+{
+public:
+    MarchingCubesHoppe(const Points::PointKernel&, Mesh::MeshObject&);
+    /** \brief Set the number of k nearest neighbors to use for the normal estimation.
+      * \param[in] k the number of k-nearest neighbors
+      */
+    void perform(int ksearch=5);
+    /** \brief Pass the normals to the points given in the constructor.
+      * \param[in] normals the normals to the given points.
+      */
+    void perform(const std::vector<Base::Vector3f>& normals);
+
+private:
+    const Points::PointKernel& myPoints;
+    Mesh::MeshObject& myMesh;
+};
+
 } // namespace Reen
 
 #endif // REEN_SURFACETRIANGULATION_H