一種基于局部平均法向變形的網(wǎng)格參數(shù)化方法

焦 沖, 蘇科華, 吳博文, 任術(shù)波, 辛 寧

(1. 武漢大學(xué) 計(jì)算機(jī)學(xué)院, 武漢 430072; 2. 中國(guó)空間技術(shù)研究院 通信與導(dǎo)航衛(wèi)星總體部, 北京100094)

0 引 言

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

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

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

能量?jī)?yōu)化方法可應(yīng)用于球面拓?fù)涞哪Ｐ? 但對(duì)于高曲率的區(qū)域目標(biāo)能量可能會(huì)導(dǎo)致翻轉(zhuǎn)[18]. 一種解決方式是通過(guò)一個(gè)動(dòng)態(tài)調(diào)整的參數(shù)生成雙射的球面參數(shù)化[19]; 另一種方式是通過(guò)最小化調(diào)和能量將0虧格的流形變分為全局共形參數(shù)化[20-21]. 如文獻(xiàn)[22]使用一些中間的參考三角形定義能量, 并生成具有低等距扭曲的無(wú)折疊參數(shù)化.

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

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

圖1 Torse模型的參數(shù)化過(guò)程Fig.1 Parameterization process of Torse model

1 算法描述

類似于MCF方法, 本文的變形方法通過(guò)將頂點(diǎn)移向其鄰居頂點(diǎn)實(shí)現(xiàn), 但不直接計(jì)算頂點(diǎn)的鄰域平均位置, 而是通過(guò)面法向量場(chǎng)平均實(shí)現(xiàn)頂點(diǎn)位置的計(jì)算. 因此, 本文算法共分為兩個(gè)階段.

1) 局部平均法向旋轉(zhuǎn)： 計(jì)算每個(gè)三角形k-階鄰居域的平均法向量, 并以該法向?yàn)槟繕?biāo)法向, 以笛卡爾坐標(biāo)系下的原點(diǎn)作為旋轉(zhuǎn)中心, 為每個(gè)三角形執(zhí)行旋轉(zhuǎn)操作；

2) Poisson表面形變： 結(jié)合Poisson系統(tǒng)“縫合”網(wǎng)格, 通過(guò)優(yōu)化一種拉伸能量計(jì)算頂點(diǎn)的新位置.

交替執(zhí)行上述兩個(gè)階段, 直到網(wǎng)格收斂, 即網(wǎng)格上每個(gè)點(diǎn)的平均曲率都相等. 本文用M={V,F}表示三角網(wǎng)格, 其中V={vi,i=1,2,…,Nv}為頂點(diǎn)集合,F={fi,i=1,2,…,Nf}為三角形面集合.算法的實(shí)現(xiàn)過(guò)程中采用半邊結(jié)構(gòu)存儲(chǔ)網(wǎng)格信息.

1.1 局部平均法向旋轉(zhuǎn)

本文以k-階鄰居域的平均法向作為目標(biāo)法向量, 對(duì)每個(gè)三角形應(yīng)用旋轉(zhuǎn)變換, 并記錄旋轉(zhuǎn)矩陣.

1.1.1 鄰居面法向平均

(1)

1.1.2k-階鄰居面選取

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

1.2 Poisson表面形變

通過(guò)局部的平均法向分別對(duì)每個(gè)三角形執(zhí)行旋轉(zhuǎn)操作會(huì)導(dǎo)致整個(gè)網(wǎng)格的拓?fù)浒l(fā)生改變, 頂點(diǎn)被分裂. 因此, 算法基于Poisson系統(tǒng)重建變形的三角網(wǎng)格, 保證每個(gè)面都滿足當(dāng)前指定的目標(biāo)法向量. 由于Poisson系統(tǒng)只能逼近當(dāng)前目標(biāo)法向, 采用迭代式的過(guò)程可引導(dǎo)面法向不斷靠近最終的目標(biāo)法向. 本質(zhì)上, 迭代過(guò)程可視為是網(wǎng)格表面不斷變形的過(guò)程, 常見(jiàn)網(wǎng)格變形算法的目的包括兩方面: 一是保持原始網(wǎng)格的局部特征; 二是盡可能保持網(wǎng)格的度量, 即保證最小化邊長(zhǎng)變化. 由于參數(shù)化算法需要使扭曲盡可能小, 因此本文在重建過(guò)程中應(yīng)盡量保持網(wǎng)格的度量.

由文獻(xiàn)[33]可知, 網(wǎng)格表面形變能量可顯式地分為兩部分： 拉伸能量和彎曲能量. 拉伸能量與度量相關(guān), 而彎曲能量與曲率相關(guān). 本文算法以拉伸能量為基礎(chǔ), 并以每個(gè)三角形的局部平均法向?yàn)橄拗? 建立拉伸能量方程：

(2)

對(duì)于每次迭代,Rij可通過(guò)三角形的當(dāng)前法向量和目標(biāo)法向量計(jì)算, 因此Es可視為是只關(guān)于頂點(diǎn)坐標(biāo)的二次函數(shù).通過(guò)對(duì)Es求導(dǎo)可知, 在梯度為0時(shí), 其解為最優(yōu)解.通過(guò)化簡(jiǎn), 可建立關(guān)于頂點(diǎn)坐標(biāo)的一個(gè)線性方程組為

(3)

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

(4)

因此, 可將方程組(3)簡(jiǎn)寫(xiě)為

Lv′=b.

(5)

求解方程(5)即可得到網(wǎng)格頂點(diǎn)在本次迭代中的新坐標(biāo).迭代操作的終止條件是所有頂點(diǎn)的平均曲率相等, 由于計(jì)算機(jī)精度的原因, 可將此條件放寬至頂點(diǎn)的平均曲率相差在一個(gè)較小的閾值內(nèi).由于L是對(duì)稱正定的, 因此, 可利用Cholesky因式分解求解線性方程組, 比能量?jī)?yōu)化方法更簡(jiǎn)單高效.在每次迭代中為每個(gè)三角形重新計(jì)算旋轉(zhuǎn)矩陣, 用新的旋轉(zhuǎn)矩陣與頂點(diǎn)坐標(biāo)更新L和b.對(duì)于不同拓?fù)湟约耙?guī)模大小不同的網(wǎng)格, 所需的迭代次數(shù)一般不同.對(duì)于小型網(wǎng)格, 一般僅需幾十次迭代即可達(dá)到收斂.

1.3 算法的加速策略

實(shí)驗(yàn)發(fā)現(xiàn), 僅使用k-階鄰居平均法向量時(shí), 鄰居域范圍越大, 頂點(diǎn)移動(dòng)范圍越大, 收斂速度越快.但在某些區(qū)域, 尤其是高曲率區(qū)域, 過(guò)大的鄰居域可能會(huì)導(dǎo)致三角形翻轉(zhuǎn).本文通過(guò)動(dòng)態(tài)調(diào)整鄰居域的平均法向量權(quán)重, 在保證三角形不翻轉(zhuǎn)的同時(shí)加快算法收斂.即在引進(jìn)全局平均法向量的同時(shí), 設(shè)置懲罰函數(shù)調(diào)整其對(duì)目標(biāo)法向量的影響, 保證三角形不發(fā)生翻轉(zhuǎn).

(6)

由于在形變過(guò)程中, 網(wǎng)格的整體形狀會(huì)發(fā)生變化,n*也是不斷變化的, 固定權(quán)重μ通常無(wú)法取得理想結(jié)果, 因此本文定義懲罰函數(shù)動(dòng)態(tài)地改變?chǔ)讨?對(duì)于平面參數(shù)化, 懲罰函數(shù)的選擇需滿足以下條件: 對(duì)于面fi, 如果n*與ni之間為鈍角, 則應(yīng)減小權(quán)重μ; 由于翻轉(zhuǎn)的非邊界點(diǎn)高斯曲率絕對(duì)值收斂于2π, 所以在高曲率區(qū)域使用n*加速收斂也會(huì)導(dǎo)致翻轉(zhuǎn), 因此當(dāng)頂點(diǎn)的曲率較大時(shí), 應(yīng)減小權(quán)重μ.基于上述條件, 可設(shè)計(jì)懲罰函數(shù)為

(7)

(8)

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

(9)

(10)

2 實(shí)驗(yàn)與評(píng)估

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

2.1 扭曲度量

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

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

2.2 數(shù)據(jù)集測(cè)試

為驗(yàn)證本文提出的基于局部平均法向變形參數(shù)化方法的有效性和魯棒性, 在2 063個(gè)模型上測(cè)試該算法, 所有模型均來(lái)自文獻(xiàn)[22]中數(shù)據(jù)集.對(duì)每個(gè)模型的參數(shù)化結(jié)果, 先計(jì)算其扭曲并記錄運(yùn)行時(shí)間, 再將所有結(jié)果繪制成頻率直方圖, 如圖3所示, 其中max,min,ave和std分別表示對(duì)應(yīng)扭曲或時(shí)間的最大值、 最小值、 平均值和標(biāo)準(zhǔn)偏差.由圖3可見(jiàn), 本文算法能高效生成低扭曲的參數(shù)化.

圖3 2 063個(gè)模型的參數(shù)化結(jié)果直方圖Fig.3 Histogram of parameterization results with 2 063 models

對(duì)具有多邊界的開(kāi)網(wǎng)格及高曲率模型, 本文方法也能計(jì)算出高質(zhì)量的參數(shù)化, 圖4和圖5分別為這兩種類型網(wǎng)格的參數(shù)化結(jié)果.

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

圖5 具有高曲率區(qū)域的網(wǎng)格模型及其球面參數(shù)化結(jié)果Fig.5 Spherical parameterization results on mesh model with high curvature area

2.3 算法對(duì)比測(cè)試

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

表1 不同算法平面參數(shù)化方法的對(duì)比結(jié)果

表2 不同算法球面參數(shù)化方法的對(duì)比結(jié)果

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

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

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

2.4 算法分析

2.4.1 算法效率

通過(guò)統(tǒng)計(jì)算法在數(shù)據(jù)集上的運(yùn)行效率, 可得出以下結(jié)論: 算法的收斂速度與網(wǎng)格的規(guī)格大小密切相關(guān), 這是因?yàn)閷?duì)于規(guī)模較大的網(wǎng)格, 針對(duì)每個(gè)頂點(diǎn)與三角形, 都需要進(jìn)行一系列的計(jì)算, 因此在每次迭代過(guò)程中, 計(jì)算規(guī)模隨著網(wǎng)格規(guī)模增大而增加, 總運(yùn)算時(shí)間也會(huì)隨之延長(zhǎng). 但無(wú)論在平面參數(shù)化還是在球面參數(shù)化中, 本文算法的效率均具有明顯優(yōu)勢(shì).

2.4.2 對(duì)不同網(wǎng)格剖分的魯棒性

為測(cè)試本文算法對(duì)不同網(wǎng)格剖分的敏感性, 對(duì)部分網(wǎng)格模型進(jìn)行重新處理, 例如隨機(jī)增加頂點(diǎn)數(shù), 隨機(jī)減少頂點(diǎn)數(shù), 或使頂點(diǎn)分布均勻等. 將本文算法應(yīng)用于不同剖分的網(wǎng)格模型并統(tǒng)計(jì)其運(yùn)行時(shí)間及度量扭曲. 圖8為本文算法在不同剖分網(wǎng)格上的測(cè)試結(jié)果. 通過(guò)將棋盤(pán)格重新貼回模型中發(fā)現(xiàn), 同一網(wǎng)格不同剖分下的參數(shù)化扭曲基本一致; 對(duì)比不同剖分模型的運(yùn)行時(shí)間, 即使模型剖分變大, 運(yùn)行時(shí)間卻僅有很小的變化. 實(shí)驗(yàn)結(jié)果表明, 本文算法在不同剖分的網(wǎng)格中都能獲得高質(zhì)量的參數(shù)化, 具有較強(qiáng)的魯棒性.

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

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