四面体网格边界恢复改进算法

曾寅家，陈建军，傅珂杰

(浙江大学 航空航天学院，杭州 310027)

0 引言

Delaunay 三角化（delaunay triangulation,DT）是生成四面体网格的主流方法之一，具有运算快速、网格总体质量较高等优点。给定表面网格，DT 方法通过不断插入网格点形成一个由表面网格点和内部点组成的初始四面体网格，该网格满足Delaunay 准则[1]，但不能保证所有边界约束（在这指输入的表面网格边和网格面）都存在于其中，因此需要调用边界恢复算法修改初始体网格，获得与边界约束兼容一致的体网格。

三维边界恢复问题的困难在于，仅通过拓扑变换而不插入额外的网格点（以下简称Steiner 点）无法保证能恢复给定的边界约束，典型例子如Schönhardt 多面体[2]与Chazelle 多面体[3]。因此，所有鲁棒的边界恢复算法都需要考虑如何插入Steiner 点，以及消除因Steiner 点增加引起的负面效应：1）影响网格可用性。在分区体网格生成等应用场景中，常要求输出边界与输入边界一致。然而，在边界上插入的Steiner 点会改变输入边界的拓扑结构，导致网格不可用。2）损害算法鲁棒性。Steiner 点坐标通常采用浮点数存储，其计算过程必然会伴随舍入误差，复杂输入情形下舍入误差的累计会导致算法失效。3）降低网格质量。增加的Steiner 点通常位于低质量面表面单元附近，其附近连接了一定数量的低质量、零体积单元。因边界约束限制，后续质量优化也很难大幅提升这些单元的质量。

鉴于Steiner 点的负面效应，好的边界恢复算法需要尽量减少Steiner 点插入的数量，以及优化Steiner点的插入位置。然而，Ruppert 和Seidel 证明了确定任意多面体是否可以在不加入Steiner 点的情况下进行三角化是一个NP 完全问题[4]。因此，如何同时保证边界恢复算法的时间效率也是难题之一。

完整的边界恢复流程可以分为前处理、Steiner 点插入、后处理三个方面。前处理指的是仅通过拓扑变换恢复部分甚至全部边界，后处理则为Steiner 点的移动和删除。其中前处理和后处理负责实现约束边界恢复，Steiner 点插入则保证了算法的鲁棒性和完备性。

陈建军等近期提出了一类代表性的前处理算法[5]。该算法通过引入一类可递归运算的拓扑变换操作（递归壳变换），在无需增加Steiner 点的情形下尽量多地减少边界约束与当前网格的相交频次，以最小化边界恢复所需Steiner 点的数量。

本文算法是对文献[5]算法的改进。改进算法先通过基于递归壳变换的前处理步骤恢复大部分遗失边界，再增加Steiner 点恢复所有遗失边界。与文献[5]算法相比，新算法的创新点包括：1）引入模拟退火思想，使新算法更具有全局最优性；2）实现了两类Steiner 点插入方案，减少边界上Steiner 点数量；3）优化了边界恢复流程，充分利用每一个插入的Steiner点。改进算法旨在减少Steiner 点数量和减轻Steiner点的负面效应，并提升了算法的鲁棒性。为了展现算法效果，本文利用Thingi10k 数据集[6]进行测试。

1 相关工作

1.1 边界恢复相关研究

国内外关于前处理、Steiner 点插入、后处理已有许多研究工作。

关于前处理，Radon[7]提出四种三维空间中的基本变换：1-4 变换、4-1 变换、2-3 变换、3-2 变换（如图1 所示）；Joe[8]则通过多种变换的组合，实现了更多复杂变换。Shewchuk[9]的研究指出，大多组合变换可以通过两次n-m 变换（n-to-m flip）实现，n-m 变换指的是一种通用的2-3、3-2 变换的组合。Si[10]将n-m变换加以推广，并系统性地实现在网格生成程序TetGen 中。Liu等[11]则提出了一种称为小多面体重连（small polyhedron reconnection,SPR）的方法，相比于前述变换方法找到局部最优解的效率更高，许多后继研究都对其进行了改进[12-13]，但由于该方法为NP 复杂度，因此建议输入多面体面数低于40[11]。陈建军等[5]则提出了一种递归壳变换边界恢复策略，鉴于该算法与本文的密切联系，下段先对其进行简要回顾。

图1 1-4、4-1 变换和3-2、2-3 变换示意图Fig.1 Schematic diagram of transformation 1-4,4-1,3-2 and 2-3

壳指的是共享同一条边的所有单元构成的集合，如图2[14]所示；而壳变换则是对裙多边形进行重三角化的网格局部重连技术，如图3[14]所示。针对壳变换问题，Shewchuk[9]提出了一种动态规划算法，可以高效地找到裙多边形三角化的最优解。陈建军等[5]则对该算法进行了拓展，以适应裙多边形部分三角化的情形，并提出壳变换的递归策略，大幅增加了其成功率。文献[5]进一步提出利用递归壳变换移除与边界约束相交的网格边/面，从而恢复边界的方法，并为了保证算法的收敛性与效率，引入了多个约束条件，其中包括：1）不对已经恢复的边界边执行壳变换；2）输出网格与边界相交数不大于输入网格。

图2 壳的示意图[14]Fig.2 Schematic diagram of the shell[14]

图3 利用壳变换重三角化裙多边形示意图[14]Fig.3 Schematic diagram of the shell transformation for skirt polygon re-triangulation[14]

递归壳变换算法相较于n-m 变换等一般拓扑变换算法更具侵略性，更容易找到局部最优解，且相比SPR 等局部重连方法，其优化范围更广。因此，本文基于递归壳变换，提出改进的约束边界恢复方法。

对于Steiner 点的插入步骤，其关键问题在于插入点的位置。George等[15]提出了一系列启发式规则，在内部直接加入Steiner 点以恢复边界，然而Liu 和Baida[16]证明该算法在理论上不收敛。另一类思路则基于Shewchuk[17]提出的约束Delaunay 四面体化（constrained Delaunay tetrahedralization）算法，这类算法一般在边界边或边界面的中点附近加入Steiner点，并理论保证了恢复后网格的质量，但缺点在于一般会加入过多Steiner 点。Weathrill 和Hassan[18]最早提出了在遗失边界和当前四面体网格交点处插入Steiner 点的思路，作者所在课题组也针对该算法进行过研究，验证了其收敛性，并实现了一套完备的基于交点处插点的一致边界恢复算法[12,19-20]，但该思路缺点在于常产生低质量单元，并且易受浮点误差影响导致算法失效。针对Steiner 点的数量和位置问题，Goerigk和Si 对一系列典型多面体进行了系统性研究[21-23]，但其成果尚未得到广泛的实际工程应用。

关于后处理，目前一般采用的是“点分解”策略。该策略的收敛性已有理论证明，George等[24]、Du 和Wang[25]、Si[10]、陈建军[12,19]等都研究或实现过该算法。研究表明，“点分解”策略会引入“薄元”等单元质量较低的局部[26]，这些差单元需要通过后续的Steiner 点删除、网格质量优化等算法加以解决，但这些算法由于受到已恢复的边界约束限制，效果常无法达到预期。此外，后处理也无法完全消除前述Steiner 点产生的负面效应。为此，本文也旨在提升前处理的边界恢复效果，从而减少后处理的使用。

1.2 模拟退火算法

本文引入模拟退火算法，对递归壳变换算法[5]加以改进。

模拟退火（simulated annealing,SA）算法是一种全局优化方法，由Metropolis等[27]提出。该算法受到金属物质退火时的统计规律启发。研究表明[28]，在温度T时，粒子处于状态x的概率满足玻尔兹曼分布：

其中：E(x) 为状态x的能量；kB为玻尔兹曼常数；Z(T)为概率归一化因子。

SA 算法的思想是：将优化目标视作E(x)；温度T由高至低迭代，每一T下，通过状态转移使得稳态分布趋近于Px,T(x) ；T=0 时，E(x)达到最优解。

从任意初始状态出发，要使最终状态满足P(x)分布，转移过程需要满足如下的马尔科夫链细致平衡方程。其中x、x′分别为当前状态和新状态，f为状态转移概率函数，a为接受概率函数。

在f(x→x′)=f(x′→x),∀x,x′条件下，Metropolis等[27]提出了以下接受准则，称为Metropolis 准则：

SA 算法具有在高温时倾向于跳出局部最优，而在低温时收敛于局部最优的性质，因此该方法是改良前述递归壳变换算法的有力工具。

2 本文方法

2.1 总体流程

本文涉及的Delaunay 四面体网格生成流程如图4 所示。首先，基于输入网格的外包围盒构建初始Delaunay 四面体网格；然后，将输入网格点采用Bowyer-Watson[29-30]算法插入四面体；其次，执行本文提出的改进的边界边/面恢复算法，恢复所有输入边界约束；最后，删除外部点和外部单元，输出有效四面体网格。四面体网格生成的其他常见流程，诸如Steiner 点移除、网格加密、网格质量优化等，不在本文赘述。

图4 网格生成流程图Fig.4 Flowchart of mesh generation process

边界边恢复整体流程如图5 所示。流程表现为“两层+两轮”的迭代过程，其中“两层”指的是外层的边恢复迭代以及内层的Steiner 点插入迭代，“两轮”指的是Steiner 点插入迭代分两轮进行。

图5 边界边恢复整体流程图Fig.5 Overall flowchart of boundary edge recovery

流程说明如下：边恢复迭代负责恢复当前的所有边界边，每次迭代首先执行2.2 节提出的改进约束边界恢复策略，以减少整体的遗失边界数量；对于无法约束恢复的边，再分两轮迭代插入Steiner 点，每轮遍历所有遗失边界，第一轮只采用体内Steiner 点插入策略，第二轮再开启边界Steiner 点插入，以减少Steiner 点对于输入边界的破坏；在每次插点恢复一边前，考虑到其他边插入的Steiner 点可能为该边的约束恢复创造有利条件，为了充分降低Steiner 点数量，尝试使用递归壳变换策略[5]恢复该边；边界插点后，需要更新边界拓扑，以维护新产生的子边信息，子边由下一次边恢复迭代进行恢复。

边界面恢复与边界边恢复流程类似，但更加简单。实践中注意到恢复边界边后遗失面的数量一般较少，因此可以直接使用递归壳变换策略[5]与面心加点方法进行恢复。

2.2 改进的约束边界恢复算法

本节利用模拟退火算法，对文献[5]提出的递归壳变换算法进行优化，提出一种改进的约束边界恢复算法，并介绍其中多个关键方法。

2.2.1 结合模拟退火的约束边界恢复算法

给定点集S∈R2，当三维单纯复形 T满足以下条件，称 T 为S的四面体化[31]：1）T 的点集等于S；2）T的覆盖空间等于S的凸包 conv(S)。

本文将S的四面体化 T视为状态，状态空间即为所有四面体化的集合 Ω，将 T与输入边界边集合E的交点数量N(T) 作为状态 T的能量。本节所描述的算法以边界边恢复为例，其目标是从Delaunay 四面体化 D ∈Ω 出发，找到 T ∈Ω，使得N(T)尽可能低。该优化目标与1.1 节所述递归壳变换边界恢复算法的约束条件2 思路一致。

本文提出一种结合模拟退火的约束边界恢复方法。算法流程如下所示：

该算法的说明如下：

先求出四面体化 T 与边界边集合Eg的所有相交实体insts。简要方法为，遍历Eg中每一边e，由e其一端点作为起始位置在 T中开展深度优先搜索，并将路径中的相交实体加入insts，直到搜索至另一端点。

模拟退火过程，由预设的初始温度 T0起始，在固定温度下执行多次“取相交实体—递归壳变换”迭代，至交点数N(T)收敛（判断方式为：最后50 次迭代与倒数第51 至100 次迭代的交点数量的平均值之差绝对值 ≤0.5 ），再降低温度。当温度低于阈值Tthreshold（本文取0.1），结束退火过程。

每次迭代从insts随机选择一个实体entity，利用模拟退火改进的递归壳变换例程RecursiveShell-TransformSA 对其进行移除操作。该例程与文献[5]的recursiveST 例程大致类似，即壳变换的递归执行，单次壳变换步骤为：先利用Shewchuk[9]提出的动态规划算法，产生最优三角化子问题的解矩阵；再将矩阵视作有向图，遍历回路找到最优解 Tnew；最后根据1.1 节所述若干约束条件比较 Tnew与原四面体化Told，决定是否接收 Tnew。但本文例程有两点不同：1）允许已经恢复的边界被壳变换；2）每次产生并比较可行状态 T′时（包括子问题解间的比较、回路间的比较以及最优解与原四面体化间的比较），本算法首先通过2.2.2 节所述方法计算交点数N(T′)，这里不再严格通过比较交点数大小来判定解的优劣，而是通过式（4）的Metropolis准则计算出接收概率 0，再比较0-1 随机数r是否小于 0 来决定是否转移至 T′。

式（4）中的 T为退火温度，随迭代过程逐渐降温。对此本文采用了非常快速模拟退火（very fast simulated annealing,VFSA）[32]算法的降温计划：

其中：k为迭代次数；α与 T0为预设参数，后文实验中将给出参考默认值，也可根据具体案例自行设定。

对网格进行拓扑变换后，需要更新相交实体数据insts。该过程与第一步求insts 的方法类似，但只遍历那些相交单元被拓扑变换修改的边界边。

考虑到递归壳变换状态转移过程的复杂性，本文忽略了前述Metropolis 准则的前置条件f(T →T′)=f(T′→T)，并假设其成立。

改进的约束边界恢复算法的总体流程如算法2所示。首先，执行一次常规的递归壳变换[5]例程，使得网格达到局部最优解；然后，从局部最优出发，再进行结合SA 的边界恢复例程，可加快收敛速度；由于SA 算法可能迭代次数不足以收敛到最优解，最后再执行一次递归壳变换，以达到局部最优状态。

2.2.2 交点数的计算与查询

要计算任意四面体化 T′与边界约束的交点数N(T′)，本文方法是在拓扑变换过程中进行增量维护。其中的关键算法是获取任意网格边/面与所有边界约束的交点数量。为此，本文利用了多个关键方法保证其鲁棒、高效地执行。

该算法外层是对于边界约束的遍历，内层则是“网格实体—边界实体求交”算法。首先，程序引入了AABB（axis-aligned bounding boxes）树[33]对外层遍历进行加速。具体方法是：在初始化阶段中，将所有边界约束实体加入AABB 树，再在求交时通过AABB树快速过滤不相交的边界约束。在包围盒重叠程度较低的情况下，该方法可将遍历复杂度由 O(M)降为O(logM)，其中M为边界实体数量。

其次，为了避免浮点误差导致的求交算法错误，本算法利用了Guigue等[34]实现的“三角形-三角形相交”检测程序，其中包含几何精确的“线面相交”函数。

最后，使用了一个哈希表存储已计算出网格实体与边界的交点数，以减少二次计算带来的性能损耗。图6 展示了网格面交点数哈希表的C++数据结构定义，查询逻辑伪代码如算法 3 所示。哈希表利用二维变长数组储存，利用其内存局部性，可加速哈希碰撞时的搜索时间；对于任意网格面p1p2p3，本算法在HashSort 函数中对其进行结果固定的乱序排序，取首位为哈希值，后两位为哈希碰撞时的识别标识；并利用了C++Union 关键字，将两位的标识等价为单个长整型，将两次相等判断缩减到了一次，进一步加快了查询时间。

图6 三角面哈希表数据结构Fig.6 Hash table data structure for triangular polygons

2.3 Steiner 点插入策略

仅通过拓扑变换不能解决所有边界恢复问题，此时需要考虑插入Steiner 点。由于体内Steiner 点不改变输入边界约束，首先尝试在四面体内部插点，进一步地，再利用插点打断边界的方法，以恢复所有边界。

本文受到TetGen 研究成果[35]的启发，实现了一种针对特殊情形的体内Steiner 点插入方法，图7 展示了其过程示例。所谓特殊情形，指的是遗失边ab与四面体网格 T的所有相交单元 Tab具有共边cd。换而言之，Tab由边 Δd所在壳上的连续单元p0p1cd,p1p2cd,...,pn-1pncd构成，如图7（a）所示，pi(i=0,1,···,n) 为壳的裙点，且p0=a、pn=b。方法具体步骤如下：

图7 体内Steiner 点插入原理示例图Fig.7 Principle illustration of the interior Steiner point insertion

1）将 Tab的覆盖空间记作空腔C，首先从C中分离四面体abcd，得到空腔C′，如图7（b）所示。

2）从ab中点出发，通过点光滑化例程得到点 6，使其满足与C′的所有表面三角形所构成的四面体具有正体积。

3）若上述点e不存在，则失败；否则将点e与C′的所有表面三角形相连接，得到C′的四面体化，如图7（c）所示。

由此通过插入点 6，得到了空腔C的一种四面体化，其包含所需恢复的遗失边ab。

针对体内插点无法解决的情况，采用了Si 和Gärtner[36]提出的一种保质Steiner 点插入策略对边界边进行分裂。该方法根据遗失边端点周围的面网格性质，分三类情况计算Steiner 点的位置，并且在计算时额外考虑了遗失边周围网格点的分布以及端点的局部特征尺寸[37]。

对于遗失边界面，则考虑直接往面心加点。

3 实验及结果分析

所有实验均在AMD Ryzen 7 5800H、3.2GHz CPU、16GB RAM 平台上进行。测试集取自Thingi10k 数据集[6]，该数据集中包含1 万个三角网格模型，其中多数样本存在自相交、重节点等问题。鉴于本文重点非面网格修复，实际测试中采用了几何正确的面网格的4 303 例。

3.1 模拟退火参数的评估与选取

就模拟退火算法而言，参数的选择至关重要。本文方法一共有三个待定参数，列于表1。实验前首先为每个待定参数设计了多个候选值；再对所有候选值进行组合，共产生210 种策略；最后从4 303 例样本中随机抽取200 例子样本集，对每一策略进行测试，并统计其平均耗时与平均Steiner 点数量。测试结果见图8。

表1 模拟退火算法待定参数Table 1 Pending parameters for the simulated annealing algorithm

图8 210 种策略测试结果统计图Fig.8 Statistical chart of 210 strategy test results

不同参数下的测试结果表明，要使算法效果更佳，即Steiner 点数量更低，一般要牺牲时间效率，反之亦然。为了在尽可能提升时间性能同时保证算法效果，本文选取了图8 中箭头所指的参数组合，分别为lmax=7,T0=1.1,α=0.86。为了表现本文方法的普适性，后文实验过程均固定取上述参数值，不再手动调整。但在实际算例中，用户仍可以根据具体情况调整参数。

3.2 与其他方法的比较

本节将本文方法与文献[5]提出的边界恢复程序以及TetGen1.6 版本（下文以TetGen 指代）进行比较，分别统计各程序在边界恢复中插入的Steiner 点数量。为了公平性，关闭了所有程序的面网格优化、Steiner 点后处理功能。

所有4 303 个样本的运行结果分类情况如表2 所示。对比文献[5]的程序，本文方法的无Steiner 点边界恢复数量大幅提高。其程序另有45 例网格生成失败的样本，经研究，所有失败样本的原因在于使用了直接在网格与遗失边界交点处插入Steiner 点的方法[20]。当面网格质量较差，局部Steiner 点数量较多时，该方法易受浮点误差累积影响从而导致程序崩溃。图9展示了该程序报错的代号238 420 的样本，该样本在局部聚集有较多最小角极低的针状单元。本文所实现的程序由于降低了Steiner 点数量，并且采用了保质的Steiner 点插入方法，可以处理诸如此类的例子。后续实验将剔除上述失败的45 例样本，使用剩余的4 258 例。

表2 4 303 例样本的运行结果统计数据Table 2 Statistical data of running results with 4 303 samples

图9 文献[5]提出方法处理的面网格结果Fig.9 Surface mesh processed by the method in Ref.[5]

对于2 829 例TetGen 生成的Steiner 点数量大于0 的样本，使用本文方法能够将其中的495 例转化为零Steiner 点恢复，转化率达到17.50%。另外，在2 402例本文程序含Steiner 点的样本中，TetGen 仅68 例零Steiner 点恢复。说明本文方法较TetGen 也具有更高的零Steiner 点恢复率。

后文将进一步将本文方法与另两程序进行横向对比。为了更加直观地比较各程序的Steiner 点数量，在4 258 例上述程序都能成功生成的样本中，剔除了1 259例所有程序都无Steiner 点的样本，其余2 999 例的Steiner 点数量分布如图10 所示，统计数据如表3 所示。从数据上看，本文方法较TetGen 将2 999 例样本的合计Steiner 点数量降低了29.86%，并且生成结果中Steiner 点数量≥5 的样本数量也最少，说明本文方法总体效果显著。另外注意到文献[5]方法的Steiner 点数量偏多，原因在于该程序采用了单次“递归壳变换边界恢复—插入Steiner 点”的线性流程，导致加入冗余Steiner 点。而本文方法通过拓扑变换与插点步骤的穿插执行，即“两层+两轮”迭代，充分利用每个Steiner 点，有效降低了Steiner 点数量。

表3 2 999 例样本Steiner 点数量统计数据Table 3 Statistical data for the Steiner point number in 2 999 samples

图10 2 999 例样本Steiner 点数量分布图Fig.10 Distribution of the Steiner point number in 2999 samples

图11 展示了2 999 例样本中的3 例典型网格。表4 展示了其Steiner 点数据。由此可见，本文算法能以较少Steiner 点数量处理不同形态的输入网格。

表4 3 例典型样本的Steiner 点数量Table 4 Steiner point numbers for 3 typical samples

图11 3 例典型样本网格Fig.11 Three typical test meshes

表5 展示了上述程序都能通过的4 258 例样本的边界恢复耗时统计数据。数据表明，针对Thingi10k数据集，本文方法和文献[5]提出的方法均与TetGen有一个数量级左右的性能差距。这是由于递归壳变换本身复杂度高，并且对于大规模的复杂网格，较大的递归深度带来了更高的耗时；另一方面，模拟退火算法需要通过降温步骤从局部最优向全局最优过渡，为本文方法带来了额外耗时。

表5 4 258 例样本耗时统计数据Table 5 Time consumption statistics of 4 258 samples

相较于其他程序，本文方法以一定时间为代价，提高了零Steiner 点恢复成功率，并降低Steiner 点数量。尤其相对于文献[5]提出的方法，本文方法效果提升显著，并且只增加了少量耗时。上述提升对于分区体网格生成等需要尽可能降低Steiner 点数量的应用场景具有重要意义，并且为复杂模型的约束边界恢复创造了新的可能。

上述时间性能差距也可以通过多种方法优化，如利用更加简单的拓扑变换方法预先处理网格、优化模拟退火策略与参数等。文献[5]指出，在实际网格生成应用中，边界恢复步骤往往只占据较小一部分耗时，并在其文章中展示了一个400 万单元的四面体网格生成案例，其边界恢复耗时只占0.23%。

为了评估本文方法耗时对于实际网格生成耗时的影响程度，使用伦敦塔桥模型（如图12 所示）开展进一步实验。实验中，首先开启TetGen 网格优化功能（参数为“ -q1.1/10”），测试其完整的网格生成流程，最终输出体网格单元数达到千万量级，基本达到工程标准。再利用本文方法测试该模型的边界恢复过程，实验数据如表6 所示。数据表明，理想情况下，本例中若以本文方法替代TetGen1.6 原有边界恢复算法，能起到降低Steiner 点数量的作用，同时网格生成总耗时仅增加4.86%，说明本文方法具有实用价值。

表6 伦敦塔桥模型网格生成统计数据Table 6 Mesh generation statistics for the Tower Bridge model

图12 伦敦塔桥模型面网格Fig.12 Surface mesh of Tower Bridge model

4 结论

本文提出了一种改进的Delaunay 四面体网格生成的边界恢复流程。该流程主要由改进的约束边界恢复步骤以及体内、边界Steiner 点插入步骤构成，表现为“两层+两轮”的迭代流程，旨在降低Steiner 点数量并减少Steiner 点的负面效应。其中改进的约束边界恢复方法基于递归壳变换算法，并结合模拟退火算法，利用Metropolis 准则使递归壳变换脱离局部最优解，进一步减少遗失边界数量。同时，针对“网格实体—边界实体求交”的效率与鲁棒性问题，引入和提出了多个关键方法进行优化。对于拓扑变换无法恢复的边界，实现了体内、边界两种插入Steiner 点方法。

Thingi10k 数据集的测试结果表明，本文方法相较于文献[5]提出的程序以及TetGen1.6，可以零Steiner点恢复更多例子，并且普遍降低了边界恢复所需的Steiner 点数量。这对于复杂例子的约束边界恢复、提升程序鲁棒性以及网格质量有着重要意义。

本文方法在时间性能上尚有不足。后续计划就拓扑变换方法、模拟退火策略与参数等方面开展优化工作。