From 93b53f16629f81f1fa53fdca6cf8a1755f991054 Mon Sep 17 00:00:00 2001 From: wmayer Date: Sat, 5 Dec 2015 16:19:18 +0100 Subject: [PATCH] + add several surface reconstruction methods from pcl to Reen module --- .../App/AppReverseEngineering.cpp | 274 ++++++++- .../App/SurfaceTriangulation.cpp | 522 +++++++++++++++++- .../App/SurfaceTriangulation.h | 103 +++- 3 files changed, 871 insertions(+), 28 deletions(-) diff --git a/src/Mod/ReverseEngineering/App/AppReverseEngineering.cpp b/src/Mod/ReverseEngineering/App/AppReverseEngineering.cpp index ff2f32ce3..d04e9531c 100644 --- a/src/Mod/ReverseEngineering/App/AppReverseEngineering.cpp +++ b/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(pcObj); - Points::PointKernel* points = pPoints->getPointKernelPtr(); - + Points::PointKernel* points = static_cast(pts)->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 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(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}; - if (!PyArg_ParseTupleAndKeywords(args.ptr(), kwds.ptr(), "O!|iiid", kwds_poisson, - &(Points::PointsPy::Type), &pcObj, - &ksearch, &octreeDepth, &solverDivide, &samplesPerNode)) + "SamplesPerNode", "Normals", NULL}; + if (!PyArg_ParseTupleAndKeywords(args.ptr(), kwds.ptr(), "O!|iiidO", kwds_poisson, + &(Points::PointsPy::Type), &pts, + &ksearch, &octreeDepth, &solverDivide, &samplesPerNode, &vec)) throw Py::Exception(); - Points::PointsPy* pPoints = static_cast(pcObj); - Points::PointKernel* points = pPoints->getPointKernelPtr(); + Points::PointKernel* points = static_cast(pts)->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 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(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(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(pts)->getPointKernelPtr(); + + Mesh::MeshObject* mesh = new Mesh::MeshObject(); + GridReconstruction tria(*points, *mesh); + if (vec) { + Py::Sequence list(vec); + std::vector 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(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(pts)->getPointKernelPtr(); + + Mesh::MeshObject* mesh = new Mesh::MeshObject(); + MarchingCubesRBF tria(*points, *mesh); + if (vec) { + Py::Sequence list(vec); + std::vector 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(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(pts)->getPointKernelPtr(); + + Mesh::MeshObject* mesh = new Mesh::MeshObject(); + MarchingCubesHoppe tria(*points, *mesh); + if (vec) { + Py::Sequence list(vec); + std::vector 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(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, °ree, &refinement, &iterations, &interiorSmoothness, &interiorWeight, &boundarySmoothness, &boundaryWeight)) throw Py::Exception(); - Points::PointsPy* pPoints = static_cast(pcObj); - Points::PointKernel* points = pPoints->getPointKernelPtr(); + Points::PointKernel* points = static_cast(pts)->getPointKernelPtr(); BSplineFitting fit(points->getBasicPoints()); fit.setOrder(degree+1); diff --git a/src/Mod/ReverseEngineering/App/SurfaceTriangulation.cpp b/src/Mod/ReverseEngineering/App/SurfaceTriangulation.cpp index 96af58437..4526e3f39 100644 --- a/src/Mod/ReverseEngineering/App/SurfaceTriangulation.cpp +++ b/src/Mod/ReverseEngineering/App/SurfaceTriangulation.cpp @@ -26,8 +26,10 @@ #include "SurfaceTriangulation.h" #include #include +#include #include #include +#include // http://svn.pointclouds.org/pcl/tags/pcl-1.5.1/test/ #if defined(HAVE_PCL_SURFACE) @@ -42,6 +44,8 @@ //#include //#include #include +#include +#include #include #include #include @@ -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::Ptr cloud (new PointCloud); PointCloud::Ptr cloud_with_normals (new PointCloud); search::KdTree::Ptr tree; search::KdTree::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); 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 parts = gp3.getPartIDs(); + //std::vector states = gp3.getPointStates(); +} + +void SurfaceTriangulation::perform(const std::vector& normals) +{ + if (myPoints.size() != normals.size()) + throw Base::RuntimeError("Number of points doesn't match with number of normals"); + + PointCloud::Ptr cloud_with_normals (new PointCloud); + search::KdTree::Ptr tree; + + cloud_with_normals->reserve(myPoints.size()); + std::size_t num_points = myPoints.size(); + const std::vector& points = myPoints.getBasicPoints(); + for (std::size_t index=0; indexpush_back(pn); + } + + // Create search tree + tree.reset (new search::KdTree); + tree->setInputCloud (cloud_with_normals); + + // Init objects + GreedyProjectionTriangulation 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::Ptr tree; search::KdTree::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); tree2->setInputCloud (cloud_with_normals); @@ -180,6 +247,449 @@ void PoissonReconstruction::perform(int ksearch) MeshConversion::convert(mesh, myMesh); } +void PoissonReconstruction::perform(const std::vector& normals) +{ + if (myPoints.size() != normals.size()) + throw Base::RuntimeError("Number of points doesn't match with number of normals"); + + PointCloud::Ptr cloud_with_normals (new PointCloud); + search::KdTree::Ptr tree; + + cloud_with_normals->reserve(myPoints.size()); + std::size_t num_points = myPoints.size(); + const std::vector& points = myPoints.getBasicPoints(); + for (std::size_t index=0; indexpush_back(pn); + } + + // Create search tree + tree.reset (new search::KdTree); + tree->setInputCloud (cloud_with_normals); + + // Init objects + Poisson 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::Ptr cloud (new PointCloud); + PointCloud::Ptr cloud_with_normals (new PointCloud); + search::KdTree::Ptr tree; + search::KdTree::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 (false)); + tree->setInputCloud (cloud); + + // Normal estimation + NormalEstimation n; + PointCloud::Ptr normals (new PointCloud ()); + 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); + tree2->setInputCloud (cloud_with_normals); + + // Init objects + GridProjection 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& normals) +{ + if (myPoints.size() != normals.size()) + throw Base::RuntimeError("Number of points doesn't match with number of normals"); + + PointCloud::Ptr cloud_with_normals (new PointCloud); + search::KdTree::Ptr tree; + + cloud_with_normals->reserve(myPoints.size()); + std::size_t num_points = myPoints.size(); + const std::vector& points = myPoints.getBasicPoints(); + for (std::size_t index=0; indexpush_back(pn); + } + + // Create search tree + tree.reset (new search::KdTree); + tree->setInputCloud (cloud_with_normals); + + // Init objects + GridProjection 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::Ptr cloud_organized (new pcl::PointCloud ()); + cloud_organized->width = width; + cloud_organized->height = height; + cloud_organized->points.resize (cloud_organized->width * cloud_organized->height); + + const std::vector& 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 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::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 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 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::Ptr cloud (new PointCloud); + PointCloud::Ptr cloud_with_normals (new PointCloud); + search::KdTree::Ptr tree; + search::KdTree::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 (false)); + tree->setInputCloud (cloud); + + // Normal estimation + NormalEstimation n; + PointCloud::Ptr normals (new PointCloud ()); + 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); + tree2->setInputCloud (cloud_with_normals); + + // Init objects + pcl::MarchingCubesRBF 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& normals) +{ + if (myPoints.size() != normals.size()) + throw Base::RuntimeError("Number of points doesn't match with number of normals"); + + PointCloud::Ptr cloud_with_normals (new PointCloud); + search::KdTree::Ptr tree; + + cloud_with_normals->reserve(myPoints.size()); + std::size_t num_points = myPoints.size(); + const std::vector& points = myPoints.getBasicPoints(); + for (std::size_t index=0; indexpush_back(pn); + } + + // Create search tree + tree.reset (new search::KdTree); + tree->setInputCloud (cloud_with_normals); + + + // Init objects + pcl::MarchingCubesRBF 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::Ptr cloud (new PointCloud); + PointCloud::Ptr cloud_with_normals (new PointCloud); + search::KdTree::Ptr tree; + search::KdTree::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 (false)); + tree->setInputCloud (cloud); + + // Normal estimation + NormalEstimation n; + PointCloud::Ptr normals (new PointCloud ()); + 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); + tree2->setInputCloud (cloud_with_normals); + + // Init objects + pcl::MarchingCubesHoppe 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& normals) +{ + if (myPoints.size() != normals.size()) + throw Base::RuntimeError("Number of points doesn't match with number of normals"); + + PointCloud::Ptr cloud_with_normals (new PointCloud); + search::KdTree::Ptr tree; + + cloud_with_normals->reserve(myPoints.size()); + std::size_t num_points = myPoints.size(); + const std::vector& points = myPoints.getBasicPoints(); + for (std::size_t index=0; indexpush_back(pn); + } + + // Create search tree + tree.reset (new search::KdTree); + tree->setInputCloud (cloud_with_normals); + + + // Init objects + pcl::MarchingCubesHoppe 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) diff --git a/src/Mod/ReverseEngineering/App/SurfaceTriangulation.h b/src/Mod/ReverseEngineering/App/SurfaceTriangulation.h index 0e0fdf200..e075ae9c8 100644 --- a/src/Mod/ReverseEngineering/App/SurfaceTriangulation.h +++ b/src/Mod/ReverseEngineering/App/SurfaceTriangulation.h @@ -24,6 +24,9 @@ #ifndef REEN_SURFACETRIANGULATION_H #define REEN_SURFACETRIANGULATION_H +#include +#include + 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& 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& 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& 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& 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& normals); + +private: + const Points::PointKernel& myPoints; + Mesh::MeshObject& myMesh; +}; + } // namespace Reen #endif // REEN_SURFACETRIANGULATION_H