一種基于局部平均法向變形的網格參數化方法

焦 沖, 蘇科華, 吳博文, 任術波, 辛 寧

(1. 武漢大學 計算機學院, 武漢 430072; 2. 中國空間技術研究院 通信與導航衛星總體部, 北京100094)

0 引 言

三維網格模型在工業制造、 輔助醫療、 3D打印、 文化娛樂和建筑設計等領域應用廣泛. 由于物體的幾何形狀通常較復雜, 而這些復雜的幾何形狀使后續的網格處理很困難, 因此, 通常需要將三維網格參數化到一個簡單的參數域, 例如平面參數域、 球面參數域等. 作為計算機圖形學、 計算機輔助幾何設計和數字幾何處理中的一個重要工具, 網格參數化在紋理映射[1]、 網格變形[2-3]、 形狀建模[4]及網格優化[5-6]等網格處理中都具有重要作用.

網格參數化問題一般可描述為: 給定一個二維流形的三維網格和參數域, 尋找從參數域的點到三維網格點的一一映射, 使參數域的網格與原網格拓撲同構, 且保證三角形之間互不重疊. 理想狀態下, 三維網格到參數域之間的參數化映射是等距的, 即原始網格的邊長與夾角在映射后應保持不變. 但除可展曲面外, 一般曲面都達不到這一理想條件. 一種參數化方法只能盡可能保持角度或者面積不產生扭曲, 但不能同時消除角度扭曲和面積扭曲. 因此, 為提高算法的適用性, 參數化的核心任務之一就是盡可能降低某種類型的扭曲, 例如保角參數化能最小化三角網格的角度扭曲, 從而有效地保留網格的形狀信息.

早期的平面參數化算法一般采用凸組合的方式, 通過固定邊界求解線性方程得到平面參數化[7-9]. 這類方法雖然能保證一一映射, 但都具有很大的扭曲. 因此, 一些方法采用自由邊界或分割展平的方式直接控制扭曲[10-12]. 優化能量的方法使用無翻轉的參數化作為初始, 并優化某種能量以降低扭曲實現參數化[13-17]. 由于扭曲的上界很難預先確定, 因此該類方法需要不斷嘗試尋找一個最佳上界, 同時由于扭曲度量通常是高度非線性的, 這些基于能量優化的方法計算代價較高.

能量優化方法可應用于球面拓撲的模型, 但對于高曲率的區域目標能量可能會導致翻轉[18]. 一種解決方式是通過一個動態調整的參數生成雙射的球面參數化[19]; 另一種方式是通過最小化調和能量將0虧格的流形變分為全局共形參數化[20-21]. 如文獻[22]使用一些中間的參考三角形定義能量, 并生成具有低等距扭曲的無折疊參數化.

基于幾何流的方法成功地實現了在同一框架下實現多種參數域的參數化. 平均曲率流(MCF)[23]是用于演化曲面幾何的最基本流之一, 其可等價地表示為最小化曲面嵌入的梯度或最小化曲面面積的流. 但奇點的存在導致MCF計算過程不穩定. Bobenko等[24]提出了Willmore流, 但Willmore流需要曲面自身接近于球面, 并且依賴于高階導數； Zhao等[25]提出了單位法向流的定義, 并將其應用于網格參數化, 但其計算效率較低并且對于某些復雜的模型效果不佳. 基于Monge-Brenier理論, 最優傳輸理論(OMT)也可用于計算保面積參數化[26]. 對于poly-annulus曲面, Su等[27]結合OMT和Ricci流計算保面積的參數化. 基于離散Calabi流, Zhao等[28]提出了一種保角的參數化算法. 為提高算法性能, Su等[29]為離散Calabi流引進了邊的翻轉、 擬牛頓法、 最優步長、 優先級嵌入和邊界自定義等操作； Su等[30]利用離散的李導數流計算圓盤和球面的保面積參數化.

基于此, 本文針對單/多邊界的開網格以及虧格為0的封閉網格, 提出一種基于局部平均法向變形的網格參數化方法. 首先通過計算局部平均面法向, 以該平均法向為引導, 交替地將三角形旋轉到新位置, 然后利用Poisson系統解算旋轉拉伸能量方程, 得到分散的三角形頂點新位置, 在視覺上如同將分散的三角形重新縫合. 通過不斷迭代, 將頂點不斷向鄰居位置推移, 使重建網格達到目標曲率. 該方法將局部平均法向的思想引入到MCF的框架中, 并基于該思想用動態控制平均法向權重的方式避免三角形翻轉, 極大減少了算法的迭代次數, 提高了求解速度, 解決了傳統MCF算法求解不穩定、 迭代速度過慢的問題. 同時, 優化策略也盡量保證了最小化參數化扭曲. Torse模型的原始網格通過迭代變形為最終球面參數的參數化過程如圖1所示.

圖1 Torse模型的參數化過程Fig.1 Parameterization process of Torse model

1 算法描述

類似于MCF方法, 本文的變形方法通過將頂點移向其鄰居頂點實現, 但不直接計算頂點的鄰域平均位置, 而是通過面法向量場平均實現頂點位置的計算. 因此, 本文算法共分為兩個階段.

1) 局部平均法向旋轉： 計算每個三角形k-階鄰居域的平均法向量, 并以該法向為目標法向, 以笛卡爾坐標系下的原點作為旋轉中心, 為每個三角形執行旋轉操作；

2) Poisson表面形變： 結合Poisson系統“縫合”網格, 通過優化一種拉伸能量計算頂點的新位置.

交替執行上述兩個階段, 直到網格收斂, 即網格上每個點的平均曲率都相等. 本文用M={V,F}表示三角網格, 其中V={vi,i=1,2,…,Nv}為頂點集合,F={fi,i=1,2,…,Nf}為三角形面集合.算法的實現過程中采用半邊結構存儲網格信息.

1.1 局部平均法向旋轉

本文以k-階鄰居域的平均法向作為目標法向量, 對每個三角形應用旋轉變換, 并記錄旋轉矩陣.

1.1.1 鄰居面法向平均

(1)

1.1.2k-階鄰居面選取

圖2 不同k值下不同算法的收斂速度對比Fig.2 Comparison of convergence speed of different algorithms under different k values

1.2 Poisson表面形變

通過局部的平均法向分別對每個三角形執行旋轉操作會導致整個網格的拓撲發生改變, 頂點被分裂. 因此, 算法基于Poisson系統重建變形的三角網格, 保證每個面都滿足當前指定的目標法向量. 由于Poisson系統只能逼近當前目標法向, 采用迭代式的過程可引導面法向不斷靠近最終的目標法向. 本質上, 迭代過程可視為是網格表面不斷變形的過程, 常見網格變形算法的目的包括兩方面: 一是保持原始網格的局部特征; 二是盡可能保持網格的度量, 即保證最小化邊長變化. 由于參數化算法需要使扭曲盡可能小, 因此本文在重建過程中應盡量保持網格的度量.

由文獻[33]可知, 網格表面形變能量可顯式地分為兩部分： 拉伸能量和彎曲能量. 拉伸能量與度量相關, 而彎曲能量與曲率相關. 本文算法以拉伸能量為基礎, 并以每個三角形的局部平均法向為限制, 建立拉伸能量方程：

(2)

對于每次迭代,Rij可通過三角形的當前法向量和目標法向量計算, 因此Es可視為是只關于頂點坐標的二次函數.通過對Es求導可知, 在梯度為0時, 其解為最優解.通過化簡, 可建立關于頂點坐標的一個線性方程組為

(3)

其解為變形后網格頂點的新坐標, 其中vj～vi表示頂點vj與頂點vi相鄰接.在線性方程組(3)中, 需已知網格中每個三角形面的旋轉矩陣Rji.雖然預先未知迭代變形后頂點的確切位置, 但通過上述計算可知每個三角形在此次迭代過程中的目標法向量, 因此可為每個面計算旋轉矩陣.由于式(3)右邊為已知量, 因此可將其簡寫為三維向量bi.對于方程組(3)的系數, 定義為

(4)

因此, 可將方程組(3)簡寫為

Lv′=b.

(5)

求解方程(5)即可得到網格頂點在本次迭代中的新坐標.迭代操作的終止條件是所有頂點的平均曲率相等, 由于計算機精度的原因, 可將此條件放寬至頂點的平均曲率相差在一個較小的閾值內.由于L是對稱正定的, 因此, 可利用Cholesky因式分解求解線性方程組, 比能量優化方法更簡單高效.在每次迭代中為每個三角形重新計算旋轉矩陣, 用新的旋轉矩陣與頂點坐標更新L和b.對于不同拓撲以及規模大小不同的網格, 所需的迭代次數一般不同.對于小型網格, 一般僅需幾十次迭代即可達到收斂.

1.3 算法的加速策略

實驗發現, 僅使用k-階鄰居平均法向量時, 鄰居域范圍越大, 頂點移動范圍越大, 收斂速度越快.但在某些區域, 尤其是高曲率區域, 過大的鄰居域可能會導致三角形翻轉.本文通過動態調整鄰居域的平均法向量權重, 在保證三角形不翻轉的同時加快算法收斂.即在引進全局平均法向量的同時, 設置懲罰函數調整其對目標法向量的影響, 保證三角形不發生翻轉.

(6)

由于在形變過程中, 網格的整體形狀會發生變化,n*也是不斷變化的, 固定權重μ通常無法取得理想結果, 因此本文定義懲罰函數動態地改變μ值.對于平面參數化, 懲罰函數的選擇需滿足以下條件: 對于面fi, 如果n*與ni之間為鈍角, 則應減小權重μ; 由于翻轉的非邊界點高斯曲率絕對值收斂于2π, 所以在高曲率區域使用n*加速收斂也會導致翻轉, 因此當頂點的曲率較大時, 應減小權重μ.基于上述條件, 可設計懲罰函數為

(7)

(8)

對于平面參數域, 本文的基本思想是通過預測每個面的法向方向加速收斂.對于球面參數域, 該策略同樣可行, 其懲罰函數需滿足以下條件: 對于面fi, 如果n*與ni之間為鈍角, 則應減小權重μ; 如果網格M是非凸的, 則應減小權重μ.基于上述條件, 可設計懲罰函數為

(9)

(10)

2 實驗與評估

本文以C++實現該算法, 并在AHSP[22]提供的模型上進行測試, 同時與其他算法進行比較. 首先在一個模型集合中對本文的參數化算法進行測試, 包括各種拓撲模型, 驗證其有效性和高效性; 然后與目前其他先進的網格參數化方法進行對比; 最后利用具有不同剖分的網格證明其對網格剖分不敏感. 實驗平臺為配置AMD Ryzen-7-1700處理器, GTX 1060顯卡, 16 GB內存, Windows 10操作系統的臺式電腦； 圖像處理庫為OpenMesh 6.3, 矩陣處理庫為Eigen 3.3.3; 輸入的網格文件為obj格式.

2.1 扭曲度量

為評價參數化效果, 目前已有很多種類參數化扭曲的度量方式. 假設σ1和σ2分別是從原三角形到參數化三角形轉換的Jacobi矩陣的最大和最小特征值, 用ρi表示fi的面積與網格總面積之比, 扭曲的度量方式如下：

文獻[34]研究表明,Darea和Dangle的值越接近于2, 其面積扭曲或角度扭曲越小.此外, 為進一步比較算法性能, 用等距能量EARAP、 保角能量EASAP和Green-Lagrange能量(EGL)[1]評估參數化效果.

2.2 數據集測試

為驗證本文提出的基于局部平均法向變形參數化方法的有效性和魯棒性, 在2 063個模型上測試該算法, 所有模型均來自文獻[22]中數據集.對每個模型的參數化結果, 先計算其扭曲并記錄運行時間, 再將所有結果繪制成頻率直方圖, 如圖3所示, 其中max,min,ave和std分別表示對應扭曲或時間的最大值、 最小值、 平均值和標準偏差.由圖3可見, 本文算法能高效生成低扭曲的參數化.

圖3 2 063個模型的參數化結果直方圖Fig.3 Histogram of parameterization results with 2 063 models

對具有多邊界的開網格及高曲率模型, 本文方法也能計算出高質量的參數化, 圖4和圖5分別為這兩種類型網格的參數化結果.

圖4 具有多邊界的汽車模型平面參數化結果Fig.4 Plane parameterization results of car model with multiple boundaries

圖5 具有高曲率區域的網格模型及其球面參數化結果Fig.5 Spherical parameterization results on mesh model with high curvature area

2.3 算法對比測試

為驗證本文算法的實用性和高效性, 將本文算法與其他參數化方法進行比較, 包括單位法向流(UNF)[25]、 可伸縮局部內射映射(SLIM)[17]和多層次球面參數化(AHSP)[22]. 對于UNF, 遵循其默認參數, 其他算法均采用其源代碼, 所有實驗均在同一臺電腦上運行. 由于篇幅所限, 本文僅選取部分網格的參數化度量結果列于表1和表2. 圖6和圖7分別為Bear和VaseLion模型的參數化視覺效果.

表1 不同算法平面參數化方法的對比結果

表2 不同算法球面參數化方法的對比結果

由表1、 表2及圖6、 圖7可見： 基于局部平均法向思想的加速策略減少了迭代時間, 相比其他算法有更高的計算效率； 在平面參數化中, 雖然保面積的效果稍遜于其他方法, 但卻擁有較好的保形效果, 并能很好地降低拉伸扭曲; 在球面參數化中, 本文方法盡管在降低面積扭曲方面略差于AHSP方法, 但在其他度量中都取得了較好結果. 因此, 本文方法可在同一個框架下高效地實現高質量的平面參數化和球面參數化.

圖6 不同算法對Bear模型平面參數化方法的比較Fig.6 Comparison of different algorithms for Bear model by plane parameterization methods

圖7 不同算法對VaseLion模型球面參數化方法的比較Fig.7 Comparison of different algorithms for VaseLion model by spherical parameterization methods

2.4 算法分析

2.4.1 算法效率

通過統計算法在數據集上的運行效率, 可得出以下結論: 算法的收斂速度與網格的規格大小密切相關, 這是因為對于規模較大的網格, 針對每個頂點與三角形, 都需要進行一系列的計算, 因此在每次迭代過程中, 計算規模隨著網格規模增大而增加, 總運算時間也會隨之延長. 但無論在平面參數化還是在球面參數化中, 本文算法的效率均具有明顯優勢.

2.4.2 對不同網格剖分的魯棒性

為測試本文算法對不同網格剖分的敏感性, 對部分網格模型進行重新處理, 例如隨機增加頂點數, 隨機減少頂點數, 或使頂點分布均勻等. 將本文算法應用于不同剖分的網格模型并統計其運行時間及度量扭曲. 圖8為本文算法在不同剖分網格上的測試結果. 通過將棋盤格重新貼回模型中發現, 同一網格不同剖分下的參數化扭曲基本一致; 對比不同剖分模型的運行時間, 即使模型剖分變大, 運行時間卻僅有很小的變化. 實驗結果表明, 本文算法在不同剖分的網格中都能獲得高質量的參數化, 具有較強的魯棒性.

圖8 Buddha模型在4種網格剖分下的參數化扭曲魯棒性測試Fig.8 Robustness test of parametric distortion under four mesh subdivisions of Buddha model

綜上所述, 本文針對具有單/多邊界的開網格以及虧格為0的封閉網格設計了一種基于局部平均法向變形的網格參數化方法. 通過交替執行兩階段的操作將頂點不斷向鄰居域移動, 每次迭代通過求解稀疏線性系統使網格變形, 使網格逐步收斂為一個常平均曲率曲面, 從而生成低扭曲的參數化結果. 一方面, 使用k-階鄰居面平均法向的思想極大加快了算法的運行效率； 另一方面, 通過懲罰函數動態調整每次迭代的目標法向, 有效保證了網格無翻轉情況的發生. 通過在2 063個模型上的測試, 驗證了本文算法的有效性和魯棒性. 由于局部平均法向思想僅受曲率限制, 因此本文算法對于三角剖分不敏感, 實用性較強. 此外, 與其他參數化算法相比, 本文算法能高效生成低扭曲的參數化結果.