《点云处理》 点云平滑

前言

通常情况下，使用3D相机得到的点云在深度方向上（z方向）是上下起伏的，这种波动有时候不利于点云真实特征的提取。因此，在进行特征提取之前，可以先对输入的点云进行平滑。

环境：
Windows11 + VS2019 + PCL1.11.1

1.高斯平滑

平滑相关的算子有很多，无论是2D还是3D，都可以使用高斯滤波进行平滑处理。PCL库中也集成了相关的实现方法。

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

代码如下：

/// <summary>
/// 使用高斯滤波进行点云平滑
/// </summary>
/// <param name="cloud">输入点云</param>
/// <param name="cloudout">输出点云</param>
/// <param name="kernels">高斯核数量</param>
/// <param name="cigma">方差</param>
/// <param name="radius">搜索半径</param>
/// <returns>return true表示成功，false表示失败</returns>
bool GaussianSmooth(const pcl::PointCloud<pcl::PointXYZ>::Ptr& cloud, pcl::shared_ptr<pcl::PointCloud<pcl::PointXYZ>>& cloudout, const int& kernels, 
	const double& cigma, const double& radius)
{
	if (cloud == nullptr) return false;
	if (cloud->points.size() < 10) return false;

	pcl::shared_ptr<pcl::filters::GaussianKernel<pcl::PointXYZ, pcl::PointXYZ>> kernel(new pcl::filters::GaussianKernel<pcl::PointXYZ, pcl::PointXYZ>);
	kernel->setSigma(cigma);
	kernel->setThresholdRelativeToSigma(4);
	kernel->setThreshold(kernels);
	pcl::shared_ptr<pcl::search::KdTree<pcl::PointXYZ>> kdtree(new pcl::search::KdTree<pcl::PointXYZ>);
	kdtree->setInputCloud(cloud);
	pcl::filters::Convolution3D <pcl::PointXYZ, pcl::PointXYZ, pcl::filters::GaussianKernel<pcl::PointXYZ, pcl::PointXYZ> > convolution;
	convolution.setKernel(*kernel);
	convolution.setInputCloud(cloud);
	convolution.setSearchMethod(kdtree);
	convolution.setRadiusSearch(radius);
	convolution.setNumberOfThreads(10);            //important! Set Thread number for openMP
	convolution.convolve(*cloudout);

	return true;
}

对线扫相机所获取的一帧点云进行平滑处理，效果如下，绿色点是原始点，红色为平滑后输出点云：

运行时间为0.009s，由此可见对于单帧点云平滑速度和效果还是比较显著的。输入参数是kernels = 5，cigma = 1.4，搜索半径为0.1（这个跟点间距有关，设置为2倍点间距）。

对点云轮廓进行平滑，从下述结果来看，高斯滤波可以起到明显的平滑作用，且并不会破坏掉点云转折处的特征。

扩大输入点云的数量，发现即便输入100w级别的点云，平滑时间也仅需2.4s。

2.MLS平滑

MLS，Moving Least Square移动最小二乘法，也是比较常见的点云滤波算法，该算法还能够计算点云的法向和曲率。

头文件：
#include <pcl/surface/mls.h>

代码：

/// <summary>
/// 使用MLS进行平滑
/// </summary>
/// <param name="cloud">输入点云</param>
/// <param name="cloudout">输出点云</param>
/// <param name="radius">搜索半径</param>
/// <returns>return true表示成功，false表示失败</returns>
bool MLSSmooth(const pcl::PointCloud<pcl::PointXYZ>::Ptr& cloud, pcl::shared_ptr<pcl::PointCloud<pcl::PointXYZ>>& cloudout, const double& radius)
{
	if (cloud == nullptr) return false;
	if (cloud->points.size() < 10) return false;

	pcl::search::KdTree<pcl::PointXYZ>::Ptr tree(new pcl::search::KdTree<pcl::PointXYZ>);
	pcl::PointCloud<pcl::PointNormal>::Ptr mls_points_normal(new pcl::PointCloud<pcl::PointNormal>);
	pcl::MovingLeastSquares<pcl::PointXYZ, pcl::PointXYZ> mls;
	mls.setInputCloud(cloud);
	mls.setComputeNormals(false);					// 是否计算法线，设置为ture则计算法线
	mls.setPolynomialFit(false);					    // 设置为true则在平滑过程中采用多项式拟合来提高精度
	mls.setPolynomialOrder(2);						// 设置MLS拟合的阶数，默认是2
	mls.setSearchMethod(tree);						// 邻域点搜索的方式
	mls.setSearchRadius(radius);				    // 邻域搜索半径
	mls.setNumberOfThreads(10);						// 设置多线程加速的线程数
	mls.process(*cloudout);				            // 曲面重建

	return true;
}

MLS与高斯滤波一样，需要设置一个搜索半径，当搜索半径设置为0.1时，两种平滑方法节拍差不多，MLS甚至要略快：

增大搜素半径为0.3时，可以发现搜索半径增加，两种算法运行时间会变长，但是相对而言，MLS时间消耗更加显著：

注意，下述两个方法如果无必要，可以设置为false，设置为true会更消耗时间：

mls.setComputeNormals(false);					// 是否计算法线，设置为ture则计算法线
mls.setPolynomialFit(false);					// 设置为true则在平滑过程中采用多项式拟合来提高精度

这里存在一个问题：MLS平滑点云，可以获取平滑后的点云，还可以获取点云的法向和曲率，但是比较耗时。前面的测试发现，如果搜索半径比较小，设置为0.1时，那么算法运行的时间比较快。此外，还可以通过体素对点云进行稀释，然后再进行MLS平滑。但是使用体素滤波稀释后的点云，点间距必然增大，搜索半径肯定也要增加，最起码为点距的两倍，那么搜索半径的增加必然会导致运行时间的增加。所以到底是设置一个小的搜索半径直接对大数量级的点云进行处理快，还是稀释之后再进行处理快。接下来要对其做一个测试。

对于一个1000w以上的庞大点云，直接使用MLS和稀释之后再使用MLS，区别很大。运行时间上有明显的差距。

3.中值滤波

PCL中封装的中值滤波用于有序点云，不适用于无序点云。可以通过cloud->isOrganized()查看是否是有序点云。

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

代码：

/// <summary>
/// 中值滤波
/// </summary>
/// <param name="cloud">输入点云</param>
/// <param name="cloudout">输出点云</param>
/// <param name="size">窗口大小</param>
/// <returns>return true表示成功，false表示失败</returns>
bool MedianSmooth(const pcl::PointCloud<pcl::PointXYZ>::Ptr& cloud, pcl::shared_ptr<pcl::PointCloud<pcl::PointXYZ>>& cloudout, const int& size)
{
	if (cloud == nullptr) return false;
	if (cloud->points.size() < 10) return false;
	
	pcl::MedianFilter<pcl::PointXYZ> median_filter;
	median_filter.setInputCloud(cloud);
	median_filter.setWindowSize(size);
	median_filter.filter(*cloudout);

	return true;
}