An error occurred while loading the file. Please try again.
areafromconcavehull.h 13.54 KiB
#ifndef AREAFROMCONCAVEHULL_H
//PCL
#include <pcl/io/pcd_io.h>
#include <pcl/point_types.h>
#include <pcl/sample_consensus/method_types.h>
#include <pcl/sample_consensus/model_types.h>
#include <pcl/sample_consensus/sac_model_line.h>
#include <pcl/segmentation/sac_segmentation.h>
#include <pcl/filters/voxel_grid.h>
#include <pcl/kdtree/kdtree_flann.h>
#include <pcl/common/pca.h>
#include <pcl/filters/extract_indices.h>
#include <pcl/filters/statistical_outlier_removal.h>
#include <pcl/surface/mls.h>
#include <pcl/surface/poisson.h>
#include <pcl/surface/gp3.h>
#include <pcl/features/normal_3d.h>
#include <pcl/visualization/pcl_visualizer.h>
#include <pcl/io/vtk_io.h>
// STL
#include <tgmath.h>
// own
#include "removeconcavehullframe.h"
#include "../h/inrange.h"
#define AREAFROMCONCAVEHULL_H
template <typename PointInT>
class AreaFromConcaveHull: public  pcl::PCLBase<PointInT>
protected:
    typedef typename pcl::PointCloud<PointInT>::Ptr PointCloudPtr;
    typedef typename pcl::PointCloud<PointInT>::ConstPtr PointCloudConstPtr;
public:
    typedef boost::shared_ptr< AreaFromConcaveHull<PointInT> > Ptr;
    typedef boost::shared_ptr< const AreaFromConcaveHull<PointInT> > ConstPtr;
    AreaFromConcaveHull () :
          lineDistThreshold_ (0.1),
          pointIncrement_ (0.1),
          area_ (0.0),
          kdRadius_ (1),
          coefficients_ (new pcl::ModelCoefficients),
          gridLine_ (new pcl::PointCloud<PointInT>),
          extractedLine_ (new pcl::PointCloud<PointInT>)
          { };
  //  inline void setPointIncrement(float inc) {pointIncrement_ = inc;}
    inline void setKdTreeRadius(float radius) {kdRadius_ = radius;}
    inline void setLineDistThreshold(float threshold) {lineDistThreshold_ = threshold;}
    void greedyArea(float &area);
    void detectLines(typename pcl::PointCloud<PointInT>::Ptr cloud_out);
    void calculateAreaFromHull(typename pcl::PointCloud<PointInT>::Ptr cloudOut);
void makeGridLine(); virtual ~AreaFromConcaveHull () { } protected: float kdRadius_; float area_; /** \brief The grid spacing in x,yz, direction. Smaller values lead to to larger and denser point clouds (default: 0.25). */ const float pointIncrement_; float lineDistThreshold_; /** \brief An input point cloud describing the surface that is to be used for nearest neighbors estimation. */ using pcl::PCLBase<PointInT>::input_; /** \brief The subtracted (inverted) point cloud. */ typename pcl::PointCloud<PointInT>::Ptr extractedLine_; typename pcl::PointCloud<PointInT>::Ptr gridLine_; pcl::ModelCoefficients::Ptr coefficients_; }; template <typename PointInT> void AreaFromConcaveHull<PointInT>::greedyArea(float &area){ // read in parameter file std::vector<float> params; readParamsFromFile("../cfg/GreedyArea.cfg", params, true); typename pcl::PointCloud<pcl::PointNormal>::Ptr cloudWithNormals (new pcl::PointCloud<pcl::PointNormal>); typename pcl::PointCloud<pcl::PointNormal>::Ptr cloudWithNormalsMLS (new pcl::PointCloud<pcl::PointNormal>); typename pcl::PointCloud<pcl::PointXYZ>::Ptr inputXYZ (new pcl::PointCloud<pcl::PointXYZ>); typename pcl::PointCloud<pcl::PointXYZ>::Ptr filteredXYZ (new pcl::PointCloud<pcl::PointXYZ>); // manually copy all points over to pointXYZ data type pcl::PointXYZ p; for (int i = 0; i < input_->points.size(); i++){ p.x= input_->points[i].x; p.y= input_->points[i].y; p.z= input_->points[i].z; inputXYZ->points.push_back(p); } inputXYZ->width = inputXYZ->points.size(); inputXYZ->height = 1; // Create the filtering object pcl::StatisticalOutlierRemoval<pcl::PointXYZ> sor; sor.setInputCloud (inputXYZ); sor.setMeanK (50); sor.setStddevMulThresh (0.75); sor.filter (*inputXYZ); pcl::PointIndices::Ptr inliers (new pcl::PointIndices); // Create the segmentation object pcl::SACSegmentation<pcl::PointXYZ> seg; // Optional seg.setOptimizeCoefficients (true); // Mandatory seg.setModelType (pcl::SACMODEL_PLANE); seg.setMethodType (pcl::SAC_RANSAC); seg.setDistanceThreshold (params[0]); seg.setInputCloud (inputXYZ); seg.segment (*inliers, *coefficients_);
pcl::ExtractIndices<pcl::PointXYZ> extract; extract.setInputCloud (inputXYZ); extract.setIndices (inliers); extract.setNegative (false); // extract.filter (*filteredXYZ); extract.filter (*inputXYZ); // estimate normals typename pcl::NormalEstimation<pcl::PointXYZ, pcl::PointNormal> ne; ne.setInputCloud (inputXYZ); typename pcl::search::KdTree<pcl::PointXYZ>::Ptr tree (new pcl::search::KdTree<pcl::PointXYZ> ()); ne.setSearchMethod (tree); // Use all neighbors in a sphere of radius 3cm ne.setRadiusSearch (params[1]); ne.compute (*cloudWithNormals); // concaternate the normals and the input xyz data pcl::concatenateFields(*inputXYZ, *cloudWithNormals, *cloudWithNormals); /* typename pcl::search::KdTree<pcl::PointXYZ>::Ptr treemls (new pcl::search::KdTree<pcl::PointXYZ> ()); // Init object (second point type is for the normals, even if unused) //viewer<pcl::PointXYZ> (inputXYZ); pcl::MovingLeastSquares<pcl::PointXYZ, pcl::PointNormal> mls; mls.setInputCloud (inputXYZ); mls.setComputeNormals(true); mls.setSearchRadius (0.1); mls.setSearchMethod(treemls); mls.setPolynomialFit (false); //mls.setPolynomialOrder (2); //mls.setUpsamplingMethod (pcl::MovingLeastSquares<pcl::PointXYZ, pcl::PointXYZ>::RANDOM_UNIFORM_DENSITY); //mls.setUpsamplingMethod (pcl::MovingLeastSquares<pcl::PointXYZ, pcl::PointXYZ>::VOXEL_GRID_DILATION); //mls.setDilationVoxelSize (0.1); //mls.setDilationIterations(5); //mls.setUpsamplingRadius (0.5); //mls.setUpsamplingStepSize (0.3); // Reconstruct std::cout << "starting MovingLeastSqr Normal Estimation and Smoothing..."; mls.process (*cloudWithNormalsMLS); std::cout << "Process completed." << std::endl; // viewer<pcl::PointXYZ> (*cloudWithNormalsMLS); */ typename pcl::search::KdTree<pcl::PointNormal>::Ptr tree1 (new pcl::search::KdTree<pcl::PointNormal> ()); pcl::GreedyProjectionTriangulation<pcl::PointNormal> gp; pcl::PolygonMesh mesh; // GreedyTriangulation gp.setInputCloud(cloudWithNormals); gp.setSearchMethod(tree1); /*gp.setMu(2.20); gp.setSearchRadius(10.40); gp.setMaximumNearestNeighbors(500); gp.setMaximumSurfaceAngle(M_PI/4); gp.setMaximumAngle(3*M_PI/4); gp.setMinimumAngle(M_PI/18); */ gp.setMu(params[2]);
gp.setSearchRadius(params[3]); gp.setMaximumNearestNeighbors(params[4]); gp.setMaximumSurfaceAngle(params[5]); gp.setMaximumAngle(params[6]); gp.setMinimumAngle(params[7]); gp.setNormalConsistency(bool(params[8])); std::cout << "\nstarting reconstructing surface mesh..." << std::endl; std::cout << " \n"; gp.reconstruct(mesh); std::cout << "Process completed." << std::endl; /**/ /* pcl::Poisson<pcl::PointNormal> poisson; poisson.setDepth(9); poisson.setInputCloud(cloudWithNormals); poisson.reconstruct(mesh); */ cout<<"The size of polygons is "<<mesh.polygons.size()<<endl; pcl::io::saveVTKFile("potting-reconstruction.vtk",mesh); //pcl::io::savePLYFile("potting-reconstruction.ply",mesh); //calculate area; pcl::PointCloud<pcl::PointXYZ> cloud1; cloud1=*inputXYZ; int index_p1,index_p2,index_p3; float x1,x2,x3,y1,y2,y3,z1,z2,z3,a,b,c,q; for(int i=0;i<mesh.polygons.size();++i) { index_p1=mesh.polygons[i].vertices[0]; index_p2=mesh.polygons[i].vertices[1]; index_p3=mesh.polygons[i].vertices[2]; x1=cloud1.points[index_p1].x; y1=cloud1.points[index_p1].y; z1=cloud1.points[index_p1].z; x2=cloud1.points[index_p2].x; y2=cloud1.points[index_p2].y; z2=cloud1.points[index_p2].z; x3=cloud1.points[index_p3].x; y3=cloud1.points[index_p3].y; z3=cloud1.points[index_p3].z; //Heron's formula; a=sqrt(pow((x1-x2),2)+pow((y1-y2),2)+pow((z1-z2),2)); b=sqrt(pow((x1-x3),2)+pow((y1-y3),2)+pow((z1-z3),2)); c=sqrt(pow((x3-x2),2)+pow((y3-y2),2)+pow((z3-z2),2)); q=(a+b+c)/2; float radizierer = q*(q-a)*(q-b)*(q-c); if (radizierer < 0) std::cout << "negativ of sqrt not allowed" << std::endl; //area=area+sqrt(fabs(radizierer)); else area=area+sqrt(radizierer); } std::cout << area << std::endl; area_ = area; boost::shared_ptr<pcl::visualization::PCLVisualizer> viewer (new pcl::visualization::PCLVisualizer ("3D Viewer Greedy Mesh Reconstruction")); viewer->setBackgroundColor (0, 0, 0); viewer->addPolygonMesh(mesh,"meshes",0); viewer->addCoordinateSystem (1.0);
viewer->initCameraParameters (); while (!viewer->wasStopped ()){ viewer->spinOnce (100); boost::this_thread::sleep (boost::posix_time::microseconds (100000)); } } template <typename PointInT> void AreaFromConcaveHull<PointInT>::makeGridLine(){ PointInT startPt, endPt, minPt, maxPt; pcl::getMinMax3D(*extractedLine_, minPt, maxPt); size_t nPts = extractedLine_->points.size(); startPt = extractedLine_->points[0]; endPt = extractedLine_->points[nPts-1]; ////////// detect plane //////////////////////////////// pcl::PointIndices::Ptr inliers (new pcl::PointIndices); // Create the segmentation object pcl::SACSegmentation<PointInT> seg; // Optional seg.setOptimizeCoefficients (true); // Mandatory seg.setModelType (pcl::SACMODEL_PLANE); seg.setMethodType (pcl::SAC_RANSAC); seg.setDistanceThreshold (lineDistThreshold_); seg.setInputCloud (input_); seg.segment (*inliers, *coefficients_); std::cout << "Normalized coeffs of Pt-Normal Eq. (Ax+By+Cz=D): " << coefficients_->values[0] << " " << coefficients_->values[1] << " " << coefficients_->values[2] << " " << coefficients_->values[3] << std::endl; float z = (startPt.z + endPt.z) / 2.0; PointInT point; // construct new points on the plane by rearranging the point normal equation for one of the coordinates. for (float x=minPt.x; x<=maxPt.x; x+=pointIncrement_) { for (float y=minPt.y; y<= maxPt.y; y+=pointIncrement_) { point.x = x; if (coefficients_->values[1] != 0) // to avoid division by zero point.y = -1.0 * (1.0 * coefficients_->values[3] + coefficients_->values[0] * x + coefficients_->values[2] * z ) / coefficients_->values[1] ; else // divisor is set to 0.0000001 if the coefficient is zero point.y = -1.0 * (1.0 * coefficients_->values[3] + coefficients_->values[0] * x + coefficients_->values[2] * z ) / 0.000001 ; point.z = z; gridLine_->points.push_back (point); } } std::cerr << "GridLine_ has: " << gridLine_->points.size () << " points. " << std::endl; viewer<PointInT> (gridLine_); } template <typename PointInT> void AreaFromConcaveHull<PointInT>::detectLines(typename pcl::PointCloud<PointInT>::Ptr cloudOut){ typename pcl::PointCloud<PointInT>::Ptr cloudFiltered (new pcl::PointCloud<PointInT>);
pcl::copyPointCloud(*input_, *cloudFiltered); pcl::ExtractIndices<PointInT> extract; pcl::PointIndices::Ptr inliers (new pcl::PointIndices); // Create the segmentation object pcl::SACSegmentation<PointInT> seg; // Optional seg.setOptimizeCoefficients (true); // Mandatory //seg.setModelType (pcl::SACMODEL_LINE); seg.setModelType(pcl::SACMODEL_LINE); seg.setMethodType (pcl::SAC_RANSAC); seg.setDistanceThreshold (lineDistThreshold_); seg.setInputCloud (cloudFiltered); seg.segment (*inliers, *coefficients_); std::cerr << "PointCloud after LINE segmentation has: " << inliers->indices.size () << " inliers." << std::endl; // Extract the inliers extract.setInputCloud (cloudFiltered); extract.setIndices (inliers); extract.setNegative (false); extract.filter (*extractedLine_); viewer<PointInT> (extractedLine_); Eigen::VectorXf coeffs_; std::vector<int> indx; std::vector<double> distToModel; indx.push_back(*inliers->indices.begin()); indx.push_back(*inliers->indices.end()); pcl::SampleConsensusModelLine<PointInT> segLine(input_, false); segLine.setInputCloud(input_); segLine.computeModelCoefficients(indx,coeffs_); segLine.getDistancesToModel(coeffs_, distToModel); std::cout << coeffs_ << std::endl; makeGridLine(); pcl::copyPointCloud(*gridLine_, *cloudOut); } #endif // AREAFROMCONCAVEHULL_H