四面體網(wǎng)格邊界恢復(fù)改進算法

曾寅家，陳建軍，傅珂杰

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

0 引言

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

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

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

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

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

本文算法是對文獻[5]算法的改進。改進算法先通過基于遞歸殼變換的前處理步驟恢復(fù)大部分遺失邊界，再增加Steiner 點恢復(fù)所有遺失邊界。與文獻[5]算法相比，新算法的創(chuàng)新點包括：1）引入模擬退火思想，使新算法更具有全局最優(yōu)性；2）實現(xiàn)了兩類Steiner 點插入方案，減少邊界上Steiner 點數(shù)量；3）優(yōu)化了邊界恢復(fù)流程，充分利用每一個插入的Steiner點。改進算法旨在減少Steiner 點數(shù)量和減輕Steiner點的負面效應(yīng)，并提升了算法的魯棒性。為了展現(xiàn)算法效果，本文利用Thingi10k 數(shù)據(jù)集[6]進行測試。

1 相關(guān)工作

1.1 邊界恢復(fù)相關(guān)研究

國內(nèi)外關(guān)于前處理、Steiner 點插入、后處理已有許多研究工作。

關(guān)于前處理，Radon[7]提出四種三維空間中的基本變換：1-4 變換、4-1 變換、2-3 變換、3-2 變換（如圖1 所示）；Joe[8]則通過多種變換的組合，實現(xiàn)了更多復(fù)雜變換。Shewchuk[9]的研究指出，大多組合變換可以通過兩次n-m 變換（n-to-m flip）實現(xiàn)，n-m 變換指的是一種通用的2-3、3-2 變換的組合。Si[10]將n-m變換加以推廣，并系統(tǒng)性地實現(xiàn)在網(wǎng)格生成程序TetGen 中。Liu等[11]則提出了一種稱為小多面體重連（small polyhedron reconnection,SPR）的方法，相比于前述變換方法找到局部最優(yōu)解的效率更高，許多后繼研究都對其進行了改進[12-13]，但由于該方法為NP 復(fù)雜度，因此建議輸入多面體面數(shù)低于40[11]。陳建軍等[5]則提出了一種遞歸殼變換邊界恢復(fù)策略，鑒于該算法與本文的密切聯(lián)系，下段先對其進行簡要回顧。

圖1 1-4、4-1 變換和3-2、2-3 變換示意圖Fig.1 Schematic diagram of transformation 1-4,4-1,3-2 and 2-3

殼指的是共享同一條邊的所有單元構(gòu)成的集合，如圖2[14]所示；而殼變換則是對裙多邊形進行重三角化的網(wǎng)格局部重連技術(shù)，如圖3[14]所示。針對殼變換問題，Shewchuk[9]提出了一種動態(tài)規(guī)劃算法，可以高效地找到裙多邊形三角化的最優(yōu)解。陳建軍等[5]則對該算法進行了拓展，以適應(yīng)裙多邊形部分三角化的情形，并提出殼變換的遞歸策略，大幅增加了其成功率。文獻[5]進一步提出利用遞歸殼變換移除與邊界約束相交的網(wǎng)格邊/面，從而恢復(fù)邊界的方法，并為了保證算法的收斂性與效率，引入了多個約束條件，其中包括：1）不對已經(jīng)恢復(fù)的邊界邊執(zhí)行殼變換；2）輸出網(wǎng)格與邊界相交數(shù)不大于輸入網(wǎng)格。

圖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)化范圍更廣。因此，本文基于遞歸殼變換，提出改進的約束邊界恢復(fù)方法。

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

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

1.2 模擬退火算法

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

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

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

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

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

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

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

2 本文方法

2.1 總體流程

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

圖4 網(wǎng)格生成流程圖Fig.4 Flowchart of mesh generation process

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

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

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

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

2.2 改進的約束邊界恢復(fù)算法

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

2.2.1 結(jié)合模擬退火的約束邊界恢復(fù)算法

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

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

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

該算法的說明如下：

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

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

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

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

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

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

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

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

2.2.2 交點數(shù)的計算與查詢

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

該算法外層是對于邊界約束的遍歷，內(nèi)層則是“網(wǎng)格實體—邊界實體求交”算法。首先，程序引入了AABB（axis-aligned bounding boxes）樹[33]對外層遍歷進行加速。具體方法是：在初始化階段中，將所有邊界約束實體加入AABB 樹，再在求交時通過AABB樹快速過濾不相交的邊界約束。在包圍盒重疊程度較低的情況下，該方法可將遍歷復(fù)雜度由 O(M)降為O(logM)，其中M為邊界實體數(shù)量。

其次，為了避免浮點誤差導(dǎo)致的求交算法錯誤，本算法利用了Guigue等[34]實現(xiàn)的“三角形-三角形相交”檢測程序，其中包含幾何精確的“線面相交”函數(shù)。

最后，使用了一個哈希表存儲已計算出網(wǎng)格實體與邊界的交點數(shù)，以減少二次計算帶來的性能損耗。圖6 展示了網(wǎng)格面交點數(shù)哈希表的C++數(shù)據(jù)結(jié)構(gòu)定義，查詢邏輯偽代碼如算法 3 所示。哈希表利用二維變長數(shù)組儲存，利用其內(nèi)存局部性，可加速哈希碰撞時的搜索時間；對于任意網(wǎng)格面p1p2p3，本算法在HashSort 函數(shù)中對其進行結(jié)果固定的亂序排序，取首位為哈希值，后兩位為哈希碰撞時的識別標識；并利用了C++Union 關(guān)鍵字，將兩位的標識等價為單個長整型，將兩次相等判斷縮減到了一次，進一步加快了查詢時間。

圖6 三角面哈希表數(shù)據(jù)結(jié)構(gòu)Fig.6 Hash table data structure for triangular polygons

2.3 Steiner 點插入策略

僅通過拓撲變換不能解決所有邊界恢復(fù)問題，此時需要考慮插入Steiner 點。由于體內(nèi)Steiner 點不改變輸入邊界約束，首先嘗試在四面體內(nèi)部插點，進一步地，再利用插點打斷邊界的方法，以恢復(fù)所有邊界。

本文受到TetGen 研究成果[35]的啟發(fā)，實現(xiàn)了一種針對特殊情形的體內(nèi)Steiner 點插入方法，圖7 展示了其過程示例。所謂特殊情形，指的是遺失邊ab與四面體網(wǎng)格 T的所有相交單元 Tab具有共邊cd。換而言之，Tab由邊 Δd所在殼上的連續(xù)單元p0p1cd,p1p2cd,...,pn-1pncd構(gòu)成，如圖7（a）所示，pi(i=0,1,···,n) 為殼的裙點，且p0=a、pn=b。方法具體步驟如下：

圖7 體內(nèi)Steiner 點插入原理示例圖Fig.7 Principle illustration of the interior Steiner point insertion

1）將 Tab的覆蓋空間記作空腔C，首先從C中分離四面體abcd，得到空腔C′，如圖7（b）所示。

2）從ab中點出發(fā)，通過點光滑化例程得到點 6，使其滿足與C′的所有表面三角形所構(gòu)成的四面體具有正體積。

3）若上述點e不存在，則失敗；否則將點e與C′的所有表面三角形相連接，得到C′的四面體化，如圖7（c）所示。

由此通過插入點 6，得到了空腔C的一種四面體化，其包含所需恢復(fù)的遺失邊ab。

針對體內(nèi)插點無法解決的情況，采用了Si 和G?rtner[36]提出的一種保質(zhì)Steiner 點插入策略對邊界邊進行分裂。該方法根據(jù)遺失邊端點周圍的面網(wǎng)格性質(zhì)，分三類情況計算Steiner 點的位置，并且在計算時額外考慮了遺失邊周圍網(wǎng)格點的分布以及端點的局部特征尺寸[37]。

對于遺失邊界面，則考慮直接往面心加點。

3 實驗及結(jié)果分析

所有實驗均在AMD Ryzen 7 5800H、3.2GHz CPU、16GB RAM 平臺上進行。測試集取自Thingi10k 數(shù)據(jù)集[6]，該數(shù)據(jù)集中包含1 萬個三角網(wǎng)格模型，其中多數(shù)樣本存在自相交、重節(jié)點等問題。鑒于本文重點非面網(wǎng)格修復(fù)，實際測試中采用了幾何正確的面網(wǎng)格的4 303 例。

3.1 模擬退火參數(shù)的評估與選取

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

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

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

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

3.2 與其他方法的比較

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

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

表2 4 303 例樣本的運行結(jié)果統(tǒng)計數(shù)據(jù)Table 2 Statistical data of running results with 4 303 samples

圖9 文獻[5]提出方法處理的面網(wǎng)格結(jié)果Fig.9 Surface mesh processed by the method in Ref.[5]

對于2 829 例TetGen 生成的Steiner 點數(shù)量大于0 的樣本，使用本文方法能夠?qū)⑵渲械?95 例轉(zhuǎn)化為零Steiner 點恢復(fù)，轉(zhuǎn)化率達到17.50%。另外，在2 402例本文程序含Steiner 點的樣本中，TetGen 僅68 例零Steiner 點恢復(fù)。說明本文方法較TetGen 也具有更高的零Steiner 點恢復(fù)率。

后文將進一步將本文方法與另兩程序進行橫向?qū)Ρ取榱烁又庇^地比較各程序的Steiner 點數(shù)量，在4 258 例上述程序都能成功生成的樣本中，剔除了1 259例所有程序都無Steiner 點的樣本，其余2 999 例的Steiner 點數(shù)量分布如圖10 所示，統(tǒng)計數(shù)據(jù)如表3 所示。從數(shù)據(jù)上看，本文方法較TetGen 將2 999 例樣本的合計Steiner 點數(shù)量降低了29.86%，并且生成結(jié)果中Steiner 點數(shù)量≥5 的樣本數(shù)量也最少，說明本文方法總體效果顯著。另外注意到文獻[5]方法的Steiner 點數(shù)量偏多，原因在于該程序采用了單次“遞歸殼變換邊界恢復(fù)—插入Steiner 點”的線性流程，導(dǎo)致加入冗余Steiner 點。而本文方法通過拓撲變換與插點步驟的穿插執(zhí)行，即“兩層+兩輪”迭代，充分利用每個Steiner 點，有效降低了Steiner 點數(shù)量。

表3 2 999 例樣本Steiner 點數(shù)量統(tǒng)計數(shù)據(jù)Table 3 Statistical data for the Steiner point number in 2 999 samples

圖10 2 999 例樣本Steiner 點數(shù)量分布圖Fig.10 Distribution of the Steiner point number in 2999 samples

圖11 展示了2 999 例樣本中的3 例典型網(wǎng)格。表4 展示了其Steiner 點數(shù)據(jù)。由此可見，本文算法能以較少Steiner 點數(shù)量處理不同形態(tài)的輸入網(wǎng)格。

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

圖11 3 例典型樣本網(wǎng)格Fig.11 Three typical test meshes

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

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

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

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

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

表6 倫敦塔橋模型網(wǎng)格生成統(tǒng)計數(shù)據(jù)Table 6 Mesh generation statistics for the Tower Bridge model

圖12 倫敦塔橋模型面網(wǎng)格Fig.12 Surface mesh of Tower Bridge model

4 結(jié)論

本文提出了一種改進的Delaunay 四面體網(wǎng)格生成的邊界恢復(fù)流程。該流程主要由改進的約束邊界恢復(fù)步驟以及體內(nèi)、邊界Steiner 點插入步驟構(gòu)成，表現(xiàn)為“兩層+兩輪”的迭代流程，旨在降低Steiner 點數(shù)量并減少Steiner 點的負面效應(yīng)。其中改進的約束邊界恢復(fù)方法基于遞歸殼變換算法，并結(jié)合模擬退火算法，利用Metropolis 準則使遞歸殼變換脫離局部最優(yōu)解，進一步減少遺失邊界數(shù)量。同時，針對“網(wǎng)格實體—邊界實體求交”的效率與魯棒性問題，引入和提出了多個關(guān)鍵方法進行優(yōu)化。對于拓撲變換無法恢復(fù)的邊界，實現(xiàn)了體內(nèi)、邊界兩種插入Steiner 點方法。

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

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