PCL贪婪投影三角化算法

一、原理

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

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

这里我们设置如下函数及参数：

• 函数SetMaximumNearestNeighbors(int num)和SetMu(double mu)，这两个函数的作用是控制搜索邻域大小。前者定义了可搜索的邻域数num，后者设置了样本点到最近邻域距离的乘积系数mu来获得每个样本点的最大搜索距离（是为了适应点云密度的变化）。
• 函数SetSearchRadius(double radius)，该函数设置了搜索半径radius，来确定三角化时K-邻近的球半径。
• 函数SetMinimumAngle(double min)和SetMaximumAngle(double max)，这两个函数是三角化后每个三角形的最大角和最小角。两者至少要符合一个。典型值分别是10和120度（弧度）。
• 函数SetMaximumSurfaceAgle(double angle)和SetNormalConsistency(bool con)，这两个函数是为了处理边缘或者角很尖锐以及一个表面的两边非常靠近的情况。为了处理这些特殊情况，函数SetMaximumSurfaceAgle(double angle)规定如果某点法线方向相对于采样点的法线偏离超过指定角度（注：大多数表面法线估计方法可以估计出连续变化的表面法线方向，即使在尖锐的边缘条件下），连接时不考虑该点。该角度是通过计算法向线段（忽略法线方向）之间的角度。函数SetNormalConsistency(bool con)保证法线朝向，如果设置为true，会值得算法保持法线方向一致，如果为false算法则不会进行发现一致性检查。

二、代码实现

我用的配置环境是：

    PCL1.8.0+VS2013环境配置

#include <pcl/point_types.h>  
#include <pcl/io/pcd_io.h>  
#include <pcl/io/ply_io.h>  
#include <pcl/io/obj_io.h>
#include <pcl/PolygonMesh.h>
#include <pcl/io/vtk_lib_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>  

typedef pcl::PointXYZ PointT;
typedef pcl::PointCloud<PointT> PointCloudT;
int main(int argc, char** argv)
{
	//打开文件  
	//PointCloudT::Ptr cloud(new PointCloudT);
	//if (pcl::io::loadPCDFile<pcl::PointXYZ>("Bunny.pcd", *cloud) == -1)
	//{
	//	PCL_ERROR("Couldn't read file1 \n");
	//	return (-1);
	//}
	pcl::PolygonMesh mesh;
	pcl::io::loadPolygonFile("cat.obj", mesh);

	pcl::PointCloud<pcl::PointXYZ>::Ptr cloud(new pcl::PointCloud<pcl::PointXYZ>);
	pcl::fromPCLPointCloud2(mesh.cloud, *cloud);
	// 估计法向量  
	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中  
	//* normals 不能同时包含点的法向量和表面的曲率  

	//将点云和法线放到一起  
	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  

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

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

	//设置GreedyProjectionTriangulation对象的参数  
	//第一个参数影响很大  
	gp3.setSearchRadius(200.0f);           //设置连接点之间的最大距离（最大边长）用于确定k近邻的球半径【默认值 0】  
	gp3.setMu(2.5f);                       //设置最近邻距离的乘子，以得到每个点的最终搜索半径【默认值 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  

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

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

	//保存网格图  
	pcl::io::saveOBJFile("result.obj", triangles);

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

	// 显示结果图  
	boost::shared_ptr<pcl::visualization::PCLVisualizer> viewer(new pcl::visualization::PCLVisualizer("3D Viewer"));
	viewer->setBackgroundColor(0, 0, 0);                  //设置背景  
	viewer->addPolygonMesh(triangles, "my");              //设置显示的网格
	//设置网格模型显示模式
	//viewer->setRepresentationToSurfaceForAllActors();   //网格模型以面片形式显示
	//viewer->setRepresentationToPointsForAllActors();    //网格模型以点形式显示
	viewer->setRepresentationToWireframeForAllActors();   //网格模型以线框图模式显示
	viewer->addCoordinateSystem(1.0);                     //设置坐标系  
	viewer->initCameraParameters();
	while (!viewer->wasStopped()){
		viewer->spinOnce(100);
		boost::this_thread::sleep(boost::posix_time::microseconds(100000));
	}

	return (0);
}

效果图如下：

注意：根据每个模型的点的大小不同，设置搜索半径要相对应改变，不然出不来效果图。