PCL点云边界特征检测 （附完整代码 C++）

一、概述

点云特征在定义上（以我个人理解）大致可以分为两大类：

一类是类似于深度学习的featrue map意义，通过计算一些算子来描述点云局部，这种描述只是一种标识符，并没有实际的几何意义，比如 PFH或者 FPFH 之类的，它们只是通过对每个点的局部邻域计算一个 125 125 125 维或者 33 33 33 维的向量来描述当前点，这跟机器学习中的 f e a t u r e feature feature 是一样的，通过这类特征可以用来做点云的配准（其中某些特征可以进一步处理【模式识别】得到几何特征，因为某一类的几何特征在这种 f e a t u r e feature feature的表现形式上通常都会有某种固定的模式）；

而另一类则是有实际的几何意义的，那么这类特征就是叫做几何特征了，通常几何特征包括 { 边界特征，折痕特征，尖锐点 } 等，其中 折痕特征 又包含了 { 凹特征【谷线】和 凸特征【脊线】}。对于一个模型，如果有了几何特征，那么我们就可以得到这个模型的轮廓或者框架，这对模型的压缩表达是非常有意义的。而这类几何特征只要通过合理的处理当然也可以用在特征匹配上。

二、边界特征检测

目前，在点云的边界特征检测（网格模型的边界特征检测已经是一个确定性问题了，见 网格模型边界检测）方面，PCL中有一个针对点云边界的可以称作为是 s t a t e − o f − t h e − a r t state-of-the-art state−of−the−art的方法，这个方法出自 D e t e c t i n g   H o l e s   i n   P o i n t   S e t   S u r f a c e s Detecting\space Holes\space in\space Point\space Set\space Surfaces Detecting Holes in Point Set Surfaces这篇论文，叫做 A n g l e   C r i t e r i o n Angle\ Criterion Angle Criterion，简称 A C AC AC。这篇论文中还提出了另一种检测边界的方法是 H a l f d i s c   C r i t e r i o n Halfdisc\ Criterion Halfdisc Criterion， H d C HdC HdC。目前PCL中应该只集成了 A C AC AC，因为这个方法确实比 H d C HdC HdC好，已经够用了。这两种方法的思路都非常简单，但是却非常有效，而往往流传下来的经典方法都是这种简单有效的方法。本文只介绍 A C AC AC方法。

1. AC方法思路

• 假设当前点为 p p p，对其邻域 N p N_p Np​向其切平面投影，注意经过投影后 p p p的位置并没有改变；
• 然后将投影后的近邻点以 p p p为角点，按照逆时针顺序（当然顺时针也可），两两与 p p p连接形成一系列夹角 Θ = { θ 1 , θ 2 , … , θ n } ， 其 中 n = ∥ N p ∥ − 1 \Theta=\{\theta_1,\theta_2, \dots,\theta_n\}，其中n = \|N_p\| - 1 Θ={θ1​,θ2​,…,θn​}，其中n=∥Np​∥−1（因为 p ∈ N p p\in N_p p∈Np​）；
• 找到 Θ \Theta Θ中最大的夹角 θ m a x = m a x ( Θ ) \theta_{max}=max(\Theta) θmax​=max(Θ)；
• 当 θ m a x \theta_{max} θmax​越大，那么 p p p就越可能为边界点。可以通过设置一个角度阈值 ξ \xi ξ，当 θ m a x > ξ \theta_{max}>\xi θmax​>ξ 时，判定 p p p为边界点。
  
  如图所示，第一行为非边界点邻域经过投影后，近邻点与当前点 p p p形成一系列夹角的示意图，其中红色点为 p p p。第二行为边界点邻域经投影后形成夹角的示意图。可以看到当 p p p为边界点时，其近邻与之形成的最大夹角比 p p p为非边界点时形成的最大夹角要大很多，所以可以通过这个性质来对 p p p进行判定。

2. 代码

	#include <pcl/features/boundary.h>
	#include <pcl/visualization/cloud_viewer.h>
	#include <pcl/io/ply_io.h>
	
	int main()
	{
		/*输入点云和法线*/
		pcl::PointCloud<pcl::PointXYZ>::Ptr cloud(new pcl::PointCloud<pcl::PointXYZ>);
		pcl::PointCloud<pcl::Normal>::Ptr normal(new pcl::PointCloud<pcl::Normal>);
		pcl::PLYReader reader;
		reader.read("*.ply", *cloud);
		reader.read("*.ply", *normal);
		
		/*pcl计算边界*/
		pcl::PointCloud<pcl::Boundary>::Ptr boundaries(new pcl::PointCloud<pcl::Boundary>); //声明一个boundary类指针，作为返回值
		boundaries->resize(cloud->size()); //初始化大小
		pcl::BoundaryEstimation<pcl::PointXYZ, pcl::Normal, pcl::Boundary> boundary_estimation; //声明一个BoundaryEstimation类
		boundary_estimation.setInputCloud(cloud); //设置输入点云
		boundary_estimation.setInputNormals(normal); //设置输入法线
		pcl::search::KdTree<pcl::PointXYZ>::Ptr kdtree_ptr(new pcl::search::KdTree<pcl::PointXYZ>); 
		boundary_estimation.setSearchMethod(kdtree_ptr); //设置搜寻k近邻的方式
		boundary_estimation.setKSearch(30); //设置k近邻数量
		boundary_estimation.setAngleThreshold(M_PI * 0.6); //设置角度阈值，大于阈值为边界
		boundary_estimation.compute(*boundaries); //计算点云边界，结果保存在boundaries中

		/*可视化*/
		pcl::PointCloud<pcl::PointXYZRGB>::Ptr cloud_visual(new pcl::PointCloud<pcl::PointXYZRGB>);
		cloud_visual->resize(cloud->size());
		for(size_t i = 0; i < cloud->size(); i++)
		{
			cloud_visual->points[i].x = cloud->points[i].x;
			cloud_visual->points[i].y = cloud->points[i].y;
			cloud_visual->points[i].z = cloud->points[i].z;
			if(boundaries->points[i].boundary_point != 0)
			{
				cloud_visual->points[i].r = 255;
				cloud_visual->points[i].g = 0;
				cloud_visual->points[i].b = 0;
			}
			else
			{
				cloud_visual->points[i].r = 255;
				cloud_visual->points[i].g = 255;
				cloud_visual->points[i].b = 255;
			}
		}
		pcl::visualization::CloudViewer viewer("view");
		viewer.showCloud(cloud_visual);
		system("PAUSE");
		
		return 0;
	}

3. AC边界检测结果