点云的三角化重建

对点云数据进行三角化重建大致步骤：

程序先读取点云文件；然后计算法向量，并将法向量和点云坐标放在一起；接着使用贪婪三角化投影算法进行重构，最后显示结果。

由于获取的点云数据中常常伴有杂点或噪声，影响了后续的处理，因此为了获取完整的模型，需要对点云数据进行一定的预处理，常用的方法有滤波去噪、数据精简、数据插补等。

定义并读取点云数据：

pcl::PointCloudpcl::PointXYZ::Ptr cloud_old(new pcl::PointCloudpcl::PointXYZ);

pcl::PointCloudpcl::PointXYZ::Ptr cloud_downsampled(new pcl::PointCloudpcl::PointXYZ);

pcl::PointCloudpcl::PointXYZ::Ptr cloud(new pcl::PointCloudpcl::PointXYZ);

pcl::io::loadPCDFile(“maize.pcd”, *cloud_old);//读取点云数据

voxel滤波简化：

使用体素化网格方法实现下采样，即减少点的数量，减少点云数据，并同时保持点云的形状特征，作为预处理，可以很好的提高程序的速度，在提高配准、曲面重建、形状识别等算法速度中非常实用。

PCL实现的VoxelGrid类通过输入的点云数据创建一个三维体素栅格（可把体素栅格想象为微小的空间三维立方体的集合），然后在每个体素（即，三维立方体）内，用体素中所有点的重心来近似显示体素中其他点，这样该体素就内所有点就用一个重心点最终表示，对于所有体素处理后得到过滤后的点云。这种方法比用体素中心来逼近的方法更慢，但它对于采样点对应曲面的表示更为准确。

pcl::VoxelGridpcl::PointXYZ voxelSampler;

// pcl::VoxelGrid filter 的 leaf size 为 1cm，结果储存在cloud_downsampled中

voxelSampler.setInputCloud(cloud_old);

voxelSampler.setLeafSize(0.01f, 0.01f, 0.01f);

voxelSampler.filter(*cloud_downsampled);

使用StatisticalOutlierRemoval滤波器移除离群点：

对每个点的邻域进行一个统计分析，并修剪掉那些不符合一定标准的点。稀疏离群点移除方法基于在输入数据中对点到临近点的距离分布的计算。对每个点，计算它到它的所有临近点的平均距离。假设得到的结果是一个高斯分布，其形状由均值和标准差决定，平均距离在标准范围（由全局距离平均值和方差定义）之外的点，可被定义为离群点并可从数据集中去除掉。

创建了一个pcl::StatisticalOutlierRemoval滤波器，将对每个点分析的临近点个数设为10，并将标准差倍数设为0.2，这意味着如果一个点的距离超出平均距离0.2个标准差以上，则该点被标记为离群点，并将被移除。计算后的输出结果储存在cloud中。

pcl::StatisticalOutlierRemovalpcl::PointXYZ statFilter; // 创建滤波器对象

statFilter.setInputCloud(cloud_downsampled); //设置待滤波的点云

statFilter.setMeanK(10); //设置在进行统计时考虑查询点邻近点数

statFilter.setStddevMulThresh(0.2); //设置判断是否为离群点的阈值

statFilter.filter(*cloud); //执行滤波处理保存内点到cloud_filtered

法向量估计：

使用的三角化算法输入必须为有向点云，所以需要使用PCL中的法线估计方法预先估计出数据中每个点的法线。

表面法线是几何体表面的重要属性，在很多领域都有大量应用。对于一个已知的几何体表面，根据垂直于点表面的矢量，因此推断表面某一点的法线方向通常比较简单。然而，由于我们获取的点云数据集在真实物体的表面表现为一组定点样本，这样就会有两种解决方法：

（1）使用曲面重建技术，从获取的点云数据集中得到采样点对应的曲面，然后从曲面模型中计算表面法线；

（2）直接从点云数据集中近似推断表面法线。

理论基础：确定表面一点法线的问题近似于估计表面的一个相切面法线的问题，因此转变成一个最小二乘法平面拟合估计问题。

pcl::NormalEstimation<pcl::PointXYZ, pcl::Normal> ne; //创建法向估计对象

ne.setInputCloud(cloud); //为法向估计对象输入点云

pcl::search::KdTreepcl::PointXYZ::Ptr tree(new pcl::search::KdTreepcl::PointXYZ());

ne.setSearchMethod(tree); //设置法向估计时采用的搜索方式为kdtree

pcl::PointCloudpcl::Normal::Ptr normals(new pcl::PointCloudpcl::Normal);

ne.setRadiusSearch(0.01); // Use all neighbors in a sphere of radius 1cm

ne.compute(*normals); //计算法向量

std::cerr << “法线计算完成” << std::endl;

三角化：

由于XYZ坐标字段和法线字段需要在相同PointCloud对象中，所以创建一个新的PointNormal类型的点云来存储坐标字段和法线连接后的点云。

以下代码是对三角化对象相关变量进行定义

// Concatenate the XYZ and normal fields*

pcl::PointCloudpcl::PointNormal::Ptr cloud_with_normals(new pcl::PointCloudpcl::PointNormal);

pcl::concatenateFields(*cloud, *normals, *cloud_with_normals);

//* cloud_with_normals = cloud + normals //将点云和法线放到一起

// Create search tree* //创建搜索树

pcl::search::KdTreepcl::PointNormal::Ptr tree2(new pcl::search::KdTreepcl::PointNormal);

tree2->setInputCloud(cloud_with_normals); //利用点云构建搜索树

// Initialize objects //初始化GreedyProjectionTriangulation对象，并设置参数

pcl::GreedyProjectionTriangulationpcl::PointNormal gp3;

pcl::PolygonMesh triangles; //创建多变形网格，用于存储结果

以下代码是设置参数特征值和实际三角化的过程。

// Set the maximum distance between connected points (maximum edge length)

//设置GreedyProjectionTriangulation对象的参数 //第一个参数影响很大

gp3.setSearchRadius(0.025); //设置连接点之间的最大距离（三角形最大边长）用于确定k近邻的球半径【默认值 0】

gp3.setMu(2.5); //设置最近邻距离的乘子，以得到每个点的最终搜索半径【默认值 0】

gp3.setMaximumNearestNeighbors(100); //设置搜索的最近邻点的最大数量

gp3.setMaximumSurfaceAngle(M_PI / 4); // 45 degrees（pi）最大平面角

gp3.setMinimumAngle(M_PI / 18); // 10 degrees三角形内角最小角度

gp3.setMaximumAngle(2 * M_PI / 3); // 120 degrees三角形内角最大角度

gp3.setNormalConsistency(false); //设置该参数保证法线朝向一致，如果法向量一致，设置为true

// Get result

//设置搜索方法和输入点云

gp3.setInputCloud(cloud_with_normals); //设置输入点云为有向点云cloud_with_normals

gp3.setSearchMethod(tree2); //设置搜索方式为tree2

gp3.reconstruct(triangles); //执行重构，结果保存在triangles中

//保存网格图

pcl::io::savePLYFile(“result.ply”, triangles);

// Additional vertex information

std::vector parts = gp3.getPartIDs();

std::vector states = gp3.getPointStates();

// Finish结束

可视化：

boost::shared_ptrpcl::visualization::PCLVisualizer viewer(new pcl::visualization::PCLVisualizer(“3D Viewer”));

viewer->setBackgroundColor(0, 0, 0); //设置背景

viewer->addPolygonMesh(triangles, “my”); //设置显示的网格

// viewer->setPointCloudRenderingProperties(pcl::visualization::PCL_VISUALIZER_POINT_SIZE, 3, “sample cloud”);

viewer->addCoordinateSystem(1.0); //设置坐标系

viewer->initCameraParameters();

while (!viewer->wasStopped())

{

viewer->spinOnce(100);

boost::this_thread::sleep(boost::posix_time::microseconds(100000));

}

示例

代码如下：

#include <pcl/point_types.h>
#include <pcl/io/pcd_io.h>
#include <pcl/kdtree/kdtree_flann.h>
#include <pcl/features/normal_3d.h>
#include <pcl/surface/gp3.h>
#include <pcl/visualization/pcl_visualizer.h>
#include <boost/thread/thread.hpp>

int
main(int argc, char** argv)
{
	// Load input file into a PointCloud<T> with an appropriate type
	pcl::PointCloud<pcl::PointXYZ>::Ptr cloud(new pcl::PointCloud<pcl::PointXYZ>);
	pcl::PCLPointCloud2 cloud_blob;
	pcl::io::loadPCDFile("F://PCL online data//thirteen//3 greedy_projection//source//table_scene_lms400_downsampled.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>);//定义kd树指针
	tree->setInputCloud(cloud);//用cloud构造tree对象
	n.setInputCloud(cloud);//为法线估计对象设置输入点云
	n.setSearchMethod(tree);//设置搜索方法
	n.setKSearch(30);//设置k邻域搜素的搜索范围
	n.compute(*normals);//估计法线存储结果到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_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);//利用有向点云构造tree

	// Initialize objects
	pcl::GreedyProjectionTriangulation<pcl::PointNormal> gp3;//定义三角化对象
	pcl::PolygonMesh triangles;//存储最终三角化的网络模型

	// Set the maximum distance between connected points (maximum edge length)
	gp3.setSearchRadius(0.025);         //设置搜索半径radius，来确定三角化时k一邻近的球半径。

	// Set typical values for the parameters
	gp3.setMu(2.5);                     //设置样本点到最近邻域距离的乘积系数 mu 来获得每个样本点的最大搜索距离，这样使得算法自适应点云密度的变化
	gp3.setMaximumNearestNeighbors(100);//设置样本点最多可以搜索的邻域数目100 。
	gp3.setMaximumSurfaceAngle(M_PI / 4);  //45 degrees，设置连接时的最大角度 eps_angle ，当某点法线相对于采样点的法线偏离角度超过该最大角度时，连接时就不考虑该点。
	gp3.setMinimumAngle(M_PI / 18);        //10 degrees，设置三角化后三角形的最小角，参数 minimum_angle 为最小角的值。
	gp3.setMaximumAngle(2 * M_PI / 3);    //120 degrees，设置三角化后三角形的最大角，参数 maximum_angle 为最大角的值。
	gp3.setNormalConsistency(false);     //设置一个标志 consistent ，来保证法线朝向一致，如果设置为 true 则会使得算法保持法线方向一致，如果为 false 算法则不会进行法线一致性检查。

	// Get result
	gp3.setInputCloud(cloud_with_normals);//设置输入点云为有向点云
	gp3.setSearchMethod(tree2);           //设置搜索方式tree2
	gp3.reconstruct(triangles);           //重建提取三角化
   // std::cout << triangles;
	// Additional vertex information
	std::vector<int> parts = gp3.getPartIDs();//获得重建后每点的 ID, Parts 从 0 开始编号， a-1 表示未连接的点。
	/*
	获得重建后每点的状态，取值为 FREE 、 FRINGE 、 BOUNDARY 、 COMPLETED 、 NONE 常量，
	其中 NONE 表示未定义，
	FREE 表示该点没有在 三 角化后的拓扑内，为自由点，
	COMPLETED 表示该点在三角化后的拓扑内，并且邻域都有拓扑点，
	BOUNDARY 表示该点在三角化后的拓扑边缘，
	FRINGE 表示该点在 三 角化后的拓扑内，其连接会产生重叠边。
	*/
	std::vector<int> states = gp3.getPointStates();


	//三角网格化的可视化
	boost::shared_ptr<pcl::visualization::PCLVisualizer> viewer(new pcl::visualization::PCLVisualizer("3D Viewer"));
	viewer->setBackgroundColor(255, 255, 0);     //设置背景
	viewer->addPolygonMesh(triangles, "my");

	viewer->addCoordinateSystem(1.0);
	viewer->initCameraParameters();
	// 主循环
	while (!viewer->wasStopped())
	{
		viewer->spinOnce(100);
		boost::this_thread::sleep(boost::posix_time::microseconds(100000));
	}
	// Finish
	return (0);
}

运行显示：

局部区域放大：

谢谢！