金 耀,宋 丹 ,俞成海,馬文娟 ,宋 瀅 ,何利力

1(浙江理工大學(xué) 信息學(xué)院,浙江 杭州 310018)

2(天津大學(xué) 電氣自動化與信息工程學(xué)院,天津 300072)

在網(wǎng)格曲面上進行曲線設(shè)計是進行幾何處理與編輯的一項重要內(nèi)容,且有著廣泛的應(yīng)用背景,如網(wǎng)格處理中的部件分割[1,2]、虛擬手術(shù)中的模擬切割[3,4]、制造業(yè)中的圖案設(shè)計[5-7]以及流形上的計算幾何[8-11]等.這種曲線有多種設(shè)計方式,其中常見的有插值與光順兩種:用戶在網(wǎng)格曲面上輸入幾個控制點或一條初始曲線,系統(tǒng)便生成一條附著于曲面的光滑曲線,使之經(jīng)過這些輸入點(插值)或逼近于初始曲線(光順).換言之,在網(wǎng)格曲面上設(shè)計曲線,需同時滿足流形約束(曲線位于曲面上)、光滑約束、插值點約束(若有)等,而其中的流形約束使其較之于歐式空間的曲線設(shè)計問題更為復(fù)雜,難度更大.

現(xiàn)有針對網(wǎng)格曲面的曲線設(shè)計方法大致可以分為3 類:投影法[12-14]、光順法[1,2,4]以及參數(shù)化法[6,15].投影法通常將光滑的空間曲線迭代地投影到網(wǎng)格流形上,方法簡單、高效,但其迭代過程緩慢且投影步驟往往不太魯棒.光順法即對給定曲線在保持流形約束下進行光順,方法通常較為魯棒,但其效率較低.參數(shù)化法即將局部曲面參數(shù)化到歐式平面,在平面上設(shè)計曲線并將結(jié)果映射回曲面,方法往往魯棒、高效且靈活,但其受限于局部區(qū)域,難以進行大范圍的曲線設(shè)計.

針對上述不足,本文提出了一種高效且實用的曲線設(shè)計方法.該方法的思想是:將流形約束轉(zhuǎn)化為距離約束,用切平面逼近局部曲面,將距離約束松弛為點到切平面的距離,并與光滑約束、插值點約束共同作為優(yōu)化目標(biāo)進行求解.求解時采用“整體-局部”交替迭代的思想,并借助Gauss-Newton 法控制其收斂行為.該方法不僅提高了魯棒性與計算效率,而且具有形狀可控、適用范圍更廣等特點,能夠應(yīng)用于多種場合.

1 相關(guān)工作

曲線設(shè)計是計算機圖形學(xué)與計算機輔助設(shè)計的一個充滿活力且具有較長歷史的研究話題,現(xiàn)有的大部分研究工作著眼于歐氏空間的曲線設(shè)計,其相關(guān)的理論與算法被深入研究并已趨于成熟[16,17].隨著數(shù)字幾何的廣泛應(yīng)用以及工業(yè)界與學(xué)術(shù)界需求的增加,在流形空間,尤其是在曲面上進行曲線設(shè)計,越來越多地引起人們的關(guān)注.但是,針對這類方法的理論與算法未及前者成熟.根據(jù)曲面表示形式的不同,這些方法可分為兩大類:針對參數(shù)曲面的方法與針對網(wǎng)格曲面的方法.前者通常采用參數(shù)化法:由于參數(shù)曲面具有一個天然的參數(shù)域,因此通?？稍趨?shù)域上設(shè)計曲線并將結(jié)果映射至原曲面空間[18-20].但是絕大多數(shù)方法僅考慮具有單個參數(shù)域的參數(shù)曲面,對于形狀與拓?fù)鋸?fù)雜的曲面則不適用.由于在實際應(yīng)用中網(wǎng)格曲面更為廣泛,因此本文僅討論與研究內(nèi)容相關(guān)的針對網(wǎng)格曲面的工作.

網(wǎng)格曲面上的曲線由于其所在的曲面具有分段線性特點,通常表示為一條分段線性的折線段.其中,具有廣泛應(yīng)用背景的離散測地線便是其中一種,且一直是圖形學(xué)相關(guān)領(lǐng)域的研究熱點[9-11].但測地線旨在計算“局部最短”的曲線,而本文所研究的是離散曲面上的“光滑”曲線,因此本文僅考慮網(wǎng)格曲面上自由曲線的設(shè)計.現(xiàn)有方法大致分為以下3 類.

· 投影法

這類方法將流形約束進行松弛,首先生成經(jīng)過插值點的光滑曲線,然后將其迭代地投影到流形上.Hofer 等人[12,13]提出了一種基于能量極小化的樣條曲線插值方法,他們將曲線離散成折線,把問題描述成一個帶流形約束的優(yōu)化問題,其中,目標(biāo)函數(shù)是樣條能量與測地能量的組合;并采用投影梯度法進行求解,即將松弛流形約束得到的空間曲線迭代地投影到曲面上.該方法與流形的表示形式無關(guān),可用于諸如隱式曲面、參數(shù)曲面與網(wǎng)格曲面等,但其投影步驟并不能保證一定能找到理想的投影點,魯棒性較差.同時,由于投影采用了移動最小二乘法(MLS)局部擬合法,效率較低.Pottmann[14]隨后提出了高階樣條能量,并將該方法進行推廣在曲面上擬合曲線.Wallner 和Pottmann[21]分別提出兩種方法,將歐氏空間中設(shè)計細(xì)分曲線的線性細(xì)分法推廣至流形曲面:第1種提出用“測地平均法”替代傳統(tǒng)的線性平均法,但需計算大量的測地線,較為耗時;第2 種是在外圍空間設(shè)計細(xì)分曲線,然后將結(jié)果投影到曲面上,但是投影算法對于復(fù)雜曲面并不魯棒.Morera 等人[22]與Sarlabous 等人[23]均采用測地細(xì)分策略,分別提出測地Bézier 曲線與測地圓錐Bézier 曲線的細(xì)分曲線生成算法.該方法雖然魯棒,但涉及到測地線的計算,時間復(fù)雜度較高.劉斌等人[5]提出了一種測地B 樣條插值方法,通過插值點反求控制點,并運用最近點投影法將其映射到網(wǎng)格曲面上,然后采用插值誤差反向補償法使曲線逐漸逼近曲面.但該方法同樣耗時且投影算法可能存在多解,從而無法處理復(fù)雜曲面.Bock 等人[24]為生成曲邊網(wǎng)格,在線性網(wǎng)格邊上采樣點,將其投影到曲面上,并對投影點運用高階樣條擬合,其中,投影使用了光滑的Phong 法向場,但生成曲線并不要求嚴(yán)格地位于曲面.Panozzo 等人[25]借助Frechet 均值定義,將歐式空間的加權(quán)平均法推廣到流形空間,并用于設(shè)計樣條曲線.他們把網(wǎng)格曲面嵌入到高維空間并在該空間進行加權(quán)平均,然后采用Phong 法向場將其投影到嵌入網(wǎng)格.該方法無需迭代且所得曲線光滑性好,但其高維嵌入步驟極為耗時.綜上所述,投影法直觀且易于實現(xiàn),但在計算效率或魯棒性方面表現(xiàn)較差.

· 光順法

這類方法通常被用于優(yōu)化特征邊,往往首先松弛光滑約束,然后在光滑能量驅(qū)動下對初始曲線進行光順,并保持曲線始終在曲面上.Jung 等人[26]將圖像中的主動輪廓模型用于網(wǎng)格曲面,對分割線進行光順,通過優(yōu)化表示曲線形狀的內(nèi)部能量和網(wǎng)格特征的外部能量迭代地演化曲線.Benhabiles 等人[27]改進了該方法中的外部能量,對其進行平均和歸一化處理,所得到的曲線能夠更好地與特征邊對齊.Kaplansky 等人[2]采用測地流方程光順封閉曲線,通過修改擴散準(zhǔn)則(方程)以適應(yīng)不同需求.Wu 等人[28]與Zhang 等人[29]基于測地曲率流方程演化曲線形狀,得到了較好的結(jié)果,但其效率難以達到交互要求.由于這些方法所生成的曲線均沿著網(wǎng)格邊,雖然對于提取特征邊有益,但難以設(shè)計一般的光滑曲線.Ji 等人[1]也采用了主動輪廓模型,通過構(gòu)造新的蛇形能量,并結(jié)合曲線節(jié)點分裂、移動和移除等拓?fù)洳僮魇骨€能夠穿越網(wǎng)格面片,所構(gòu)造的曲線更為光滑.Lawonn 等人[4]提出了一種網(wǎng)格域曲線的拉普拉斯光順?biāo)阕?通過降低曲線的測地曲率以平滑曲線.該方法能夠控制曲線的光滑度以及與初始曲線的貼近度,但是由于采用了局部光順策略,需要較多的迭代次數(shù)且無法插值控制點.總而言之,這類方法魯棒性較好,所設(shè)計的曲線通常具有良好的光滑性,但是由于在嚴(yán)格的流形約束下進行光順,導(dǎo)致其效率低下,且難以用于設(shè)計插值曲線.

· 參數(shù)化法

這類方法運用平面局部參數(shù)化方法,在參數(shù)域進行曲線設(shè)計,并將結(jié)果映射到原曲面.Lee 等人[15]提出了“幾何蛇形”法演變曲線形狀:他們對初始曲線所在的局部區(qū)域進行參數(shù)化,并在參數(shù)域上運用蛇形能量演化曲線,并將結(jié)果映射回網(wǎng)格空間.進而基于極小化原則與局部顯著性提出“智能剪刀”法,并用于網(wǎng)格切割[30].朱文明等人[31]采用保角映射構(gòu)造參數(shù)域,并在參數(shù)域上生成細(xì)分曲線,最終將結(jié)果映射回網(wǎng)格曲面.劉斌等人提出了基于指數(shù)映射的測地B 樣條曲線設(shè)計方法[6,32],他們將源曲線控制頂點映射到切空間,根據(jù)曲線遷移前后控制頂點法保持坐標(biāo)不變的原則,建立曲線遷移前后控制頂點的對應(yīng)關(guān)系,實現(xiàn)曲線的幾何變換與樣條曲線陣列的構(gòu)造.但該方法僅限于局部區(qū)域,且需頻繁更新參數(shù)域.這類方法可照搬歐式空間的曲線設(shè)計方法,簡單且靈活,然而現(xiàn)有方法大多考慮局部參數(shù)化方法,難以設(shè)計大范圍的曲線,且未曾考慮參數(shù)化形變誤差對結(jié)果的影響.

本文提出的方法屬于投影法,但它將流形約束松弛為點到切平面的距離約束進行求解,克服了上述3 類方法存在的缺點:相比于傳統(tǒng)的投影法,該方法同時優(yōu)化光滑與流形能量,使曲線能夠逐步光順的同時,穩(wěn)步地貼近曲面,從而使投影步驟的魯棒性和可靠性大為提高;相比于光順法,該方法由于松弛了流形約束,優(yōu)化自由度更大,且計算效率更高;相比于參數(shù)化法,該方法能夠處理局部參數(shù)化法難以適用的具有復(fù)雜拓?fù)涞木W(wǎng)格曲面.

2 曲線設(shè)計方法

2.1 問題描述

給定一個嵌入于三維歐氏空間R3的光滑曲面S以及一堆有序點陣{pj}?S,要求一條位于該曲面上且插值于{pj}的光滑曲線c∈S.若采用參數(shù)方程表示曲線:c(t)→R3,并運用3 次樣條能量度量曲線的光滑度[12,15],則求解曲面上的光滑插值曲線問題可轉(zhuǎn)化為計算如下位置與流形約束下的優(yōu)化問題:

其中,x(t)∈S表示位于曲面S的曲線,x″(t)為曲線關(guān)于位置參數(shù)t的二階導(dǎo),其幾何意義與曲線的曲率正相關(guān).

當(dāng)S為網(wǎng)格曲面時,其上的曲線可表示為分段線性的離散曲線c=〈V,E〉,其中,V={q1,q2,…,qn}?S為曲線上的頂點集,E為形成曲線的邊集.

·當(dāng)曲線為開曲線時,E==1,2,...,n-1}?S;

·當(dāng)其為閉曲線時,E==1,2,...,n}?S(%為取模運算).

由于曲線是離散的,公式(1)中曲線的二階導(dǎo)數(shù)x″(t)并不存在,因此需首先對該公式進行離散化,并采用數(shù)值方法進行求解.

2.2 流形約束轉(zhuǎn)化

對于離散網(wǎng)格曲面S,由于其沒有解析表達式,導(dǎo)致公式(1)中的流形約束難以顯式表示,從而給求解該帶約束的優(yōu)化問題帶來較大的困難.為此,本文引入一個距離函數(shù)D:R3×S→R,用于度量空間中任意點到流形S的距離,從而將流形約束轉(zhuǎn)化為距離約束,即等價地表示為曲線上任意點q∈x(t)到曲面S的距離為0:

而距離函數(shù)D則可定義為點到局部曲面的最短距離:

其中,d為點到點的歐式距離函數(shù),B(q,r)為q的半徑為r的球形鄰域.

2.3 離散化

公式(1)中帶約束的優(yōu)化問題可用數(shù)值方法求解,但首先需進行離散化.為方便對二階微分表示的3 次樣條能量進行離散化,類似文獻[12],本文采用弦長參數(shù)化方法計算采樣密度,即對相鄰兩個插值點所形成的線段pi pi+1,確定其采樣點個數(shù)Ni為

其中,D為離散步長,亦即每條線段的長度(默認(rèn)情況下設(shè)置為網(wǎng)格平均邊長的五分之一),用于控制曲線的采樣密度.由此,該曲線可離散成折線段,其有序采樣點為{qj},其中,=pi(ki為插值點在采樣點中的序號).由于該離散方式使得曲線上的采樣點能近似均勻分布,因此公式(1)的目標(biāo)函數(shù)可近似地離散化表示為均權(quán)拉普拉斯能量[33]:

其中,對于開曲線,M表示折線段中除卻首尾端點的離散點個數(shù);而對于閉曲線,M表示所有離散點個數(shù).于是,公式(1)可改寫為如下離散表達式:

為求解該問題,本文運用罰函數(shù)法,結(jié)合距離公式(3),并將公式(5)中硬約束表示的插值點約束與流形約束分別作成軟約束,則公式(5)進一步轉(zhuǎn)化為如下優(yōu)化問題:

其中,

· 參數(shù)ω默認(rèn)取作大的數(shù)值(108),使之近似滿足插值約束;

·λ由用戶控制,用于平衡光滑性與流形約束的滿足程度.

2.4 數(shù)值求解

由公式(7.1)、公式(7.2)所定義的優(yōu)化問題可知:每一對qj與j′q都是彼此相關(guān)的,其非線性程度較高,因此直接求解難度較大.為此,本文計算相鄰插值點之間的Dijsktra 路徑作為初始曲線,根據(jù)公式(4)所得的采樣點數(shù)量在曲線上均勻采樣離散點{qj},并設(shè)置{q′j}={qj};然后借鑒局部/整體交替迭代的優(yōu)化策略[34]進行求解,即將該兩組優(yōu)化變量進行分離,并把每次迭代分解為整體階段和局部階段,直至前后兩次計算得到的曲線差異很小為止,其計算過程如下.

整體階段.固定{q′j},求解{qj}.此時即計算能量函數(shù)(7.1)的極小值.由于該式是關(guān)于{qj}的二次函數(shù),因此可轉(zhuǎn)化為矩陣方程組AQ=B,其中,A為常系數(shù)矩陣.

局部階段.固定{qj},求解{q′j}.此時對每一個qj,求解問題(7.2),即在qj的球形鄰域與網(wǎng)格曲面的交集中搜索最近點.由于計算球體與網(wǎng)格的交集以及在該交集處的最近點較為耗時,因此可采用拓?fù)溧徲蛱娲蝮w幾何鄰域,即將前一次迭代得到的投影點(n為迭代次數(shù))所在三角形面片f0的r-環(huán)拓?fù)溧徲虻拿嫫螰作為交集.在該交集中,搜索qj到該局部區(qū)域的最近點距離:

其中,點到三角形的距離d(qj,fi)可采用文獻[35]中介紹的方法.具體地,以f0作為種子點,以廣度優(yōu)先的方式遍歷F,對每個面片計算該點在面片上的法向投影:若投影點在其內(nèi)部,則將其作為投影點并結(jié)束搜索;否則,計算點到面片的距離并更新最近點,直至遍歷完所有面片.當(dāng)局部區(qū)域非凸時,有可能存在多個最近點,此時可選取使局部拉普拉斯能量,即最小的最近點作為候選點.

易證,上述交替迭代的優(yōu)化策略可使能量函數(shù)(7.1)單調(diào)下降,從而收斂.具體地,在每一次迭代中,其整體階段求解能量函數(shù)極小值能獲得比當(dāng)前能量更低的值;而在局部階段,由于固定{qj},能量函數(shù)(7.1)的前2 項為常數(shù),第3 項由于D(qj,S)≤||qj-q′j||,因此也能得到更低的值.從而對于每一次迭代,該策略均將使能量值下降.

但由于q′j由前一階段曲線采樣點計算得到,而非當(dāng)前qj的投影點,直接用歐式距離||qj-q′j||度量點到曲面的距離并不精確.因此,上述算法雖能收斂,但收斂速度較為緩慢.本文采用該方法在駱駝頭模型上設(shè)計一條經(jīng)過3 個插值點的閉合曲線,如圖1 所示.隨著迭代的進行,曲線變得越來越光滑,但其變化過程緩慢,迭代次數(shù)過多.雖然該算法在理論上是收斂的,且其整體階段僅需高效地回代計算更新曲線坐標(biāo),但由于其局部階段的投影計算耗時相對較多,若迭代次數(shù)過多將嚴(yán)重影響計算效率.在本實驗中,達到最大迭代次數(shù)(1 000)依然未能收斂.

Fig.1 Iteration results and its energy graph with point-to-point distance constraints (λ=0.1)圖1 基于點到點距離約束的迭代結(jié)果與能量圖(λ=0.1)

為此,本文借鑒用于曲面配準(zhǔn)的迭代最近點(ICP)[36]的算法思想,用對應(yīng)點處的切平面逼近局部曲面,并采用點到切平面的距離近似度量點到曲面的距離.設(shè)q′j處的切平面為Tj,其對應(yīng)的法向為N(q′j),則qj到曲面S的距離可近似表示如下:

其示意圖如圖2 所示.

Fig.2 Illustration of the distance from a point to its tangent plane圖2 點到切平面距離的示意圖

從而,公式(7.1)可近似轉(zhuǎn)化為如下優(yōu)化問題:

相對于點到點的距離約束,點到切平面的距離約束由于約束曲線采樣點在近切平面上移動,使得一方面給予其足夠的自由度進行光順,另一方面又能貼近網(wǎng)格曲面,因此能夠較快地實現(xiàn)曲線的光順.

類似基于點到點距離約束的數(shù)值方法,本文同樣采用“整體-局部”交替迭代的策略求解公式(10).然而,若直接采用該策略,將可能導(dǎo)致曲線在整體階段“過度”光順而偏離曲面,且不能保證算法收斂——雖然能量函數(shù)值在整體階段得到下降,但在局部投影階段可能上升.為此,本文借助Gauss-Newton 法[37]的思想,采用局部線搜索控制其收斂行為.利用公式(10)的Hessian 矩陣H(E(q))計算當(dāng)前曲線頂點的迭代方向,并結(jié)合線性回溯搜索法控制迭代步長.即:當(dāng)能量函數(shù)值大于當(dāng)前值時,采用對步長折半的方式保證能量函數(shù)在迭代過程中單調(diào)下降,從而使算法能夠穩(wěn)定收斂.由于公式(10)確定的能量函數(shù)為平方和形式,其對應(yīng)的Hessian 矩陣是對稱正定的,因此其每次迭代的“整體”階段均可轉(zhuǎn)化為凸問題進行求解,無需進行正則化.其迭代求解的每一步通過如下公式更新曲線頂點位置:

其中,ω∈(0,1]為迭代步長,通過折半線性回溯法確定大小.

為了提高計算效率,本文分別從整體和局部階段兩方面進行考慮.整體階段需求解一個線性系統(tǒng),利用公式(11)中Hessian 矩陣的結(jié)構(gòu)在迭代中始終保持不變的性質(zhì),通過對矩陣結(jié)構(gòu)進行預(yù)分解以加速.局部階段需搜索最近點,本文用切平面法向量與局部鄰域面片的交點進行近似,一旦找到交點則結(jié)束,不必遍歷剩余面片.由于交點不一定存在,此時可采用最近點替代.

由上述算法得到的曲線,其每一段離散的折線將非常接近曲面但往往不嚴(yán)格位于曲面上(對位于不同面片的相鄰兩個點形成的線段).為此,本文采用文獻[38]中介紹的投影法將其映射到曲面上:對于每一條未處于網(wǎng)格表面的線段,由一個端點(起始點)向另一個端點(終止點)逐步構(gòu)造切割平面,與相關(guān)的網(wǎng)格邊進行求交,并將交點更新為起始點,逐步往前傳播.如此循環(huán),直至起始點與終止點重疊.具體細(xì)節(jié)可參考文獻[38].

綜上所述,本文提出的基于距離約束的曲線優(yōu)化算法可描述如下.

算法1.基于距離約束的曲線優(yōu)化算法.

3 實驗結(jié)果與討論

本文用C++語言實現(xiàn)了上述算法,并在8 核Intel i7-7700U 處理器、8G 內(nèi)存的Linux 虛擬機環(huán)境下運行.其中,算法所涉及的線性方程組的求解采用Eigen 庫[39]提供的稀疏Cholesky 分解法——LLT 算法.下面,本文通過一系列實驗展示該算法的特點.

3.1 算法的收斂性實驗

本文在人臉模型上設(shè)計了一條經(jīng)過3 個插值點的開曲線.該曲線在迭代20 次后收斂,圖3 給出了迭代的中間結(jié)果.由圖可見:隨著迭代次數(shù)的增加,曲線在曲面上滑動而趨于光滑,最終收斂并達到穩(wěn)定狀態(tài).大量實驗結(jié)果表明:本算法在步長控制策略的作用下均能穩(wěn)定地收斂,且在絕大多數(shù)情況下均在梯度條件下收斂.

Fig.3 Sampled curves on the face model generated during the iterations (λ=0.1)圖3 人臉模型上插值曲線的迭代序列(λ=0.1)

圖4 所示為3 個迭代曲線圖,分別為能量函數(shù)值(公式(10))、對應(yīng)的梯度模長的平方(施加了log 函數(shù))以及平均距離的變化.隨著迭代的進行,能量函數(shù)值逐步下降;其梯度值盡管在局部跳動,但總體亦呈現(xiàn)下降趨勢;曲線到網(wǎng)格的平均距離值也較為穩(wěn)定地下降,且當(dāng)算法收斂時,其平均距離也達到較小值.此時利用投影法將其映射到曲面上,對曲線光滑度的影響較小.

本文算法較為耗時部分為迭代中局部階段的投影點計算.該算法在前次迭代的投影點所在面片的r-環(huán)拓?fù)溧徲蛩阉魍队包c.為加速起見,本文采用法向與鄰域的交點替代最近點以過濾部分面片方式,并比較了兩種策略(r=2),實驗結(jié)果如圖5(a)所示.由圖可見:兩者所設(shè)計的曲線形狀較為相似,幾乎重合,但法向交點策略迭代88 次,耗時104ms,而最近點策略迭代101 次,耗時431ms.因此,前者能在不影響結(jié)果的前提下有效地提升效率.本文同時比較了不同鄰域大小對結(jié)果的影響,如圖5(b)所示(其中的紅點為插值點,下同).由圖可見:1-環(huán)鄰域結(jié)果與其他結(jié)果有一定差異,而2-環(huán)與3-環(huán)鄰域的結(jié)果幾無區(qū)別,三者的迭代次數(shù)分別為126、95 與88 次.由此說明,r對算法收斂快慢有一定影響.這是由于r較小時,曲線移動的步長相應(yīng)減少,因此需要更多的迭代次數(shù)到達極小值點;而當(dāng)r≥2 時,其結(jié)果差異較小,雖然能在一定程度上加快收斂,但同時也意味著將增加局部搜索次數(shù).綜合權(quán)衡算法效率,本文將r默認(rèn)設(shè)置為2.

Fig.4 Generated curves during the iteration on the face model圖4 人臉模型上曲線的迭代變化過程

Fig.5 Comparison results with different strategies of local neighboring search圖5 采用不同局部鄰域搜索策略的結(jié)果比較

3.2 可控參數(shù)λ

在公式(10)中,參數(shù)λ用于平衡曲線的光滑度與流形約束的滿足程度.當(dāng)λ過大時,流形約束能量將占主導(dǎo).此時,算法行為將與光滑法(例如文獻[4])類似.曲線在迭代時易陷入局部極小值,即曲線得不到充分光順,同時也將增大方程條件數(shù),使得數(shù)值性質(zhì)變差.當(dāng)λ過小時,光滑約束能量將占主導(dǎo),此時,算法行為將與傳統(tǒng)的投影法(如文獻[12])相似.雖然曲線得到了充分的光順,但是距離曲面較遠,投影將破壞曲線的光滑度.當(dāng)λ取值適宜時,算法將兼具兩類方法的優(yōu)勢.圖6 展示了插值點相同但λ取值不同情況下,人臉模型上的插值曲線.

Fig.6 Interpolatory curves on the face model with different values of λ圖6 λ取不同值時,人臉模型上的插值曲線

由圖可見:當(dāng)λ取值為1 和10 時,曲線偏離初始曲線較近,未能得到充分光順;當(dāng)其取值為0.1 和0.01 時,曲線偏離初始曲線較遠,但光滑度得到了改善.為了定量比較各曲線的光滑度,本文采用表示曲線曲率的離散拉普拉斯能量,即公式(5)進行度量.λ從0.01 變化到10,各曲線的光滑度分別為[0.0200,0.0133,0.0237,0.0247].通過大量實驗發(fā)現(xiàn),λ=0.1 時獲得的曲線能夠取得較好的光滑度.因此,本文實驗中將參數(shù)λ默認(rèn)取作0.1.

3.3 比較實驗

相比于經(jīng)典的基于能量極小化的投影法[12],本文算法具有眾多優(yōu)點.本文實現(xiàn)了文獻[12]中描述的算法(部分參數(shù)原文未給出,其結(jié)果可能與原文有差異),并與本文方法進行了比較.

· 首先,本文方法效率較高.文獻[12]中的方法首先松弛流形約束,然后逐步將其投影到流形上,其中,投影步驟采用了耗時的移動最小二乘法(MLS)估算局部曲率,使得效率大為降低.圖7(a)所示為兩種方法在雕像上的曲線設(shè)計結(jié)果,其中,本文方法(右圖)耗時117.2ms,文獻[12]中的方法(左圖)耗時281.9ms(迭代22 次收斂);

· 其次,兩者采用同樣的離散拉普拉斯能量作為優(yōu)化目標(biāo),但本文方法得到的結(jié)果通常比文獻[12]的結(jié)果更為光滑.文獻[12]在逐步投影時并未考慮光滑約束,且最終還需將隱式曲面上的曲線投影到網(wǎng)格上,因此在一定程度上影響了其光滑性;而本文方法中光滑步驟與投影步驟交替進行,且由于曲面上大部分采樣點落在網(wǎng)格面上,很多時候?qū)嶋H投影點即為該點本身,因此能夠有效地保持光滑性.圖7(a)所示結(jié)果的光滑性在視覺上的差異并不明顯,為此本文采用公式(5)來度量光滑度并對其進行比較:左圖的光滑度值為2.49-6,右圖的光滑度值為2.76-6;

· 再次,本文方法更為魯棒.文獻[12]在具有尖銳邊特征邊的模型上設(shè)計曲線時,可能因法曲率估算不夠準(zhǔn)確使其步長估算不正確,從而導(dǎo)致投影錯誤,使曲線出現(xiàn)較明顯的鋸齒狀(如圖7(b)左所示),而本文方法不受特征邊影響(如圖7(b)右所示).

Fig.7 Comparison results of the algorithm[12] and our method圖7 文獻[12]算法與本文算法的結(jié)果比較

本文同時比較了文獻[24]中提出的基于Phong 投影的B 樣條設(shè)計方法.

該方法首先運用多維尺度變換法(multidimensional scaling)將網(wǎng)格曲面嵌入到八維歐式空間,然后在該空間設(shè)計B 樣條曲線,并運用Phong 投影法將其映射到嵌入網(wǎng)格上,最后將所得曲線映射回原網(wǎng)格.文獻[24]僅給出了B 樣條的正向設(shè)計方法,即根據(jù)控制點設(shè)計曲線.為方便比較,本文實現(xiàn)了B 樣條的逆向算法設(shè)計插值曲線(采用夾持端點條件)[40],并與本文方法進行比較,其結(jié)果如圖8 所示.由圖8 可見:兩者結(jié)果雖然在曲線外形上存在一定差異,但視覺上均較為光滑.就曲線設(shè)計步驟而言,文獻[24]由于無需迭代,效率較高,在該實例中(曲線的離散點數(shù)目相同)僅耗時1ms;而本文算法耗時6ms.但是該方法需計算點點測地距離與網(wǎng)格的高維嵌入,計算復(fù)雜性很高,在實際應(yīng)用中可能會帶來不便.

相比于最新的基于局部拉普拉斯光順的曲線設(shè)計算法[4],本文算法亦具有較大優(yōu)勢.文獻[4]的算法有一個參數(shù)t∈(0,1],用于控制與初始曲線的逼近程度.為使其具有近似插值效果,將t統(tǒng)一取為1,與本文算法進行比較.首先,本文算法可設(shè)計插值曲線與閉合曲線,而文獻[4]算法則不能,其比較結(jié)果如圖9 與圖10 所示.圖9 比較了兩種方法在瓶子與梨模型上所生成的開曲線.前者算法得到的曲線沒有經(jīng)過所有控制點,而本文算法則能插值于控制點;圖10 比較了兩種方法在腳與狗頭模型上所生成的閉曲線.前者算法由于采用了局部拉普拉斯光順?biāo)惴?無法處理閉合曲線,因此所生成的曲線在首尾連接處出現(xiàn)尖角,而本文算法得到的曲線則在視覺上處處光滑.同時,前者算法所生成曲線的采樣點均在網(wǎng)格邊上,無法做到自適應(yīng)細(xì)分,因此其光滑度受到網(wǎng)格分辨率的影響,如圖9(b)、圖10(b)所示,曲線上的“折角”較為明顯.此外,在多數(shù)情況下,本文算法效率較高.對于圖9、圖10 中的4 個例子,本文算法的運行時間分別為27.8、25.4、28.4 與44.0(ms),而文獻[4]算法則分別需要123、27、68、14(ms).

Fig.8 Comparison results of the algorithm[24] and our method on the wings of the gargoyle model圖8 文獻[24]中算法與本文算法在怪獸翅膀模型的結(jié)果比較

Fig.9 Comarison results for open curves of the algorithm[4] and ours圖9 文獻[4]中算法與本文算法所設(shè)計的開曲線的比較結(jié)果

Fig.10 Comarison results for close curves of the algorithm[4] and ours圖10 文獻[4]中算法與本文算法所設(shè)計的閉曲線的比較結(jié)果

相比于基于參數(shù)化的曲線設(shè)計方法[6,15],本文算法適用范圍廣泛.本文實現(xiàn)了文獻[15]中提出的基于參數(shù)化的幾何蛇形算法,該方法通過優(yōu)化內(nèi)部能量和外部能量的組合演化給定曲線.其中,內(nèi)部能量表示為曲線一階導(dǎo)與二階導(dǎo)能量的線性組合,外部能量與網(wǎng)格特征相關(guān);而其二階導(dǎo)能量即為本文所采用的能量函數(shù)(見公式(5)),區(qū)別在于該算法在參數(shù)域,而本文算法直接在曲面上計算曲線.為在同等條件下進行比較,本文對該算法進行簡化,僅運用二階導(dǎo)能量優(yōu)化曲線,采用最小二乘共形映射(LSCM)[41]在初始曲線的5-環(huán)鄰域計算局部參數(shù)域.由于該優(yōu)化方程沒有外部能量,因此可簡化為線性方程組進行高效求解.圖11 給出了具有2 個虧格的“8 字”模型上所設(shè)計的曲線的比較結(jié)果.由圖11(a)與圖11(b)的結(jié)果可以看出,兩者均能在曲面上設(shè)計光滑曲線.文獻[15]中算法由于無需迭代,在效率上比本文算法更有優(yōu)勢,其僅耗時9ms,而本文算法耗時85.3ms.但是文獻[15]無法處理圖11(c)所示曲線——其局部鄰域需經(jīng)過切割才能進行參數(shù)化,而在切割線處難以保證曲線的連續(xù)性.而本文算法可在任意虧格的曲面上自由地設(shè)計曲線,與模型的拓?fù)鋸?fù)雜性無關(guān).

Fig.11 Interpolatory curve on the eight model圖11 “8 字”模型上的插值曲線

3.4 應(yīng) 用

本文方法不僅能夠設(shè)計插值曲線,且能光順給定的曲線.該問題即等價于設(shè)計逼近于初始曲線c0∈S的光滑曲線.因此僅需對插值曲線模型稍加改造,即將公式(10)中的插值點約束修改成數(shù)據(jù)擬合項,并以軟約束的形式加入到目標(biāo)函數(shù),使得曲線在得到光順的同時,盡可能逼近初始曲線,即.其中,為采樣自c0與qj對應(yīng)的點;α可根據(jù)實際需求(逼近程度)確定,默認(rèn)設(shè)置為0.01.圖12(a)所示為鴨子脖子部位的初始分割線,由于其分割邊界由網(wǎng)格邊依次相連而成,其曲線形狀的鋸齒狀較為明顯;而經(jīng)過本文算法光順后,可得到圖12(b)所示的曲線,視覺效果上變得更為光滑.

Fig.12 Curve smoothing result on the duck model圖12 鴨子模型分割線的光順結(jié)果

本文曲線設(shè)計算法可用于多種場合,圖13 給出了幾個應(yīng)用實例.其中,圖13(a)展示了在魚身模型上交互勾勒五角星圖案的例子;圖13(b)給出了在牛模型的身體部位切割孔洞的例子,可用于虛擬手術(shù),觀察器官內(nèi)部環(huán)境;圖13(c)所示為部件切割的例子,可用于交互分割模型部件.

Fig.13 Applications of our curve design method圖13 曲線設(shè)計的應(yīng)用實例

本文算法最為耗時的部分為迭代中所涉及的線性方程組的求解與局部投影計算.由于該算法在整體與局部階段均采用了加速策略,因此求解效率相對較高.表1列出了在各個模型上設(shè)計曲線的運行時間.由表1可見:算法效率與網(wǎng)格模型的規(guī)模相關(guān)度較小,而主要取決于曲線的離散采樣點的數(shù)目,且總體效率可滿足交互應(yīng)用的要求.

Table 1 Cost time of generating curve on various models表1 各個模型上曲線的生成時間

4 結(jié) 語

本文提出了一種針對網(wǎng)格曲面的高效且實用的曲線設(shè)計方法,該方法將流形約束松弛為點到切平面距離的約束,并將所有約束共同描述成一個優(yōu)化問題進行求解.該方法魯棒、高效,適用范圍廣,可用于任意復(fù)雜拓?fù)?如多虧格)的網(wǎng)格曲面.此外,相比于前沿算法,本文方法不僅能夠設(shè)計閉合與開曲線、插值與光順曲線,而且在效率上占有明顯優(yōu)勢,能夠滿足交互響應(yīng)要求.

本文方法通過采用局部線搜索策略控制迭代步長,能夠保證算法收斂,但是該策略由于對步長進行了限制,有時未能充分松弛曲線形狀,從而可能陷入局部最優(yōu)解.在今后,我們將探索更好的收斂策略,在保證收斂性的同時,更好地控制曲線形狀.此外,本文算法的復(fù)雜度與曲線采樣點的個數(shù)相關(guān),對于設(shè)計采樣密度較大的曲線,其效率將成為一個瓶頸.因此,這也是未來需要研究的工作.