-
随着激光扫描技术的高速发展，三维激光扫描技术作为测绘领域继GPS技术之后的新一轮技术革命，广泛应用于测绘行业的各个领域[1-2]。利用三维激光扫描仪获取的点云数据通常是散乱无序的，由于设备精度、操作者经验、视线遮挡、偶然误差等原因，获取的点云数据往往会缺失部分采样点数据信息，进而导致出现孔洞的情况[3-5]，这将严重影响后续点云曲面重建的精度与效率，因此如何高效地对点云孔洞进行边界检测以及修补也逐渐成为点云数据处理的重要课题之一。
针对散乱点云的空洞修补问题，算法主要可以分为以下三种：基于体素的修补方法[6]、基于三角剖分的方法[7]以及构建隐式曲面的修补方法[8]，但这些方法在修补复杂曲面的点云数据时，修补结果往往不够理想。近年来许多学者进行了大量研究，如文献[9]利用最小二乘支持向量机建立曲面以达到点云孔洞修补，此方法能有效修补曲面平滑的点云，但由于缺少特征信息而不能修补曲率变化较大的点云；文献[10]提出一种全局拟合算法，利用泊松方程拟合结构复杂的孔洞三维模型的曲面孔洞，再利用三角剖分与孔洞进行拼接缝合，能有效对模型进行全局修复，但对于面积较大且具有尖锐特征的孔洞还原性较差。文献[11]提出了利用三次B样条曲线进行孔洞修补，该方法对大面积孔洞具有较好的修补效果，但修补局部孔洞时易受噪点影响导致修补精度得不到保障；文献[12]基于运动结构对点云孔洞进行修补，该方法对纹理复杂的点云孔洞具有良好修补效果，反之修补效果较差；文献[13]提出原始点云与孔洞特征匹配的孔洞修补方法，利用点云的特征信息提高匹配程度以提高孔洞修补精度，但该方法计算量大，导致效率较低。对于点云曲率变化较小、密度较为均匀或者不具有复杂空间特征的曲面来说，上述方法修补效果都较为良好，但对于包含多种曲面的复杂孔洞而言，上述方法的修补效果往往不够理想。
近年来，有学者将切片的概念融入到点云数据缺失处理中，如文献[14]提出在保留点云特征信息的前提下，进行自适应切片的孔洞修补算法。文献[15]在单向切片的基础上，提出双向切片概念，综合多方向的拟合结果对点云孔洞进行修补。为获取吻合程度高且表面光滑的修补点云常引入最小二乘法[16-18]，如文献[19]采用传统的最小二乘法进行平面拟合，该算法虽然易于实现，但受噪声点以及离散群点的影响较大，而且点云数据存在异方差性，所以使用传统最小二乘法计算曲率会存在一定的误差。文献[20]提出一种基于加权最小二乘法曲率计算的点云精简算法，该算法在相对平坦的区域精简得比较好，但是对于不平坦的区域却是保留了更多的特征点，因此该方法不适用于曲率较大的点云。
文中提出一种附加增值条件的移动最小二乘法的孔洞修补方法，尤其对于密度较大且空间特征明显的点云孔洞具有更优的修补效果。首先利用最小二乘微切面[21]检测孔洞边界，其次利用迭代切片法平均化点云密度，最后通过附加增值条件移动最小二乘法达到恢复点云孔洞特征以及修补点云孔洞的目的。文中的修补误差相较于三角格网法、定宽切片法、文献[7]方法以及文献[10]方法都有不同程度的降低，且修补后点云孔洞的曲率分布以及变化与原始点云模型契合度较高。
-
文中算法总体流程如图1所示，首先对点云数据模型构建最小二乘微切面检测封闭的孔洞边界并删除影响修补准确性的外边界点；其次利用迭代切片法从散乱点云中获取二维轮廓截面信息，并达到平均化点云密度的效果，可确保修补点云过渡自然；最后利用附加增值条件的移动最小二乘法进行增值填补点云，并通过曲率约束制约增值点云，可更大程度地恢复模型的空间特征，以达到契合原模型曲率变化的迭代增值，最终完成模型孔洞的修补。
图 1 总流程图
Figure 1. General flow chart
-
通过采样点及其K临近点集构建最小二乘微切面，完成孔洞边界检测与提取。
首先构建KdTree建立点云数据集间的拓扑关系，再利用最小二乘法构建对应的微切平面，该微切平面的法向量即为该点的法向量，并对点云法向量的方向进行一致化检测，保证其指向一致。然后将样本点的最近邻点集投影到对应的微切平面上，计算向量与基准向量间的夹角并排序，保证提取封闭的边界。外边界点对后续进行孔洞修补存在极大的影响，利用单坐标值搜索法对其进行删除。
-
（1）构建最小包围盒。最小包围盒是指三维坐标系中将所有散乱点云数据都包含在内的最小长方体，该长方体的长宽高分别与笛卡尔坐标系各坐标轴平行。沿着各坐标轴方向将坐标值快速排序，确定各方向上的最大最小值，得到最小包围盒的长宽高
${d_x}$ 、${d_y}$ 、${d_z}$ ，其中${d_x} = {x_{\max }}- {x_{\min }}$ ，$ {d_y} = {y_{\max }}- {y_{\min }}$ ，$ {d_z} = {z_{\max }} - {z_{\min }}$ 。（2）均匀分割。对最小包围盒进行初步栅格化，使用长宽高均为
$ L $ 的立方体对最小包围盒$Box$ 进行初步分割，得到一系列大小相同的小立方体，公式为：$$ L = \sqrt[3]{{a \cdot \varepsilon \cdot {d_x} \cdot {{d_y}} \cdot {d_z}/b}} $$ (1) 式中：
$ a $ 是一个标量；$\varepsilon $ 是邻域点个数；$b$ 是总点云数。通过初步栅格化得到的一系列小立方体是建立点云切片的基础，为降低计算量并保证小立方体个数充足，取$a = 3$ 。（3）迭代分割。在进行初步栅格化时未考虑到散乱点云可能存在密度不均、局部区域曲率变化大等问题，需对密度过大的立方体进行迭代分割。步骤如下：
1）确定切片方向。比较
$ \left| {{X_{\text{max}}} - {X_{\min }}} \right| $ 、$ \left| {{Y_{\text{max}} }- {Y_{\min }}} \right| $ 、$ \left| {{Z_{\text{max}} }- {Z_{\min }}} \right| $ 三者大小，选取最大值作为切片方向，再将切片投影后逐一拟合。2）确定切片宽度。对初次分割得到的立方体进行密度分析，若密度大于阈值，则需要进行循环迭代切割，该阈值取小立方中最大密度的0.2，将密度过大的立方体视为新的包围盒
$Box'$ ，密度分析结果决定分割块数量，依次迭代分割直至所有切片的密度相近即完成分割。密度分析公式为：$$ \rho = \frac{{\displaystyle\sum\limits_{i = 1}^{{b}} {\displaystyle\sum\limits_{j = 1}^\varepsilon {d_i^j} } }}{{\varepsilon \times b}} $$ (2) 切片宽度的计算公式为：
$$ \delta = \beta \times \rho $$ (3) 式中：
$\delta $ 为切片宽度；$\rho $ 为点云密度；$\beta $ 为常数，为避免切片宽度过小导致误差堆积，以及切片过大导致密度增大而加重计算量，当$\beta $ 取0.4~0.8时效果较优。 -
移动最小二乘法（moving least square, MLS）相较于其他拟合方法具有完备性和连续性，且精度较高。该算法是在点云原始数据点的基础上通过对周围数据点进行高阶多项式插值来重建表面缺失部分，这不仅可以保持点云原始样本不变，而且可以使填补的孔洞区域与原始表面较为契合。
-
设
$ U $ 为局部拟合域的一个子域，并构建拟合函数$f(x)$ 为：$$ {{f}}(x) = \sum\limits_{i = 1}^m {{\alpha _i}(x){p_i}(x) = {p^{\rm T}}(x)} \alpha (x) $$ (4) 式中：
$\alpha (x) = {[{a_1}(x),{a_2}(x), \cdots ,{a_m}(x)]^{\rm T}}$ 为待求系数；$p(x) = {[{p_1}(x),{p_2}(x), \ldots ,{p_m}(x)]^{\rm T}}$ 为基函数；$m$ 为为函数的项数。残差的离散加权
${L_2}$ 范式为：$$ \begin{split} J = & {\sum\limits_{j = 1}^N {w(x - {x_j})[f(x) - {y_j}]} ^2} =\\ & \sum\limits_{j = 1}^N {w(x - {x_j}){{[{p^{\rm T}}({x_j})\alpha (x) - {y_j}]}^2}} \end{split}$$ (5) $$ J = {[P\alpha (x) - Y]^{\rm T}}W(x)[P\alpha (x) - Y] $$ (6) 当局部近似值
$f({x_j})$ 和节点值${y_j}$ 之差平方带权越小，局部近似值精确的度越高，式中，$N$ 是求解区域内节点的数目，$c$ 为拟合函数，$w(x - {x_j})$ 为节点${x_j}$ 的权函数。权函数只与节点和拟合点的距离有关，权函数
$w(x - {x_j})$ 非负而且随着距离的增加而单调递减。常用的权函数为三次样条权函数，其表达式为：$$ w(\bar s) = \left\{ {\begin{array}{*{20}{c}} {\dfrac{2}{3} - 4{s^{ - 2}} + 4{s^{ - 3}}}&{(\bar s) \leqslant \dfrac{1}{2}}\\ {\dfrac{4}{3} - 4{s^{ - 1}} + 4{s^{ - 2}} - \dfrac{4}{3}{s^{ - 3}}}&{(\dfrac{1}{2} < \bar s \leqslant 1)}\\ {\rm{0}}&{(\bar s > 1)} \end{array}} \right. $$ (7) 式中：
$\bar s $ 为距离相对量；$s = x - {x_j}$ 为点到节点的距离。系数${\alpha _i}(x)$ 的选取是使拟合函数$f(x)$ 在计算点$x$ 的邻域内求得最佳近似值。为确定系数$\alpha (x)$ ，令$J$ 取最小值，对于任意函数$h(x)$ 和$g(x)$ ，有：$$ (h,g) = \sum\limits_{i = 1}^n {w(x - {x_i})h({x_i})} g({x_i}) $$ (8) $$ \begin{split} &{\alpha _1}(x)({p_i},{p_1}) + {\alpha _2}(x)({p_i},{p_2}) + \cdots + {\alpha _m}(x)({p_i},{p_m}) =\\ &\qquad\qquad ({p_i},{y_i}),i = 1,2, \cdots m \end{split} $$ (9) 以矩阵形式求解方程组得：
$$ \alpha = {\left( {{p^{\rm T}}wp} \right)^{ - 1}}{p^{\rm T}}wy $$ (10) 代入公式(9)可得MLS的拟合函数：
$$ f(x) = \sum\limits_{i = 1}^n {{\phi _i}(x)} {y_i} $$ (11) 其中形函数为：
$$ {\phi _i}(x) = {\sum\limits_{j = 1}^m {{p_j}(x)\left[ {{{({p^{\rm T}}wp)}^{ - 1}}{p^{\rm T}}w} \right]} _{ji}} $$ (12) -
由于引入近支概念中拟合点受邻域节点的影响，则拟合曲线经过两个间距较大的节点时，如果直接使用传统的移动最小二乘法进行拟合，则易出现偏差，所以文中提出了一种附加增值的移动最小二乘法来减小该偏差对孔洞修补效果的不良影响，具体步骤如下：
（1）在获得密度均匀的点云切片后，对切片进行逐一的投影，获取对应的二维散点图，依据传统的移动最小二乘法原理构建拟合函数
$f(x)$ 并求出初始节点点集$H$ ，拟合得到初始曲面，用曲面外的点云定义离散群点，并计算离散群点对应的法线；（2）对初始节点点集
$H$ 进行重采样，得到样本点并定义搜索邻域半径${r_v}$ ，由迭代切片获得的分割块密度均匀，则以样本点所在的分割块的宽度$ \delta $ 为基准，取${r_v} = {\delta \mathord{\left/ {\vphantom {\delta 3}} \right. } 3}$ ，以样本点为圆心的平面邻域内画圆，再将圆内的邻域离群点投影至初始拟合曲面上，使得拟合曲面包含充足的点云（见图2），避免进行孔洞修补时由于节点数较少或距离较大导致拟合出现较大偏差；
图 2 节点邻域离散群点投影示意图
Figure 2. Projection diagram of the discrete points in the node neighbourhood
（3）由于节点
$H$ 为拟合曲线必经之点，为避免节点间出现间距较大或过于稀疏的情况，将初始拟合曲面上所有点再次进行拟合，得到的节点点集$A$ ，不仅保证有足够的边界邻域节点为基础进行孔洞修补，也能更有效地拟合出更为光滑的曲面；（4）对点云孔洞而言，当孔洞较大或空间特征明显时，只基于孔洞边界邻域的点进行修补时，可类比于节点间距离过大的情况，该情况下通常无法恢复原始模型的特征，对此，文中提出一种以增值条件新增点云的方法，先对边界点构建最小包围盒并进行栅格化处理获取边长最大值
${l_{\max }}$ ，再搜索边界邻域的节点${Q_k}\left( {{x_k},{y_k}} \right)$ ，${Q_k} \in A$ ，并定义搜索邻域半径为${{{r}}_k} = {{{l_{\max }}} \mathord{\left/ {\vphantom {{{l_{\max }}} 6}} \right. } 6}$ ，以节点${Q_k}$ 为圆心画圆并统计圆内包含的点云数$q$ ，$q = 1,2, \cdots ,n$ ，点云孔洞邻域处的节点通常位于孔洞边界处，增值的点云数应在契合原模型密度的同时，也要避免新孔洞的出现，所以定义增值半径${{3{{{r}}_k}} \mathord{\left/ {\vphantom {{3{{t{r}}_k}} 2}} \right. } 2}{\text{ \gt }} {{{r}}_s} \gt {{\text{r}}_k}$ ，并使得新增点云符合增值条件${G_s}({x_s},{y_s})$ ，$s = 1,2, \cdots ,t$ ，$t \leqslant n$ ，可达到在某一方向上新增点云，则附加增值的移动最小二乘法拟合的曲线表示为：$$ y = f(x) - \sum\limits_{s = 1}^t {{h_s}(x)} {\gamma _s} $$ (13) 式中：
${h_s}(x) = \displaystyle\prod\limits_{j = 1,s \ne j}^t {\dfrac{{(x - {x_j})}}{{({x_s} - {x_j})}}}$ ，$s = 1,2, \cdots ,t$ ，${\gamma _s} = f({x_s}) - {y_s}$ 。（5）为确定点云的增值方向，计算边界点的质心坐标
$g\left( {{x_l},{y_l},{z_l}} \right)$ ，$g$ 点也为点云孔洞的质心，其中边界点质心坐标为：$$ g = \frac{1}{k}\left( {\sum\limits_{l = 0}^k {{x_l}} ,\sum\limits_{l = 0}^k {{y_l}} ,\sum\limits_{l = 0}^k {{z_l}} } \right) $$ (14) 将质心
$g\left( {{x_l},{y_l},{z_l}} \right)$ 投影至二次拟合曲面上获得质心二维坐标，以孔洞边界邻域内的最终节点${Q_k} \left( {{x_k},{y_k}} \right)$ 为初始点，以增值半径${{{r}}_s}$ 向$g$ 质心方向进行增值；（6）为确保增值点云更契合原始模型，则需对新增点云进行曲率约束，需涉及平均曲率与高斯曲率，构建的拟合曲面为
$z = \lambda (x,y)$ ，则拟合曲线的单位法向量为$\overrightarrow {n_z} $ ：$$ {\overrightarrow{n}}_{z}=\frac{\dfrac{\partial \lambda }{\partial x} \cdot \dfrac{\partial \lambda }{\partial y}}{\left|\dfrac{\partial \lambda }{\partial x} \cdot \dfrac{\partial \lambda }{\partial y}\right|} $$ (15) 曲面平均曲率
$W$ 与高斯曲率$V$ 表示为：$$ \left\{ \begin{gathered} W = \frac{{EN - 2FM + TU}}{{2\left( {ET - {F^2}} \right)}} \\ V = \frac{{UN - {M^2}}}{{ET - {F^2}}} \\ \end{gathered} \right. $$ (16) 式中：
$E = \dfrac{\partial \lambda }{\partial x} \cdot \dfrac{\partial \lambda }{\partial x}$ ；$F = \dfrac{\partial \lambda }{\partial x} \cdot \dfrac{\partial \lambda }{\partial y}$ ；$T = \dfrac{\partial \lambda }{\partial y} \cdot \dfrac{\partial \lambda }{\partial y}$ ；$U = \dfrac{\partial \lambda }{\partial x\partial x} \cdot \overrightarrow{{n}_{z}}$ ；$ M=\dfrac{\partial \lambda }{\partial x\partial y} \cdot \overrightarrow{{n}_{z}} $ ；$ N=\dfrac{\partial \lambda }{\partial y\partial y} \cdot \overrightarrow{{n}_{z}} $ 。计算整体点云模型的平均曲率
$W$ ，以及孔洞边界点的邻域高斯曲率${V_q}$ ，其中，孔洞边界点邻域以分割块长度为单位，距离依次增加三个分割块，从远到近依次对比$W$ 与${V_q}$ 的大小关系，构建函数$f(q) = W - {V_q}$ ，$q = 1,2, \cdots ,n$ ，当$f(q)$ 的变化值突增或突减，则结束邻域获取（达到$n$ 值），如若$f(q)$ 为变化率较大的增函数，则孔洞的空间特征应为凸区域；若$f(q)$ 为变化率较大的减函数，则孔洞的空间特征应为凹区域；若$f(q)$ 为变化率较小的函数，则孔洞的空间特征较为平缓。以$f(q)$ 的变化规律制约新增点云的曲率大小，即利用$f(q)$ 作为新增点云曲率${V_o}$ 的预判函数，新增点云的曲率值应处于曲率预判函数上或${V_{q - }} \leqslant 2{V_o} \leqslant {V_{q + }}$ 。为避免孔洞边界处曲率变化较大，即选择${V_{q - }}$ 与${V_{q - }}$ 是除去$q = 1,2,3$ 以外的孔洞边界邻域内曲率最小值与最大值。（7）返回步骤（3）将新增点云进行拟合后获取新的边界邻域节点，若增值圆内密度大于或等于邻域值，则说明已增值至质心处，孔洞修补结束，反之继续向
$g$ 质心方向进行迭代增值。 -
为验证文中修补方法的有效性，以石碑模型与三个斯坦福大学3D点云库模型作为实验数据，并在点云模型上人为制造不同类型的点云孔洞。
通过模型孔洞分别在五种不同算法修补后与原始模型之间的完整率、准确率、配准均方根误差以及平均曲率的均方根误差进行精度对比分析。完整率和准确率是通过求修补模型与原始模型的交集比（大于50%）构建对应空间关系，再利用公式(17)以点云模型的单个对应点为基元评价修补效果，公式为：
$$ \left\{ \begin{gathered} comp = \frac{{TP}}{{TP + FN}} \\ corr = \frac{{TP}}{{TP + FP}} \\ \end{gathered} \right. $$ (17) 式中：
$TP$ 表示正确修补的基元数，由修补点云在原始模型邻域内的点数决定；$FN$ 表示没有交集的基元数，代表与原始模型邻域无交集且平滑成面的点云数；$FP$ 表示错误的基元数，为远离拟合面的杂乱噪声点。均方根误差是利用修补模型与原始模型进行ICP配准，通过已知的原始模型点集
$P$ 与修补模型点集$Q$ 按照约束条件搜索最邻近点$\left( {{p_i},{q_i}} \right)$ ，计算出误差函数最小的最优匹配参数$R$ 与$t$ ，表达式为：$$ E(R,t) = \frac{1}{u}{\sum\limits_{i = 1}^u {\left\| {R{p_i} + t - {q_i}} \right\|} ^2} $$ (18) 式中：
$u$ 为最邻近点对应的个数；${p_i}$ 为$P$ 中的点；${q_i}$ 为$Q$ 中的点，由点沿坐标轴的旋转角度可计算得旋转矩阵$R$ ；平移矩阵$t$ 由点沿坐标轴的平移量计算得到。再对模型进行曲率分析与可视化，并对修补的孔洞点云进行高斯平均曲率计算，利用修补后模型点云的平均曲率
${H_i}$ 与原模型的平均曲率${F_i}$ 求差的平方和实验次数$e$ 比值的平方根，计算得修补后模型的平均曲率的均方根误差，公式为：$$ Rmse = \sqrt {\frac{1}{e}{{\sum\limits_{i = 1}^e {\left( {{H_i} - {F_i}} \right)^2} }}} $$ (19) -
文中方法是在Intel Core i9-12900 k 3.2 GHz的CPU、64 GB内存的PC机，利用PCL点云库与Visual studio2017环境下的C++语言实现。
实验1数据为2021年6月通过RIEGLVZ-1000型号高精度地面三维激光扫描仪获取的景观石碑点云模型。为了便于观察效果，在石碑表面规则截取包含多种曲面且密度不均的点云作为实验数据。选取石碑点云模型中密度分布不均匀且含有多种曲面的部分并制造规则孔洞，红色点集为提取的孔洞边界信息。
从图3中可知，TG生成的点云分布散乱，未能较好修补孔洞；FWS能以距离分布均匀的点云较好地修补孔洞，MLS能均匀地修补孔洞，文献[10]方法相较于前三种方法能更接近模型密度修补孔洞，但都未能修复模型的曲面特征信息；而文中方法能有效地修补点云孔洞并具有原模型的特征信息以及更契合的密度。
图 3 石碑模型孔洞修补结果
Figure 3. Results of the repair of the hole in the stela model
由表1可知文中方法对模型密度的适应性较强，修补后的点云数最接近于原始孔洞的点云数，由于经历了二次曲面拟合通过附加增值条件进行修补，修补完整率与准确率都高于其余方法，配准均方根误差较低，且计算效率较高，则文中方法更契合原始模型特征。
表 1 石碑点云模型孔洞修补结果对比
Table 1. Comparison of hole repair results of stela model point cloud
Model point
cloud countHole point
cloud numberBox parameters/m Method Repair point
cloud numberIntegrity rate Accuracy rate ICP-error/mm t/s 2890 216 dx=0.316
dy=1.278
dz=0.980TG 307 92.32% 88.46% 0.105 20.24 FWS 126 96.51% 91.07% 0.064 13.41 MLS 105 90.43% 86.62% 0.142 17.58 Literature [10] 169 98.01% 94.74% 0.025 9.83 Proposed 254 99.84% 97.38% 0.008 3.16 实验2选用马首模型，在进行曲率分析后选择在多种曲率变化的位置造不规则的人为孔洞，该孔洞包含橙、黄、绿和蓝四种颜色，则说明该部分点云拥有四种明显的曲率变化，使得马首模型表面缺失多曲率变化的不规则曲面域。
由图4可知，TG不能适应原始模型的密度，甚至修补后出现新的孔洞；FWS与MLS虽能得到密度分布均匀的修补点云，但修补后的点云过于稀疏也不能适应原始模型的密度；文献[10]的修补方法虽能适应原始模型的密度，但未能恢复模型的特征信息；而文中方法能较好地适应原始模型的密度，且与孔洞边界点云进行了更好的融合，以确保能具有原模型的特征信息。
图 4 马首模型孔洞修补结果
Figure 4. Hole repair results of horse head model
图5为马首模型的修补曲率分布对比图，模型曲率变化从大至小是由红至蓝表示，其中，TG修补结果存在曲率过高、不符合原始模型的噪声点以及出现新孔洞的问题；FWS、MLS以及文献[10]的修补结果过于稀疏且曲率过度不自然；而文中方法能有效地适应原始模型的曲率变化且密度分布均匀，产生的噪声点较少，能较好地体现出原始模型的形状。
由表2可知，在孔洞曲率增大时，修补完整率变化不大，但前三者准确率下降幅度较大且计算效率明显降低，由于文中方法在增值点云时具有曲率约束，使得文中方法的配准均方根误差与平均曲率均方根误差得到一定的控制。
图 5 马首模型曲率分布对比图
Figure 5. Comparison chart of curvature distribution of horse head model
表 2 马首点云模型孔洞修补结果对比
Table 2. Comparison of point cloud hole repair results of horse head model
Model point
cloud countHole point
cloud numberBox parameters/m Method Repair point
cloud numberIntegrity rate Accuracy rate Rmse/mm ICP-error/mm t/s 9695 302 dx=0.068
dy=0.084
dz=0.044TG 562 91.38% 86.14% 0.556 1.342 47.35 FWS 127 95.75% 85.72% 0.507 1.914 39.48 MLS 101 88.90% 83.48% 0.574 1.275 44.19 Literature [10] 176 97.04% 91.37% 0.172 0.851 34.54 Proposed 279 98.96% 94.39% 0.029 0.258 13.26 实验3选用Bunny模型，该模型体积最大，且不规则人为孔洞面积也远大于其余模型，使得Bunny点云模型表面缺失不规则且面积较大的曲面域。Bunny 模型孔洞修补结果如图6所示。
图 6 Bunny模型孔洞修补结果
Figure 6. Hole repair results of Bunny model
由图7可知，TG修补结果密度过大导致(图7(b))①处掩盖了一部分边界点，②与③处出现了新的孔洞；FWS的修补结果过于稀疏导致(图7(c))①处绿色点云面积较大但数量少，(图7(c))②处点云的稀疏程度较大；MLS的修补结果中(图7(d))①和(d)②处绿色点云分布均大于原始模型；文献[10]的修补结果在曲率变化最大的位置仍未能适应原始曲率；而文中方法虽出现了部分修补点云走向与原始点云走向不同，但该部分点云走向契合边界点，并能较好地填充孔洞且曲率相似度高，说明已大致恢复原始模型的空间特征。由表3可知，当孔洞面积增大时，文中方法由于通过迭代增值点云，能更好地适应模型的点云密度并在增值时更为迅速，使得修补完整率仍较高，配准均方根误差变大，但平均曲率均方根误差仍较小。
图 7 Bunny模型曲率分布图对比图
Figure 7. Comparison chart of curvature distribution of Bunny model
表 3 Bunny模型点云孔洞修补结果对比
Table 3. Comparison of point cloud hole repair results of Bunny model
Model point
cloud countHole point
cloud numberBox parameters/m Method Repair point
cloud numberIntegrity rate Accuracy rate Rmse/mm ICP-error/mm t/s 35947 620 dx=15.570
dy=15.433
dz=12.067TG 1188 92.12% 83.76% 3.785 0.303 85.59 FWS 253 90.68% 83.04% 3.506 0.099 76.60 MLS 225 89.76% 80.49% 3.668 0.282 80.31 Literature [10] 356 93.71% 89.72% 2.336 0.094 63.39 Proposed 539 99.05% 93.87% 1.051 0.003 28.49 实验4选用dragon模型，该模型密度与空间特征的复杂程度较大，选择弯曲程度最大的部分制造相对该模型面积较大的不规则人为孔洞，使得dragon孔洞点云具有曲率变化大、密度大且空间特征复杂的特点，则需恢复原模型特征以及空间结构的难度与复杂程度更高。dragon 模型孔洞修补过程以及原始模型与修补结果的对比如图8所示。图8(e)中红色点云为原始点云，黄色点云为文中修补结果。
图 8 dragon模型孔洞修补过程以及原始模型与修补结果的对比
Figure 8. Hole repair process of the dragon model and the comparison between the original model and the repair result
图 9 dragon模型孔洞修补结果
Figure 9. Hole repair results of the dragon model
dragon模型点云孔洞修补的实验结果(图9)与曲率分布对比结果(图10)是最显而易见的，TG、FWS、MLS和文献[10]不能适应原始模型的曲率与密度，未恢复空间结构与原模型特征，出现了多处粘连且无序的噪声点云，导致修补效果较差，而文中通过附加增值条件的移动最小二乘法对于弯曲程度高的点云孔洞依旧能较好地恢复该区域的空间结构与特征，使得修补后的点云不存在错误的粘连且平滑程度较好。
由表4可知，对于模型空间结构含有弯曲程度较大且密度较大的区域，五种对比方法的修补误差增大，说明无论是曲率与模型特征都未能得到良好恢复，但文中方法能较好的恢复孔洞原有的空间结构，且各误差值相较于其余方法较小，计算效率较高，则表明文中附加增值条件的移动最小二乘法能生成更为契合的修补点云。
图 10 dragon模型曲率分布对比图
Figure 10. Comparison of the curvature distribution of the dragon model
表 4 dragon模型点云孔洞修补结果对比
Table 4. Comparison of point cloud hole repair results of dragon model
Model point
cloud countHole point
cloud numberBox parameters/m Method Repair number Integrity rate Accuracy rate Rmse/mm ICP-error/mm t/s 47449 2158 dx=0.205
dy=0.144
dz=0.092TG 1003 83.68% 79.75% 2.857 4.186 163.40 FWS 395 81.14% 73.46% 3.394 5.800 146.17 MLS 368 72.62% 65.20% 5.166 6.572 98.46 Literature [10] 647 83.23% 80.09% 2.825 4.311 112.53 Proposed 1515 98.67% 91.95% 0.772 0.735 57.32 -
针对空间复杂度不同的点云模型孔洞问题，提出了一种利用附加增值条件的移动最小二乘法，通过迭代切片保留模型特征，再对离散点构建二次拟合曲面减小修补的累积误差，最后利用附加增值条件的最小二乘法实现契合原始模型的孔洞修补。利用四种不同密度、不同大小以及不同复杂程度的点云模型为数据源，对比四种已有的点云孔洞修补方法可知，文中方法在精度上有较大的优势，其中完整率提高1.83%~35.87%，准确率提高2.64%~41.03%，配准均方根误差与平均曲率均方根分别降低了68%、69.68%以上，耗费的计算时间均低于其余方法，可较好地平均化点云密度、契合原模型曲率变化，并保持了点云模型的局部信息特征，对点云模型具有较强适用性。该方法主要面向三维点云模型，如何运用到大场景是下一步的研究重点。
Additional value-added conditional moving least squares method for point cloud hole repair
-
摘要: 由于扫描设备局限或模型结构复杂等因素导致点云模型出现孔洞，这严重影响模型的后续处理。针对点云孔洞的修补问题，文中提出了一种附加增值条件移动最小二乘法的点云孔洞修补方法。首先提取封闭的孔洞边界，通过密度分析进行迭代切片，不仅削弱点云分布不均的影响，还提高模型细节特征的保留程度；再将离散群点投影至拟合曲面，投影点集二次拟合以获取拟合面节点，保证有足够的边界邻域节点为基础进行孔洞修补；最后利用附加增值条件移动最小二乘法对孔洞进行迭代修补，并对增值点云进行曲率约束，从而达到契合原始模型空间特征的重建。实验采用人为在四个点云模型上制造不同类型的孔洞，并与现有的四种方法进行对比，验证所提方法的有效性，结果表明，文中方法相较于现有的四种方法，完整率、准确率提高了1.83%以上，配准均方根误差与平均曲率均方根降低了68%以上，对比证明了文中方法对于点云模型孔洞具有较强的适用性，可为重建三维点云模型提供可靠信息。Abstract: Due to the limitations of scanning equipment or the complexity of the model structure, holes appear in the point cloud model, which seriously affects the subsequent processing of the model. To address the problem of point cloud hole repair, this paper proposes a point cloud hole repair method with an additional value-added conditional moving least squares method. Firstly, the closed hole boundary is extracted and iteratively sliced through density analysis, which not only reduces the impact of the uneven distribution of the point cloud, but also improves the retention standard of detailed features of the model; Besides, the discrete group of points is projected onto the fitted surface, and the projected point set is fitted twice to obtain the nodes of the fitted surface to ensure that there are enough boundary neighborhood nodes for hole repair; Finally, the holes are repaired by using the additional value-added conditional moving least squares method. Meanwhile, the curvature of the value-added point cloud is constrained, so as to achieve the reconstruction that fits the spatial characteristics of the original model. In the test, different types of holes are artificially made on four point cloud models, and the effectiveness of proposed method is verified by comparison with the existing four methods. The results show that, compared with the four existing methods, the completeness and accuracy of this method are improved by more than 1.83%, and the root mean square error of the alignment and the root mean square of the curvature are reduced by more than 68%, which proves the applicability of this method for point cloud model holes, which can provide reliable information for the reconstruction of 3D point cloud models.
-
图 2 节点邻域离散群点投影示意图
Figure 2. Projection diagram of the discrete points in the node neighbourhood
图 5 马首模型曲率分布对比图
Figure 5. Comparison chart of curvature distribution of horse head model
图 7 Bunny模型曲率分布图对比图
Figure 7. Comparison chart of curvature distribution of Bunny model
图 8 dragon模型孔洞修补过程以及原始模型与修补结果的对比
Figure 8. Hole repair process of the dragon model and the comparison between the original model and the repair result
图 10 dragon模型曲率分布对比图
Figure 10. Comparison of the curvature distribution of the dragon model
表 1 石碑点云模型孔洞修补结果对比
Table 1. Comparison of hole repair results of stela model point cloud
Model point
cloud countHole point
cloud numberBox parameters/m Method Repair point
cloud numberIntegrity rate Accuracy rate ICP-error/mm t/s 2890 216 dx=0.316
dy=1.278
dz=0.980TG 307 92.32% 88.46% 0.105 20.24 FWS 126 96.51% 91.07% 0.064 13.41 MLS 105 90.43% 86.62% 0.142 17.58 Literature [10] 169 98.01% 94.74% 0.025 9.83 Proposed 254 99.84% 97.38% 0.008 3.16 
下载: 导出CSV
表 2 马首点云模型孔洞修补结果对比
Table 2. Comparison of point cloud hole repair results of horse head model
Model point
cloud countHole point
cloud numberBox parameters/m Method Repair point
cloud numberIntegrity rate Accuracy rate Rmse/mm ICP-error/mm t/s 9695 302 dx=0.068
dy=0.084
dz=0.044TG 562 91.38% 86.14% 0.556 1.342 47.35 FWS 127 95.75% 85.72% 0.507 1.914 39.48 MLS 101 88.90% 83.48% 0.574 1.275 44.19 Literature [10] 176 97.04% 91.37% 0.172 0.851 34.54 Proposed 279 98.96% 94.39% 0.029 0.258 13.26
下载: 导出CSV
表 3 Bunny模型点云孔洞修补结果对比
Table 3. Comparison of point cloud hole repair results of Bunny model
Model point
cloud countHole point
cloud numberBox parameters/m Method Repair point
cloud numberIntegrity rate Accuracy rate Rmse/mm ICP-error/mm t/s 35947 620 dx=15.570
dy=15.433
dz=12.067TG 1188 92.12% 83.76% 3.785 0.303 85.59 FWS 253 90.68% 83.04% 3.506 0.099 76.60 MLS 225 89.76% 80.49% 3.668 0.282 80.31 Literature [10] 356 93.71% 89.72% 2.336 0.094 63.39 Proposed 539 99.05% 93.87% 1.051 0.003 28.49
下载: 导出CSV
表 4 dragon模型点云孔洞修补结果对比
Table 4. Comparison of point cloud hole repair results of dragon model
Model point
cloud countHole point
cloud numberBox parameters/m Method Repair number Integrity rate Accuracy rate Rmse/mm ICP-error/mm t/s 47449 2158 dx=0.205
dy=0.144
dz=0.092TG 1003 83.68% 79.75% 2.857 4.186 163.40 FWS 395 81.14% 73.46% 3.394 5.800 146.17 MLS 368 72.62% 65.20% 5.166 6.572 98.46 Literature [10] 647 83.23% 80.09% 2.825 4.311 112.53 Proposed 1515 98.67% 91.95% 0.772 0.735 57.32
下载: 导出CSV
-
[1] Shao X J, Pan S, Song S, et al. Research on point cloud splicing for inner surface inspection of deep hole [J]. Infrared and Laser Engineering, 2021, 50(12): 20210210. (in Chinese) doi: 10.3788/IRLA20210210 [2] Zhao F Q, Zhou M Q. Hierarchical point cloud denoising algorithm [J]. Optics and Precision Engineering, 2020, 28(7): 1618-1625. (in Chinese) doi: 10.37188/OPE.20202807.1618 [3] Chang B T, Chen C F, Guo J J, et al. Interpolation-based filtering with segmentation for airborne LiDAR point clouds [J]. Infrared and Laser Engineering, 2021, 50(9): 20200369. (in Chinese) doi: 10.3788/IRLA20200369 [4] Pang C, Zhong X, Hu H, et al. Adaptive obstacle detection for mobile robots in urban environments using downward-looking 2D LiDAR [J]. Sensors, 2018, 18(6): 1749-1767. doi: 10.3390/s18061749 [5] Wu Q H, Cai Q J S, Lai C A, et al. Registration of losing point cloud based on clustering extended Gaussian image [J]. Optics and Precision Engineering, 2021, 29(5): 1199-1206. (in Chinese) doi: 10.37188/OPE.20212905.1199 [6] Huang X, Xie C, Fang X, et al. Combining pixel and object-based machine learning for identification of water-body types from urban high-resolution remote-sensing imagery [J]. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 2015, 8(5): 2097-2110. doi: 10.1109/JSTARS.2015.2420713 [7] Han J, Jiang B C, Xia L, et al. A B-spline curve fitting algorithm based on contour key points [J]. Applied Mathematics and Mechanics, 2015, 36(4): 423-431. (in Chinese) doi: 10.3879/j.issn.1000-0887.2015.04.010 [8] Hu Y Y, Wang J J, Fan Y Y, et al. LiDAR-based three-dimensional modeling and volume calculation for space objects [J]. Chinese Journal of Lasers, 2020, 47(5): 0510001. (in Chinese) doi: 10.3788/CJL202047.0510001 [9] Yang Y Q, Li S H. Hole repairing algorithm for point cloud data based on least square support vector machine [J]. Journal of Jilin University (Science Edition), 2018, 56(3): 692-696. (in Chinese) [10] Li Y W, Geng G H, Wei X R. Hole-filling algorithm based on poisson equation [J]. Computer Engineering, 2017, 480(10): 209-221. (in Chinese) [11] He J D, Niu J Y, Zhang Z R, et al. Repairing method of missing area of dairy cows’point cloud based on improved cubic B-spline curve [J]. Journal of Agricultural Machinery, 2018, 49(6): 225-231. (in Chinese) [12] Gai S, Da F P, Zeng L L, et al. Research on a hole filling algorithm of a point cloud based on structure from motion [J]. Journal of the Optical Society of America A, 2019, 36(2): 39-46. doi: 10.1364/JOSAA.36.000A39 [13] Harary G, Tal A, Grinspun E. Context-based coherent surface completion [J]. ACM Transactions on Graphics, 2014, 33(1): 1-12. doi: https://doi.org/10.1145/2532548 [14] He G Z. Hole patching of adaptive slicing based on feature-data segmentation [J]. Journal of East China Jiaotong University, 2014, 31(4): 95-99. (in Chinese) [15] Wang Y G, Xu Ya J. Research about point cloud hole-repairing based on bidirectional slice method [J]. Geomatics & Spatial Information Technology, 2015, 38(10): 218-220. (in Chinese) [16] Salkauskas P L. Surfaces generated by moving least squares methods [J]. Mathematics of Computation, 1981, 37(155): 141-158. doi: 10.1090/S0025-5718-1981-0616367-1 [17] Hui Z Y, Li N, Hu H Y, et al. Multi-scale progressive digital terrain model construction method based on backpack LiDAR point cloud [J]. Chinese Journal of Lasers, 2022, 49(4): 0410001. (in Chinese) doi: 10.3788/CJL202249.0410001 [18] Xu Y T, Li G Y, Qiu C X, et al. Single photon laser data processing technology based on Terrain Correlation and least square curve fitting [J]. Infrared and Laser Engineering, 2019, 48(12): 1205004. (in Chinese) doi: 10.3788/IRLA201948.1205004 [19] Ma D, Li G, Gui J, et al. Research on redundant data reduction algorithm of laser scanning point cloud based on curvature and normal [J]. Boletin Tecnico/Technical Bulletin, 2017, 55(6): 272-278. [20] Tang Z Y, Gao B L, Dou M L. Point clouds simplification algorithm based on weighted least squares surface fittiong for curvature computation [J]. Computer Engineering and Design, 2019, 40(6): 1606-1610, 1659. (in Chinese) [21] Sun D Z, Fan Z X, Li Y R. Automatic extraction of boundary characteristic from scatter data [J]. Journal of Huazhong University of Science and Technology (Natural Science Edition), 2008, 297(8): 82-84. (in Chinese) -
点击查看大图
计量
- 文章访问数: 166
- HTML全文浏览量: 60
- PDF下载量: 43
- 被引次数: 0
