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

一、概述

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

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

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

二、边界特征检测

目前，在点云的边界特征检测（网格模型的边界特征检测已经是一个确定性问题了，见 网格模型边界检测）方面，PCL中有一个针对点云边界的可以称作为是$state-of-the-art$的方法，这个方法出自 $DetectingHolesinPointSetSurfaces$这篇论文，叫做 $AngleCriterion$，简称 $AC$。这篇论文中还提出了另一种检测边界的方法是$HalfdiscCriterion$，$HdC$。目前PCL中应该只集成了$AC$，因为这个方法确实比$HdC$好，已经够用了。这两种方法的思路都非常简单，但是却非常有效，而往往流传下来的经典方法都是这种简单有效的方法。本文只介绍$AC$方法。

1. AC方法思路

• 假设当前点为$p$，对其邻域$N_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\in N_p$）；
• 找到$\Theta$中最大的夹角${\theta }_{max}=max(\Theta )$；
• 当${\theta }_{max}$越大，那么$p$就越可能为边界点。可以通过设置一个角度阈值 $\xi$，当 ${\theta }_{max}>\xi$ 时，判定$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边界检测结果