0 引　言

由于建模误差、格式转换等一系列问题，模型往往存在法向错误、穿插、狭缝、贴合等一系列问题，我们将这类无法直接应用于数值计算的几何模型称为脏几何^[1]。因此，在传统的网格生成流程中，网格生成前往往需要对几何进行清理操作，用户需要识别所有的几何错误，并依次修复^[2]。这一修复过程往往需要大量人工交互，并借助复杂的图形界面，如CADfix^[3]等专业的几何处理软件，且十分依赖于操作人员的经验。对于大量部件组成的复杂模型，人工修复的成本是难以忍受的。

Shephard等提出了一种基于八叉树的全自动网格生成方法^[4]，基于八叉树的网格生成方法可以一定程度上容忍几何错误，但仍存在一些固有的缺陷。模型的几何特征并非初始构造，该方法需要向基网格嵌入特征，使网格逼近真实几何^[5]，基网格的质量严重影响几何特征嵌入的效果。在嵌入过程中，往往需要引入大量的单元来捕捉细小的几何特征，导致单元数目的大幅上升。

刘等人基于混合曲面表征^{[2, 6-7]}实现了装配体的“曲面嵌入”，可以解决多个部件之间可能存在的相交和重叠问题。其方法是在离散曲面上计算边界相交图，得到离散嵌入结果，之后根据连续-离散的映射关系完成连续曲面边界表征上的曲面自动嵌入^[3]。该算法具有较高的时间效率，也保证了计算结果的几何精度，但其适用范围有限，鲁棒性不足，针对模型贯穿的问题，由于贯穿处缺少边界表征，曲面嵌入失效。

Hu 等通过将三角形-三角形求交升维至四面体-三角形求交，在初始三角形表面约束下生成纯四面体网格，并在有理数范围内实现了精确求交^[8]。这种方法在一定程度上可以解决脏几何模型中出现的相交问题，但其局限性亦存在。一是相交计算十分耗时，二是此类方法仅能设定全局容差值，当缝隙间距大于其允许的容差值时，使用此方法不仅无法消除几何误差，反而可能急剧影响性能和网格质量，而容差值过大时，几何特征不能被完全保留。刘^[9]也采用四面体化的方法来完成脏几何的网格生成，并创新性地引入了边界恢复算法^[10]，从而大幅减少相交计算量。两者的四面体化方法中均要求输入边界法向正确，对于存在法向错误的几何，仍然需要大量的手工交互调整输入模型的法向。

本文针对脏几何中可能出现的法向错误、穿插、狭缝、贴合问题，通过四面体化方法自动修复脏几何，并生成高质量的曲面网格。其核心思想是一个有效的四面体网格的表面一定是一个有效的曲面网格。本文在输入模型的离散表征的基础上直接生成四面体网格，通过离散边界切割四面体解决了模型中存在的穿插问题，基于四面体优化解决狭缝和贴合问题，此时该四面体网格的表面便是有效的曲面网格，提取该表面并应用曲面网格重构提升表面网格质量，得到最终的曲面网格。

相较于八叉树方法，本方法具有更高的几何精度，可以生成贴体的曲面网格；相较于基于混合曲面表征的曲面嵌入方法，本方法具有更强的普适性和健壮性；相较于Hu等的四面体化方法，本方法的主要优点有：

1）应用基于视图采样的方法完成法向修复，允许输入的脏几何中存在法向错误；

2）实现了基于容差场的局部修复与特征保持，避免四面体优化导致的特征丢失或修复不完全；

3）实现了基于多级Hausdorff距离^[11]的网格重构技术，可以在修复的曲面网格的基础上进行高精度的网格重构，避免网格重构导致的表面自相交。

1 算法流程

图1和图2给出了本文所提出的基于四面体化方法的面向脏几何曲面网格生成算法流程，该算法的输入为脏几何模型，输出为保形的高质量曲面网格。该算法主要包含以下步骤：

1）构建混合曲面表征。建立模型的离散-连续表征数据结构^[2]，采用B-rep表示模型的连续表征，采用三角网格表示模型的离散表征，同时维护连续模型和离散模型之间的映射关系。得到离散表征后，我们采用基于视图采样的方法，完成法向调整。

2）初始四面体化。本文采用Bowyer-Watson方法^{[6, 12]}完成四面体化，我们提取离散表征中的所有点，作为B-W算法的输入，构造一个Delaunay三角化。单次B-W增量插点的过程为：（1）在当前Delaunay三角化中找到包含待插入点的基单元；（2）利用单元的相邻关系在基单元附近快速搜索所有打破外接圆准则的单元；（3）删除所有（2）中的单元，形成空腔，连接待插入点与空腔边界顶点，形成新的Delaunay三角化。

3）边界恢复与相交计算。初始四面体化中并不保证包含所有的原始曲面边界，不能包含的边界单元可分为两类。第一类是由于Delaunay三角化只能保证点存在于三角化中，不能保证网格边和网格面也都存在于三角化中，这类边界可以通过边界恢复^{[9-10, 13-14]}的策略来找回原始边界，第二类是由于脏几何模型本身存在相交或重叠的情形，得到的三角化必然无法包含此类边界，如图1（c）所示蓝色区域，这类边界需要通过相交计算对四面体进行裁剪以找回边界。

4）四面体网格质量优化。相交计算仅能解决部件相交的问题，但模型中仍可能存在狭缝等不必要的特征。我们采用四面体网格质量优化的方法，在容差范围内抹除此类特征。

5）曲面网格重构。提取上述四面体网格的表面得到修复后的表面网格，修复后的表面网格解决了原始模型的一系列错误，但网格质量与网格密度均难以达到数值计算的要求，因此需要对4）中得到的表面网格进行重构。

2 算法实现 2.1 法向调整

本文算法的核心思想是生成有效的四面体网格，提取其表面得到最终有效的曲面网格。如图1所示，我们在输入模型的外部加入包围盒，在整个包围盒区域内生成四面体网格，我们计算每个四面体单元相对输入表面的绕卷数^[15]，筛选得到内部的四面体单元，从而方便地提取出有效的表面。绕卷数计算正确的前提是输入表面的法向是正确的，此外，在曲面网格重构中，也需要保持每个面上单元的法向一致。因此在完成混合曲面表征的构造后，我们需要对离散表征的面法向进行修复。

若离散表征中所有的三角单元都是相连的，即两个三角单元之间存在一条共享边，我们可以通过BFS算法^[16]完成面法向调整，即锁定某一个三角单元的法向，以该单元为起点，不断搜索其邻居，使与其相连的单元的法向调整为与该起点单元同向。但考虑最坏情况下，若几乎所有的面片均不存在与其他面片的拓扑连接，此时我们无法通过BFS算法快速地完成法向调整。为了尽可能地解决脏几何中可能出现的法向问题，我们采用基于视图采样的法向调整策略。通过绘制该离散模型在各个方向上的视图，可以计算该模型中任意一个面片法向的可能性。理想情况下，我们在任意一个观测点 $\boldsymbol{v}$ 观察该模型，对于观察到的某一个三角形中的某一点 $\boldsymbol{p}$ ，视线方向 $(\boldsymbol{v}-\boldsymbol{p})$ 与该三角形的法向 $\boldsymbol{n}$ 应满足 $(\boldsymbol{v}-\boldsymbol{p})\cdot \boldsymbol{n}\leqslant 0$ 。

定义三角面片 $ {t}_{i} $ 的法向正确的概率为 $ P\left({t}_{i}\right) $ 。我们对该模型建立一个包围球，从包围球上均匀选取 $ N $ 个观测点对该模型进行绘制，得到 $ N $ 张 $ 1\;024\times 1\;024 $ 的图像。定义 $ {O}_{+}\left({t}_{i}\right) $ 为三角面片 $ {t}_{i} $ 在所有视图中法向与视线夹角大于90°的像素点数， $ {O}_{-}\left({t}_{i}\right) $ 为三角面片 $ {t}_{i} $ 在所有视图中法向与视线夹角小于90°的像素点数。最终可计算得到三角面片 $ {t}_{i} $ 法向正确的概率为：

$ P\left({t}_{i}\right) = \frac{{O}_+\left({t}_{i}\right)}{{O}_+\left({t}_{i}\right)+{O}_-\left({t}_{i}\right)} $

(1)

若 $ P\left({t}_{i}\right) ＞ 0.5 $ ，将保持该面片的法向，否则将该三角面片的法向翻转。在实际应用中，首先采用BFS算法尽可能地将面片组合，然后采用上述视图采样法向调整策略完成法向的全局调整。对于单张视图，我们采用扫描线Z缓冲器算法^[17]，完成视图的快速绘制。图3为法向调整前后的对比，其中绿色表示该三角面片的法向指向体的内部。

2.2 曲面求交

为了恢复模型的所有原始边界，我们使用原始边界的三角面片所在的平面来切割四面体。非退化情况下四面体单元被原始表面单元所在平面切割如图4所示，被切割的四面体单元的侧面将会被截断为一个四边形和一个三角形，该四面体单元将形成一个四面体单元和一个五面体单元。首先对侧面的三个四边形分别进行Delaunay三角化，再对该五面体单元进行分割，在其重心处插入新点，与侧面的新三角形顶点连接后形成分裂后的新四面体单元。贴合问题将产生零体积的四面体单元，本算法中的三角形求交计算全部采用基于GMP的高精度有理数计算^[18]，退化情况下舍弃零体积的新单元即可。基于精确的求交计算，我们可以保证输入模型中的所有自相交问题可以得到修复。

图5展示了利用原始边界切割四面体解决相交问题，图5（a）中红色四面体与蓝色的原始边界出现相交，图5（b）为该四面体的底面。利用原始边界切割四面体，其结果如图5（c）和图5（d）所示。

2.3 四面体网格优化

上述四面体-平面求交可以在有理数的范围内修复脏几何的自相交问题，但模型仍可能存在狭缝等不必要的细小特征。由于前述步骤已经完成了输入模型原始表面的恢复，那么在狭缝处必然存在狭窄甚至是退化的四面体单元，因此优化此处的四面体单元质量，即可抹除此类特征。图6展示了通过四面体优化解决狭缝的详细过程。图示两圆柱之间存在狭缝需要被抹除，经过初始四面体化后，如图6（b）所示，狭缝处产生了图示红色的四面体单元，这类四面体单元的质量较差，经过四面体质量提升后，这类单元被抹去，得到的表面如图6（c）所示。

2.3.1 四面体质量提升

四面体单元质量的提升主要操作为边分裂、边折叠、面交换与点优化四种局部操作^[19]，其示意图如图7所示。本文通过循环执行这四种局部操作来逐步优化四面体单元。

本文采用Rabinovich等^{[8, 20]}提出的能量函数作为四面体单元的质量标准，该函数具有天然的各向同性，可以有效地甄别出各向异性的和负体积的四面体单元，公式如下：

$ {{ E}}_{t} = \dfrac{\mathrm{t}\mathrm{r}\left({\mathit{J}}_{t}^{{\;\rm{T}}}{\mathit{J}}_{t}\right)}{{\mathrm{det}\left({\mathit{J}}_{t}\right)}^{2/D}}$

(2)

其中， $ {\mathit{J}}_{\mathrm{t}} $ 是将单元 $ t $ 转化为正单元（三维时为正四面体）的雅可比矩阵，D为网格维度，衡量四面体单元质量时取3。当某一单元的能量函数值大于临界值时，将被视为差单元，并通过上述质量提升操作提升单元质量。

定义理想单元尺寸为 $ l $ ，上述四种操作中可以影响单元边长的操作为边分裂和边折叠操作。当四面体单元某一边长大于 $ \sqrt{2}l $ 时，将对该边执行边分裂操作；当四面体单元某一边长小于 $ \dfrac{\sqrt{2}}{2}l $ ，且执行边折叠后合并得到的单元质量提高，则对该边执行边折叠操作，边折叠操作是修复狭缝的主要手段。

四面体网格质量优化的目的是解决狭缝或重叠处产生的狭窄的或退化的四面体单元，对于空间中的四面体，其质量评判标准不同。为了加速优化效率，我们采用了多级能量阈值判别方案，即每个四面体单元优化的目标能量值不同，这里我们采用四面体单元至边界的距离来计算该单元的能量阈值。我们令原始表面附近的单元的目标能量值为 $ {E}_{ts} $ ，该三角化中允许的最大能量值为 $ k{E}_{ts} $ ，则空间中某一四面体的目标优化质量 $ {E}_{t} $ 可以通过线性插值得到，即：

$ {E}_{t} = \left[1+\frac{2{d}_{s}}{L}\left(k-1\right)\right]{E}_{ts} $

(3)

其中， $ {d}_{s} $ 为该四面体与输入表面的最短距离。靠近表面的四面体单元，其目标能量值较小，而远离所有表面的四面体单元，其目标能量值较大，我们允许在距离表面远处产生较差的四面体单元。

2.3.2 包络检测与容差场

上述四面体网格质量优化方法中，边折叠与点优化涉及网格点的位移操作，若网格点的偏移量过大，可能导致网格与输入模型之间产生较大的偏差，我们将禁止这样的边折叠与点优化操作。为了确保本算法生成的表面网格与输入模型之间的误差在可控范围内，每次边折叠与点优化的操作执行后，需要进行包络检测以确保生成保形的表面网格，包络检测的方法如下。

记待检测的三角形为 $ T $ ，允许的最大容差为 $\varepsilon$ ，对 $ T $ 进行均匀采样，采样点以边长为 $ d $ 的正三角形均匀分布，并在边界处添加采样，采样点的最大误差为 $ d/\sqrt{3} $ ，如图8所示。若 $ T $ 中的所有采样点与输入边界的最小距离均小于 $\varepsilon-{d}/{\sqrt{3}}$ ，我们则认为该三角形单元位于输入边界的包络内。根据不同的场景需求，我们可以设定不同精度的采样点，随着 $ d $ 的减小，包络检测耗费的时间越长，检测的精度也越高。

由上述采样方法可知，实际采样点允许的容差值小于 $\varepsilon$ ，倘若某一采样点计算得到的最近距离小于 $\varepsilon$ 但大于 $\varepsilon-{d}/{\sqrt{3}}$ ，将会产生错误判断。这种错误的判断将影响网格质量优化的性能与结果，因为一次边折叠或点优化的失败，该单元将继续执行边分裂操作，产生更多的单元，随着单元数量的上升，降低优化质量的同时也会降低网格优化的速度。

为解决上述问题，我们采用多级的包络检测^[8]，第 $ k $ 级包络检测的采样点间距 ${d}_{k} = d/k$ ，若第 $ k $ 级检测中所有的采样点至输入边界的最小距离均小于 $\varepsilon-d_k/ \sqrt{3}$ ，则认为该三角形单元位于输入边界的包络内。

容差 $\varepsilon$ 的设置将会影响包络检测的时间效率与精确性，同时也会影响模型缺陷的修复效果。若要修复模型中的狭缝，那么此处的容差值应当大于狭缝的宽度，但其它区域应当尽可能地保证边界不发生改变。因此，我们定义了基于离散表征的容差场，任意一点处的容差分为两部分，全局容差值 ${\varepsilon }_{0}$ 与局部容差值 ${\varepsilon}_{p}$ ，该点处的容差值取两者中较小的一个。某一部件的 ${\varepsilon}_{0}$ 为该部件的包围盒对角线长度的千分之一，局部容差值 ${\varepsilon}_{p}$ 需根据模型需求，在需要修复的狭缝区域设置较大的容差值，在必须严格保持的特征处设置较小的容差值。如图9所示，我们通过设置楔形体脊边处的容差值为一个较小的容差值，实现特征保持的效果。

曲面求交与四面体质量提升两个环节中均可能产生新的边界点，新的边界点都是由边分裂形成的，故新边界点的容差值应取其边分裂前的两端点处容差值的最小值。

2.4 曲面网格重构

本文中采用经典的网格重构策略完成曲面网格重构，其核心算法与四面体网格重构相同，对需要优化的单元依次执行边分裂、边折叠、边交换与点优化操作，直至所有的表面网格均满足质量要求。

本文中的尺寸场定义采用Liu等^[21]提出的定义在离散边界上的尺寸场，并采用其全局光滑化策略^[22]。曲面网格单元的衡量标准采用面积与边长之比的平方和，表面网格中的每个点都记录其理想边长 ${L}_{{\rm{ideal}}}$ ， ${L}_{{\rm{ideal}}}$ 的值可以通过查询离散尺寸场并线性插值得到，当网格边长大于 $\sqrt{2}{L}_{{\rm{ideal}}}$ 时，执行边分裂操作，当网格边长小于 $\dfrac{\sqrt{2}}{2}{L}_{{\rm{ideal}}}$ 时，执行边折叠操作。

为了生成保形的表面网格，我们希望在重构过程中由于边分裂产生的新点、点优化后的点都能尽可能地贴近原始表面。网格重构中的误差估计通过计算Hausdorff距离^[11-21]来衡量，我们要求所有的网格点距离输入边界的Hausdorff距离均小于 $ {h}_{\mathrm{m}\mathrm{a}\mathrm{x}} $ 。由于表面网格优化中存在点优化的操作，在狭窄的区域有可能出现网格自相交的情况。为了避免这种情况，我们采用多级误差来避免由表面网格重构引起的表面自相交。在执行完毕一次曲面网格重构后，检测所有的表面网格是否存在自相交的情况，对于存在自相交的表面，令 ${h}_{\mathrm{max}}' = {h}_{\mathrm{max}}/{2}$ ，重新进行曲面网格重构，直至表面不存在自相交的区域。

3 数值验证

本算法采用C++语言实现，借助GMP库完成基于有理数的精确计算。测试平台为Intel Core i7-6700，3.4 GHz，内存24 GB。

3.1 网格生成效果

本文选取了靶球、格栅鳍、F-16战斗机、AIM-54导弹、Ford发动机和某压气机六个典型的脏几何模型作为实验对象，模型网格如图10所示。

为验证本文算法的有效性，这里详细展示三个复杂脏几何的网格生成结果，分别为：靶球模型、格栅鳍模型和F-16战斗机模型。

图11展示了靶球模型的几何及其生成结果，该模型包含115个体部件，存在15177处相交，图11（a）为曲面离散后的局部放大结果，可以明显地看出该处存在部件之间的相交问题，经过本文的四面体化方法修复的表面网格如图11（b）所示，部件之间的相交问题已经得到处理。

图12展示了格栅鳍模型的几何及其网格生成结果，该模型的各片格栅鳍之间存在大量的贯穿问题，如图12（a）所示，对于此类模型，曲面嵌入的算法由于贯穿处缺少对应的边界表征，故无法支持部件之间互相贯穿的情形。而本文算法由于不依赖原始的边界表征，可以很好地解决这一问题，网格生成结果如图12（b）所示。

图13展示了F-16战斗机模型的几何及其网格生成结果，该模型由8个部件组装而成，部件之间存在相交与狭缝等问题。如图13（a）展示了该模型存在的部件之间的相交问题，粉色尾鳍与淡蓝色机身之间存在部件相交问题，图13（b）为该处网格的生成结果。图13（c）展示了该模型部件之间存在的狭缝问题，深蓝色为该模型机翼，淡蓝色区域为该模型机身，机翼与机身之间存在狭缝，图13（d）为该处网格生成结果，可以看出狭缝已经被抹除。

F-16战斗机模型中存在一些特殊的几何特征在网格生成中需要得到保留，如图14所示机翼后缘，我们希望在此处能够保留直线特征，后缘处应用2.3.2节所述的局部容差可以实现该处的特征保留，在机翼后缘处设置较小的容差值，实现了该处特征保持的效果，TetWild^[8]在此处生成的网格如图14（a）所示，无法保留此类特征。

此外，由于本文算法的网格重构过程中引入了尺寸场^[21]，可以根据曲率特征实现网格的局部加密，即在曲率较大的区域自动生成较密的网格，图14（b）中的机翼前后缘区域的网格受曲率特征的控制，网格尺寸较小，这一特征也使得该网格更能满足后续数值分析的需要。

AIM-54导弹、Ford发动机和某压气机的部分模型错误与修复结果如图15所示。图15（a）中AIM-54的机翼与机身之间，图15（c）中Ford发动机红圈处存在贴合错误，其修复结果如图15（b）和15（d）所示。图15（e）中压气机存在狭缝，修复结果如图15（f）所示。

3.2 对比与数据统计

常用的网格生成算法均不具备表面修复的能力，Hu等完成的TetWild^[8]旨在任意离散表面上生成四面体网格，具备一定的表面修复能力，3.1节中的结果表明，本文算法具有更好的修复效果，并且由于本文加入了曲面网格重构的策略，可以在保形的基础上尽可能地生成高质量的表面网格。

我们将前述6个典型的脏几何模型作为实验对象，与TetWild^[8]进行网格质量与时间效率对比。由于TetWild只能接受离散输入，且无法处理法向错误的情况，我们将混合曲面表征中已经完成法向调整的离散表征结果作为TetWild的输入。

我们选用平均最小角度和平均最大边长比作为网格质量对比项，平均最小角度为该网格所有单元中每个单元的最小角度的平均值，平均最大边长比为该网格所有单元中每个单元的最大边长比的平均值。前述6个模型分别采用TetWild与本文算法处理所得网格的质量统计与对比结果如表1所示。对比表明，本文算法生成的网格的平均最小角度与平均最大边长比均优于TetWild。

前述6个模型分别采用TetWild与本文算法处理所耗时间的数据统计与对比结果如表2所示。对比表明，TetWild所耗时间是本文算法的3倍及以上，本文算法在时间效率上领先TetWild。

4 结　论

本文针对脏几何存在的法向错误、穿插、狭缝、贴合问题，应用基于四面体化的曲面网格修复方法，开发了适用于脏几何的自动网格生成算法。本算法实现了以往四面体化方法无法完成的法向修复；应用边界恢复算法，可以大幅减少参与求交的三角形数量，提高求交性能；应用分级的误差检测，可以大幅提高包络检测精度，从而保证生成保形的曲面网格；引入容差场，实现了特征保持与特征消除；应用多级误差估计，避免曲面网格重构产生自相交，最终生成了高质量的表面网格。上述对比实验也展示了本算法在修复效果、时间效率与网格质量上的优势。

本文算法需要手动设置容差，从而达到好的修复效果，未来将持续探索更加简便与自动化的容差设置方式。此外，本文算法针对存在破洞的几何尚无法完成修复，需要进一步研究补面的算法。