PCL使用贪婪三角形算法曲面重构

内容介绍

贪婪投影三角化算法是一种对原始点云进行快速三角化的算法，该算法假设曲面光滑，点云密度变化均匀，不能在三角化的同时对曲面进行平滑和孔洞修复。

方法：

（1）将三维点通过法线投影到某一平面

（2）对投影得到的点云作平面内的三角化

（3）根据平面内三位点的拓扑连接关系获得一个三角网格曲面模型

在平面区域的三角化过程中用到了基于Delaunay的空间区域增长算法，该方法通过选取一个样本三角片作为初始曲面，不断扩张曲面边界，最后形成一张完整的三角网格曲面，最后根据投影点云的连接关系确定各原始三维点间的拓扑连接，所得的三角网格即为重建得到的曲面模型。

该算法适用于采样点云来自表面连续光滑的曲面且点云的密度变化比较均匀的情况。

/*
* GreedyProjection是根据点云进行三角化，而 poisson 则是对water-tight的模型进行重建，
* 所以形成了封闭mesh和很多冗余信息，需要对poisson的重建进行修剪才能得到相对正确的模型
*
*/

#include <pcl/point_types.h>
#include <pcl/io/pcd_io.h>
#include <pcl/io/ply_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>
#include <fstream>
#include <iostream>
#include <stdio.h>
#include <string.h>
#include <string>
#include<pcl/filters/voxel_grid.h>
#include <pcl/surface/mls.h>

int main(int argc, char** argv)
{

	//读取文件
	pcl::PointCloud<pcl::PointXYZ>::Ptr cloud2(new pcl::PointCloud<pcl::PointXYZ>); //创建点云对象指针，用于存储输入
	if (pcl::io::loadPCDFile("sparse_points.pcd", *cloud2) == -1) {
		PCL_ERROR("Could not read pcd file!\n");
		return -1;
	}

	// 体素滤波
	pcl::VoxelGrid<pcl::PointXYZ> vg;
	vg.setInputCloud(cloud2);
	vg.setLeafSize(30.0f, 30.0f, 30.0f);  // 设置体素大小
	pcl::PointCloud<pcl::PointXYZ>::Ptr cloud2_smoothed1(new pcl::PointCloud<pcl::PointXYZ>);
	vg.filter(*cloud2_smoothed1);
	//最小二乘法平滑
	pcl::search::KdTree<pcl::PointXYZ>::Ptr treeSampling(new pcl::search::KdTree<pcl::PointXYZ>);
	pcl::PointCloud<pcl::PointXYZ>::Ptr cloud2_smoothed(new pcl::PointCloud<pcl::PointXYZ>);
	pcl::MovingLeastSquares<pcl::PointXYZ, pcl::PointXYZ> mls;
	mls.setComputeNormals(false);
	mls.setInputCloud(cloud2_smoothed1);
	mls.setPolynomialOrder(3);
	//mls.setPolynomialFit(false);
	mls.setSearchMethod(treeSampling);
	mls.setSearchRadius(150);
	mls.process(*cloud2_smoothed);
	// 估计法向量
	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(cloud2_smoothed);
	n.setInputCloud(cloud2_smoothed);
	n.setSearchMethod(tree);
	n.setKSearch(10);
	n.compute(*normals); //计算法线，结果存储在normals中
	//* normals 不能同时包含点的法向量和表面的曲率

	//将点云和法线放到一起
	pcl::PointCloud<pcl::PointNormal>::Ptr cloud2_with_normals(new pcl::PointCloud<pcl::PointNormal>);
	pcl::concatenateFields(*cloud2_smoothed, *normals, *cloud2_with_normals);
	//* cloud2_with_normals = cloud2 + normals

	//创建搜索树
	pcl::search::KdTree<pcl::PointNormal>::Ptr tree2(new pcl::search::KdTree<pcl::PointNormal>);
	tree2->setInputCloud(cloud2_with_normals);

	//初始化GreedyProjectionTriangulation对象，并设置参数
	pcl::GreedyProjectionTriangulation<pcl::PointNormal> gp3;
	//创建多变形网格，用于存储结果
	pcl::PolygonMesh triangles;

	//设置GreedyProjectionTriangulation对象的参数
	  //第一个参数影响很大
	gp3.setSearchRadius(160.0f); //设置连接点之间的最大距离（最大边长）用于确定k近邻的球半径【默认值 0】
	gp3.setMu(30.0f); //设置最近邻距离的乘子，以得到每个点的最终搜索半径【默认值 0】
	gp3.setMaximumNearestNeighbors(150); //设置搜索的最近邻点的最大数量
	gp3.setMaximumSurfaceAngle(M_PI /1); // 45 degrees（pi）最大平面角
	gp3.setMinimumAngle(M_PI / 36); // 10 degrees 每个三角的最小角度
	gp3.setMaximumAngle(M_PI/1); // 120 degrees 每个三角的最大角度
	gp3.setNormalConsistency(false); //如果法向量一致，设置为true

	//设置搜索方法和输入点云
	gp3.setInputCloud(cloud2_with_normals);
	gp3.setSearchMethod(tree2);

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

	//保存网格图
	//pcl::io::savePLYFile("result2.ply", triangles);

	// Additional vertex information
	//std::vector<int> parts = gp3.getPartIDs();
	//std::vector<int> states = gp3.getPointStates();

	  // 显示结果图
	boost::shared_ptr<pcl::visualization::PCLVisualizer> viewer1(new pcl::visualization::PCLVisualizer("original_cloud2"));
	viewer1->setBackgroundColor(0, 0, 0); //设置背景
	viewer1->addPointCloud(cloud2_smoothed);
	boost::shared_ptr<pcl::visualization::PCLVisualizer> viewer2(new pcl::visualization::PCLVisualizer("rebuild_cloud2"));
	viewer2->setBackgroundColor(0, 0, 0); //设置背景
	viewer2->addPolygonMesh(triangles, "my"); //设置显示的网格
	//viewer->addCoordinateSystem(1.0); //设置坐标系
	//viewer->initCameraParameters();
	while (!viewer1->wasStopped()) {
		viewer1->spinOnce(100);
		viewer2->spinOnce(100);
		boost::this_thread::sleep(boost::posix_time::microseconds(100000));
	}

	return (0);
}

代码实现了下采样和最小二乘法平滑，具体参数大家自己修改即可。