基于改进多边滤波的点云三维重建与分析

3D Reconstruction and Analysis of Point Cloud Based on Improved Multilateral Filtering

摘要: 为了提高三维重建过程点云噪声的滤除准确度与重建精确度，提出一种基于尺度分类的点云多边滤波算法，在双边滤波基础上引入尺度分类的自适应参数，将曲率作为第三个滤波影响因子，兼顾滤波过程的噪声滤除与点云平滑。实验结果表明，所提多边滤波算法与双边滤波算法、高斯滤波算法相比，噪声滤除准确度最高，达97.1%。且在重建过程中保证了细节与平滑度，偏离模型真实表面程度最低。

Abstract: In order to improve the filtering accuracy and reconstruction accuracy of point cloud noise in 3D reconstruction process, a multi-lateral filtering algorithm of point cloud based on scale classification was proposed. The adaptive parameters of scale classification were introduced on the basis of bilat-eral filtering, and the curvature was taken as the third filtering influence factor to give considera-tion to noise filtering and point cloud smoothing in the filtering process. The experimental results show that the proposed multi-lateral filtering algorithm has the highest noise removal accuracy of 97.1% compared with bilateral filtering algorithm and Gaussian filtering algorithm. In the recon-struction process, details and smoothness are guaranteed, and the deviation from the real surface of the model is minimal.

1. 引言

随着计算机视觉技术的快速发展，三维点云的曲面重建技术越来越多的被应用到智慧城市、文物保护、虚拟现实、工业零件三维重建 [1] [2] [3] [4] 等场景中。然而点云的实际获取过程中会因为诸多因素引入不同尺度的噪声，因此对三维点云进行滤波去噪是一项具有重要现实意义的工作。袁华等 [5] 提出改进的双边滤波算法，平衡了过度光滑与细节保持的矛盾，但算法效率有所减低。张志斌等 [6] 提出结合多尺度点云特征和随机森林模型的滤波算法，能有效分离植被点与地面点且精度较高，但是应用场景单一。Wu [7] 提出一种基于分解图滤波的自适应点云数据简化方法，大大减少了内存空间的使用，方法有效且可行，但是计算复杂度与分解框架仍有改进空间。Duan [8] 等提出一种基于主成分分析的自适应聚类方法，通过将点云降维、聚类，有效去除点云噪声并降低点云复杂度，但噪声滤除效果仍有待提高。龙家乐等 [9] 提出一种基于图像阈值分割和区域生长法的点云去噪方法，避免了三维数据的复杂计算，但只适用于结构光点云数据。Ren [10] 等提出点云数据的整体滤波算法，通过视觉判断、正态分布检验等验证了算法的可靠性。武军安等 [11] 提出了基于一种基于导向滤波器的边缘区域约束的超分辨率重建算法，解决了TOF (time of flight)相机获取的深度图像分辨率低的问题，但在深度图或灰度图纹理丰富时会抑制算法性能。Digne [12] 提出通过频率分解，基于向量场相似度进行无组织点云去噪的方法。Liu [13] 等提出高斯过程回归(Gaussian Process Regression, GPR)和鲁棒局部加权回归(Robust Locally Weighted Regression, RLWR)相结合的混合回归技术，针对点云滤波具有普适性。

以上方法部分针对点云数据类型具有局限性，部分方法步骤繁琐，本文提出一种对点云数据进行尺度分类并结合改进多边滤波进行精准滤波的方法，使后续重建过程更加精准。

2. 点云的尺度分类方法

点云的噪声一般由传感器本身精度、环境因素、操作员操作精确度、目标物本身特性引入，并会对后续的逆向工程产生干扰，影响点云配准，重建等过程的精度，因此对点云噪声进行分类并有效滤波具有重要意义。噪声按尺度一般分为大尺度噪声和小尺度噪声，大尺度噪声一般表现远离原始点云模型的个别点或点集，小尺度噪声一般存在于原始点云附近，按分布区域又可以分为平坦区域噪声和纹理丰富区域噪声。本文基于点云的欧式距离、密度、曲率三个参数对点云进行分类。

首先基于欧式距离将点云划分为离散区与密集区，然后基于密度与曲率将密集区点云再分为平坦区域纹理丰富区，对离散区点云进行大尺度噪声滤波，对平坦区和纹理丰富区部分区域进行针对性多边滤波，其流程示意图如图1所示。

Figure 1. Algorithm flow chart of this paper

图1. 本文算法流程图

2.1. 基于二叉树的k近邻搜索

二叉树是一种可应用于分割三维点云的高效树形数据结构，在点云的k近邻搜索和范围搜索性能优异。

步骤1：确定点云样本集合 $P\left[p_i|i=1,2,\cdot \cdot \cdot ,n\right]$ 中的几何中心：

$p_c=\frac{1}{n}{\displaystyle {\sum }_{i=1}^n{p}_i}$ (1)

步骤2：以 $p_c$ 为原点并作为初级父节点，根据右手定则建立三维坐标系，三个坐标轴分别记作 $X_1$ 轴， $Y_1$ 轴， $Z_1$ 轴。

步骤3：计算点云在三个坐标轴方向上的方差：

$\{\begin{array}{c}D\left(X_1\right)=\frac{{\displaystyle {\sum }_{i=1}^n{\left(x_i-\frac{1}{n}{\displaystyle {\sum }_{i=1}^n{x}_i}\right)}^2}}{n}\\ D\left(Y_1\right)=\frac{{\displaystyle {\sum }_{i=1}^n{\left(y_i-\frac{1}{n}{\displaystyle {\sum }_{i=1}^n{y}_i}\right)}^2}}{n}\\ D\left(Z_1\right)=\frac{{\displaystyle {\sum }_{i=1}^n{\left(z_i-\frac{1}{n}{\displaystyle {\sum }_{i=1}^n{z}_i}\right)}^2}}{n}\end{array}$ (2)

步骤4：依据方差取值较大两个方向轴所在平面划分左右子空间。

步骤5：在子空间内，选取方差取值最大的方向轴上坐标取值中位数作为下级节点，循环步骤3、4直至点云划分完毕。

基于上述二叉树遍历方法确定点 $p_i$ 的k邻域 ${\delta }_i$ ，其中包含了距离点 $p_i$ 欧式距离最近的k个点。

2.2. 尺度分类参数

基于欧式距离对点云密度进行评估，以便于对点云进行尺度分类，定义 $d\left(p_i,q_j\right)$ 表示任意两点p，q间的欧式距离， $d_{p_i}$ 表示距离点 $p_i$ 距离最近的点，则有：

$d_{p_i}=\min \left(\left[d\left(p_i,q_j\right)|j=1,2,\cdot \cdot \cdot ,n且i\ne j\right]\right)$ (3)

则点云的密度可定义为：

${\bar{d}}_i=\frac{{\displaystyle {\sum }_{i=1}^n{d}_{p{}_i}}}{n}$ (4)

可见，通过欧式距离进行的点云密度评估与点云实际密度呈负相关， $\bar{d}$ 越小，点云越密集。

另一个尺度分类的重要参数是表面曲率，表面曲率 ${\rho }_i$ 的定义如下：

${\rho }_i=\frac{{\lambda }_{i_1}}{{\lambda }_{i_1}+{\lambda }_{i_2}+{\lambda }_{i_3}}$ (5)

其中 ${\lambda }_{i_1},{\lambda }_{i_2},{\lambda }_{i_3}$ 为点云 $p_i$ 的三个特征值且满足 ${\lambda }_{i_1}<{\lambda }_{i_2}<{\lambda }_{i_3}$ ， ${\rho }_i$ 越小说明点云所在邻域越平坦。

2.3. 参数自适应方法

根据前文描述，可以通过两个尺度分类参数 ${\bar{d}}_i,{\rho }_i$ 对点云进行尺度分类，但是两个参数的划分阈值难以选择，下面介绍一种自动选取阈值的方法，能针对不同模型自适应的选择阈值。

设点云的尺度参数序列为 $\left\{{\lambda }_i|i=1,2,3,\cdots ,n\right\}$ ，其中 ${\lambda }_i$ 表示 ${\bar{d}}_i$ 或 ${\rho }_i$ ，选择 ${\lambda }_t$ 作为初始分离阈值 ${\bar{d}}_{th}$ 或 ${\rho }_{th}$ ，这样点云尺度参数被分为两个部分，当 $i\in \left[1,t\right]$ 时，第一部分尺度参数所占比例为：

$S_1=\frac{{\displaystyle {\sum }_{i=1}^t{\lambda }_i}}{{\displaystyle {\sum }_{i=1}^n{\lambda }_i}}$ (6)

同理，当 $i\in \left[t+1,n\right]$ 时，第二部分尺度参数所占比例为：

$S_2=\frac{{\displaystyle {\sum }_{i=t+1}^n{\lambda }_i}}{{\displaystyle {\sum }_{i=1}^n{\lambda }_i}}$ (7)

设置分离函数：

$I={\left(S_1-1\right)}^2+{\left(S_2-1\right)}^2$ (8)

这样，当 $\Delta S=\left|S_1-S_2\right|$ 越大，分离函数取值越大，对点云进行尺度分类的效果也就越好，当分类函数取值序列 $\left\{I_1,I_2,I_3,\cdot \cdot \cdot ,I_n\right\}$ 取最大值时，尺度分类效果最好。当确定欧式距离阈值 ${\bar{d}}_{th}$ 后，将点云模型划分为大、小尺度空间两个部分，确定曲率阈值 ${\rho }_{th}$ 后，再将小尺度空间点云划分为平坦区域纹理丰富区两部分，这样就将点云模型划分为两大类、三小类，以便于后续针对性滤波重建。

3. 多边滤波算法

3.1. 双边滤波算法

1997年，Tomasi和Manduchi [14] 首次提出双边滤波算法并应用于图像的去噪过程中，图像的双边滤波仅需考虑像素点灰度值改变的大小，而点云的滤波过程不仅需要考虑采样点的剔除与稠密化，还要考虑其移动的方向和距离。因此，点云的滤波过程相比图像滤波更为复杂，目前，针对散乱点云的双边滤波算法多采用沿点云法向量移动的形式，具体形式如下：

${p^{\prime }}_i=p_i+d\cdot n_i$ (9)

其中， $p_i$ 原始采样点， ${p^{\prime }}_i$ 为滤波后的输出采样点， $n$ 为点 $p_i$ 的法向量，d为 $p_i$ 沿其法向量移动的距离，文献 [15] 将d定义如下：

$d=\frac{{\displaystyle \underset{p_j\in N(p_i)}{\sum }{W}_c\left(\Vert p_j-p_i\Vert \right)\cdot W_s\left(\Vert n,p_j-p_i\Vert \right)\cdot \langle n,p_j-p_i\rangle }}{{\displaystyle \underset{p_j\in N(p)}{\sum }{W}_c\left(\Vert p_j-p_i\Vert \right)\cdot W_s\left(\Vert p_j-p_i\Vert \right)}}$ (10)

其中， $N\left(p_i\right)$ 是点 $p_i$ 的邻域， $W_c\left(x\right)=\exp \left(-x^2/2{\sigma }_c^2\right)$ ， $W_s\left(x\right)=\exp \left(-x^2/2{\sigma }_s^2\right)$ 分别为距离权重函数和特征保持权重函数， ${\sigma }_c$ 和 ${\sigma }_s$ 分别是空域权重因子和特征权重因子， ${\sigma }_c$ 反映了采样点与其邻域点之间的欧式距离对该采样点的影响，其值与所选邻域点云数量成正比，值越大，滤波后的点云模型越光滑， ${\sigma }_s$ 反映了采样点与其邻域点的欧式距离在其法向量投影对该采样点的影响，其值与滤波过程点云模型的特征保持性能成正相关。虽然在一定邻域内双边滤波算法可以通过采样点沿法向量方向的移动逼近点云模型的表面实现光滑去噪的目的。但也存在一定的弊端，第一，在特征丰富的区域容易混淆特征点与噪点，出现错误滤波的情况，导致原点云模型特征细节被过滤掉；第二，双边滤波算法以k近邻搜索为依据， ${\sigma }_c$ 和 ${\sigma }_s$ 两个固定值的滤波参数不能随各邻域实际分布特性和细节特征情况自适应调整，这样会出现部分区域过度光滑甚至出现孔洞，部分区域欠光滑。

3.2. 改进的多边滤波算法

点云的双边滤波算法虽然具有较好的特征保持能力，但在邻域选择不合适时会出现过度平滑的情况。为了平衡双边滤波算法在滤波平滑与细节特征保持之间的矛盾，在原有的双边滤波算法基础上引入曲率影响因子 ${\rho }_i$ 作为第三个滤波因子，提出法向修正的三边滤波算法，更新点 $p_i$ 处的法向量 $n_i$ 以及 $p_i$ 的坐标：

${n^{\prime }}_i=n_i+\frac{{\displaystyle \underset{p_j\in N(p_i)}{\sum }\Vert n_i-n_j\Vert W_c{W}_s\left(1+{{\rho }^{\prime }}_i\right)}}{{\displaystyle \underset{p_j\in N(p_i)}{\sum }{W}_c{W}_s\left(1+{{\rho }^{\prime }}_i\right)}}$ (11)

这样，可以将邻域内法向量夹角过大且曲率相差较大的点权重限制，避免了采样点在其它方向上的漂移，在法向量得到修正的基础上，更新 $p_i$ 的坐标：

${p^{\prime }}_i=p_i+\frac{{\displaystyle \underset{p_j\in N(p_i)}{\sum }\langle {n^{\prime }}_i,p_i-p_j\rangle W_c{W}_s\left(1+{{\rho }^{\prime }}_i\right)}}{{\displaystyle \underset{p_j\in N(p_i)}{\sum }{W}_c{W}_s\left(1+{{\rho }^{\prime }}_i\right)}}$ (12)

其中：

$W_c\left(x\right)=\exp \left(-{\Vert p_i-p_j\Vert }^2/2{\sigma }_c^2\right)$ (13)

$W_s=\{\begin{array}{l}\exp \left(-{\langle {n^{\prime }}_i,p_i-p_j\rangle }^2/2{\sigma }^2\right),n_i\cdot n_j>{V^{\prime }}_i\\ 0,n_i\cdot n_j\le {V^{\prime }}_i\end{array}$ (14)

且

${{\rho }^{\prime }}_i=\frac{{\rho }_i-{\rho }_{\min }}{{\rho }_{\max }-{\rho }_{\min }}{\rho }_i$ (15)

${\rho }_{\min }$ 和 ${\rho }_{\max }$ 分别是邻域内点曲率值的最小值和最大值， ${{\rho }^{\prime }}_i$ 意义在于曲率大、噪声尺度大的区域采样点可移动距离范围大，可以快速有效滤除大尺度噪声。 ${V^{\prime }}_i$ 为采样点与邻域内点法向量夹角的阈值，其作

用是将邻域内实际影响点限制在法向量方向相近的区域内。 ${\displaystyle \underset{p_j\in N(p_i)}{\sum }{W}_c{W}_s\left(1+{{\rho }^{\prime }}_i\right)}$ 是正则化函数，经过正则

化后，采样点在更新法向量方向移动的距离量就限定在邻域点法向量在采样点法向量方向投影的最小与最大值之间，因此不会引起过度的平滑，不会引入新的噪声。

4. 实验结果及评价

本文采用Artec 3D公司的涡轮与曲轴箱点云模型，为了充分验证本文算法以及对比算法的滤波性能，对点云数据添加随机噪声，进行对比实验。

4.1. 滤波效果对比

图2、图3分别为涡轮与曲轴箱点云经高斯滤波、双边滤波、本文算法滤波的输出结果，并对输出点云的表面密度对点云赋色，不同的颜色代表点云表面密度不同。由图可知，与高斯滤波、双边滤波算法对比，本文算法的输出结果在具有相同特征的区域点云密度均匀，变化最小，有益于重建光滑均匀的表面模型。表1则表示点云模型滤波前后的点云数量。在不考虑错误滤波与引入新噪声的情况下，高斯滤波、双边滤波、本文算法针对噪声滤除的准确度分别为92%、94.3%、97.1%，可见，本文算法由于曲率影响因子 ${\rho }_i$ 的引入，在多边滤波平滑点云的过程兼顾了大尺度噪声的滤除，使噪声滤除准确度进一步提高。

Table 1. Number of point clouds before and after filtering

表1. 滤波前后点云数量统计

(a) (b) (c) (d)

Figure 2. The filtering result of turbine (a) point cloud with noises; (b) result of gaussian filtering; (c) result of bilateral filtering; (d) result of our algorithm

图2. 涡轮滤波结果(a) 带噪声点云；(b) 高斯滤波结果；(c) 双边滤波结果；(d) 本文算法结果

(a) (b) (c) (b)

Figure 3. The filtering result of crankcase (a) point cloud with noises; (b) result of gaussian filtering; (c) result of bilateral filtering; (d) result of our algorithm

图3. 曲轴箱滤波结果(a) 带噪声点云；(b) 高斯滤波结果；(c) 双边滤波结果；(d) 本文算法结果

4.2. 算法效率对比

针对两个点云模型做10次对比实验，各算法的运行时间如图4所示。在涡轮模型的算法效率对比实验中，高斯滤波、双边滤波、本文算法的平均运行时间分别为24.504s、20.348s、21.733s。在曲轴箱模型的算法效率对比实验中，高斯滤波、双边滤波、本文算法的平均运行时间分别为339.175s、275.820s、294.536s。高斯滤波算法效率最低，其效率与算法实现过程距离统计的期望和标准差、邻域大小的选择相关。本文算法由于在双边滤波算法基础上引入了参数自适应分离函数I与曲率影响因子 ${\rho }_i$ ，算法复杂性提高，因此时间效率较双边滤波算法有所降低。

(a) (b)

Figure 4. Algorithm running time (a) turbine; (b) crankcase

图4. 算法运行时间(a) 涡轮；(b) 曲轴箱

4.3. 重建结果对比

将滤波后的点云使用屏蔽泊松算法重建其表面模型，如图5、图6所示，由于高斯滤波只能滤波点云模型明显的离群点，对距离模型真实表面的小尺度噪声滤波效果不佳，所以经高斯滤波后点云模型仍存在少数大尺度噪声和较多小尺度噪声，导致重建结果部分细节处不够光滑。双边滤波则能兼顾大尺度与小尺度噪声，但 $\sigma {}_c$ 和 ${\sigma }_s$ 两个固定值的滤波参数不能随各邻域实际分布特性和细节特征情况自适应调整，这样会出现部分区域过度光滑，部分区域欠光滑，导致重建结果出现孔洞或者细节较差。结合滤波效果对比实验结果，本文算法未引入新噪声且兼顾了噪声滤除与平滑要求，根据尺度分类对不同区域点云进行针对性滤波，可以滤除绝大部分噪声点，重建结果细节处表现良好。

(a) (b) (c)

Figure 5. Reconstruction result of turbine (a) reconstruction result of gaussian filtering; (b) reconstruction result of bilateral filtering; (c) result of our algorithm after filtering and reconstruction

图5. 涡轮点云重建结果(a) 高斯滤波重建结果；(b) 双边滤波重建结果；(c) 本文算法滤波重建结果

(a) (b) (c)

Figure 6. Reconstruction result of crankcase (a) reconstruction result of gaussian filtering; (b) reconstruction result of bilateral filtering; (c) result of our algorithm after filtering and reconstruction

图6. 曲轴箱点云重建结果(a) 高斯滤波重建结果；(b) 双边滤波重建结果；(c) 本文算法滤波重建结果

为了定量描述各滤波算法的准确性，将重建结果与原始点云模型进行三维比较，衡量重建结果偏离原始点云模型的程度，设置可接受偏差范围 $\left[-0.01mm,+0.01mm\right]$ 。如图7、图8所示，其中红色阴影部分是可接受偏差范围内点云比例。涡轮点云模型经高斯滤波、双边滤波、本文算法滤波并重建后可接受偏差比例分别为8.96%、51.60%、72.02%，曲轴箱点云模型经高斯滤波、双边滤波、本文算法滤波并重建后可接受偏差比例分别为43.73%、57.06%、63.55%。因为本文所提算法引入了曲率影响因子作为第三个滤波影响因子，可以限制滤波过程的作用区域且不会过度平滑点云。因此，本文算法滤波结果重建模型偏离真实模型程度最小。

除可接受偏差比例外，重建结果的其它误差项统计如表2、表3所示。其中本文算法重建结果在正向最大偏差、负向最大偏差项均不是最小值，因为这两项是统计结果的极值，具有偶然性。其它针对全局的误差项统计，本文算法表现最优。

(a) (b) (c)

Figure 7. Deviation statistics of turbine reconstruction result (a) gaussian filtering; (b) bilateral filtering; (c) ours

图7. 涡轮重建结果偏差统计(a) 高斯滤波；(b) 双边滤波；(c) 本文算法

(a) (b) (c)

Figure 8. Deviation statistics of crankcase reconstruction result (a) gaussian filtering; (b) bilateral filtering; (c) ours

图8. 曲轴箱重建结果偏差统计(a) 高斯滤波；(b) 双边滤波；(c) 本文算法

Table 2. Statistics of turbine reconstruction results

表2. 涡轮重建结果误差统计

Table 3. Statistics of crankcase reconstruction results

表3. 曲轴箱重建结果误差统计

5. 结论

本文从点云尺度分类的角度出发对双边滤波算法提出改进，提出一种新的点云参数自适应分类的点云多边滤波算法，提高了点云的噪声滤除准确度与重建准确度。在小尺度点云噪声分类过程中引入自适应分离参数完成尺度的自动划分。经过多组对比，本文算法的噪声滤除效果准确度对比高斯滤波、双边滤波分别提高5.54%、2.97%。且重建精确度在对比实验中表现最优。为点云的滤波与三维重建相关研究提供了参考。下一步工作将针对算法的运行效率与尺度分类的准确性进行研究，以提高算法的适用性。