《点云处理》 点云去噪

前言

通常从传感器（3D相机、雷达）中获取到的点云存在噪点（杂点、离群点、孤岛点等各种叫法）。噪点产生的原因有不同，可能是扫描到了不想要扫描的物体，可能是待测工件表面反光形成的，也可能是相机内部的原因。在进行其他算法处理之前，通常需要先去除噪声，避免带来干扰。去噪的方法有很多，效率和效果也是各不相同，应用场景也不太一样，本篇内容就是想要将不同的去噪方法进行归纳。

环境：

Windows + VS2019 + PCL1.11.1

在开始之前，先贴出用于测试的点云，可以看出来空间中漂浮了很多噪点。原始点云总数为5882482个点，x和y方向间距为0.05。

1.半径滤波

PCL中集成有半径滤波，可以用于噪点的去除。主要需要设定两个参数，一个是搜索半径，另一个是在搜索半径内近邻点的最小数量。这两个参数需要根据点云的x，y，z方向上的点间距来设定（或者说分辨率）。一般使用线扫相机或者结构光相机这类3D相机得到的点在x，y方向上间距都比较均匀。假设点间距是0.05，那么将搜索半径设置为0.1，则搜索半径内理论上在上下左右方向上分别有2个点，共8个点，所以近邻点最起码也有8个。

头文件：
#include <pcl/filters/radius_outlier_removal.h>

代码如下：

/// <summary>
/// 使用半径滤波进行点云噪点去除
/// </summary>
/// <param name="cloud">输入点云</param>
/// <param name="cloud_out">输出除噪后的点云</param>
/// <param name="radius">搜索半径</param>
/// <param name="neighbors">单个搜索内近邻点数</param>
/// <returns>return true表示去噪成功，return false表示失败</returns>
bool NoiseRemoveROR(const pcl::PointCloud<pcl::PointXYZ>::Ptr& cloud, pcl::PointCloud<pcl::PointXYZ>::Ptr& cloud_out, const float& radius,
	const int& neighbors)
{
   
	if (cloud == nullptr) return false;
	if (cloud->points.size() < 10) return false;

	pcl::RadiusOutlierRemoval<pcl::PointXYZ> rorfilter(true); // 初始化为true以获取被移除的索引
	rorfilter.setInputCloud(cloud);
	rorfilter.setRadiusSearch(radius);                 // 设置搜索半径
	rorfilter.setMinNeighborsInRadius(neighbors);      // 设置搜索半径内满足的点数
	//rorfilter.setUserFilterValue(1);                 // 将不符合要求的点更改为新的值，而不是去去除它们，与setKeepOrganized一起使用
	//rorfilter.setKeepOrganized(false);               // 保证有序点云的结构，与setUserFilterValue一起使用
	rorfilter.filter(*cloud_out);                      // 进行离群点去除
	return true;
}

一般近邻点的数量设置与搜索半径有直接关系。而经过测试发现，搜索的近邻点数量需要越多则耗时越长。较小的搜索半径只需要设置较小数量的近邻点，但是这样可能会导致有的噪点去除不干净。稍大一点的搜索半径以及近邻点数量适当调大可以得到更稳定的去噪效果。因此，效果和效率之间需要找到一个平衡点。
测试及结果：搜索半径设置为0.2，近邻点数量为20，耗时10.94s。效果如下：

针对于半径滤波耗时问题其实可以将点云使用直通滤波分成若干份，然后使用OpenMP对这几份点云并行处理进行半径滤波，将结果再合并成一个点云，这个方法可能能够起到加速的作用。

有一篇论文《Fast Radius Outlier Filter Variant for Large Point Clouds》，更快速的半径滤波，还未找到源码，若有人找到请联系我，感谢。

2.统计学滤波

统计学滤波也是一种比较常见的去噪算法。其官方解释如下：
The statistical outlier removal process is a bit more refined. First, for every point, the mean distance to its K neighbors is computed. Then, if we asume that the result is a normal (gaussian) distribution with a mean μ and a standard deviation σ, we can deem it safe to remove all points with mean distances that fall out of the global mean plus deviation. Basically, it runs a statistical analysis of the distances between neighboring points, and trims all which are not considered “normal” (you define what “normal” is with the parameters of the algorithm).

个人理解就是遍历点云，对于每个点，都先搜索与其最相近的k个点，计算这k个点与该点的距离，并得到一个平均距离和一个平均距离的标准差。然后比较距离是否是大于μ+stddev*σ，如果大于则表示是离群点。

头文件：#include <pcl/filters/statistical_outlier_removal.h>

代码如下：

/// <summary>
/// 使用PCL中集成的统计学滤波进行去噪
/// </summary>
/// <param name="cloud">输入点云</param>
/// <param name="cloud_out">输出点云</param>
/// <param name="k">搜索近邻点的个数</param>
/// <param name="stddev">平均距离标准差的乘数</param>
/// <returns>return true表示去噪成功，return false表示失败</returns>
bool NoiseRemoveSOR(const pcl::PointCloud<pcl::PointXYZ>::Ptr& cloud, pcl::PointCloud<pcl::PointXYZ>::Ptr& cloud_out, const int& k, const double& stddev)
{
   
	if (cloud == nullptr) return false;
	if (cloud->points.size() < 10) return false;
	pcl::StatisticalOutlierRemoval<pcl::PointXYZ> filter; // Filter object.
	filter.setInputCloud(cloud);
	filter.setMeanK(k);                  // Set number of neighbors to consider to K.
	filter.setStddevMulThresh(stddev);   // Points with a distance larger than stddev standard deviation of the mean distance will be outliers.
	filter.filter(*cloud_out);
	return true;
}

使用这个统计学滤波感觉还是十分耗时，但是效果确实是不错的。
测试及结果：近邻点k值的数量为20，系数设定为1.0，耗时11.07s。效果如下：

有一篇论文《Fast statistical outlier removal based method for large 3d point clouds of outdoor environments》，根据其题目理解来说是一种更快的统计学滤波方法去除离群点，论文还没看，源码也没找到，谁能找到麻烦联系我，感谢！

3.RANSAC

没错，就是RANSAC。使用RANSAC去去噪需要满足条件，那就是目标点云是具有几何特征的。如果目标点云是一个平面，那么就可以使用RANSAC拟合一个平面，并且将距离平面较远的点（外点）去除。这样当然也可以达到去噪的效果，而且速度还比较快。

头文件：
#include <pcl/kdtree/kdtree.h>
#include <pcl/segmentation/sac_segmentation.h>

/// <summary>
/// 使用PCL中集成的用于点云分割的RANSAC方法进行平面拟合
/// </summary>
/// <param name="cloud_in">输入待拟合的点云</param>
/// <param name="inliers">RANSAC拟合得到的内点</param>
/// <param name="coefficients">得到的平面方程参数</param>
/// <param name="iterations">平面拟合最大迭代次数</param>
/// <param name="threshold">RANSAC拟合算法距离阈值</param>
void SEG_RANSAC(const pcl::PointCloud<pcl::PointXYZ>::Ptr& cloud_in, pcl::PointIndices::Ptr& inliers, Eigen::VectorXf& coefficients,
	const int& iterations, const double& threshold)
{
   
	if (inliers == nullptr) inliers.reset(new pcl::PointIndices);
	pcl::ModelCoefficients::Ptr coefficients_m(new pcl::ModelCoefficients);
	pcl::shared_ptr<pcl::search::KdTree<pcl::PointXYZ>> tree(new pcl::search::KdTree<pcl::PointXYZ>);
	tree->setInputCloud(cloud_in);
	pcl::SACSegmentation<pcl::PointXYZ> seg;
	seg.setOptimizeCoefficients(true);
	seg.setModelType(pcl::SACMODEL_PLANE);
	seg.setMethodType(pcl::SAC_RANSAC);
	seg.setMaxIterations(iterations);                     // 设置最大迭代次数
	seg.setDistanceThreshold(threshold);                  // 设定阈值
	seg.setNumberOfThreads(10);                           // 设置线程数量
	seg.setSamplesMaxDist(3, tree);
	seg.setInputCloud(cloud_in);
	seg.segment(*inliers, *coefficients_m);

	coefficients.resize(4);
	coefficients[0] = coefficients_m->values[0]; coefficients[1] = coefficients_m->values[1];
	coefficients[2] = coefficients_m->values[2]; coefficients[3] = coefficients_m->values[3];
	std::cout << "SEG coefficients: " << coefficients[0] << ", " << coefficients[1] << ", " << coefficients[2] << ", " << coefficients[3] << std::endl;
}

运行时间只需要1.5s，迭代了200次。速度比半径滤波和统计学滤波要快很多。但是这个方法也是很有局限性的，只适合特定点云才能用，而且如下图所示，距离平面较近的噪点无法去除，所以是存在去除噪点不干净的情况，如果追求速度，而对去噪要求没那么高，则可以考虑使用该方法。

4.欧式聚类

欧式聚类既可以用于分割，也可以用于去噪，其实跟上述半径滤波区别不大。噪点肯定是距离想要的点云比较“远”的，设置好minSize，把想要的点聚成一个类，噪点自然就去除了。

代码如下：

/// <summary>
/// PCL中集成的欧式聚类
/// </summary>
/// <param name="cloud">输入点云</param>
/// <param name="cluster_indices">聚类索引的数组</param>
/// <param name="tolerance ">距离阈值</param>
/// <param name="MinClusterSize">单个类最少的点数</param>
/// <param name="MaxClusterSize">单个类最大的点数</param>
/// <returns>return true表示有聚类结果，return false表示聚类失败</returns>
bool NoiseRemoveEC(const pcl::PointCloud<pcl::PointXYZ>::Ptr& cloud, std::vector<pcl::PointIndices>& cluster_indices, const double& tolerance,
	const int& MinClusterSize, const int& MaxClusterSize)
{
   
	if (cloud == nullptr) return false;
	if (cloud->points.size() < 10) return false;
	int maxSize = MaxClusterSize;
	if (maxSize < 0) maxSize = cloud->points.size();
	if (MinClusterSize > maxSize) return false;

	pcl::shared_ptr<pcl::search::KdTree<pcl::PointXYZ>> tree(new pcl::search::KdTree<pcl::PointXYZ>); // Creating the KdTree object for the search method of the extraction
	tree->setInputCloud(cloud);

	cluster_indices.clear();
	pcl::EuclideanClusterExtraction<pcl::PointXYZ> ec;    // 欧式聚类对象
	ec.setClusterTolerance(tolerance);                       // 设置近邻搜索的搜索半径为  单位是m
	ec.setMinClusterSize(MinClusterSize);                 // 设置一个聚类需要的最少的点数目
	ec.setMaxClusterSize(maxSize);                        // 设置一个聚类需要的最大点数目
	ec.setSearchMethod(tree);                             // 设置点云的搜索机制
	ec.setInputCloud(cloud);                              // 输入聚类点云
	ec.extract(cluster_indices);                          // 从点云中提取聚类，并将点云索引保存在cluster_indices中

	if (cluster_indices.empty()) return false;
	return true;
}

设置的距离阈值为0.1，是两倍的点距大小（x，y方向点距都是0.05）。耗时10.9s。当距离阈值tolerance设置的比较大时，就会特别耗时，该方法就几乎不可用了。