pcl--第十节 点云曲面重建

曲面重建技术在逆向工程、数据可视化、机器视觉、虚拟现实、医疗技术等领域中得到了广泛的应用 。 例如，在汽车、航空等工业领域中，复杂外形产品的设计仍需要根据手工模型，采用逆向工程的手段建立产品的数字化模型，根据测量数据建立人体以及骨骼和器官的计算机模型，在医学、定制生产等方面都有重要意义 。

除了上述传统的行业，随着新兴的廉价 RGBD 获取设备在数字娱乐行业的病毒式扩展，使得更多人开始使用点云来处理对象并进行工程应用 。 根据重建曲面和数据点云之间的关系，可将曲面重建分为两大类：插值法和逼近法。前者得到的重建曲面完全通过原始数据点，而后者则是用分片线性曲面或其他形式的曲面来逼近原始数据点,从而使得得到的重建曲面是原始点集的一个逼近曲面。

关联知识：

Search、KdTree、Octree

基于多项式重构的平滑和法线估计¶

本教程说明如何使用移动最小二乘（MLS）曲面重构方法来平滑和重采样噪声数据。

使用统计分析很难消除某些数据不规则性（由较小的距离测量误差引起）。要创建完整的模型，必须考虑光滑的表面以及数据中的遮挡。在无法获取其他扫描的情况下，一种解决方案是使用重采样算法，该算法尝试通过周围数据点之间的高阶多项式插值来重新创建表面的缺失部分。通过执行重采样，可以纠正这些小的错误，并且可以将多个扫描记录在一起执行平滑操作合并成同一个点云。

在上图的左侧，我们在包含两个配准点云的数据集中看到了配准后的效果及表面法线估计。由于对齐错误，所产生的法线有噪声。在右侧，使用移动最小二乘法对表面法线估计进行平滑处理后，在同一数据集中看到了该效果。绘制每个点的曲率，作为重新采样前后特征值关系的度量，我们得到：

 

#include <pcl/point_types.h>
#include <pcl/io/pcd_io.h>
//#include <pcl/kdtree/kdtree_flann.h>
#include <pcl/surface/mls.h>
#include <pcl/search/kdtree.h>
int
main ()
{// Load input file into a PointCloud<T> with an appropriate typepcl::PointCloud<pcl::PointXYZ>::Ptr cloud (new pcl::PointCloud<pcl::PointXYZ> ());// Load bun0.pcd -- should be available with the PCL archive in test pcl::io::loadPCDFile ("table_scene_lms400_downsampled.pcd", *cloud);// Create a KD-Tree//pcl::search::KdTree<pcl::PointXYZ>::Ptr tree (new pcl::search::KdTree<pcl::PointXYZ>);pcl::search::KdTree<pcl::PointXYZ>::Ptr tree(new pcl::search::KdTree<pcl::PointXYZ>()); // kd树对象// Output has the PointNormal type in order to store the normals calculated by MLSpcl::PointCloud<pcl::PointNormal> mls_points;// Init object (second point type is for the normals, even if unused)pcl::MovingLeastSquares<pcl::PointXYZ, pcl::PointNormal> mls;mls.setComputeNormals (true);// Set parametersmls.setInputCloud (cloud);mls.setPolynomialOrder (2);mls.setSearchMethod (tree);mls.setSearchRadius (0.03);// Reconstructmls.process (mls_points);// Save outputpcl::io::savePCDFile ("table_scene_lms400_downsampled-mls.pcd", mls_points);
}

右边是重建后，会发现比左边更平整

在平面模型上提取凸(凹)多边形

学习如何为平面模型上的点集提取其对应的凹多边形的例子,该例首先从点云中提取平面模型,再通过该估计的平面模型系数从滤波后的点云投影一组点集形成点云，最后为投影后的点云计算其对应的二维凸(凹)多边形。

#include <pcl/ModelCoefficients.h>
#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/filters/passthrough.h>
#include <pcl/filters/project_inliers.h>
#include <pcl/segmentation/sac_segmentation.h>
#include <pcl/surface/convex_hull.h>intmain ()
{pcl::PointCloud<pcl::PointXYZ>::Ptr cloud (new pcl::PointCloud<pcl::PointXYZ>), cloud_filtered (new pcl::PointCloud<pcl::PointXYZ>), cloud_projected (new pcl::PointCloud<pcl::PointXYZ>);pcl::PCDReader reader;reader.read ("table_scene_mug_stereo_textured.pcd", *cloud);// Build a filter to remove spurious NaNs and scene backgroundpcl::PassThrough<pcl::PointXYZ> pass;pass.setInputCloud (cloud);pass.setFilterFieldName ("z");pass.setFilterLimits (0, 1.1);pass.filter (*cloud_filtered);std::cerr << "PointCloud after filtering has: " << cloud_filtered->size () << " data points." << std::endl;pcl::ModelCoefficients::Ptr coefficients (new pcl::ModelCoefficients);pcl::PointIndices::Ptr inliers (new pcl::PointIndices);// Create the segmentation objectpcl::SACSegmentation<pcl::PointXYZ> seg;// Optionalseg.setOptimizeCoefficients (true);// Mandatoryseg.setModelType (pcl::SACMODEL_PLANE);seg.setMethodType (pcl::SAC_RANSAC);seg.setDistanceThreshold (0.01);seg.setInputCloud (cloud_filtered);seg.segment (*inliers, *coefficients);// Project the model inlierspcl::ProjectInliers<pcl::PointXYZ> proj;proj.setModelType (pcl::SACMODEL_PLANE);proj.setInputCloud (cloud_filtered);proj.setIndices (inliers);proj.setModelCoefficients (coefficients);proj.filter (*cloud_projected);// Create a Convex Hull representation of the projected inlierspcl::PointCloud<pcl::PointXYZ>::Ptr cloud_hull (new pcl::PointCloud<pcl::PointXYZ>);pcl::ConvexHull<pcl::PointXYZ> chull;chull.setInputCloud (cloud_projected);chull.reconstruct (*cloud_hull);std::cerr << "Convex hull has: " << cloud_hull->size () << " data points." << std::endl;pcl::PCDWriter writer;writer.write ("table_scene_mug_stereo_textured_hull.pcd", *cloud_hull, false);return (0);
}

#include <pcl/ModelCoefficients.h>
#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/filters/passthrough.h>
#include <pcl/filters/project_inliers.h>
#include <pcl/segmentation/sac_segmentation.h>
#include <pcl/surface/concave_hull.h>int
main ()
{pcl::PointCloud<pcl::PointXYZ>::Ptr cloud (new pcl::PointCloud<pcl::PointXYZ>), cloud_filtered (new pcl::PointCloud<pcl::PointXYZ>), cloud_projected (new pcl::PointCloud<pcl::PointXYZ>);pcl::PCDReader reader;reader.read ("table_scene_mug_stereo_textured.pcd", *cloud);// Build a filter to remove spurious NaNs and scene backgroundpcl::PassThrough<pcl::PointXYZ> pass;pass.setInputCloud (cloud);pass.setFilterFieldName ("z");pass.setFilterLimits (0, 1.1);pass.filter (*cloud_filtered);std::cerr << "PointCloud after filtering has: "<< cloud_filtered->size () << " data points." << std::endl;pcl::ModelCoefficients::Ptr coefficients (new pcl::ModelCoefficients);pcl::PointIndices::Ptr inliers (new pcl::PointIndices);// Create the segmentation objectpcl::SACSegmentation<pcl::PointXYZ> seg;// Optionalseg.setOptimizeCoefficients (true);// Mandatoryseg.setModelType (pcl::SACMODEL_PLANE);seg.setMethodType (pcl::SAC_RANSAC);seg.setDistanceThreshold (0.01);seg.setInputCloud (cloud_filtered);seg.segment (*inliers, *coefficients);std::cerr << "PointCloud after segmentation has: "<< inliers->indices.size () << " inliers." << std::endl;// Project the model inlierspcl::ProjectInliers<pcl::PointXYZ> proj;proj.setModelType (pcl::SACMODEL_PLANE);// proj.setIndices (inliers);proj.setInputCloud (cloud_filtered);proj.setModelCoefficients (coefficients);proj.filter (*cloud_projected);std::cerr << "PointCloud after projection has: "<< cloud_projected->size () << " data points." << std::endl;// Create a Concave Hull representation of the projected inlierspcl::PointCloud<pcl::PointXYZ>::Ptr cloud_hull (new pcl::PointCloud<pcl::PointXYZ>);pcl::ConcaveHull<pcl::PointXYZ> chull;chull.setInputCloud (cloud_projected);chull.setAlpha (0.1);chull.reconstruct (*cloud_hull);std::cerr << "Concave hull has: " << cloud_hull->size ()<< " data points." << std::endl;pcl::PCDWriter writer;writer.write ("table_scene_mug_stereo_textured_hull.pcd", *cloud_hull, false);return (0);
}

 无序点云的快速三角化

本小节介绍了怎样使用贪婪投影三角化算法对有向点云进行三角化,具体方法是先将有向点云投影到某一局部二维坐标平面内，再在坐标平面内进行平面内的三角化,再根据平面内三位点的拓扑连接关系获得一个三角网格曲面模型。贪婪投影三角化算法原理是处理一系列可以使网格“生长扩大”的点(边缘点)延伸这些点直到所有符合几何正确性和拓扑正确性的点都被连上。该算法的优点是可以处理来自一个或者多个扫描仪扫描得到并且有多个连接处的散乱点云。但该算法也有一定的局限性,它更适用于采样点云来自于表面连续光滑的曲面且点云密度
变化比较均匀的情况。该算法的三角化过程是局部进行的,首先沿着一点的法线将该点投影到局部二维坐标平面内并连接其他悬空点,然后再进行下一点。所以这里设置如下参数:

1. 函数 SetMaximumNearestNeighbors(unsigned)和 SetMu(double),这两个函数的作用是控制搜索邻域大小。前者定义了可搜索的邻域个数,后者规定了被样本点搜索其邻近点的最远距离(是为了适应点云密度的变化),特征值一般是 50~100和2.5~3(或者15每栅格)。
2. 函数SetSearchRadius(double),该函数设置了三角化后得到的每个三角形的最大可能边长。
3. 函数SetMinimumAngle(double)和 SetMaximumAngle(double),这两个函数是三角化后每个三角形的最大角和最小角。两者至少要符合一个,典型值分别是10和120°(弧度)。
4. 函数 SetMaximumSurfaceAgle(double)和 SetNormalConsistency(bool),这两个函数是为了处理边缘或者角很尖锐以及一个表面的两边非常靠近的情况。为了处理这些特殊情况，函数SetMaximumSurfaceAgle(double)规定如果某点法线方向的偏离超过指定角度(注:大多数表面法线估计方法可以估计出连续变化的表面法线方向，即使在尖锐的边缘条件下),该点就不连接到样本点上。该角度是通过计算法向线段(忽略法线方向)之间的角度。函数SetNormalConsistency(bool)保证法线朝向，如果法线方向一致性标识没有设定,就不能保证估计出的法线都可以始终朝向一致。第一个函数特征值为45(弧度)第二个函数默认值为false。

#include <pcl/point_types.h>
#include <pcl/io/pcd_io.h>
#include <pcl/search/kdtree.h> // for KdTree
#include <pcl/features/normal_3d.h>
#include <pcl/surface/gp3.h>
#include <pcl/visualization/pcl_visualizer.h>
#include <pcl/surface/organized_fast_mesh.h>
int
main ()
{// Load input file into a PointCloud<T> with an appropriate typepcl::PointCloud<pcl::PointXYZ>::Ptr cloud (new pcl::PointCloud<pcl::PointXYZ>);pcl::PCLPointCloud2 cloud_blob;pcl::io::loadPCDFile ("bunny.pcd", cloud_blob);pcl::fromPCLPointCloud2 (cloud_blob, *cloud);//* the data should be available in cloud// Normal estimation*pcl::NormalEstimation<pcl::PointXYZ, pcl::Normal> n;pcl::PointCloud<pcl::Normal>::Ptr normals (new pcl::PointCloud<pcl::Normal>);pcl::search::KdTree<pcl::PointXYZ>::Ptr tree (new pcl::search::KdTree<pcl::PointXYZ>);tree->setInputCloud (cloud);n.setInputCloud (cloud);n.setSearchMethod (tree);n.setKSearch (20);n.compute (*normals);//* normals should not contain the point normals + surface curvatures// Concatenate the XYZ and normal fields*pcl::PointCloud<pcl::PointNormal>::Ptr cloud_with_normals (new pcl::PointCloud<pcl::PointNormal>);pcl::concatenateFields (*cloud, *normals, *cloud_with_normals);//* cloud_with_normals = cloud + normals// Create search tree*pcl::search::KdTree<pcl::PointNormal>::Ptr tree2 (new pcl::search::KdTree<pcl::PointNormal>);tree2->setInputCloud (cloud_with_normals);// Initialize objectspcl::GreedyProjectionTriangulation<pcl::PointNormal> gp3;pcl::PolygonMesh triangles;// Set the maximum distance between connected points (maximum edge length)gp3.setSearchRadius (0.02);// Set typical values for the parametersgp3.setMu (2.5);gp3.setMaximumNearestNeighbors (100);gp3.setMaximumSurfaceAngle(M_PI/4); // 45 degreesgp3.setMinimumAngle(M_PI/18); // 10 degreesgp3.setMaximumAngle(2*M_PI/3); // 120 degreesgp3.setNormalConsistency(false);// Get resultgp3.setInputCloud (cloud_with_normals);gp3.setSearchMethod (tree2);gp3.reconstruct (triangles);// Additional vertex informationstd::vector<int> parts = gp3.getPartIDs();std::vector<int> states = gp3.getPointStates();// Finishpcl::visualization::PCLVisualizer viewer("Triangulation");viewer.setBackgroundColor(0.0, 0.0, 0.0);viewer.addPointCloud<pcl::PointXYZ>(cloud, "input_cloud");viewer.addPolygonMesh(triangles, "triangles", 0);viewer.setPointCloudRenderingProperties(pcl::visualization::PCL_VISUALIZER_COLOR, 1.0, 0.0, 0.0, "input_cloud");viewer.setPointCloudRenderingProperties(pcl::visualization::PCL_VISUALIZER_COLOR, 0.0, 1.0, 0.0, "triangles");viewer.spin();return (0);
}