曾寅家，陳建軍，傅珂杰

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

0 引言

Delaunay 三角化（delaunay triangulation,DT）是生成四面體網格的主流方法之一，具有運算快速、網格總體質量較高等優(yōu)點。給定表面網格，DT 方法通過不斷插入網格點形成一個由表面網格點和內部點組成的初始四面體網格，該網格滿足Delaunay 準則[1]，但不能保證所有邊界約束（在這指輸入的表面網格邊和網格面）都存在于其中，因此需要調用邊界恢復算法修改初始體網格，獲得與邊界約束兼容一致的體網格。

三維邊界恢復問題的困難在于，僅通過拓撲變換而不插入額外的網格點（以下簡稱Steiner 點）無法保證能恢復給定的邊界約束，典型例子如Sch?nhardt 多面體[2]與Chazelle 多面體[3]。因此，所有魯棒的邊界恢復算法都需要考慮如何插入Steiner 點，以及消除因Steiner 點增加引起的負面效應：1）影響網格可用性。在分區(qū)體網格生成等應用場景中，常要求輸出邊界與輸入邊界一致。然而，在邊界上插入的Steiner 點會改變輸入邊界的拓撲結構，導致網格不可用。2）損害算法魯棒性。Steiner 點坐標通常采用浮點數存儲，其計算過程必然會伴隨舍入誤差，復雜輸入情形下舍入誤差的累計會導致算法失效。3）降低網格質量。增加的Steiner 點通常位于低質量面表面單元附近，其附近連接了一定數量的低質量、零體積單元。因邊界約束限制，后續(xù)質量優(yōu)化也很難大幅提升這些單元的質量。

鑒于Steiner 點的負面效應，好的邊界恢復算法需要盡量減少Steiner 點插入的數量，以及優(yōu)化Steiner點的插入位置。然而，Ruppert 和Seidel 證明了確定任意多面體是否可以在不加入Steiner 點的情況下進行三角化是一個NP 完全問題[4]。因此，如何同時保證邊界恢復算法的時間效率也是難題之一。

完整的邊界恢復流程可以分為前處理、Steiner 點插入、后處理三個方面。前處理指的是僅通過拓撲變換恢復部分甚至全部邊界，后處理則為Steiner 點的移動和刪除。其中前處理和后處理負責實現約束邊界恢復，Steiner 點插入則保證了算法的魯棒性和完備性。

陳建軍等近期提出了一類代表性的前處理算法[5]。該算法通過引入一類可遞歸運算的拓撲變換操作（遞歸殼變換），在無需增加Steiner 點的情形下盡量多地減少邊界約束與當前網格的相交頻次，以最小化邊界恢復所需Steiner 點的數量。

本文算法是對文獻[5]算法的改進。改進算法先通過基于遞歸殼變換的前處理步驟恢復大部分遺失邊界，再增加Steiner 點恢復所有遺失邊界。與文獻[5]算法相比，新算法的創(chuàng)新點包括：1）引入模擬退火思想，使新算法更具有全局最優(yōu)性；2）實現了兩類Steiner 點插入方案，減少邊界上Steiner 點數量；3）優(yōu)化了邊界恢復流程，充分利用每一個插入的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變換加以推廣，并系統(tǒng)性地實現在網格生成程序TetGen 中。Liu等[11]則提出了一種稱為小多面體重連（small polyhedron reconnection,SPR）的方法，相比于前述變換方法找到局部最優(yōu)解的效率更高，許多后繼研究都對其進行了改進[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]提出了一種動態(tài)規(guī)劃算法，可以高效地找到裙多邊形三角化的最優(yōu)解。陳建軍等[5]則對該算法進行了拓展，以適應裙多邊形部分三角化的情形，并提出殼變換的遞歸策略，大幅增加了其成功率。文獻[5]進一步提出利用遞歸殼變換移除與邊界約束相交的網格邊/面，從而恢復邊界的方法，并為了保證算法的收斂性與效率，引入了多個約束條件，其中包括：1）不對已經恢復的邊界邊執(zhí)行殼變換；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 變換等一般拓撲變換算法更具侵略性，更容易找到局部最優(yōu)解，且相比SPR 等局部重連方法，其優(yōu)化范圍更廣。因此，本文基于遞歸殼變換，提出改進的約束邊界恢復方法。

對于Steiner 點的插入步驟，其關鍵問題在于插入點的位置。George等[15]提出了一系列啟發(fā)式規(guī)則，在內部直接加入Steiner 點以恢復邊界，然而Liu 和Baida[16]證明該算法在理論上不收斂。另一類思路則基于Shewchuk[17]提出的約束Delaunay 四面體化（constrained Delaunay tetrahedralization）算法，這類算法一般在邊界邊或邊界面的中點附近加入Steiner點，并理論保證了恢復后網格的質量，但缺點在于一般會加入過多Steiner 點。Weathrill 和Hassan[18]最早提出了在遺失邊界和當前四面體網格交點處插入Steiner 點的思路，作者所在課題組也針對該算法進行過研究，驗證了其收斂性，并實現了一套完備的基于交點處插點的一致邊界恢復算法[12,19-20]，但該思路缺點在于常產生低質量單元，并且易受浮點誤差影響導致算法失效。針對Steiner 點的數量和位置問題，Goerigk和Si 對一系列典型多面體進行了系統(tǒng)性研究[21-23]，但其成果尚未得到廣泛的實際工程應用。

關于后處理，目前一般采用的是“點分解”策略。該策略的收斂性已有理論證明，George等[24]、Du 和Wang[25]、Si[10]、陳建軍[12,19]等都研究或實現過該算法。研究表明，“點分解”策略會引入“薄元”等單元質量較低的局部[26]，這些差單元需要通過后續(xù)的Steiner 點刪除、網格質量優(yōu)化等算法加以解決，但這些算法由于受到已恢復的邊界約束限制，效果常無法達到預期。此外，后處理也無法完全消除前述Steiner 點產生的負面效應。為此，本文也旨在提升前處理的邊界恢復效果，從而減少后處理的使用。

1.2 模擬退火算法

本文引入模擬退火算法，對遞歸殼變換算法[5]加以改進。

模擬退火（simulated annealing,SA）算法是一種全局優(yōu)化方法，由Metropolis等[27]提出。該算法受到金屬物質退火時的統(tǒng)計規(guī)律啟發(fā)。研究表明[28]，在溫度T時，粒子處于狀態(tài)x的概率滿足玻爾茲曼分布：

其中：E(x) 為狀態(tài)x的能量；kB為玻爾茲曼常數；Z(T)為概率歸一化因子。

SA 算法的思想是：將優(yōu)化目標視作E(x)；溫度T由高至低迭代，每一T下，通過狀態(tài)轉移使得穩(wěn)態(tài)分布趨近于Px,T(x) ；T=0 時，E(x)達到最優(yōu)解。

從任意初始狀態(tài)出發(fā)，要使最終狀態(tài)滿足P(x)分布，轉移過程需要滿足如下的馬爾科夫鏈細致平衡方程。其中x、x′分別為當前狀態(tài)和新狀態(tài)，f為狀態(tài)轉移概率函數，a為接受概率函數。

在f(x→x′)=f(x′→x),?x,x′條件下，Metropolis等[27]提出了以下接受準則，稱為Metropolis 準則：

SA 算法具有在高溫時傾向于跳出局部最優(yōu)，而在低溫時收斂于局部最優(yōu)的性質，因此該方法是改良前述遞歸殼變換算法的有力工具。

2 本文方法

2.1 總體流程

本文涉及的Delaunay 四面體網格生成流程如圖4 所示。首先，基于輸入網格的外包圍盒構建初始Delaunay 四面體網格；然后，將輸入網格點采用Bowyer-Watson[29-30]算法插入四面體；其次，執(zhí)行本文提出的改進的邊界邊/面恢復算法，恢復所有輸入邊界約束；最后，刪除外部點和外部單元，輸出有效四面體網格。四面體網格生成的其他常見流程，諸如Steiner 點移除、網格加密、網格質量優(yōu)化等，不在本文贅述。

圖4 網格生成流程圖Fig.4 Flowchart of mesh generation process

邊界邊恢復整體流程如圖5 所示。流程表現為“兩層+兩輪”的迭代過程，其中“兩層”指的是外層的邊恢復迭代以及內層的Steiner 點插入迭代，“兩輪”指的是Steiner 點插入迭代分兩輪進行。

圖5 邊界邊恢復整體流程圖Fig.5 Overall flowchart of boundary edge recovery

流程說明如下：邊恢復迭代負責恢復當前的所有邊界邊，每次迭代首先執(zhí)行2.2 節(jié)提出的改進約束邊界恢復策略，以減少整體的遺失邊界數量；對于無法約束恢復的邊，再分兩輪迭代插入Steiner 點，每輪遍歷所有遺失邊界，第一輪只采用體內Steiner 點插入策略，第二輪再開啟邊界Steiner 點插入，以減少Steiner 點對于輸入邊界的破壞；在每次插點恢復一邊前，考慮到其他邊插入的Steiner 點可能為該邊的約束恢復創(chuàng)造有利條件，為了充分降低Steiner 點數量，嘗試使用遞歸殼變換策略[5]恢復該邊；邊界插點后，需要更新邊界拓撲，以維護新產生的子邊信息，子邊由下一次邊恢復迭代進行恢復。

邊界面恢復與邊界邊恢復流程類似，但更加簡單。實踐中注意到恢復邊界邊后遺失面的數量一般較少，因此可以直接使用遞歸殼變換策略[5]與面心加點方法進行恢復。

2.2 改進的約束邊界恢復算法

本節(jié)利用模擬退火算法，對文獻[5]提出的遞歸殼變換算法進行優(yōu)化，提出一種改進的約束邊界恢復算法，并介紹其中多個關鍵方法。

2.2.1 結合模擬退火的約束邊界恢復算法

給定點集S∈R2，當三維單純復形 T滿足以下條件，稱 T 為S的四面體化[31]：1）T 的點集等于S；2）T的覆蓋空間等于S的凸包 conv(S)。

本文將S的四面體化 T視為狀態(tài)，狀態(tài)空間即為所有四面體化的集合 ?，將 T與輸入邊界邊集合E的交點數量N(T) 作為狀態(tài) T的能量。本節(jié)所描述的算法以邊界邊恢復為例，其目標是從Delaunay 四面體化 D ∈? 出發(fā)，找到 T ∈?，使得N(T)盡可能低。該優(yōu)化目標與1.1 節(jié)所述遞歸殼變換邊界恢復算法的約束條件2 思路一致。

本文提出一種結合模擬退火的約束邊界恢復方法。算法流程如下所示：

該算法的說明如下：

先求出四面體化 T 與邊界邊集合Eg的所有相交實體insts。簡要方法為，遍歷Eg中每一邊e，由e其一端點作為起始位置在 T中開展深度優(yōu)先搜索，并將路徑中的相交實體加入insts，直到搜索至另一端點。

模擬退火過程，由預設的初始溫度 T0起始，在固定溫度下執(zhí)行多次“取相交實體—遞歸殼變換”迭代，至交點數N(T)收斂（判斷方式為：最后50 次迭代與倒數第51 至100 次迭代的交點數量的平均值之差絕對值 ≤0.5 ），再降低溫度。當溫度低于閾值Tthreshold（本文取0.1），結束退火過程。

每次迭代從insts隨機選擇一個實體entity，利用模擬退火改進的遞歸殼變換例程RecursiveShell-TransformSA 對其進行移除操作。該例程與文獻[5]的recursiveST 例程大致類似，即殼變換的遞歸執(zhí)行，單次殼變換步驟為：先利用Shewchuk[9]提出的動態(tài)規(guī)劃算法，產生最優(yōu)三角化子問題的解矩陣；再將矩陣視作有向圖，遍歷回路找到最優(yōu)解 Tnew；最后根據1.1 節(jié)所述若干約束條件比較 Tnew與原四面體化Told，決定是否接收 Tnew。但本文例程有兩點不同：1）允許已經恢復的邊界被殼變換；2）每次產生并比較可行狀態(tài) T′時（包括子問題解間的比較、回路間的比較以及最優(yōu)解與原四面體化間的比較），本算法首先通過2.2.2 節(jié)所述方法計算交點數N(T′)，這里不再嚴格通過比較交點數大小來判定解的優(yōu)劣，而是通過式（4）的Metropolis準則計算出接收概率 0，再比較0-1 隨機數r是否小于 0 來決定是否轉移至 T′。

式（4）中的 T為退火溫度，隨迭代過程逐漸降溫。對此本文采用了非?？焖倌M退火（very fast simulated annealing,VFSA）[32]算法的降溫計劃：

其中：k為迭代次數；α與 T0為預設參數，后文實驗中將給出參考默認值，也可根據具體案例自行設定。

對網格進行拓撲變換后，需要更新相交實體數據insts。該過程與第一步求insts 的方法類似，但只遍歷那些相交單元被拓撲變換修改的邊界邊。

考慮到遞歸殼變換狀態(tài)轉移過程的復雜性，本文忽略了前述Metropolis 準則的前置條件f(T →T′)=f(T′→T)，并假設其成立。

改進的約束邊界恢復算法的總體流程如算法2所示。首先，執(zhí)行一次常規(guī)的遞歸殼變換[5]例程，使得網格達到局部最優(yōu)解；然后，從局部最優(yōu)出發(fā)，再進行結合SA 的邊界恢復例程，可加快收斂速度；由于SA 算法可能迭代次數不足以收斂到最優(yōu)解，最后再執(zhí)行一次遞歸殼變換，以達到局部最優(yōu)狀態(tài)。

2.2.2 交點數的計算與查詢

要計算任意四面體化 T′與邊界約束的交點數N(T′)，本文方法是在拓撲變換過程中進行增量維護。其中的關鍵算法是獲取任意網格邊/面與所有邊界約束的交點數量。為此，本文利用了多個關鍵方法保證其魯棒、高效地執(zhí)行。

該算法外層是對于邊界約束的遍歷，內層則是“網格實體—邊界實體求交”算法。首先，程序引入了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]的啟發(fā)，實現了一種針對特殊情形的體內Steiner 點插入方法，圖7 展示了其過程示例。所謂特殊情形，指的是遺失邊ab與四面體網格 T的所有相交單元 Tab具有共邊cd。換而言之，Tab由邊 Δd所在殼上的連續(xù)單元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中點出發(fā)，通過點光滑化例程得到點 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 萬個三角網格模型，其中多數樣本存在自相交、重節(jié)點等問題。鑒于本文重點非面網格修復，實際測試中采用了幾何正確的面網格的4 303 例。

3.1 模擬退火參數的評估與選取

就模擬退火算法而言，參數的選擇至關重要。本文方法一共有三個待定參數，列于表1。實驗前首先為每個待定參數設計了多個候選值；再對所有候選值進行組合，共產生210 種策略；最后從4 303 例樣本中隨機抽取200 例子樣本集，對每一策略進行測試，并統(tǒng)計其平均耗時與平均Steiner 點數量。測試結果見圖8。

表1 模擬退火算法待定參數Table 1 Pending parameters for the simulated annealing algorithm

圖8 210 種策略測試結果統(tǒng)計圖Fig.8 Statistical chart of 210 strategy test results

不同參數下的測試結果表明，要使算法效果更佳，即Steiner 點數量更低，一般要犧牲時間效率，反之亦然。為了在盡可能提升時間性能同時保證算法效果，本文選取了圖8 中箭頭所指的參數組合，分別為lmax=7,T0=1.1,α=0.86。為了表現本文方法的普適性，后文實驗過程均固定取上述參數值，不再手動調整。但在實際算例中，用戶仍可以根據具體情況調整參數。

3.2 與其他方法的比較

本節(jié)將本文方法與文獻[5]提出的邊界恢復程序以及TetGen1.6 版本（下文以TetGen 指代）進行比較，分別統(tǒng)計各程序在邊界恢復中插入的Steiner 點數量。為了公平性，關閉了所有程序的面網格優(yōu)化、Steiner 點后處理功能。

所有4 303 個樣本的運行結果分類情況如表2 所示。對比文獻[5]的程序，本文方法的無Steiner 點邊界恢復數量大幅提高。其程序另有45 例網格生成失敗的樣本，經研究，所有失敗樣本的原因在于使用了直接在網格與遺失邊界交點處插入Steiner 點的方法[20]。當面網格質量較差，局部Steiner 點數量較多時，該方法易受浮點誤差累積影響從而導致程序崩潰。圖9展示了該程序報錯的代號238 420 的樣本，該樣本在局部聚集有較多最小角極低的針狀單元。本文所實現的程序由于降低了Steiner 點數量，并且采用了保質的Steiner 點插入方法，可以處理諸如此類的例子。后續(xù)實驗將剔除上述失敗的45 例樣本，使用剩余的4 258 例。

表2 4 303 例樣本的運行結果統(tǒng)計數據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 所示，統(tǒng)計數據如表3 所示。從數據上看，本文方法較TetGen 將2 999 例樣本的合計Steiner 點數量降低了29.86%，并且生成結果中Steiner 點數量≥5 的樣本數量也最少，說明本文方法總體效果顯著。另外注意到文獻[5]方法的Steiner 點數量偏多，原因在于該程序采用了單次“遞歸殼變換邊界恢復—插入Steiner 點”的線性流程，導致加入冗余Steiner 點。而本文方法通過拓撲變換與插點步驟的穿插執(zhí)行，即“兩層+兩輪”迭代，充分利用每個Steiner 點，有效降低了Steiner 點數量。

表3 2 999 例樣本Steiner 點數量統(tǒng)計數據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 點數量處理不同形態(tài)的輸入網格。

表4 3 例典型樣本的Steiner 點數量Table 4 Steiner point numbers for 3 typical samples

圖11 3 例典型樣本網格Fig.11 Three typical test meshes

表5 展示了上述程序都能通過的4 258 例樣本的邊界恢復耗時統(tǒng)計數據。數據表明，針對Thingi10k數據集，本文方法和文獻[5]提出的方法均與TetGen有一個數量級左右的性能差距。這是由于遞歸殼變換本身復雜度高，并且對于大規(guī)模的復雜網格，較大的遞歸深度帶來了更高的耗時；另一方面，模擬退火算法需要通過降溫步驟從局部最優(yōu)向全局最優(yōu)過渡，為本文方法帶來了額外耗時。

表5 4 258 例樣本耗時統(tǒng)計數據Table 5 Time consumption statistics of 4 258 samples

相較于其他程序，本文方法以一定時間為代價，提高了零Steiner 點恢復成功率，并降低Steiner 點數量。尤其相對于文獻[5]提出的方法，本文方法效果提升顯著，并且只增加了少量耗時。上述提升對于分區(qū)體網格生成等需要盡可能降低Steiner 點數量的應用場景具有重要意義，并且為復雜模型的約束邊界恢復創(chuàng)造了新的可能。

上述時間性能差距也可以通過多種方法優(yōu)化，如利用更加簡單的拓撲變換方法預先處理網格、優(yōu)化模擬退火策略與參數等。文獻[5]指出，在實際網格生成應用中，邊界恢復步驟往往只占據較小一部分耗時，并在其文章中展示了一個400 萬單元的四面體網格生成案例，其邊界恢復耗時只占0.23%。

為了評估本文方法耗時對于實際網格生成耗時的影響程度，使用倫敦塔橋模型（如圖12 所示）開展進一步實驗。實驗中，首先開啟TetGen 網格優(yōu)化功能（參數為“ -q1.1/10”），測試其完整的網格生成流程，最終輸出體網格單元數達到千萬量級，基本達到工程標準。再利用本文方法測試該模型的邊界恢復過程，實驗數據如表6 所示。數據表明，理想情況下，本例中若以本文方法替代TetGen1.6 原有邊界恢復算法，能起到降低Steiner 點數量的作用，同時網格生成總耗時僅增加4.86%，說明本文方法具有實用價值。

表6 倫敦塔橋模型網格生成統(tǒng)計數據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 準則使遞歸殼變換脫離局部最優(yōu)解，進一步減少遺失邊界數量。同時，針對“網格實體—邊界實體求交”的效率與魯棒性問題，引入和提出了多個關鍵方法進行優(yōu)化。對于拓撲變換無法恢復的邊界，實現了體內、邊界兩種插入Steiner 點方法。

Thingi10k 數據集的測試結果表明，本文方法相較于文獻[5]提出的程序以及TetGen1.6，可以零Steiner點恢復更多例子，并且普遍降低了邊界恢復所需的Steiner 點數量。這對于復雜例子的約束邊界恢復、提升程序魯棒性以及網格質量有著重要意義。

本文方法在時間性能上尚有不足。后續(xù)計劃就拓撲變換方法、模擬退火策略與參數等方面開展優(yōu)化工作。