随着3维扫描技术的不断发展和进步, 人们可以用各种3维扫描设备获取不同形状的数字点云模型, 然后完成物体重建。该技术尤其在工业领域的逆向工程中被广泛应用。然而, 由于扫描物体本身具有复杂的几何形状, 扫描过程中物体遮挡, 物体本身残缺, 扫描过程的光线反射等原因, 会造成不同程度的扫描物体重建后出现孔洞, 孔洞的存在对后续模型处理比如逆向工业原件生产等造成不同程度的障碍。因此需要对孔洞完成修补, 这也是当前3维数据处理的一个难题。

对于孔洞修补, 符号距离场体积扩散的方法是最为典型的方法Davis^[1]。Argudo等人^[2]基于符号距离场提出双谐波场的体积网格修复算法, 该算法对处理岛屿和一般复杂拓扑孔洞取得了很好的修补效果。Carr和张丽艳等人^[3-5]运用径向基函数来模拟孔洞隐式曲面和通过径向基函数进行孔洞修补。Centin等人^[6]结合泊松重建和三角网格化方法也对不同类型孔洞取得理想的修补结果。Wang等人^[7]提出一种使用有序双法线和高斯映射来恢复缺失特征曲线和拐角的补洞方法, 应用于CAD模型孔洞修补中。刘振等人^[8]应用动态规划方法并根据匹配特征线信息对显著几何特征取得不错的恢复效果。Qiu等人^[9]利用均值坐标和自动控制变形程度避免碰撞的方法，把补洞应用在牙齿和牙龈的孔洞上也取得了令人满意的结果。Pérez和Guo等人^[10-11]对已有补洞算法进行了实验分析，可以发现，径向基函数能够取得更好的补洞效果^[12]。

但是对于基于径向基函数隐式曲面表示的孔洞修补算法，一个主要问题是隐式曲面对于散乱点噪声十分敏感, 从而导致修复效果呈现不连续.因此, 论文提出采用对孔洞及周围邻域进行拉普拉斯平滑处理.然后通过插值细分得到采样点.最后通过径向基函数通过隐式曲面完成约束求解, 并通过梯度下降得到符合孔洞曲率的修补曲面.通过实验分析表明，算法对于不同的扫描数据，修复效果都要优于逆向软件GeoMagic和文献[13]。

为了能够消除孔洞边界存在的突变点，在这里采用拉普拉斯对孔洞进行平滑处理。虽然径向基函数对散乱点隐式曲面表示有最好的逼近效果。但是径向基函数在逼近曲面计算中, 对参与计算的噪声散乱点十分敏感, 对有噪声的散乱点集无法很好地计算出符合曲率变化的隐式曲面。特别是有部分噪声点参与计算的隐式曲面是畸形的, 甚至无法计算的情况，无法很好地拟合原始曲面。因此在针对有噪声点的隐式曲面求解前，需要结合滤波相关算法，先对孔洞周围噪声点进行基础光顺, 以达到径向基函数求解对散乱点的基础要求。

拉普拉斯算子是 $n$ 维欧几里得空间中的一个二阶微分算子。在3维笛卡儿坐标系中表示所有非混合二阶偏导数, 拉普拉斯算子表示为

$ \Delta = {\Delta ^2} = \frac{{{\partial ^2}}}{{\partial {x^2}}} + \frac{{{\partial ^2}}}{{\partial {y^2}}} + \frac{{{\partial ^2}}}{{\partial {z^2}}} $

(1)

在3维空间模型应用中，拉普拉斯滤波可以看成是一种逐步扩散的过程，即

$ \frac{{\partial {p_i}}}{{\partial t}} = \lambda L\left( {{p_i}} \right) $

(2)

通过上述扩散过程，曲面的一些钉状物、噪声起伏会很快扩散到周围领域中, 最终使得整个曲面变得光顺。若拉普拉斯算子采用显式欧拉积分法表示，则

$ p_i^{n + 1} = \left( {1 + \lambda {\rm{d}}t \cdot L} \right)p_i^n $

(3)

拉普拉斯平滑就是找到一个顶点的邻域，然后把该顶点移动到邻域的重心位置上，即

$ L\left( {{p_i}} \right) = {p_i} + \lambda \left( {\frac{{\sum\limits_{j = 1}^k {{\omega _j}{q_j}} }}{{\sum\limits_{j = 1}^k {{\omega _j}} }} - {p_i}} \right) $

(4)

式中, $ {q_j} $表示点 $p$ 的 $k$ 邻域内的点。在网格模型中用一阶或多阶邻域得到：$ \lambda $表示控制收缩度的较小的正数；$ \omega $表示邻域点的权重系数，此处应用中可以简单地用邻域点数量的倒数表示。拉普拉斯光顺利用周围邻域点，把噪声能量扩散到局部邻域从而达到光顺的目的。拉普拉斯光顺实现最简单，光顺效果较好，能较快地收敛到计算孔洞曲面的光顺位置, 很好地符合应用的需求。

针对径向基函数模拟隐式曲面, 需要比计算隐式曲面更大的邻域数据点集来进行拉普拉斯光顺计算。例如需要孔洞周围三阶邻域数据来计算径向基函数, 在计算拉普拉斯基础平滑光顺时，可以选择五阶邻域进行拉普拉斯光顺，一方面由于本身拉普拉斯在平滑过程中有较大的收缩，另外一方面为了后续补洞结构与周围有更加自然的衔接过渡。

在完成拉普拉斯平滑后，接下来就需要根据孔洞周围数据，来构建符合孔洞曲率变化的隐式曲面。

对于孔洞采样点的隐式曲面重建, 具体问题描述为：给出R³空间的 $n$ 个散乱点$ \left\{ {{c_1}, {c_2}, \cdots, {c_n}} \right\} $, 每一个散乱点都对应一个约束值$ \left\{ {{h_1}, {h_2}, \cdots, {h_n}} \right\} $, 找到一个隐式函数方程 $f(r)$ , 使得$ f\left( {{c_i}} \right) = 0, i = 1, 2, \cdots, n $。对于隐函数 $f$ 的构造采用径向基函数，主要是在隐式曲面重建逼近中, 径向基函数有最好的逼近结果。

在3维空间中基于径向基函数的隐式曲面方程可以定义为

$ f\left( r \right) = \sum\limits_{j = 1}^n {{\omega _j}\phi \left( {r - {c_j}} \right)} + P\left( r \right) $

(5)

式中， $r$ 表示逼近隐式曲面的插值散乱点, $ {\omega _j} $表示每个插值点在径向基函数中的权重, $P(r)$ 为一个一阶多项式, 保证了隐式曲面的线性和连续性, 对每一个数据点 $(x, y, z)$ , 其表达式为

$ P\left( r \right) = {p_0} + {p_1}x + {p_2}y + {p_3}z $

(6)

本文进行散乱点插值时采用的径向基核函数为

$ \phi \left( r \right) = {\left\| r \right\|^3}。$

令$ {c_i} = \left( {c_i^x, c_i^y, c_i^z} \right) $, 为求解权重值$ {\omega _j} $和多项式系数$ {p_0}, {p_1}, {p_2}, {p_3} $可用插值条件求解，即

$ \begin{array}{*{20}{c}} {f\left( {{c_i}} \right) = \sum\limits_{j = 1}^n {{\omega _j}\phi \left( {{c_i} - {c_j}} \right) + P\left( {{c_i}} \right)} = {h_i}}\\ {i = 1,2, \cdots ,n} \end{array} $

(7)

和正交条件

$ \sum\limits_{j = 1}^n {{\omega _j}} = \sum\limits_{j = 1}^n {{\omega _j}c_j^x} = \sum\limits_{j = 1}^n {{\omega _j}c_j^y} = \sum\limits_{j = 1}^n {{\omega _j}c_j^z} = 0 $

(8)

令$ {\phi _{ij}} = \phi \left( {{c_i} - {c_j}} \right) $，则可以得到线性方程

$ \left[ {\begin{array}{*{20}{c}} {{\phi _{11}}}&{{\phi _{12}}}& \cdots &{{\phi _{1n}}}&1&{c_1^x}&{c_1^y}&{c_1^z}\\ {{\phi _{21}}}&{{\phi _{22}}}& \cdots &{{\phi _{2n}}}&1&{c_2^x}&{c_2^y}&{c_2^z}\\ \vdots&\vdots &{}& \vdots&\vdots&\vdots&\vdots&\vdots \\ {{\phi _{n1}}}&{{\phi _{n2}}}& \cdots &{{\phi _{nn}}}&1&{c_n^x}&{c_n^y}&{c_n^z}\\ 1&1& \cdots &1&0&0&0&0\\ {c_1^x}&{c_2^x}& \cdots &{c_n^x}&0&0&0&0\\ {c_1^y}&{c_2^y}& \cdots &{c_n^y}&0&0&0&0\\ {c_1^z}&{c_2^z}& \cdots &{c_n^z}&0&0&0&0 \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{\omega _1}}\\ {{\omega _2}}\\ \vdots \\ {{\omega _n}}\\ {{p_0}}\\ {{p_1}}\\ {{p_2}}\\ {{p_3}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{h_1}}\\ {{h_2}}\\ \vdots \\ {{h_n}}\\ 0\\ 0\\ 0\\ 0 \end{array}} \right] $

(9)

由于上述方程是对称和半正定的, 所以上述线性系统有唯一解。

上述用于在隐式曲面插值中, 曲面上的点约束为$ f\left( {{c_i}} \right) = 0, \;i = 1, 2, \cdots, n $, 如果直接代入会出现无意义的平凡解, 所以Turkey等人^[15]为每一个采样点增加了离面约束点, 取自该点的法线方向$ f\left( {{{c'}_i}} \right) = {d_i}, \;\;i = 1, 2, \cdots, n $。针对这两类点定义为边界约束点和法向约束点, 实际应用中取法线约束点$ {{c'}_i} = {c_i} + 0.1 * \overrightarrow {{n_i}} $在隐式曲面中取$ f\left( {{{c'}_i}} \right) = 1, \;\;i = 1, 2, \cdots, n $。

则上述构建(2×$n$+4, 2×$n$+4)隐式曲面方程的采样数据中, 其中 $n$ 为边界约束点的个数.把边界约束点和法线约束点代入方程系统式(9)可以利用LU分解或者直接求解得到散乱点权值和多项式系数$ \left\{ {{\omega _1}, {\omega _2}, \cdots, {\omega _n}, {p_0}, {p_1}, {p_2}, {p_3}} \right\} $。将得到的结果代入隐式方程式(5)得到孔洞隐式曲面方程：

$ \begin{array}{*{20}{c}} {f\left( {x,y,z} \right) = }\\ {\sum\limits_{j = 1}^n {{\omega _j}{{\left( {\sqrt {{{\left( {x - c_j^x} \right)}^2} + {{\left( {y - c_j^y} \right)}^2} + {{\left( {z - c_j^z} \right)}^2}} } \right)}^3}} + }\\ {{p_0} + {p_1}x + {p_2}y + {p_3}z} \end{array} $

(10)

由于边界约束点和法线约束点本身是包含孔洞多边形曲率变化信息的, 这样计算的径向基函数拟合孔洞曲面是和孔洞周围曲面曲率变化一致的, 调整插值点到隐式曲面后, 其修补效果很好地拟合了周围曲率的变化。由于问题需求一致, 求解边界约束点和法线约束点的过程与已有的方法基本一致。

使用径向基函数插值隐式曲面求解过程中, 当插值散乱点数量增多, 线性系统增大, 造成求解效率很低。然而, 由于补洞插值数据点仅取自孔洞多边形周围, 对于一般孔洞, 点量相对较少, 计算效率较高。

在完成孔洞隐式曲面求解之后，接下来需要再针对孔洞进行插值细分，并通过梯度调整，从而使得孔洞内部的散乱点能够符合周围曲率，维持隐式曲面的特征。

为了提高算法效率，论文对孔洞边界为凸多边形时在插值细分上采用重心法。首先求出孔洞多边形重心，然后对计算出的重心, 重新连接孔洞多边形顶点, 组成新的三角面，不断如此迭代求重心连接细分，定义孔洞多边形重心和孔洞点连接的新三角面的重心为边缘重心(图 1黑色点)，连接对第1次细分的三角面的重心构造临时孔洞多边形, 定义为收敛孔洞多边形(图 1红色多边形)，构造收敛孔洞多边形的点和原重心连接的三角面的重心定义为收敛重心(图 1红色点)，本文快速插值点就是由边缘重心和收敛重心组成。如图 1所示。

具体孔洞多边形快速插值步骤:

1) 首先计算孔洞多边形重心 $G$ , 连接 $G$ 与孔洞多边形各顶点构成新三角面求一次重心, 做第1次插值。

2) 对上述插值点与 $G$ 构成的新三角面分别求重心, 连接各三角面顶点做二次细分。

3) 重复步骤2)迭代3次即可达到所需点量, 迭代次数可以自行设定。

4) 针对第1次细分求解的三角面重心, 重新连接为新的收敛孔洞多边形, 用新的收敛空洞多边形顶点重新连接重心 $G$ 细分。

5) 重复步骤4)迭代4~6次即可基本收敛到孔洞重心位置。

经过上述步骤就可以得到一个孔洞足够的插值数量点细分孔洞。

由于快速插值点并没有考虑孔洞曲面的曲率信息, 所以需要对快速插值点进行顶点调整, 使其符合隐式曲面。本文参照文献[14]的思想，应用梯度下降法对插值点进行调整。由于梯度方向是下降最快的方向，所以可以不断迭代最快地达到所要调整的位置。则利用第2节求解孔洞隐式曲面方程为 $f(r)$ , 其梯度为

$ \nabla f = \left( {\frac{{\partial f}}{{\partial x}},\frac{{\partial f}}{{\partial y}},\frac{{\partial f}}{{\partial z}}} \right) $

(11)

$ \begin{array}{*{20}{c}} {\frac{{\partial f}}{{\partial x}} = 3\sum\limits_{j = 1}^n {\left[ {{\omega _j}\left( {x - c_j^x} \right) \times } \right.} }\\ {\left. {\left( {\sqrt {{{\left( {x - c_j^x} \right)}^2} + {{\left( {y - c_j^y} \right)}^2} + {{\left( {z - c_j^z} \right)}^2}} } \right) + {p_1}} \right]} \end{array} $

$ \begin{array}{*{20}{c}} {\frac{{\partial f}}{{\partial y}} = 3\sum\limits_{j = 1}^n {\left[ {{\omega _j}\left( {y - c_j^y} \right) \times } \right.} }\\ {\left. {\left( {\sqrt {{{\left( {x - c_j^x} \right)}^2} + {{\left( {y - c_j^y} \right)}^2} + {{\left( {z - c_j^z} \right)}^2}} } \right) + {p_2}} \right]} \end{array} $

$ \begin{array}{*{20}{c}} {\frac{{\partial f}}{{\partial z}} = 3\sum\limits_{j = 1}^n {\left[ {{\omega _j}\left( {z - c_j^z} \right) \times } \right.} }\\ {\left. {\left( {\sqrt {{{\left( {x - c_j^x} \right)}^2} + {{\left( {y - c_j^y} \right)}^2} + {{\left( {z - c_j^z} \right)}^2}} } \right) + {p_3}} \right]} \end{array} $

具体算法流程如下:

1) 首先计算 $f\left( r \right)$ 的梯度, 对一个快速插值点, 计算梯度下降逼近新点的位置

$ {r_{k + 1}} = {r_k} - \frac{{f\left( {{r_k}} \right)}}{{{{\left\| {\nabla f\left( {{r_k}} \right)} \right\|}^2}}}\nabla f\left( {{r_k}} \right) $

(12)

2) 计算新点位置和梯度调整前位置的距离$ \left\| {{r_{k + 1}} - {r_k}} \right\| $, 如果该距离小于某一误差值$ \varepsilon $ (本文实验中为0.001), 则认为梯度调整到了隐式曲面的位置, 该点的迭代结束; 否则, 以$ {r_{k + 1}} $代替$ {r_k} $继续迭代。

3) 所有插值点都梯度调整到隐式曲面上后结束。

由于边界约束点和法线约束点计算出的隐式曲面方程拟合的是孔洞周围曲面, 充分考虑了曲率信息, 所以梯度调整到隐式曲面后, 对孔洞的修补的结果符合孔洞周围曲率, 有理想的修补效果。

为了验证论文算法的有效性，在主频2.4 GHz、内存4 GB的Intel i5上对不同网格模型孔洞进行修补和实验比较。比较对象采用目前逆向工程使用最为广泛的商业软件Geomagic Studio.该系统具有强大的拓扑修复和补洞效果, 以及文献[13]中被广泛应用的波前法进行实验对比。

在实验数据上采取两类数据，一类是学术界使用较为广泛的大卫、弥勒佛数据(图 2)；另外一类是机械部件(图 3)。采用国产Hscan330 3维扫描仪得到扫描数据。首先针对大卫头像得到扫描数据并进行修复处理。从大卫头像中针对有明显曲率变化的部位进行孔洞修复，结果表明即使在有明显曲率变化的颈部也可以达到令人满意的修补效果。本文算法修补结果在有很好拟合孔洞边缘曲率变化的前提下，比Geomagic和文献[13]有更加自然的边缘过渡。在实验对比图中，给出了孔洞修复前后的局部放大和比较效果。从细节对比图中可以发现, Geomagic补洞的孔洞曲面虽然拟合了周围曲率变化，但是其边缘过度有更加明显的修补痕迹；文献[13]中的修补效果在孔洞边缘曲率拟合以及孔洞曲面边缘过度方面都无法达到令人满意的结果，而本文算法的补洞结果更为平滑自然，大大提高了逆向工程性能。

进一步地，针对弥勒佛进行扫描和修复实验.和大卫数据相比较, 弥勒佛数据细节更为丰富, 补洞难度更大。本文算法在边缘曲率拟合程度以及补洞曲面边缘过度方面也较Geomagic和文献[13]好，达到了令人满意的修补效果。

为了进一步验证算法在逆向工程应用的有效性。论文对机械扫描工件进行了孔洞修补的实验对比，包括变速轮和手刹器。可以发现，本文算法即使在有噪声干扰的情况下依旧可以达到满意的修补效果。

表 1给出了孔洞修补的约束点数量和新增顶点数量统计以及隐式曲面求解和插值顶点调整时间统计。从表 1可以看出，在点量相对较少的情况下, 求解速度是很快的, 另外在实验中拉普拉斯平滑时间几乎可以忽略。

由于径向基函数表示的隐式曲面方程对没有尖锐特征的光顺曲面模拟较好, 所以针对有尖锐特征处出现的孔洞修补效果并不是十分理想，从图 4所示棱角处修补效果，可以看出基于径向基函数的修补不能很好地模拟尖锐特征, 修补结果还是有很明显的光顺痕迹。

论文提出了一种基于细分插值和径向基函数的孔洞修补算法。该算法针对噪声孔洞的采样点问题通过拉普拉斯平滑、细分插值和梯度调整提高数据质量，从而使得重建的曲面能够更好地维持曲率特征。和商业软件Geomagic以及文献[13]算法相比较，本文算法可以针对有明显曲率变化和周围有噪声的孔洞达到理想的修补效果。目前算法已经作为某公司的3维扫描仪配套软件发布, 算法性能得到了充分肯定。

下一步需要针对孔洞边界复杂狭长，在尖锐特征处出现的孔洞进行研究, 从而使得细分插值方法更加有效, 孔洞特征恢复更好, 进而使算法能够适应各种形式的孔洞问题。