CN109859114B - 基于局域平滑性和非局域相似性的三维点云修复方法 - Google Patents

Abstract

本发明提出了一种基于局域平滑性和非局域相似性的三维点云修复方法，涉及图信号处理领域，包括的主要步骤有：点云的体素化、点云孔洞的检测、搜索相似区域、相对位置的匹配、修复缺失区域。本方法能够很好地修复复杂的几何结构，减少几何结构的损失，不产生伪影，修复保真度高。

Description

基于局域平滑性和非局域相似性的三维点云修复方法

技术领域

本发明涉及点云修复领域，其理论基础涉及图信号处理领域，具体涉及一种基于局域平滑性和非局域相似性的三维点云修复方法。

背景技术

使用三维扫描仪和深度感测采集到的点云数据会不可避免地产生数据误差：对于一些扫描不到的区域，在扫描结果中会以缺失数据的孔洞的形式呈现。缺失数据的原因包括由于扫描角度的不完全，激光扫描仪本身的局限性，以及物体几何结构过于复杂繁复等。

在孔洞检测方面，目前有两类主流算法：基于网格的方法和基于点的方法。

基于网格的方法的主要是在点云上构建网格，Emelyanov等人、Orriols等人和Floater等人通过一些贪婪策略提取孔边界，并将未被三角形链包围的点视为边界特征点，但这类方法由于需要在大量的点之间建立边和面的连接关系，所以计算复杂度很高，并且当点分布不均匀时，提取的边界不够准确。

基于点的方法主要通过分析点周围的邻接关系直接在点云上进行检测。Nguyen等人、Aldeeb等人和Chen等人首先提取三维点云的轮廓，然后通过生长函数检测孔洞，但这种方法不能将缺失区域的边缘与物体本身含有的封闭边缘分开。Bendels等人计算每个点的边界概率，并以最短的成本路径方式构造围绕孔的闭环；He等人和Wu等人检测每个点的与其K近邻之间的边的夹角是否达到阈值来确定孔洞的边界点。但这些方法无法区分稀疏区域和孔洞，并且检测到的边界点可能是不完整和不连贯的。

在点云修复方面，目前有两类主流算法：基于搜索线上数据库的方法和基于孔洞周围局域信息的方法。

基于搜索线上数据库的方法主要通过寻找外界资源中相似的数据来填充缺失区域。Sahay等人先尝试通过参考缺失区域周围的几何信息，从已有的点云库中寻找相似数据进行填充；后又提出另一种将点云投影到深度图上的方法，使用字典学习从在线的深度图数据库中搜索相似的深度图，将未知部分对应的信息粘贴到空缺区域，最后将补好的深度图投影回点云。此类方法对被搜索的数据库要求较高，需要足够的与待修补点云相似的序列，且直接填充难免在衔接处产生几何结构的不对齐，这是在最后平滑处理时难以解决的问题。此外，在深度图方法中的两次投影过程将不可避免地产生几何结构的损失。

基于孔洞周围局域信息的方法主要参照当前点云本身的数据，通过对缺失区域周围的点的分析来预测缺失区域内部的几何结构。Shankar等人通过提取孔洞周围的边界点来计算其几何特征，并由此拟合出缺失部分的曲面，通过采样插值得到填充的点；Lin等人则将孔洞分解为几个子洞后，分别利用张量投票方法提取其几何特征单独填充；Lozes等人参考孔洞周围点云的信息列出修复孔洞的优化方程，通过偏微分方法求解来进行修补；Wu等人主要针对人体数据中趋于平面的条状孔洞，采用中值算子和平滑收敛算法将边界趋势进行延伸来填补。这类方法完全依靠孔洞周围的局部信息进行修补，由于参考信息比较单一，所以其修补结果比实际内容要更趋于平面，容易失真，且几何结构复杂时，修补的边界容易产生伪影。

发明内容

针对上述问题和相关方法的缺陷，本发明提出了一种基于局域平滑性和非局域相似性的三维点云修复方法，能够很好地修复复杂的几何结构，减少几何结构的损失，不产生伪影，修复保真度高。

为解决上述技术问题，本发明采用如下技术方案：

一种基于局域平滑性和非局域相似性的三维点云修复方法，包括以下步骤：

1.点云的体素化：

(1)将输入的点云P^o坐标进行放缩，使相邻两点间距离接近1，得到P^c。

(2)平移P^c使点云与坐标原点对齐，得到P^t。

(3)将点云分割成边长为1的体素。

(4)计算各体素的坐标值和对应附加属性值，得到点云P^v。

2.点云孔洞的检测：

(1)将输入的三维点云P^v通过向量r_p投影为二维深度图。

(2)在深度图上执行广度优先搜索，标注不同的孔洞。

(3)过滤掉可能不是孔洞的孔洞。

(4)计算孔洞在投影方向r_p上的坐标范围，得到孔洞的第三维信息，即缺失区域数据。

3.搜索相似区域：

(1)搜索相似区域，将输入的点云P^v分割成互相有重叠部分的立方块{c₁,c₂,…,c_n}，作为衡量点云的几何信息的单位。

(2)选择包含缺失区域的立方块作为目标块c_t，其他具有足够点数的块及其镜面对称的块作为候选块c_c，将从中寻找与目标块最相似的块。

(3)计算目标块和所有候选块的法向量的直流分量(Direct Component,DC)。

(4)在目标块和所有候选块上应用K-NN算法对块中的点建立图。

(5)计算目标块和所有候选块的法向量的各向异性图全变分(Anisotropic GraphTotal Variation,AGTV)。

(6)计算所有候选块与目标块的DC差距和AGTV差距，得到所有候选块与目标块间的相似度。

(7)选取相似度最高的候选块作为源块c_s。

4.相对位置的匹配：

(1)计算c_s中的几何结构与c_t中几何结构的在块中的相对距离，对c_s中的点进行块内平移得到

(2)使用简化的迭代最近点算法计算

和c_t中几何结构间的旋转矩阵，对

中的点进行块内旋转，将旋转后不在块内的点去除，并对

中的点体素化得到最终的源块

用于最后的修复步骤。

5.修复缺失区域：

(1)利用最终的源块c′_s中的信息修复目标块中的缺失区域，将这个修复工作转变为一个正则项为图平滑算子的优化问题，列出优化方程。

(2)通过对优化方程求最优解得到最终的修复结果块c_r，用结果块在点云中替换目标块，得到最终的修复后的点云。

本发明方法在点云的体素化中，通过体素划分将无序杂乱的点云的每个点的坐标及对应属性规则化。在点云孔洞的检测中，通过主成分分析选择一个主要投影方向，投影到二维深度图上进行孔洞的检测。在搜索相似区域步骤中，将点云切割成固定大小的立方块，并选择出包含缺失区域的立方块作为目标块。通过用直流分量和各向异性图全变分定义两个立方块几何结构的相似度，在所有立方块中搜索出与目标块的几何结构最相似的块作为源块。然后使用了简化的迭代最近点算法实现源块与目标块内相对位置的匹配。最后应用谱图理论和对局部频率的表示，对利用源块信息修复目标块中缺失区域的优化方程进行优化，得到修复并有一定去噪效果的结果。

本发明方法取得的技术效果是：相对于已有的技术，本方法由于在修复缺失区域的优化问题中，专门设立了一个项来保持结果块中的已知区域的部分与目标块中一致，能保持原本的几何特征。由于本方法不止依靠缺失区域的邻域信息进行修复，而是搜索了自相似区域并使用其中的内容作为填充内容，以合适的内容较好地修复了较为复杂的几何结构，而不使结果趋于平面或不符合事实。由于本方法在修复缺失区域的优化问题中，专门设立了一个项为图平滑算子，使已知内容与填充内容相融合，并有一定的去噪效果，从而不产生伪影。

附图说明

图1为针对真实存在的点云缺失区域的点云修复效果示意图，其中，(a)为原始点云，(b)～(f)分别为Meshfix方法、Meshlab方法、Wang等人的方法、Lozes等人的方法和本发明的方法进行修复后的点云。

图2为针对人工合成的点云缺失区域的点云修复效果示意图，其中，(a)为人工合成缺失区域的点云，(b)～(e)分别为Meshlab方法、Wang等人的方法、Lozes等人的方法和本发明的方法进行修复后的点云，(f)为未人工合成缺失区域前的点云。

图3为一个人工合成的点云缺失区域的点云修复结果中的去噪效果示意图，其中(a)为人工合成缺失区域前的点云，(b)为使用本发明方法获得的结果。

具体实施方式

为使本发明的上述特征和优点能更明显易懂，下文特举实施例，并配合所附图作详细说明如下。

针对本发明方法，首先介绍一下谱图理论中的基本概念，然后提出了基于图节点域的特征向量频率解释，谱图理论将是发明中的方法的核心。谱图理论包括以下几个内容：

(1)图和图拉普拉斯算子：

定义一个无向图，G＝{V,E,W}。V为图上顶点的集合，|V|＝N；E为边的集合；W为带权邻接矩阵。其中W是一个N×N的实对称矩阵，w_i,j是连接顶点i与顶点j的边的权重，通常使用非负值作为权重。

图拉普拉斯算子通常被定义为L＝D-W，其中D为对角矩阵，

(2)图信号平滑算子及其频率表示：

图信号是指驻留在图的顶点上的数据，例如社交网络、交通网络、传感器网络和神经元网络。例如，如果在点云上构造一个K近邻图(K-NN图)，那么每个点的坐标可以被看作是在K-NN图上定义的图信号。

对于在图G上定义的信号z，如果满足下述条件，那么称信号z相对于G的拓扑结构是平滑的：

其中，∈是一个参数，i～j表示图中相连的两个顶点。为了满足上述公式，z_i和z_j相对于权重较大的边时必须相似，对于权重较小的边时z_i和z_j可能完全不同。因此，上述公式强制z适应G的拓扑结构，从而创造了图信号的平滑性。由于z^TLz＝∑_i～jw_i,j(z_i-z_j)²，上述公式还可以写为z^TLz<∈。

除了平滑性质之外，这个算子还有去噪功能。由特征向量的定义有：

LX＝XΛ,

其中，X＝[x₁,x₂,…,x_Q]是L的Q个特征向量{x₁,x₂,…,x_Q}组成的矩阵，Λ＝[λ₁,λ₂,…,λ_Q]是L的特征值组成的矩阵且有0≤λ₁≤λ₂≤…≤λ_Q。因为L是一个实对称矩阵，所以有X^-1＝X^T，所以有：

L＝XΛX^-1＝XΛX^T.

然后使用特征向量来定义图傅里叶变换(GFT)。对于图G上的任意信号

它的GFT

为η＝X^Tz，其逆变换为z＝Xη＝∑_iη_ix_i。所以有

则这个算子还可以写作

为了满足这个公式，较大的λ_i对应的η_i必须较小，则其对应的在逆变换中的η_ix_i项会被抑制得很小。而特征向量x_i的频率f(x_i)与其对应的特征值λ_i是正相关的(在下一小节中证明)，所以被抑制的项是z中频率较高的部分。

这个算子将在后续的点云修复步骤中作为一个正则化项使用。

(3)特征向量在图节点域的频率表示：

从节点域的角度证明了图拉普拉斯矩阵的较大特征值通常对应于具有较高频率的特征向量，并有以下定理：

对于图拉普拉斯矩阵，当频率定义为图节点域中的过零点的数量时，特征向量的频率之间与其对应的特征值是正相关的。

对于图G上的拉普拉斯矩阵的特征向量x_i＝{x_i,1,x_i,2,…}而言，其正(负)强节点域是图G的一个最大连通子图且满足其中的所有点k都有x_i,k>0(x_i,k<0)；而其正(负)弱节点域是图G的一个最大连通子图且满足其中的所有点k都有x_i,k≥0(x_i,k≤0)，且至少有一个点x_i,k≠0。由于信号中的过零点次数可以表现其频率，所以用节点域中的过零点次数表示该节点域的频率。

Briandavies等人和Xu等人已证明了节点域中过零点次数ξ(x_i)的上限和下限值：

上限：强节点域中ξ(x_i)≤i+u-1(u为特征值λ_i的重数)，弱节点域中ξ(x_i)≤i。

下限：强节点域中ξ(x_i)≤i+u-1-q-z(q为当把图G转换为森林时需要移除的边的数量，z为特征向量x_i中零值的数量)。

所以ξ(x_i)的上限和下限都与i成正比，则通常有ξ(x_i)∝i，而定义时已有λ_i∝i，f(x_i)∝ξ(x_i)，故有f(x_i)∝λ_i。

在介绍完基本的谱图理论之后，现通过下述实施例提出的一种基于局域平滑性和非局域相似性的三维点云修复方法，具体阐述本发明方法，采用如下五个主要步骤：

1.点云的体素化：

(1)对于待处理的点云P^o，计算其每个点与K个近邻点的距离的平均值r。

(2)对点云P^o坐标进行放缩，使相邻两点间距离接近1，得到P^c，放缩点云中的点的公式为：

(3)计算平移向量s＝(s_x,s_y,s_z)，其中

是P^c中的x坐标。s_y,s_z的计算方式与s_x相同。将P^c中的点用s平移到原点对齐，得到P^t。

(4)分割出边长为1的体素后，每个体素中含有一组点{q₁,q₂,…}，每个点对应其附加属性值{g₁,g₂,…}。取每个含有点的体素的中心点组成新的点云P^v，并计算每个新点p的附加属性值如下：

其中μ_i是p所在的体素中每个点q_i对应的权重：

其中σ是权重参数。

2.点云孔洞的检测：

(1)通过对点云P^v中每个点的附加属性值中的法向量使用主成分分析，从沿着{x,y,z}三个坐标轴的六个方向中选择出投影方向r_p。主成分分析是统计学上常用的方法，是指在数据多又杂时，对数据进行降维，从中提炼出一个或多个主要数据内容(参见H.Abdiand L.J.Williams,“Principal component analysis,”Wiley InterdisciplinaryReviews Computational Statistics,vol.2,no.4,pp.433–459,2010.)。本发明利用该方法，在一个点云的几十万个点云中，提炼出一个属于{x,y,z}三个坐标轴的六个方向中的方向，作为投影方向。

(2)在深度图上将有点云中的点对应的点标记为-1，而没有点对应的点即为孔洞中的点标记为0。在深度图中执行广度优先搜索，首先将遍历出的第一个标记为0的点修改标记为1，与它相连的点标记为0的点也修改标记为1，依次类推直到与标记为1的点的邻居中没有标记为0的点，则第一个洞寻找完毕。然后继续广度优先搜索到下一个标记为0的点，修改标记为2，则第二个洞的寻找方式与第一个洞相同。以此类推，第三个洞中的点的标记均为3，……，直到修改标记完所有原本标记为0的点，得到所有不同的孔洞。

(3)将由于点云中点的密度过低形成的孔洞过滤掉，即标记点少于4的孔洞(因为经过体素化中的放缩后，相邻两点之间的距离通常小于2)。

(4)将属于点云所表示的物体之外的空点所对应的孔洞去除，即将包含深度图边界处的任何点的孔洞过滤掉，并计算该孔洞在投影方向r_p上的坐标范围，得到缺失区域Ω，即确实数据的孔洞。

3.搜索相似区域：

(1)将输入的点云P^v分割成互相有重叠部分的立方块{c₁,c₂,…,c_n}，选择包含缺失区域的立方块作为目标块c_t，将块中点数少于目标块c_t中点数80％的块去除，并将剩余的块关于x,y平面进行镜面对称得到新的一部分块，将该剩余的块和镜像得到的块作为候选块c_c。

(2)计算目标块的法向量的DC值d(c_t)和所有候选块的法向量的DC值d(c_c)。一个立方块c_i包含一组点{c_i,1,c_i,2,…,c_i,m}(m是块内点的数量)，每一个点都对应一个法向量表示为{n_i,1,n_i,2,…,n_i,m}。计算一个立方块c_i的DC值d(c_i)如下：

(3)在目标块和所有候选块上应用K-NN算法对块中的点建立图。以点的坐标为图信号，基于两点间的距离远近建立谱图理论中提到的K-NN图，并将其建立为无权图，即

(4)计算目标块的法向量的AGTV值v(c_t)和所有候选块的法向量的AGTV值v(c_c)。计算一个立方块c_i的AGTV值v(c_i)如下：

(5)计算所有候选块与目标块的DC差距δ_D(c_t,c_c)和AGTV差距δ_V(c_t,c_c)如下：

δ_D(c_t,c_c)＝|<d(c_t),d(c_c)>|，

δ_V(c_t,c_c)＝|v(c_t)-v(c_c)|，

由此计算出所有候选块与目标块间的相似度如下：

δ(c_t,c_c)＝exp{-[δ_D(c_t,c_c)+δ_V(c_t,c_c)]}。

(6)选取相似度最高的候选块作为源块c_s。

4.相对位置的匹配：

(1)计算c_s中的几何结构与c_t中几何结构的在块中的相对距离t＝(t_x,t_y,t_z)如下：

其中

表示块内缺失区域的边界。t_y和t_z同理。对c_s中的点用t进行块内平移得到

(2)随机选取c_t中的三个点{v₁,v₂,v₃}，并在在

中选取与它们距离最近的三个点{u₁,u₂,u₃}。计算三个点对之间的协方差矩阵：

其中G_i是一个4×4的矩阵：

(3)计算得到C的特征向量e＝[e₀,e₁,e₂,e₃]。

(4)计算

和c_t间的旋转矩阵：

(5)用R旋转

得到最终的源块

用于最后的修复步骤。

5.修复缺失区域：

(1)利用最终的源块

中的信息修复目标块中的缺失区域，将这个修复工作转变为一个正则项为图平滑算子的优化问题，列出优化方程如下：

其中c_r是想要得到的结果块。Ω表示c_t和

内的缺失区域，

表示块内的未缺失数据的已知区域。α和β是两个权重参数。

是

上的图的拉普拉斯矩阵，其中图是与搜索步骤中的(2)中的图相同的K-NN图。

(2)对优化方程进行c_r求导后得到最优解如下：

由此得到最终的修复结果块c_r，用结果块在点云中替换目标块，得到最终的修复后的点云。

图1、图2、图3显示了采用本发明方法和现有技术方案的修复效果图，可知本发明方法取得了更优异的修复效果。

以上实施例仅用以说明本发明的技术方案而非对其进行限制，本领域的普通技术人员可以对本发明的技术方案进行修改或者等同替换，而不脱离本发明的精神和范围，本发明的保护范围应以权利要求书所述为准。

Claims (10)

1.一种基于局域平滑性和非局域相似性的三维点云修复方法，包括以下步骤：

将原始的点云P^o坐标进行放缩和平移，分割成多个体素，计算各体素的坐标值和对应附加属性值，得到点云P^v；

将点云P^v投影为二维深度图，在深度图上标注出不同的孔洞并进行过滤，得到缺失数据的孔洞并确认为缺失区域；

将点云P^v分割成互相有重叠部分的立方块，选择包含缺失区域的立方块作为目标块c_t，并找出目标块c_t的所有候选块c_c；

计算目标块c_t和所有候选块c_c的法向量的DC值和AGTV值，根据两者的DC差距和AGTV差距，得到两者之间的相似度；

选取相似度最高的候选块作为源块c_s，计算源块c_s中的几何结构与目标块c_t中几何结构在块中的相对距离，对源块c_s中的点进行块内平移得到

计算

和目标块c_t中几何结构间的旋转矩阵R，对

中的点关于该旋转矩阵R进行块内旋转，将旋转后不在块内的点去除并对

中的点体素化，得到源块

利用源块

求解基于正则项为图平滑算子的优化方程的最优解，得到结果块c_r，将结果块c_r在点云P^o中替换目标块c_t以修复缺失区域，得到修复后的点云。

2.如权利要求1所述的方法，其特征在于，所述放缩是指先计算点云P^o中每个点与K个近邻点的距离的平均值r，再将点云P^o中的坐标按照公式

进行放缩，使相邻两点间距离接近1；

所述平移是指P^c按照平移向量s平移到与原点对齐，该平移向量s＝(s_x,s_y,s_z)，其中

是P^c中的x、y、z坐标。

3.如权利要求1所述的方法，其特征在于，将点云分割成边长为1的体素，每个体素中含有一组点{q₁,q₂,…}，每个点对应其附加属性值{g₁,g₂,…}，取每个含有点的体素的中心点组成新的点云P^v，P^v中每个新点p的附加属性值h计算公式为：

其中μ_i是p所在的体素中每个点q_i对应的权重：

其中σ是权重参数。

4.如权利要求1所述的方法，其特征在于，对点云P^v中每个点的附加属性值中的法向量进行主成分分析，从沿着{x,y,z}三个坐标轴的六个方向中选择出投影方向，点云P^v通过投影方向的向量r_p投影为二维深度图。

5.如权利要求1所述的方法，其特征在于，在深度图上进行广度优先搜索以标注出不同的孔洞，广度优先搜索方法为：在深度图上将点云中的有点对应的点标记为-1，而没有点对应的点即为孔洞中的点标记为0；首先将遍历出的第一个标记为0的点修改标记为1，与它相连的点标记为0的点也修改标记为1，依此类推，直到与标记为1的点的邻居中没有标记为0的点，则第一个孔洞寻找完毕；然后搜索到下一个标记为0的点，修改标记为2，则第二个孔洞的寻找方式与第一个孔洞相同，以此类推，第三个孔洞中的点的标记均为3，……，直到修改标记完所有原本标记为0的点，得到所有不同的孔洞；

对标注出的不同孔洞进行过滤，是指将标记点少于4的孔洞过滤掉，并将包含深度图边界处的任何点的孔洞过滤掉。

6.如权利要求1所述的方法，其特征在于，将立方块中点数少于目标块c_t中点数80％的块去除，并将剩余的块关于x,y平面进行镜面对称得到新的块，将该剩余的块和镜像得到的块作为候选块c_c。

7.如权利要求1所述的方法，其特征在于，假设一个立方块c_i包含一组点{c_i,1,c_i,2,…,c_i,m}，m是块内点的数量，每一个点都对应一个法向量表示为{n_i,1,n_i,2,…,n_i,m}，则计算立方块c_i的DC值d(c_i)的公式如下：

在目标块c_t和所有候选块c_c上应用K-NN算法对块中的点建立无权图w_k,l，即

则计算立方块c_i的AGTV值v(c_i)公式如下：

计算所有候选块c_c与目标块c_t的DC差距δ_D(c_t,c_c)和AGTV差距δ_V(c_t,c_c)如下：

δ_D(c_t,c_c)＝|<d(c_t),d(c_c)>|，

δ_V(c_t,c_c)＝|v(c_t)-v(c_c)|，

则所有候选块c_c与目标块c_t之间的相似度为：

δ(c_t,c_c)＝exp{-[δ_D(c_t,c_c)+δ_V(c_t,c_s)]}。

8.如权利要求1所述的方法，其特征在于，计算源块c_s中的几何结构与目标块c_t中几何结构在块中的相对距离t＝(t_x,t_y,t_z)如下：

其中

表示块内缺失区域的边界，t_y和t_z同理。

9.如权利要求1所述的方法，其特征在于，计算旋转矩阵R的步骤如下：

随机选取目标块c_t中的三个点{v₁,v₂,v₃}，并在

中选取与它们距离最近的三个点{u₁,u₂,u₃}，计算三个点对之间的协方差矩阵：

其中G_i是一个4×4的矩阵：

计算得到C的特征向量e＝[e₀,e₁,e₂,e₃]；

计算

和目标块c_t间的旋转矩阵R：

10.如权利要求1所述的方法，其特征在于，优化方程的最优解求解如下：

优化方程为：

其中，Ω表示目标块c_t和源块

内的缺失区域，

表示块内的已知区域，α和β是两个权重参数，

是源块

上的图的拉普拉斯矩阵，其中图是K-NN图；

对上述优化方程进行c_r求导后得到最优解如下：

Images