一种基于Hermite径向基函数的三角网格补洞方法

摘要

本发明公开了一种基于Hermite径向基函数的三角网格补洞方法，该方法包括以下步骤：提取孔洞边界后，基于孔洞周围的顶点以及法线信息，利用Hermite径向基函数插值出隐式曲面；将边界投影到平面，然后对孔洞边界进行限定Delaunay三角化，并对三角化后的网格进行细分处理；最后将新增的三角形调整到隐式曲面上，完成修补。该方法具有较高的鲁棒性而且不需要人工参与，对曲率变化较大的网格孔洞具有好的修补效果。

说明书

技术领域

本发明属于几何处理和三维扫描模型修补领域，具体涉及基于Hermite径向基函数的三 角网格补洞。

背景技术

三角网格是几何模型的一种重要表示方法，广泛应用在虚拟现实、逆向工程、计算机图 形学等领域中。随着激光三维扫描技术的快速发展，由扫描重建的三维模型得到越来越广泛 的应用。但是由于测量原理(如对黑色不敏感)的限制、物体表面反射性的影响或被扫描物 体自身结构的遮挡，扫描得到的模型往往存在着孔洞。这些孔洞极大地制约着模型的后续应 用，比如3D打印要求全封闭模型，大量图形学中的形状分析算法也要求模型不存在孔洞。 因此，孔洞的修补对模型的应用具有重要意义。

针对网格补洞人们已经做了大量研究，主要可分为两类方法：基于表面网格的方法和基 于体数据的方法，这两种方法各有优势。基于表面网格的方法能较好地保持未缺失部分的细 节，同时也能更好地利用原网格的性质，适应于在细节要求较高的场景。基于体数据的方法 采用全局处理方式，先将网格表示的模型转化为体数据表示的模型，然后对体数据进行补洞 操作。该方法使用于对算法稳定要求高的场合。

基于表面网格的补洞主要有两大类：根据孔洞及其周围信息进行曲面拟合和直接对缺失 区域进行预测和修补。在拟合方法中，Branch在孔洞提取出来后，用径向基函数(Radial Basis  Function，RBF)对孔洞缺失的部分进行拟合，然后再对孔洞区域进行采样得到补洞后模型。 由于径向基函数本身的性质，该方法在某些情况下不能拟合出符合孔洞周围网格的曲面。 Wang在提取了孔洞边界后，利用组成分分析法(Principal Component Analysis，PCA)得到 孔洞上边界点的拟合平面，将孔洞边界投影到拟合平面上，计算这些点到拟合平面的距离作 为高度场。将高度场带入移动二乘法(Moving Least Square，MLS)中得到拟合曲面。对拟 合平面上的孔洞边界进行采样，将采样点带入到MLS所确定的方程中得到最后的补洞模型。 基于MLS的补洞算法对非平面孔洞具有一定的补洞效果，但是对于高度弯曲的孔洞修补效 果不理想。

上述拟合方法对于曲率变化较大的孔洞区域修补效果不理想，而且有些方法在某些情况 下完全失效。最近，Macedo提出Hermite径向基函数(Hemite Radial basis Fucntion,HRBF) 并将其应用到点云重建上，这也成为了现阶段研究的热点。HRBF从理论上解决了RBF存在 的由于不合适约束点选取而导致的错误重建曲面。同时将法线信息考虑到重建系统中，使得 重建曲面更能反映约束点性质。对于高曲率变化的孔洞网格能够忠实地重建网格能够忠实地 反应曲率的变化。

在直接预测和修补方法中，Liepa以角度和面积作为约束，采用动态规划策略对孔洞进 行三角划分；然后根据孔洞周围网格的密度细分粗修复的网格，在细分过程中使用边交换保 证Delaunay性质得到满足；最后采用伞状平滑算子对网格进行光滑。该方法对孔洞周围曲 率变化较大的模型修复效果不好，而且还可能产生自交叉。Xie提出一种基于差分进化算法 的补洞方法：通过最小化点对之间的距离，对孔洞边界的采样点进行配对。然后根据配对结 果生成交叉平面，结合该平面与原孔洞边界及其一环邻域交点的切线，在交叉平面上采用差 分进化算法预测可能的采样点，最后将采样点投影到平面进行三角化。

在基于体数据的方法中，Ju采用八叉树对模型进行空间划分，在八叉树的结构中找到与 边界边相交的体素，然后采用类似Bresenham算法的思想，在三维空间生成一个连接这些边 界边的平面，对孔洞进行修补。最后用dual contouring重建模型。Nooruddin将网格表示的 模型转化为均匀空间体素，利用校验数(Parity Count)和光线穿透(Ray Stabbing)处理在 转化过程中产生的不同状态。然后将形态学算法中的开(open)和闭(close)算子运用到三 维空间完成空间体素的孔洞修补。最后采用Marching Cubes算法提取出等值面。相对于基于 网格修补的算法，体数据方法具有鲁棒性高、保证流形等优点。但是采用了体素的中间表示， 导致修补后的模型大量丢失原模型的细节。对于曲率变化较大的模型没有很好的修补效果。

综上所述，基于表面网格的方法在处理平面孔洞上较好的效果，但是对于曲率变化较大 的网格不能有较好的处理，而且有些方法还可能崩溃。基于体数据的方法能处理平面孔洞， 但是对于曲率变化较大的孔洞网格不能产生直观的补洞效果。而且由于体数据的方法采用全 局的处理方式，对于未缺失部分的网格细节丢失较多。本发明能够对曲率变化大的孔洞产生 符合孔洞周围网格曲率变化趋势的模型，同时保证修补后的模型在密度上和原模型保持一 致。

发明内容

本发明要解决的技术问题是：克服了现有技术对高曲率变化孔洞修补的不足，为产生能 反映孔洞周围网格变化的完整模型，提供一种三角网格补洞方法。该方法主要利用Hermite 径向基函数产生符合孔洞周围网格曲率变化的拟合曲面，并在孔洞区域惊醒采样，然后将采 样点映射到拟合曲面上。

本发明解决上述技术问题的技术方案为：基于Hermite径向基函数的三角网格补洞方法， 包括如下步骤：

1)生成拟合曲面

步骤(1)、在网格模型中检测出孔洞边界。半边数据结构的特点是将网格中的每条边分 解为带方向的两条半边(half-edge)，这两条半边的方向相反。每条半边记录与其相邻的三角 形和两个端点，每个三角形记录一条组成它的半边，每个顶点记录一条以它为出顶点 (outgoing vertex)的半边。半边数据结构保存了网格的拓扑信息，利用这些拓扑信息可以 快速检索到每个元素(点、线、面)的邻接信息。定义只与一个三角形相邻的边为边界边， 在该结构中遍历网格的边数据，找到一条种子边界边，然后沿着这个种子边界边逆时针找到 下一条与其相邻的其他边界边直到重新遍历回种子边界。这些边界边组成的环形即为孔洞边 界。

步骤(2)、根据孔洞边界及其邻域，利用Hermite径向基函数得到隐式曲面。得到孔洞 边界后，由于孔洞上的顶点不具有完整的一环领域，所以不能对其进行准确的法线计算。选 取孔洞边界的三环邻域，以三环邻域中顶点的位置和法线信息作为Hermite径向基函数输入， 重建出符合孔洞周围网格曲率变化的隐式曲面。

2)孔洞边界的三角划分

步骤(3)、利用PCA方法生成为孔洞边界上顶点的拟合平面，并将其投影到该平面上。 直接在三维空间进行三角化是一件很困难的事，而且还可能引入一些其他缺陷。而在平面上 进行三角化是一件已经被深入的研究的问题，所以第一步为利用PCA方法得到孔洞边界上 顶点的拟合平面，然后将孔洞边界投影到该平面上。PCA方法得到的平面实际是孔洞边界上 顶点的最小二乘拟合。

步骤(4)、若存在自交叉的边界，还需优化PCA方法的最小特征值，将孔洞分成小孔 洞。对于具有高度弯曲孔洞的模型，步骤(3)中的投影可能在拟合平面上产生自交叉边界， 所以还需要分割。将孔洞边界分成两部分，分别计算每部分PCA方法的最小特征值λ_i3和λ_j3， 假设λ_i3≥λ_j3。取当λ_i3为最小时的分割将孔洞分解为小孔洞。保证每个小孔洞投影到拟合平 面上的边界不会产生自交叉。

步骤(5)、在平面上进行基于扫描线的限制性Delaunay三角化，并将结果反投影回三 维空间。得到拟合平面上孔洞的投影后，采用基于扫描线的限制性Delaunay三角化 (Constrained Delaunay Triangulation，CDT)对投影多边形进行三角剖分，CDT方法和传统 的Delaunay方法相比在速度上有明显的优势。然后将三角划分的结果反投影回原模型。

3)对孔洞区域进行采样

步骤(6)、对新增三角形进行满足密度属性的细分操作。为了保证新增三角形的密度和 孔洞周围网格密度保持一致，同时也为了能将新生成的三角形映射到步骤(2)中所生成的 曲面上。为孔洞边界上的每个顶点定义一个密度属性：该顶点相邻的所有边长度的平均值。 对步骤(5)中已经得到的新增三角形进行细分操作。细分的准则为以新增三角形重心为新 增点，插值得到该新增点的密度属性，如果该密度属性和新增点到原来三角形的三个点距离 满足一定的比例关系，则对该三角形进行细分。在细分过程中利用边交换保证Delaunay性 质得到保证。这样就完成了对孔洞区域的采样。

4)将新增点映射到隐式曲面上

步骤(7)、利用梯度下降法将细分的新增顶点映射到隐式曲面上。将步骤(6)中已经 得到的采样点映射到步骤(2)中生成的隐式曲面上。梯度下降法是一种求最优化值的方法， 该方法不断迭代调整步长以到达函数的最优值。将HRBF的梯度和每个新增点带入梯度下降 的迭代公式中，得到最终的补洞网格。

本发明与现有技术相比优点在于：本发明对于周围网格曲率变化高的孔洞有较好的修补 效果，而且具有更高的稳定性。本发明主要有两点贡献：第一，采用Hermite径向基函数对 孔洞进行拟合，得到符合周围网格曲率变化的拟合曲面。第二，给出了一种新的孔洞区域采 样方法，该方法采用平面三角化和基于密度属性的细分操作。

附图说明

图1为算法整体流程图；

图2为半边数据结构示意图；

图3为孔洞邻域示意图；(a)一环邻域，(b)二环邻域；

图4为投影产生自交叉孔洞的示意图；

图5为CDT初始化示意图；

图6为CDT顶点事件示意图；(a)当前顶点在波前上投影得到的边，(b)将该边的两个端 点与当前顶点连接组成三角形，(c)若不满足Delaunay性质，则进行边交换；

图7为CDT边事件示意图；(a)与限定变相交的三角形，(b)将这些三角形删除，并根据 限定变分成上下两部分，(c)分别对这两部分进行三角化；

图8为本发明中间结果和最后结果效果图；(a)带孔洞的模型，(b)根据孔洞及邻域用 Hemite径向基函数拟合得到的曲面，(c)孔洞上顶点的拟合平面以及其在平面的投影，(d)平 面三角化完成后反投影回三维空间，(e)对新增网格进行细分得到采样点，(f)将采样点映射到 曲面上；

图9为本发明和传统RBF方法的对比效果图。从左到右：(a)带孔洞的兔子模型，基于 RBF方法后修补的模型，本方法修补后的模型，(b)带孔洞的机械零件模型，RBF方法修补 后的模型，本方法修补后的模型。

具体实施方式

下面结合附图以及本发明的具体实施方式进一步说明本发明。

步骤(1)在网格模型中检测出孔洞边界。首先用半边数据结构组织网格中的点、线、 面等元素。如图2所示，该结构的特点是将网格中的每条边分解为方向相反的两条半边 (half-edge)。在半边数据结构中每个顶点存放一条以此作为起始端点的半边(图中a)，每 个三角形存放组成它的一条半边(图中b)，每条半边存放如下内容：该半边的终止点(图中 c)，该半边属于的三角形(图中d)，该半边逆时针遍历的下条半边(图中e)，与该半边相 反方向的半边(图中f)。半边数据结构保存了网格的拓扑信息，利用这些拓扑信息可以快速 检索到每个元素(点、线、面)的邻接信息。定义只与一个三角形相邻的边为边界边，在半 边结构中遍历网格的边数据，找到一条种子边界边，然后沿着这个种子边界边逆时针找到下 一条与其相邻的其他边界边直到重新遍历回种子边界。这些边界边组成的环形即为孔洞边 界。

步骤(2)根据孔洞边界及其邻域，利用Hermite径向基函数得到隐式曲面。在得到孔洞 边界后，由于孔洞上的顶点不具有完整的一环邻域，所以不能对其进行准确的法线计算。孔 洞的邻域定义如下：

假设S是孔洞的边界边所组成的单纯形集合，M是模型网格的组合流形。N(S,M)是由至 少和S存在一个顶点相邻的三角形组成的单纯形集合，称N(S,M)是S的单纯形邻域。定义孔 洞的k阶(k-order)单纯形邻域为N(N(…N(S,M)…),M),M)，表达式有k-1层嵌套，简称k阶邻 域，记为N_K(S,M)。令孔洞边界点的集合为v_h，孔洞k阶邻域中的顶点集合为v_k，则v_k-v_h为孔洞的k环(k-ring)顶点。如图3所示，a是一环邻域，b是二环邻域。选取孔洞边界的三 环邻域，以三环邻域中顶点的位置和法线信息作为Hermite径向基函数输入，重建出符合孔 洞周围网格曲率变化的隐式曲面。

Hermite径向基函数插值问题可以描述如下：对于三维空间中给定的散乱点集 在一个函数空间H找出一个隐式函数f:R³→R，使其满足 f(xⁱ)＝0,在逼近理论中，该问题被称为一阶Hermite插值问题。假设有一个 Hilbert空间H其对偶空间为H^*。给定一个数据集其中γ_i是连续线性泛函，c_i是 实数，即γ_i∈H^*，c_i∈R。目标是找到一个函数f∈H，使得r_i(f)＝c_i，这样的函数称为H 空间的广义插值。对于H空间，这样的广义插值有无数个,但是满足H范数最小的只有一个。 这个函数可以由里斯表示vⁱ∈H的线性组合唯一确定：而且对于每一个 u∈H，γ_i(u)＝<vⁱ,u>。将f^*和γ_i联立，得到线性系统Aα＝c，其中矩阵A的元素为 (A)_i,j＝γ_i(v^j)＝<vⁱ,v^j>_H。在Hermite插值问题中其中D^λ是一个微分算子， λ∈Nⁿ，λ_j表示u中的第j个变量被微分的次数。构建这样的H空间的方法如下。

给出一个正定的径向基函数φ:R₊→R，同时定义ψ＝φ(||·||)。存在唯一的Hilbert空间， 该空间的里斯表示为$v^i\left(x\right)=-\left(D^{{\lambda }^i}\psi \right)(x-x^i),$内积为${<v^i,v^j>}_H=-\left(D^{{\lambda }^i}{D}^{{\lambda }^j}\psi \right)(x^i-x^j).$针对一 阶Hermite插值，当选用全局支撑函数时其表示形式如下。

$f\left(x\right)=\underset{j=1}{\overset{N}{\Sigma }}{\alpha }_j\psi (x-x^j)-{\Sigma }_{j=1}^N<{\beta }^j,▿\psi (x-x^j)>+P\left(x\right)---\left(1\right)$

其中ψ(x)＝φ(||x||)，||·||表示欧氏距离，对于三维重建一般取φ(x)＝x³。 P(x)＝a+bx+by+cz，(x,y,z)是点x的坐标。将插值条件带入上述方程中可以求出系数 α_j∈R和β^j∈R³：

$\underset{j=1}{\overset{N}{\Sigma }}{\alpha }_j\psi (x-x^j)-{\Sigma }_{j=1}^N<{\beta }^j,▿\psi (x-x^j)>+P\left(x\right)=0---\left(2\right)$

$\underset{j=1}{\overset{N}{\Sigma }}{\alpha }_j▿\psi (x-x^j)-{\Sigma }_{j=1}^N\mathrm{H\psi }(x-x^j){\beta }^j+▿P\left(x\right)=n^j---\left(3\right)$

H是hessian算子，定义为将上述插值条件重新组合得到适合计算形式， 对于每一个(xⁱ,nⁱ)有：

${\Sigma }_{j=1}^N\left(\begin{array}{cc}\psi (x^i-x^j)& -▿\psi {(x^i-x^j)}^T\\ ▿\psi (x^i-x^j)& \mathrm{H\psi }(x^i-x^j)\end{array}\right)\left(\begin{array}{c}{\alpha }_j\\ {\beta }^i\end{array}\right)+\left(\begin{array}{c}p\left(x\right)\\ ▿P\left(x\right)\end{array}\right)=\left(\begin{array}{c}0\\ n^j\end{array}\right)---\left(4\right)$

${\Sigma }_{j=1}^N\left(\begin{array}{cc}p\left(x\right)& ▿p{\left(x\right)}^T\end{array}\right)\left(\begin{array}{c}{\alpha }_j\\ {\beta }^j\end{array}\right)=0---\left(5\right)$

由上式可知，对于三维空间的每个(xⁱ,nⁱ)一次可计算出关于ψ的梯度和hessian矩阵， 得到一个4×4矩阵，然后与多项式项相加。所有点经上式计算后的矩阵组合得到完整的线性 系统。

Hermite径向基函数是传统径向基函数的改进，它从理论上解决了传统径向基函数在选 取约束点(off-surface points)上的主观随意性，极大地增加了基于该函数的重建鲁棒性。除 此以外，由于重建的线性系统将法线作为插值条件，使得重建的模型能更好地反映插值点的 性质。得到的中间结果如图8b所示。

步骤(3)利用PCA方法生成孔洞边界顶点的拟合平面。因为直接在三维空间进行三角 化是一件很困难的事，而且还可能引入一些其他缺陷，而在平面上进行三角化是一件已经被 深入的研究的问题，所以第一步为利用PCA方法得到孔洞边界上顶点的拟合平面，然后将 孔洞边界投影到该平面上。PCA方法得到的平面是孔洞边界上顶点的最小二乘拟合，主要步 骤如下：

a)计算孔洞上边界顶点的平均值其中p_i是孔洞上的边界顶点，n是边界 顶点的个数，是所求平面上的一个点。将该点作为拟合平面上的一个点。

b)得到点后建立矩阵M：M_i＝(x_i-x₀,y_i-y₀,z_i-z₀)，其中M_i是矩阵M的第i行， (x_i,y_i,z_i)是p_i坐标。M具有如下形式：

$M=\left(\begin{array}{ccc}{x}_1-x_0& y_1-y_0& z_1-z_0\\ x_2-x_0& y_2-y_0& z_2-z_0\\ .& .& .\\ .& .& .\\ .& .& .\\ x_n-x_0& y_n-y_0& z_n-z_0\end{array}\right)---\left(6\right)$

c)用SVD计算出M^TM的三个特征值λ₁≥λ₂≥λ₃，以λ₃对应的特征向量作为平面的 法向量。根据和确定平面，将孔洞的边界边投影到该平面上。

步骤(4)若存在自交叉的边界，还需优化PCA方法的最小特征值，将孔洞分成小孔洞。 对于高度折叠的复杂孔洞，如图4，上述投影方法会导致孔洞边界产生自交叉。需要将这类 孔洞分割为两部分，分别投影到两个平面再进行后续操作，在这步中如何获得孔洞的最优分 割是关键问题。在通过PCA获得三维点集投影平面的方法中，最小特征值的值决定了这些 点对平面的拟合程度：λ₃相对于λ₁越小，表示平面和点集的拟合度越高。因此，最小化分割 后两个点集确定的最小特征值中较大的那个得到最优的分割，即解下面的优化问题：

$\underset{1\le i,j\le n}{\min }\max \{{\lambda }_{3,{\mathrm{ps}}_1},{\lambda }_{3,{\mathrm{ps}}_2}\}---\left(7\right)$

其中，PS₁＝{p_j,...p_n,...,p_i-1},PS₂＝{p_i,...p₁,...,p_j-1},分别是PS₁,PS₂通过PCA方法 得到的最小特征值。

将孔洞分解为小孔洞后，能保证每个小孔洞投影到拟合平面上的边界不会产生自交叉。 孔洞边界投影到平面后，为了方便后续的计算的直观性，将该平面旋转到于x-y平面平行的 位置，并记录空间中的点与平面上点的对应关系。得到的中间结果如图8c所示。

步骤(5)对步骤(3)中生成的拟合平面上孔洞的投影采用基于扫描线的CDT进行三 角剖分，CDT方法和传统的Delaunay方法相比在速度上有明显的优势。CDT的步骤如下：

a)初始化。所有点按照y值从小到大排序，对于y相同的点，按照x值排序。为每个 点添加附加信息，标明该点是否为约束边的上端顶点。为了避免过多复杂情况的讨论，添加 两个虚拟点P_-1，P_-2，这两个点分别在所有输入顶点的左下方和右下方,如图5所示。将P_-1， P_-2，P₀连接构成初始三角形，将P₀P_-1,P₀P_-2作为初始扫描波前。

b)扫描。根据y值大小，从下往上移动扫描线。每当扫描线遇到一个顶点便停下来， 若该顶点不是一条约束边的上端顶点则处理顶点事件，否则处理边事件。顶点事件：将P_i垂 直投影到波前得到P_j，确定P_j所在波前的边P_LP_R(图6a)，连接P_i，P_L，P_L组成新的三角形 (图6b)，若该三角形不满足空圆性质，则对边进行调整，该步被称为边交换(图6c)；边 事件：当前顶点若为一条边的上端顶点时，就形成边事件。先确定第一个与该边相交的三角 形，然后将沿着这个三角形遍历到其他与该边相交的三角形(图7a)。遍历完成后，删除与 所有与该边相交的三角形，并将它们的顶点依据与边e的位置进行排序得到位于e下方的Π_L和位于e的上方Π_U(图7b)，最后分别对这两部分进行常规的基于增量算法的CDT(图7c)。

c)后处理。将所有与虚拟点P_-1，P_-2相连接的三角形删除，并添加边界三角形完成三 角化。

当边界边组成的多边形的三角化完成后，依据上步中的对应关系将二维平面上的点反投 影回三维空间。得到的中间结果如图8d所示。

步骤(6)为了保证新增三角形的密度和孔洞周围网格密度保持一致，同时也为了能将 新生成的三角形映射到步骤(2)中所生成的曲面上，需要对新增的三角进行细分操作。具 体步骤如下：

a)为孔洞上的每个边界点v_i都增加一个密度属性θ(x_i)，其定义为：与v_i相邻的边长度 的平均值。

b)对每个新增三角形(v_m,v_n,v_l)，计算其质心v_c和它的密度属性 θ(v_c):＝(θ(v_n)+θ(v_m)+θ(v_l))/3。对每一个i＝m,n,l：如果α||v_c-v_i||＞θ(v_c),α||v_c-v_i||＞θ(v_i)， 那么就将三角形(v_m,v_n,v_l)细分为(v_c,v_n,v_l),(v_m,v_c,v_l),(v_m,v_n,v_c)，同时检验 (v_m,v_n),(v_m,v_l),(v_l,v_n)是否需要进行边交换。

c)如果b中没有新的三角形产生跳到d，否则跳到b。

d)交换所有新增边中的需要交换的边。

在实现中，选取得到的中间结果如图8e所示。

步骤(7)当采样结束后，需要将采样点调整到HRBF拟合的曲面上。本发明采用梯度 下降法完成这一步骤。HRBF的梯度计算方法为：

$f\left(x\right)={\Sigma }_{j=1}^N{\alpha }_j▿\psi (x^i-x^j)-{\Sigma }_{j=1}^N\mathrm{H\psi }(x^i-x^j){\beta }^j+▿P\left(x\right)---\left(8\right)$

H如上述定义为Hessian矩阵。得到梯度后，对于修补网格中每一个新增顶点r₀，采用 下面的迭代公式将其映射到曲面上：

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

满足||r_k+1-r_k||≤ε条件时，迭代停止，ε为指定的误差范围。

本发明的算法用visual studio 2010和OpenGL实现，所用的计算机配置为：主频2.53GH， Xeon E5630处理器和8G内存。在解HRBF所确定的线性系统时采用Eigen库中的LU分解。 方法最终效果图如图8f所示。图9为本方法和传统基于RBF方法的对比图，从左到右分别 是带孔洞的原模型，经RBF方法补洞后模型，经HRFB方法补洞后模型。