吳 璇，張舉勇

(中國(guó)科學(xué)技術(shù)大學(xué) 數(shù)學(xué)科學(xué)學(xué)院，安徽 合肥 230026）

0 引言

三維模型是一種使用三維曲面來(lái)表述物體的三維數(shù)據(jù)，網(wǎng)格是三維模型中一種應(yīng)用廣泛的表達(dá)方式。隨著數(shù)字幾何處理技術(shù)的發(fā)展以及掃描技術(shù)的進(jìn)步，網(wǎng)格模型得以廣泛應(yīng)用于動(dòng)畫、游戲、建筑、醫(yī)療、工業(yè)設(shè)計(jì)等行業(yè)。網(wǎng)格曲面參數(shù)化是流形曲面和參數(shù)域之間的一一映射，是網(wǎng)格處理領(lǐng)域中不可或缺的基礎(chǔ)工具，在網(wǎng)格變形、紋理映射、網(wǎng)格壓縮中都發(fā)揮著重要作用。通常網(wǎng)格是在3D空間中的二維曲面，直接對(duì)于3D模型進(jìn)行網(wǎng)格處理非常復(fù)雜，通過(guò)一一映射到簡(jiǎn)單的參數(shù)域，得到的參數(shù)化結(jié)果與原始網(wǎng)格有相同的拓?fù)浣Y(jié)構(gòu)以及盡可能小的失真，然后在參數(shù)域上進(jìn)行網(wǎng)格處理，極大地降低了處理難度。

一個(gè)高質(zhì)量的參數(shù)化映射f有以下性質(zhì)：無(wú)翻轉(zhuǎn)、低失真度量。無(wú)翻轉(zhuǎn)意味著 det J(f）>0，這里 J(f）是f的雅各比矩陣。理想中的映射是在映射后網(wǎng)格與初始網(wǎng)格之間沒(méi)有形變，但這只是理想情況，一個(gè)高質(zhì)量的網(wǎng)格需要盡量減少形變，而失真度量就是用于衡量映射形變的數(shù)值。

經(jīng)典的參數(shù)化方法主要分為線性方法與非線性方法兩種。線性方法計(jì)算簡(jiǎn)單，可擴(kuò)展性強(qiáng)，因?yàn)榫€性方法通過(guò)計(jì)算一個(gè)線性系統(tǒng)來(lái)得到參數(shù)化結(jié)果。雖然線性方法在計(jì)算效率上占據(jù)優(yōu)勢(shì)，但是有許多方法都必須固定邊界，無(wú)法獲得自由邊界的參數(shù)化結(jié)果，比如針對(duì)拓?fù)鋱A盤，F(xiàn)LOATER M[1-2]通過(guò)把邊界固定到一個(gè)凸多邊形上，同時(shí)所有權(quán)重都保證為正數(shù)，得到一個(gè)無(wú)翻轉(zhuǎn)的參數(shù)化結(jié)果。自由邊界的方法可以通過(guò)虛擬邊界、增添線性方程來(lái)實(shí)現(xiàn)。自由邊界方法通?？梢詼p少固定邊界造成大的變形扭曲，卻不一定確保得到的映射是無(wú)翻轉(zhuǎn)的。非線性方法通常構(gòu)造出一個(gè)以變形能量為目標(biāo)式，包含無(wú)翻轉(zhuǎn)硬約束的全局優(yōu)化問(wèn)題[3-4]，使用牛頓法、高斯牛頓法等優(yōu)化算法降低參數(shù)化網(wǎng)格的變形能量，這些能量函數(shù)描述了參數(shù)化映射后網(wǎng)格的變形、失真程度，通常是高度非線性、非凸的，所以這些方法計(jì)算效率低，而且在處理大型網(wǎng)格時(shí)，非線性方法通常會(huì)隨著所處理網(wǎng)格的增大，收斂速度極大地降低。

為了解決以往參數(shù)化方法運(yùn)算消耗大、運(yùn)算效率低、非并行、可擴(kuò)展性差的缺陷，本文提出了一種可并行、可擴(kuò)展計(jì)算無(wú)翻轉(zhuǎn)、高質(zhì)量參數(shù)化網(wǎng)格的算法。不同于以往算法構(gòu)造出一個(gè)無(wú)法并行的全局優(yōu)化問(wèn)題，本文算法通過(guò)引入輔助變量，把參數(shù)化問(wèn)題分解為每個(gè)面上，每條內(nèi)邊上的局部子問(wèn)題。該算法的空間復(fù)雜度與網(wǎng)格模型規(guī)模成線性關(guān)系，也就是 4N+2|εint|， 其中 N 是網(wǎng)格的 面數(shù) ，|εint|是網(wǎng)格內(nèi)邊條數(shù)。相比于現(xiàn)存算法不可并行性，本文算法最大創(chuàng)新點(diǎn)在于每次迭代都可以并行處理N個(gè)關(guān)于三角面片上映射的子問(wèn)題以及|εint|個(gè)關(guān)于內(nèi)邊相容性約束的子問(wèn)題。實(shí)驗(yàn)顯示相比于現(xiàn)存算法，本文算法最終得到相同甚至更好質(zhì)量網(wǎng)格所需運(yùn)算時(shí)間縮短了至少百倍以上。隨著掃描技術(shù)的飛快發(fā)展，3D網(wǎng)格模型的規(guī)模越來(lái)越大，可擴(kuò)展的網(wǎng)格參數(shù)化算法意義重大。但計(jì)算大規(guī)模網(wǎng)格的無(wú)翻轉(zhuǎn)映射是一個(gè)具有挑戰(zhàn)性的難題，該算法可擴(kuò)展，長(zhǎng)于處理大型網(wǎng)格模型。

1 相關(guān)工作

1.1 無(wú)翻轉(zhuǎn)映射

良好的參數(shù)化結(jié)果需要具有無(wú)翻轉(zhuǎn)的性質(zhì)。TUTTE嵌入及其變體在理論上保證參數(shù)化映射無(wú)翻轉(zhuǎn)，但基本上都會(huì)引入大的失真變形?；诒硎镜姆椒╗4-5]無(wú)法確保消除所有翻轉(zhuǎn)，固定邊界方法總是可以產(chǎn)生無(wú)翻轉(zhuǎn)參數(shù)化映射，但在邊界不適宜的情況下往往帶來(lái)大的失真變形?；诰S護(hù)的方法[6-7]開(kāi)始優(yōu)化失真能量的同時(shí)保持無(wú)翻轉(zhuǎn)。然而，這些數(shù)值算法往往需要解決一個(gè)全局優(yōu)化問(wèn)題，以至于不可并行以及可擴(kuò)展性差。比如HORMANN K等人提出了MIPS方法[7]，可以獲得自由邊界的參數(shù)化結(jié)果，但是MIPS能量計(jì)算耗時(shí)久且容易陷入局部最小點(diǎn)。有界失真方法，比如LIPMAN Y等提出的方法[8-9]直接限制了失真度量的大小，但計(jì)算效率低，而且如何設(shè)置一個(gè)保證無(wú)翻轉(zhuǎn)的合適邊界成為難題。不同于以上方法，采用了一種局部推進(jìn)策略來(lái)解決復(fù)雜的網(wǎng)格參數(shù)化問(wèn)題，算法是可并行和可擴(kuò)展的。為了產(chǎn)生無(wú)翻轉(zhuǎn)映射，一個(gè)直接的方法是從一個(gè)無(wú)翻轉(zhuǎn)的初始值開(kāi)始優(yōu)化，然后使用一個(gè)防翻轉(zhuǎn)能量，通常防翻轉(zhuǎn)能量在翻轉(zhuǎn)出現(xiàn)時(shí)將會(huì)變?yōu)闊o(wú)窮大。這樣優(yōu)化問(wèn)題的結(jié)果會(huì)保證屬于可行域。比如之前的一些研究工作[10-13]，采用的方法中也使用一個(gè)無(wú)翻轉(zhuǎn)的映射作為初始值[12]，如果存在翻轉(zhuǎn)，則會(huì)將其映射到無(wú)翻轉(zhuǎn)的可行域中。在使用ADMM算法優(yōu)化的過(guò)程中，用矩陣Ai來(lái)表示三角曲面i的映射矩陣，在ADMM算法的A子問(wèn)題中，使用牛頓法來(lái)迭代A，在迭代中使用線搜索來(lái)保證每個(gè)三角曲面上的映射矩陣的行列式大于0。這樣保證了算法結(jié)果始終處于無(wú)翻轉(zhuǎn)的可行域中。

1.2 ADMM算法

ADMM是一種強(qiáng)大的優(yōu)化算法，因?yàn)槠溥m合處理大型問(wèn)題以及可以高度并行化，通過(guò)與GPU結(jié)合，在信號(hào)處理、計(jì)算機(jī)視覺(jué)、圖像處理、自動(dòng)控制等領(lǐng)域有著廣泛的應(yīng)用[14-18]。ADMM算法一開(kāi)始用來(lái)處理凸問(wèn)題，實(shí)際上，它在處理非凸問(wèn)題上也有很好的表現(xiàn)[19-21]。該算法把一個(gè)全局的優(yōu)化問(wèn)題分解為小的局部問(wèn)題。通常處理以下形式的問(wèn)題：

就像SA算法，在分散網(wǎng)格后，必須保證網(wǎng)格可以拼接起來(lái)，以及如果想要處理網(wǎng)格變形，還得保證滿足控制點(diǎn)約束。在本文中，通過(guò)引入輔助變量[22]，全局的優(yōu)化問(wèn)題被分解成了每個(gè)三角曲面上包含兩個(gè)子問(wèn)題的優(yōu)化問(wèn)題。這兩個(gè)子問(wèn)題分別對(duì)應(yīng)著網(wǎng)格拼接，能量?jī)?yōu)化兩個(gè)方面。

2 可并行參數(shù)化算法

2.1 問(wèn)題描述

考慮 3D網(wǎng)格 M 上的參數(shù)化映射 M={A1，A2，…，AN}。映射的原始網(wǎng)格是3D三角網(wǎng)格M，該三角網(wǎng)格包含 Nv個(gè)頂 點(diǎn){vi，i=1，… ，Nv}，Ne條 邊{ei，i=1，…，Ne}，N 個(gè)三 角面片{si，i=1，… ，N}，εint是包含三角網(wǎng)格M所有內(nèi)邊的集合。在每一個(gè)面片si上構(gòu)建局部坐標(biāo)系，面片上的點(diǎn)在局部坐標(biāo)系中表示為 v=(x1，x2）∈R2，參數(shù)域上的點(diǎn)表示為 u=(u1，u2）∈R2。 用{Ai，Ti}來(lái)表示面片 si上的仿射變換 Ai，其中Ai∈R2×2是線性變換矩陣，Ti∈R2是平移向量。 面片si上的點(diǎn) v映射后的結(jié)果u為 Aiv+Ti。三角網(wǎng)格M上的映射是無(wú)翻轉(zhuǎn)的當(dāng)且只當(dāng) det(Ai）>0，i=1，…，N。使用MIPS能量來(lái)度量映射誤差。對(duì)于每個(gè)三角形曲面，映射被定義成了一個(gè)仿射變換：

對(duì)于每個(gè)三角曲面，MIPS能量有以下表示：

其中 σ1、σ2是線性變換矩陣 A有符號(hào)的奇異值。優(yōu)化問(wèn)題被構(gòu)造成以下形式：

2.2 ADMM算法

對(duì)于該優(yōu)化問(wèn)題使用ADMM算法來(lái)解決。首先，把所有三角曲面上的變換矩陣Ai連接成一個(gè)矩陣 A∈R2×2N。 然 后 引 入 輔 助 變 量 ：Y∈R2×4|εint|，它 對(duì)應(yīng)于相鄰三角曲面上的線性映射的相容性約束常見(jiàn)符號(hào)如表 1 所示。

表1 常見(jiàn)符號(hào)表

可以把該問(wèn)題重寫為以下形式：

這里 σ1是相容性約束的示性函數(shù)。 SY∈R2N×4|εint|是稀疏的選擇矩陣，則有：

其中λY是對(duì)偶變量，⊙是系數(shù)乘法，μ是懲罰系數(shù)。算法的主流程見(jiàn)圖2。

本算法的初始參數(shù)化映射M0通過(guò)使用簡(jiǎn)單的參數(shù)化算法 Linear ABF[12]獲得，相對(duì)應(yīng)的初始參數(shù)化網(wǎng)格為 M0：=M0(M）。

把每個(gè)面片上的仿射變換作為優(yōu)化變量，構(gòu)造優(yōu)化問(wèn)題，引入輔助變量，應(yīng)用ADMM算法處理優(yōu)化問(wèn)題，獲取最終的參數(shù)化映射Mfinal。

等到算法收斂后，從Mfinal中恢復(fù)出相對(duì)應(yīng)的參數(shù)化網(wǎng)格 Mfinal。

2.2.1 固定 A，更新 Y

圖1 相容性約束示意圖

圖2 算法流程圖

2.2.2 固定 Y，更新 A

其 中 β=μ(|y(i）|），V1=(a，b，c，d）T是 將 Ai的 元 素 列成一列，V2=(d，-c，-b，a）T，sign 是符號(hào)函數(shù)，I是一個(gè) 4×4的單位矩陣，I′有以下表示形式：

使用牛頓法來(lái)解決這個(gè)問(wèn)題，更新流程如下，對(duì)于三角曲面i上的變換矩陣有：

令 α0=1，αk用 αk-1初始化，如果 αk使目標(biāo)式fs(Ai）上升，則乘以0.5來(lái)更新步長(zhǎng)，直到目標(biāo)式下降。為了確保無(wú)翻轉(zhuǎn)，在線搜索的過(guò)程中必須保證det(Ai）>0，在 線 搜 索 中 目 標(biāo) 式 fs(Ai）在det(Ai）＜0時(shí) 設(shè) 置 無(wú) 窮大 ，det(Ai）>0 時(shí) fs(Ai）=f(Ai）。

2.2.3 更新對(duì)偶變量

用該公式更新 λY：

2.2.4 初始化與終止條件

上述算法需要給出每個(gè)三角曲面上的初始變換矩陣。使用一些參數(shù)化方法來(lái)計(jì)算原始網(wǎng)格的參數(shù)化結(jié)果來(lái)得到初始的參數(shù)化結(jié)果，本次實(shí)驗(yàn)中主要采用的方法是Linear ABF參數(shù)化方法，除此以外還可以使用TUTTE嵌入。認(rèn)定當(dāng)||rprimal||2小于0.000 3且||rprimal||2小于1e-5的時(shí)候，ADMM算法收斂。

2.2.5 網(wǎng)格恢復(fù)

在優(yōu)化算法收斂后，將得到每個(gè)三角曲面上的變換矩陣Ai，而網(wǎng)格頂點(diǎn)在參數(shù)域上的坐標(biāo)可以從變換矩陣中恢復(fù)出來(lái)。為進(jìn)行這一系列過(guò)程，從頂點(diǎn)v0出發(fā)，進(jìn)行廣度優(yōu)先遍歷構(gòu)造以下點(diǎn)集：

其中 N(·）是一環(huán)相鄰頂點(diǎn)的集合，D0表示包含頂點(diǎn)v0的集合，Di表示距離D0中任意一個(gè)頂點(diǎn)最短路徑正好包含i-1條邊的頂點(diǎn)的集合。所有的集合都根據(jù)從D0中頂點(diǎn)出發(fā)的廣度優(yōu)先搜索來(lái)決定。構(gòu)造過(guò)程為從D0中取出點(diǎn)v0作為根節(jié)點(diǎn)，找到它未遍歷過(guò)的相鄰點(diǎn)，放入 D1中，成為v0的子節(jié)點(diǎn)，然后依次取出D1中的頂點(diǎn)，遍歷未訪問(wèn)的頂點(diǎn)，構(gòu)造D2，這樣依次構(gòu)造點(diǎn)集，直到網(wǎng)格中所有的頂點(diǎn)都被遍歷過(guò)。

依次計(jì)算D0，D1…中頂點(diǎn)的在參數(shù)域中的位置，當(dāng)已經(jīng)計(jì)算出來(lái)Di中頂點(diǎn)va在參數(shù)域中的位置時(shí) u(va），依照以下公式恢復(fù)出 Di+1中 va的子節(jié)點(diǎn)——頂點(diǎn)vb在參數(shù)域中的位置u(vb）：

其中A*是連接兩個(gè)頂點(diǎn) va、vb的邊任一相鄰三角面片上的變換矩陣。當(dāng)Di中所有頂點(diǎn)都被計(jì)算出參數(shù)域中位置之后，再依次計(jì)算Di+1中頂點(diǎn)在參數(shù)域中的位置。

3 實(shí)驗(yàn)結(jié)果

網(wǎng)格參數(shù)化上的實(shí)驗(yàn)環(huán)境是一臺(tái)2.50 GHz的Intel Core i5的筆記本電腦，RAM為 8 GB，GPU是GeForce 940MX。利用CUDA加速時(shí)，使用的是cuda 9.0。核函數(shù)中每個(gè)塊中有64個(gè)線程數(shù)，塊的個(gè)數(shù)則由該子問(wèn)題總共有多少并行計(jì)算的元素來(lái)決定的，例如在 A子問(wèn)題中，三角曲面的個(gè)數(shù)N將分散到?N/64」個(gè)塊上。牛頓法中，線搜索中α=0.5，同時(shí)當(dāng)梯度的范數(shù)▽f小于1e-3時(shí)，終止牛頓法。為了更加公正地對(duì)比不同方法得到的映射的質(zhì)量，使用文獻(xiàn)[24]中的方法來(lái)統(tǒng)一質(zhì)量度量。對(duì)于每一個(gè)三角曲面sj，采用保角損失度量其中 σj，1、σj，2是三角曲面 j上的 映射 Aj的奇異 值。用來(lái)表示所有三角曲面上的保角損失度量的平均值。

3.1 網(wǎng)格參數(shù)化

圖3展示了三角網(wǎng)格參數(shù)化結(jié)果。圖中原始網(wǎng)格之后緊跟的是該網(wǎng)格的參數(shù)化結(jié)果，第一個(gè)是SA算法的結(jié)果，第二個(gè)是本文算法的結(jié)果，使用的初始值都是由Linear ABF算法得到的。在參數(shù)化結(jié)果下面標(biāo)注著算法計(jì)算時(shí)間以及得到的參數(shù)化結(jié)果的平均保角度損失度量。可以看出，本文的算法在bunny與bimba這兩個(gè)例子上展現(xiàn)了極大的優(yōu)勢(shì)，這是因?yàn)檫@兩個(gè)網(wǎng)格相比于male，disk是大型網(wǎng)格（bunny有 1 251 046個(gè)頂點(diǎn)，bimba有 30 268個(gè)頂點(diǎn)），可以看到本文的算法因?yàn)楦叨炔⑿谢木壒剩L(zhǎng)于處理大型網(wǎng)格，在計(jì)算速度上占有極大的優(yōu)勢(shì)。同時(shí)，將本文算法同 SLIM[23]、AMIPS[24]算法進(jìn)行比較，結(jié)果展示在表2中，可以看出本文算法得到高質(zhì)量網(wǎng)格的同時(shí)，運(yùn)算時(shí)間基本縮短了至少百倍以上，比如tooth，SLIM算法所用時(shí)間是本文算法的23 726倍。圖4展示了網(wǎng)格大小遞進(jìn)的gargoyle模型MIPS能量下降趨勢(shì)，證明本文算法的可擴(kuò)展性，大多數(shù)算法在模型點(diǎn)數(shù)超過(guò)150k時(shí)，收斂會(huì)變慢，本文算法卻在網(wǎng)格模型變大為320k個(gè)面時(shí)，保持著近乎一致的收斂速度。

圖3 網(wǎng)格參數(shù)化結(jié)果展示

圖4 不同大小網(wǎng)格迭代趨勢(shì)圖

3.2 不同懲罰系數(shù)

圖5為不同懲罰系數(shù)下網(wǎng)格參數(shù)化的結(jié)果。對(duì)于camel，在懲罰系數(shù)為90 000時(shí)的最終收斂結(jié)果的平均保角度損失度量為1.32，收斂時(shí)間為0.025 s，在懲罰系數(shù)為30 000時(shí)的最終收斂結(jié)果平均保角度損失度量為 1.25，收斂時(shí)間為 0.064 s；對(duì)于 isis，在懲罰系數(shù)為50 000時(shí)的最終平均保角度損失度量為1.118，收斂時(shí)間為 0.003 5 s，在懲罰系數(shù)為 10 000時(shí)的最終平均保角度損失度量為1.099，收斂時(shí)間為0.058 s。注意當(dāng)懲罰系數(shù)過(guò)低的時(shí)候，將出現(xiàn)不收斂的情況。但參考表3，懲罰系數(shù)整體來(lái)說(shuō)有相當(dāng)大的范圍都會(huì)收斂，比如對(duì)tooth進(jìn)行保角映射，在懲罰系數(shù)1 000到50 000都是收斂的。對(duì)于MIPS能量，通?？梢赃x擇10 000。對(duì)于對(duì)稱狄利克雷能量，懲罰系數(shù)通?？梢赃x擇500。

表2 參數(shù)化算法運(yùn)行時(shí)間比較

圖5 不同懲罰系數(shù)參數(shù)化結(jié)果圖

表3 不同懲罰系數(shù)比較

4 結(jié)論

隨著參數(shù)化領(lǐng)域研究的發(fā)展，各種參數(shù)化算法被陸續(xù)提出。但快速產(chǎn)出無(wú)翻轉(zhuǎn)、高質(zhì)量的參數(shù)化映射一直是一個(gè)強(qiáng)烈需求。在本文中，提出一種針對(duì)于三角網(wǎng)格，計(jì)算三角網(wǎng)格的無(wú)翻轉(zhuǎn)、高質(zhì)量參數(shù)化結(jié)果的可并行、可擴(kuò)展方法。該網(wǎng)格參數(shù)化ADMM算法可以高效并行地計(jì)算無(wú)翻轉(zhuǎn)、高質(zhì)量的三角網(wǎng)格映射，在處理大型網(wǎng)格上尤其占有優(yōu)勢(shì)。本文從不同角度、不同例子展示了該算法的優(yōu)勢(shì)。在處理大型網(wǎng)格或者希望快速處理網(wǎng)格時(shí)，都可以從該算法中得到益處。