基于KDTree的点云数据双边滤波去噪的方法

技术领域

本发明涉及到测绘科学与技术、逆向工程、古建数字保护领域，具体涉及利用地面激光扫描技术进行对象表面重建时，基于KDTree的点云数据双边滤波去噪的方法。 

背景技术

地面三维激光扫描技术是三维激光扫描技术的非常重要组成部分，在测绘领域有广泛的应用。TLS获取的原始点云数据是密集、离散分布的，又是海量数据的集合，同时存在大量噪声与冗余数据。针对TLS的海量点云数据的去噪，需要有一套行之有效的优化方法，才能达到精化、压缩与平滑点云数据的目的。 

对TLS采集的点云数据可以删除肉眼能识别的噪声点，但物体点云的表面仍然有“浮点”、“坏点”无法去除，需要进行滤波。其中研究比较成熟的方法如中值、均值、高斯及平均曲率流滤波对于海量点云数据的滤波来说，消耗的时间难以承受，效率低；拉普拉斯、维纳与卡尔曼滤波会使模型产生不同程度的收缩，导致模型变形；小波去噪、均值漂移算法效果好，但算法复杂，实现较难。相比较而言，本发明采用的双边滤波算法很好地保持了点云数据所表达的实体表面特征信息，方法简单有效，并且运算速度快。为了克服不能处理大范围噪声的局限，可先通过交互去噪处理大量的噪声，在此基础上，结合KDTree建立点与点拓扑关系，快速搜索最近邻域点，运用双边滤波法去噪，加速滤波速度和效果。 

发明内容

本发明所要解决的技术问题是提供一种将双边滤波方法用于三维点云数据的滤波，主要是处理点云数据中的噪声，即基于KDTree的点云数据双边滤波去噪的方法。 

本发明采用的技术方案是： 

利用地面激光扫描技术进行测量对象表面重建时基于KDTree的点云数据双边滤波去噪的方法，其特征在于，其具体步骤为： 

①在地面三维激光扫描技术中的三维点云数据的滤波中，定义 

p′＝p+λn(1) 

上式中p′为数据点p滤波后的新数据点，λ为双边滤波权因子，n为数据点p的法向，双边滤波权因子λ定义如下 

$\lambda =\frac{\underset{k_{\mathrm{ij}}\in N\left(p_i\right)}{\Sigma }{H}_C\left(\right||p_i-k_{\mathrm{ij}}|\left|\right)H_S(<p_i-k_{\mathrm{ij}}>)(n_i,q_i-k_{\mathrm{ij}})}{\underset{k_{\mathrm{ij}}\in N\left(p_i\right)}{\Sigma }{H}_C\left(\right||p_i-k_{\mathrm{ij}}|\left|\right)H_S(<p_i-k_{\mathrm{ij}}>)}---\left(2\right)$

式中N(p_i)是数据点p_i的邻居点，光顺滤波是标准高斯滤波，表示为 

$H_C\left(x\right)=e^{-x^2/{2\sigma }_C^2}---\left(3\right)$

特征保持权重函数类似于光顺滤波，可定义为 

$H_S\left(y\right)=e^{-y^2/{2\sigma }_S^2}---\left(4\right)$

其中，参数σ_c是数据点p_i到邻居点的距离对该点的影响因子；参数σ_s是数据点p_i到邻近点的距离向量，在该点法向n_i上的投影对数据点q_i的影响因子； 

②散乱数据点的法矢量估计，可以利用数据点领域内的点集进行最小二乘平面拟合，以拟合平面的法矢量作为该点的法矢量；设P为散乱点云中的一个数据点，其领域点集记为k-N(P)，即包含与P距离最近的k个数据点；令拟合平面的法矢量为N，构造目标函数D表示k-N(P)中的点Q到拟合平面的距离： 

$\min D=\min \underset{i=0,Q_i\in k-N\left(P\right)}{\overset{k-1}{\Sigma }}{\left(\right||Q_i-P||·N)}^2---\left(5\right)$

式(5)的求解可转化求下面矩阵M的特征值问题： 

$M=\underset{i=0}{\overset{k-1}{\Sigma }}(Q_i-P){(Q_i-P)}^T---\left(6\right)$

利用数值方法可以得到矩阵M的最小特征值，最小特征值对应的特征向量单位化后就是拟合平面的法矢量N，即散乱数据点P的单位法矢量；它们分别表示数据点切平面方向和法向方向上的高斯滤波常量系数，反映了对任一数据点p实施双边滤波操作时的切向和法向影响范围；可以设置相邻两个数据点的N₁·N₂≈±1的条件，克服法矢量二义性的问题(即法矢量的指向可能不是曲面的同侧)； 

③读取点云数据，通过KDTree找到每一个数据点p_i的m个近邻点k_ij，j＝0，…，m-1； 

④按散乱点法向估计的方法计算每个数据点p_i的法向量n_i； 

⑤计算每个邻近点k_ij的滤波函数参数x＝||p_i-k_ij||，表示点p_i到邻近点 k_ij的距离；求特征保持权重函数的参数y＝<n_i，p_i-k_ij>，表示点p_i与邻近点的距离向量p_i-k_ij与该点法向的内积； 

⑥按照式(3)和(4)求出滤波高斯函数H_C(x )和H_S(y)； 

⑦计算双边滤波权因子λ，将H_C(x )和H_S(y)代入式(2)求解； 

⑧得到滤波后的新数据点p′_i＝p_i+λn_i； 

⑨当所有的数据点都经过更新后，双边滤波结束，把得到的优质数据源提供给三维重建技术，用于测量对象表面的三维重建。 

所述基于KDTree的点云数据双边滤波去噪的方法，其特征在于，所述的KDTree为二叉树，每个节点对应一个长方体盒和一个分割平面，分割平面与一个坐标轴平行，分割平面将长方体盒分割成两个子长方体盒，构成KDTree的两个子节点；基于这样的结构，KDTree非常适合与存储空间中的位置和大小信息，将离散的点云数据存储到KDTree中，可以利用KDTree的这个性质，更好的组织点云的存储结构，达到建立点的拓扑关系，快速遍历领域点集的目的。 

所述基于KDTree的点云数据双边滤波去噪的方法，其特征在于，求得矩阵M的最小特征值的数值方法是逆迭代法。 

本发明的有益效果在于： 

本发明的双边滤波方法用于滤除点云数据的噪声，方法简单有效，并且运算速度快，能够在保持特征的同时去除噪声，适用于处理小范围的噪声。 

附图说明

图1为本发明的点云数据双边滤波去噪的流程图。 

图2为某历史古建筑单站扫描原始点云图。 

图3为交互方法去噪后剩下的噪声点的点云图。 

图4为按本发明的方法实现滤波去噪后的最终点云图。 

具体实施方式

如图1所示，利用地面激光扫描技术进行测量对象表面重建时基于KDTree的点云数据双边滤波去噪的方法，具体步骤为： 

①在地面三维激光扫描技术中的三维点云数据的滤波中，定义 

p′＝p+λn(1) 

上式中p′为数据点p滤波后的新数据点，λ为双边滤波权因子，n为数据点p的法向，双边滤波权因子λ定义如下 

$\lambda =\frac{\underset{k_{\mathrm{ij}}\in N\left(p_i\right)}{\Sigma }{H}_C\left(\right||p_i-k_{\mathrm{ij}}|\left|\right)H_S(<p_i-k_{\mathrm{ij}}>)(n_i,q_i-k_{\mathrm{ij}})}{\underset{k_{\mathrm{ij}}\in N\left(p_i\right)}{\Sigma }{H}_C\left(\right||p_i-k_{\mathrm{ij}}|\left|\right)H_S(<p_i-k_{\mathrm{ij}}>)}---\left(2\right)$

式中N(p_i)是数据点p_i的邻居点，光顺滤波是标准高斯滤波，表示为 

$H_C\left(x\right)=e^{-x^2/{2\sigma }_C^2}---\left(3\right)$

特征保持权重函数类似于光顺滤波，可定义为 

$H_S\left(y\right)=e^{-y^2/{2\sigma }_S^2}---\left(4\right)$

其中，参数σ_c是数据点p_i到邻居点的距离对该点的影响因子；参数σ_s是数据点p_i到邻近点的距离向量，在该点法向n_i上的投影对数据点q_i的影响因子； 

②散乱数据点的法矢量估计，可以利用数据点领域内的点集进行最小二乘平面拟合，以拟合平面的法矢量作为该点的法矢量；设P为散乱点云中的一个数据点，其领域点集记为k-N(P)，即包含与P距离最近的k个数据点；令拟合平面的法矢量为N，构造目标函数D表示k-N(P)中的点Q到拟合平面的距离： 

$\min D=\min \underset{i=0,Q_i\in k-N\left(P\right)}{\overset{k-1}{\Sigma }}{\left(\right||Q_i-P||·N)}^2---\left(5\right)$

式(5)的求解可转化求下面矩阵M的特征值问题： 

$M=\underset{i=0}{\overset{k-1}{\Sigma }}(Q_i-P){(Q_i-P)}^T---\left(6\right)$

利用数值方法可以得到矩阵M的最小特征值，最小特征值对应的特征向量单位化后就是拟合平面的法矢量N，即散乱数据点P的单位法矢量；它们分别表示数据点切平面方向和法向方向上的高斯滤波常量系数，反映了对任一数据点p实施双边滤波操作时的切向和法向影响范围；可以设置相邻两个数据点的N₁·N₂≈±1的条件，克服法矢量二义性的问题(即法矢量的指向可能不是曲面的同侧)； 

③读取点云数据，通过KDTree找到每一个数据点p_i的m个近邻点k_ij，j＝0，…，m-1； 

④按散乱点法向估计的方法计算每个数据点p_i的法向量n_i； 

⑤计算每个邻近点k_ij的滤波函数参数x＝||p_i-k_ij||，表示点p_i到邻近点k_ij的距离；求特征保持权重函数的参数y＝<n_i，p_i-k_ij>，表示点p_i与邻近点的距离向量p_i-k_ij与该点法向的内积； 

⑥按照式(3)和(4)求出滤波高斯函数H_C(x)和H_S(y)； 

⑦计算双边滤波权因子λ，将H_C(x)和H_S(y)代入式(2)求解； 

⑧得到滤波后的新数据点p′_i＝p_i+λn_i； 

⑨当所有的数据点都经过更新后，双边滤波结束，把得到的优质数据源提供给三维重建技术。 

所述的KDTree为二叉树，每个节点对应一个长方体盒和一个分割平面，分割平面与一个坐标轴平行，分割平面将长方体盒分割成两个子长方体盒，构成KDTree的两个子节点；基于这样的结构，KDTree非常适合与存储空间中的位置和大小信息，将离散的点云数据存储到KDTree中，可以利用KDTree的这个性质，更好的组织点云的存储结构，达到建立点的拓扑关系，快速遍历领域点集的目的。 

求得矩阵M的最小特征值的数值方法是逆迭代法。