布料仿真中的自適應各向異性應變限制研究①

尤祖寰,董蘭芳

(中國科學技術大學 計算機科學與技術學院,合肥 230027)

布料仿真一直是計算機圖形學與虛擬現(xiàn)實領域研究的熱門課題.布料仿真技術,旨在通過計算機構(gòu)建虛擬的布料模型,模擬逼真的褶皺及折痕,從而獲得高度的真實感與動態(tài)效果.近年來,隨著虛擬現(xiàn)實的日益火熱,人們更加注重真實感與代入感.同時伴隨著硬件性能的進步,各種大型虛擬現(xiàn)實應用成為可能.因此布料仿真技術獲得了學者們的廣泛關注和進一步研究.

在進行布料仿真的過程中,布料的物理性質(zhì)通過其本構(gòu)模型得以體現(xiàn).本構(gòu)模型旨在盡可能真實地反映現(xiàn)實世界中布料所受應力與應變之間的關系.通常,將布料視為彈性材料,雖然能夠較好地仿真布料的形變特性,但容易造成過度拉伸或壓縮,即“超彈性”現(xiàn)象.應變限制算法一般指對布料質(zhì)點應變情況進行一定限制,在增強真實感的同時防止其發(fā)生“超彈性”現(xiàn)象.

另外,在布料仿真的一些實際應用中,例如虛擬試衣場合,用戶可能對服裝進行局部拉伸操作,如果不加以限制,服裝有可能因所受應力過大而發(fā)生非真實的過度形變.因此,需要對布料合理地施加應變限制約束,增強布料仿真時的非線性塑性特征,防止“超彈性”現(xiàn)象的發(fā)生.

1 仿真流程概述

布料的基本仿真流程如圖1所示.

1)讀入模型:讀入布料模型與障礙物模型,并對讀取后發(fā)生穿透的部分進行簡單的分離處理;

2)添加物理屬性:根據(jù)布料的本構(gòu)模型及相應力學性質(zhì)為布料模型添加物理屬性;

3)模型運動:根據(jù)運動時各節(jié)點受力情況計算各個節(jié)點的速度與位置變化,并更新其狀態(tài);

4)碰撞檢測與響應:根據(jù)模型包圍盒的距離、位置等信息來計算是否發(fā)生碰撞,并進行響應處理.

圖1 布料仿真流程

一般情況下,只要實現(xiàn)了布料的基本仿真流程,即使不做任何優(yōu)化,也可以得到具有一定真實感的仿真布料,其仿真精度與讀入布料的面片數(shù)量、所使用的物理模型、碰撞檢測精度與響應方法均密切相關.

除此之外,在上述所有條件都已經(jīng)確定了的情況下,還有一些可進一步豐富布料細節(jié),增強布料真實感的仿真步驟,如網(wǎng)格重劃、應變限制、塑性/脆性形變處理等等.

2 相關工作

2.1 網(wǎng)格重劃

網(wǎng)格的分辨率大小調(diào)節(jié)在布料仿真中相當重要,因為逼真的布料褶皺細節(jié)需要通過高分辨率的網(wǎng)格進行仿真,然而仿真高分辨率的網(wǎng)格又會大大降低計算效率.因此在布料仿真中,網(wǎng)格重劃無疑是兼顧仿真的實時性和真實感的有效方法.通過網(wǎng)格重劃,仿真系統(tǒng)可以自適應地改變布料模型的局部分辨率.在仿真過程中,表現(xiàn)真實的褶皺細節(jié)時就需要用小而密的三角網(wǎng)格區(qū)域,而在布料較為平坦之處即可使用低分辨率的大塊三角面片.

網(wǎng)格重劃大大提高了仿真效率,同時又能最大程度地保證了應有的布料細節(jié).目前已經(jīng)存在幾種不同的方法,通過使用自適應地網(wǎng)格重劃來進行布料仿真.Hutchinson等人是最早提出網(wǎng)格重劃算法的,其通過曲率作為質(zhì)點-彈簧模型的網(wǎng)格重劃標準[1].Villard等人將類似的方法推廣到了常規(guī)四邊形網(wǎng)格的質(zhì)點彈簧模型中[2].二者當時并未考慮布料褶皺的形成.Simnett等人在前者基礎上增加了壓縮和碰撞作為額外細化標準,進一步地提高了細化效率[3].Lee等人則是結(jié)合質(zhì)點彈簧模型與Loop細分的方法對三角網(wǎng)格模型進行網(wǎng)格重劃,通過預先計算細分步驟,得到多分辨率的層次結(jié)構(gòu),從而自適應地降低線性方程組的維數(shù),更好地進行隱式積分求解[4].Narain等人則通過動態(tài)張量場重劃網(wǎng)格,使用邊的分割,翻轉(zhuǎn)以及塌陷這三種基本操作,已經(jīng)能夠較好地展現(xiàn)出布料的褶皺效果[5].

總而言之,網(wǎng)格重劃主要根據(jù)節(jié)點法向、曲率等參數(shù),對碰撞響應結(jié)束之后的網(wǎng)格進行重劃分,以滿足展現(xiàn)布料微觀形變的網(wǎng)格精度要求.

2.2 應變限制

應變限制是對布料應變量的約束,旨在防止過度拉伸或壓縮的情況發(fā)生.應變限制算法目前可分為直接約束方法與額外約束方法.

早期的研究者們普遍通過使用直接約束方法對應變進行限制.Provot最早提出將應變限制作為構(gòu)建剛性彈簧的方法,對每個環(huán)節(jié)施加最大最小允許應變限制約束[6].Hauth等人通過增大布料剛度,從而盡可能地減少過度拉伸[7],但由于布料剛度的增大,布料數(shù)值求解的精度將無法保證.

綜上所述,當布料受力較小時,采用直接約束處理的方法限制應變可以得到較為精確的結(jié)果,而當布料受力較大時,直接約束方法的效率將大幅降低.因此,有學者提出采用額外約束的方法對應變進行限制.Thomaszewski等人提出基于連續(xù)介質(zhì)的方法,將應變限制應用到有限元方法之中,獨立地約束應變張量的三個部分[8].Narain等人提出上述方法均不適用于各向異性的高度不均勻網(wǎng)格,通過將應變限制視為一個不等式約束非線性規(guī)劃問題,對三角網(wǎng)格進行各向同性的應變限制約束[5].Ma等人在Narain等人的工作基礎上,直接使用布料平面坐標與世界坐標的映射關系提取出經(jīng)緯方向的應變,以此近似奇異值分解而得到的兩個主方向的應變,得到了各向異性的應變限制約束[9].

3 應變限制

3.1 應變限制定義

假設布料厚度及其交織模式相比其彈性行為對布料形變的影響要小得多,因此本文將布料視為一個二維連續(xù)體.將布料離散為三角形網(wǎng)格并進行分析,設三角形頂點i的平面坐標為(ui,vi),世界坐標為xi,如圖2所示.

圖2 三角形網(wǎng)格示意圖

不妨設w(u,v)為將平面坐標映射到世界坐標系的一個線性函數(shù),這里我們并不關心w(u,v)的具體形式.考慮其偏導數(shù),其大小決定了在對應方向上頂點的拉壓程度.根據(jù)上述假設,給定任意一個三角形元素,關于其三個頂點i,j,k,定義Δx1=xj-xi,Δx2=xk-xi,Δu1=uj-ui,Δu2=uk-ui,Δv1=vj-vi,Δv2=vk-vi,則可得到Δx1=wuΔu1-wvΔv1,Δx2=wuΔu2-wvΔv2,即:

式(1)可簡寫為:

其中,X=[xixjxk],Δ為 3×2 的矩陣差分算子,β-1為該三角形有限元基礎矩陣.利用矩陣乘法結(jié)合律,我們可以得到:

一般情況下,平面坐標(ui,vi)是不變的,僅在網(wǎng)格重劃之后發(fā)生改變.因此,可以將視為世界坐標xi的函數(shù).

對式(4)進行奇異值分解:

其中,U的列向量即為兩個應變主方向,Σ中的兩個對角線元素η1、η2即為其主應變率.通過對η1、η2進行限制,即可控制布料模型的應變情況.

3.2 約束非線性規(guī)劃[3]

傳統(tǒng)的求解應變限制問題的方法是通過雅可比迭代或高斯賽德爾迭代進行求解的,但當方程組較大時收斂慢是不可避免的一大問題.另外針對高度不規(guī)則與各向異性網(wǎng)格時,這種求解方法很難保證其穩(wěn)定性.因此,將應變限制問題轉(zhuǎn)化為約束非線性規(guī)劃問題,可以更為有效地處理上述問題,同時提高計算效率.定義最優(yōu)化模型如下:

其中,mn為頂點n的質(zhì)量(即其所在三角形的平均值),為未受約束時所處的位置(即頂點的當前世界坐標).關于兩個主應變率ηj(j=1,2)的約束條件c(x)為:

其中,a為材質(zhì)平面空間下三角形的面積,tj,min為壓縮應變閾值,tj,max為拉伸應變閾值.又ηj關于xn的偏導數(shù)為:

則xi方向上的應變梯度為:

3.3 基于經(jīng)緯方向的應變限制[6]

Ma等人借鑒了如下拉壓能量狀態(tài)函數(shù)[7]

通過將3.2節(jié)中由SVD得到的兩個主應變率η1、η2分別用‖wu(x)‖、‖wv(x)‖進行近似,則可以得到關于經(jīng)緯方向的最優(yōu)化約束條件

又定義剪切向量為ws=wu+wv,類似地,則可以增加兩個關于剪切向量的約束條件

該方法一方面直接使用經(jīng)緯方向的應變率近似SVD得到的主方向應變率,從而避免了SVD的開銷,提高了計算效率,另一方面通過定義剪切向量保留了剪切變形,限制了拉伸應變,從而使得布料能夠不表現(xiàn)出各向同性的行為.

通過該方法施加的應變限制約束,一方面在形變后經(jīng)緯方向無法保證正交,另一方面在每次網(wǎng)格重劃后,三角形網(wǎng)格的經(jīng)緯方向都將發(fā)生變化,如圖3所示.這樣得到的應變限制雖然大體上是各向異性的,但是局部約束方向并不一定總能沿著所希望的方向.為此,本文提出一種方法,即通過內(nèi)積進行映射,使得布料每個三角形網(wǎng)格的局部應變限制在網(wǎng)格重劃之后仍然能夠和全局應變限制保持一致,從而提高了布料仿真的準確性和真實感.

4 內(nèi)積倒數(shù)映射

考慮到布料的針織顯然是具有方向性的,因此經(jīng)向與緯向的應變一般并不相同,然而在每一幀對網(wǎng)格進行重劃分之后,重劃分的網(wǎng)格中每個三角形的主方向均發(fā)生了變化.Ma等人的工作雖然實現(xiàn)了各向異性的應變限制,但無法較好地處理上述情景,因此針對網(wǎng)格中每個三角面片的主方向不固定的特點,需要定義一種新的應變限制來處理這些問題.為此我們提出了內(nèi)積倒數(shù)映射的概念.

為了滿足需求,考慮滿足如下性質(zhì)的一個映射關系:

1)若三角形面片的主方向與布料全局應變限制方向相同,則約束值保持不變;

基于上述性質(zhì),定義映射函數(shù):

其中,ei表示布料的全局應變方向,eu表示三角面片內(nèi)的局部應變單位方向,εi表示ei方向上的全局應變,εu,i表示由ei方向上的全局應變εi映射到三角面片eu方向上的局部應變.

下面確定求解約束非線性規(guī)劃問題中所需的約束條件c(x)及其關于xn的梯度向量?xnc,分別就關于SVD主應變方向與關于布料經(jīng)緯方向的應變限制約束進行討論.

4.1 約束條件的確定

首先由式(13),我們可以得到三角面片應變閾值為

其中,εmin為壓縮應變閾值,εmax為拉伸應變閾值.結(jié)合式(7)和式(14)可得

4.2 梯度向量的確定

要求?xnc,則需要先確定?xnηj與?xneu.由式(9)可得:

為確定?xneu,首先需要確定eu與ei的關系,基于式(5),我們可以得到如下關系式:

其中,V的列向量即是全局應變方向到局部應變主方向的旋轉(zhuǎn)矩陣.

血栓形成可以導致嚴重心腦血管方面疾病，其中肺栓塞、大面積腦梗死可以直接致死和致殘，此類患者往往預后很差。靜脈血栓方面疾病具有一定的發(fā)病率，在住院患者中尤其危重病患者亦有較高的發(fā)病率，可增加死亡率[1] 。深靜脈血栓形成（deep venous thromboembolism，DVT）是最常見的靜脈血栓疾病，在大手術后及危重病患者人群中，其發(fā)病率較高，往往可能導致較嚴重的并發(fā)癥。有研究顯示，可導致肺動脈高壓[2] 等并發(fā)癥。甚至出現(xiàn)意外死亡和嚴重傷殘，這對該病的治療及預防提出了巨大挑戰(zhàn)。下面詳細闡述國內(nèi)外研究進展情況。

因此,從ei到eu的旋轉(zhuǎn)映射為:

借鑒Papadopoulo等人[10,11]的工作,由SVD得到的矩陣V對原矩陣元素wij的偏導數(shù)為:

由鏈式法則,可以得到:

5 改進乘子法

經(jīng)過上一節(jié)的詳細推導,我們得到了本文提出的自適應各向異性應變限制算法的目標函數(shù)、約束條件及其關于自變量的梯度向量.下面則需要對算法得到的約束非線性規(guī)劃問題進行求解.

求解約束非線性規(guī)劃問題,首先想到的是罰函數(shù)法和拉格朗日乘子法.考慮到二種求解算法各有各的優(yōu)點,因此本文將二者相結(jié)合,提出通過改進乘子法來處理該最優(yōu)化規(guī)劃問題.結(jié)合后的改進方法相當于在目標函數(shù)中同時引入了二次懲罰項與線性逼近項.迭代過程中仍是從小到大地依次增大懲罰因子序列,這樣在算法過程早期當懲罰因子很小時,拉格朗日乘子的更新可以很大,即通過引入二次懲罰項使得算法的收斂速度更快,同時又因為拉格朗日乘子的加速更新,一般可以在懲罰因子趨于無窮之前就收斂到最優(yōu)解,從而避免了之前提到的目標函數(shù)海塞矩陣隨懲罰因子增大而失效的情況.

下面具體介紹求解方法,結(jié)合式(6)得到如下約束非線性規(guī)劃問題:

記可行域D={x|x∈Rn,hi(x)≤0,i=1,2,...,m}.若存在不等式約束gj(x)≤0,可加入松弛變量將其轉(zhuǎn)換為等式約束,即:

下面將罰函數(shù)法和拉格朗日乘子法的思想分別引入,即通過拉格朗日乘子和懲罰函數(shù)的組合,將約束條件加入到目標函數(shù)當中,構(gòu)造函數(shù):

其中,拉格朗日乘子λ∈Rm,懲罰因子M∈R.整個過程分兩步:

1)對于固定的λ和M,通過最小化LA更新x.

2)使用前文提及的方法更新拉格朗日乘子:

與懲罰因子M.

對于足夠大的M,迭代一定會收斂,且收斂時λ→λ?,x→x?.

6 實驗結(jié)果

本文的實驗仿真是在Intel(R)Core(TM)i3-2120 CPU @ 3.30 GHz,4 GB內(nèi)存的機器下進行的,采用C++語言進行編寫,并通過OpenGL對所有模型進行渲染與顯示.

為了驗證算法,本文設計了兩個對比實驗,即分別模擬兩端懸掛布料以及仿真服裝,通過觀察每組實驗對象的形變情況,進行算法比較與分析.

6.1 兩端懸掛布料實驗

實驗分為三組,模擬兩端懸掛布料在穩(wěn)定狀態(tài)下的形變情況,分別實現(xiàn)了基于SVD主應變方向的應變限制、基于布料經(jīng)緯方向的應變限制以及本文提出的基于內(nèi)積倒數(shù)映射的應變限制三種方法,具體實驗構(gòu)成參數(shù)如表1所示.其中,圖3(a)(b)(c)垂直方向閾值為[0.99,1.01],水平方向閾值為[0.95,1.05],顯然此時垂直方向要比水平方向具有更高的剛性;圖3(d)、圖3(e)、圖3(f)則與圖3(a)、圖3(b)、圖3(c)相反,垂直方向閾值為[0.95,1.05],水平方向閾值為[0.99,1.01],此時水平方向較垂直方向更具有剛性.

圖3(a)、圖3(d)施加的是基于SVD主應變方向的各向同性應變限制,即使規(guī)定的兩個主應變方向下的限制不同,展現(xiàn)出來的整體布料形態(tài)仍十分相似,僅因為閾值不同而使得細微褶皺上有所區(qū)別,即閾值越小,布料應變越小,其褶皺也越少.

表1 布料實驗參數(shù)

圖3 布料對比實驗結(jié)果圖

圖3(b)、圖3(e)展示的是基于經(jīng)緯方向的應變限制,因為是各向異性的應變限制,所以當經(jīng)向與緯向的應變限制約束閾值不同時,影響的不僅僅是細節(jié)的褶皺部分,其展現(xiàn)出來的整體布料形態(tài)也不相同,即布料應變方向上的形變量受到其閾值的限制.

圖3(c)、圖3(f)施加的是本文所定義的應變限制,與圖3(b)、圖3(e)展示的整體布料形態(tài)比較相似,因為二者均受到各向異性的應變限制.但仔細觀察兩種方法的實驗結(jié)果可以發(fā)現(xiàn),在相同限制閾值的情況下,圖3c、圖3(f)相對于圖3(b)、圖3(e)而言展現(xiàn)了更加豐富的褶皺,以及更大的形變狀態(tài).究其原因在于本文所定義的應變限制考慮到整體閾值到局部的映射,因此圖3(c)、圖3(f)在垂直方向與水平方向上的效果差異要比前面兩種定義更加明顯.

6.2 仿真服裝實驗

將應變限制運用到服裝仿真中,通過三組實驗進一步地驗證了本文方法較基于經(jīng)緯方向的方法限制效果更加準確,實驗參數(shù)如表2所示.如圖4展示了通過三種不同的應變限制約束方法下仿真人偶在靜止狀態(tài)下其身著T恤的形變情況.其中,前兩組水平、垂直方向閾值均為±5%,不同之處在于第一組施加的是基于布料經(jīng)緯方向的應變限制,第二組施加的是本文提出的基于內(nèi)積倒數(shù)映射的應變限制;為進一步地控制實驗參數(shù)變量,第三組施加基于布料經(jīng)緯方向的應變限制,并將其水平、垂直方向閾值設置為±3%.

表2 服裝實驗參數(shù)

顯然地,細微褶皺應變的變化程度較為劇烈,而平坦或較大范圍褶皺區(qū)域的應變變化程度則較緩慢.因此在限制了T恤的應變后,最終應該得到較為平整且僅保留較大范圍褶皺的結(jié)果.可以看到,第一組使用本文方法,在應變限制閾值為±5%的情況下得到了預期結(jié)果,即細微褶皺被如期限制,如圖4(a)所示.

而第二組使用基于經(jīng)緯方向的應變限制方法,在相同的應變限制閾值(±5%)下并未取得預期結(jié)果.如圖4(b)所示,其T恤的胸口、領口處仍存在理應受到限制的細微褶皺.

第三組表明若在施加基于經(jīng)緯方向的應變限制的情況下希望得到上述效果,則需要將應變限制閾值進一步縮小.如圖4(c)所示,在閾值限定為±3%的情況下,基于經(jīng)緯方向的應變限制方法才能較好地保證網(wǎng)格中每個三角形的應變都得到了如期的限制.

圖4 服裝對比實驗結(jié)果圖

綜上所述,同為自適應各向異性應變限制方法,本文所定義的內(nèi)積倒數(shù)映射充分考慮到了整體閾值到局部閾值的映射,使得布料在預期方向上的應變限制效果差異要比基于經(jīng)緯方向的方法更加準確.

7 結(jié)語

本文在基于網(wǎng)格重劃的布料仿真流程中,對布料模型的應變限制添加了一種新的約束,解決了網(wǎng)格重劃后各個三角形形變主軸的局部應變限制與宏觀上規(guī)定的全局應變限制不一致的情況,并使用改進乘子法對約束非線性規(guī)劃問題進行求解.同時,利用OpenGL圖形庫進行顯示與交互,成功地完成了該布料仿真實驗,取得了較其他工作更為科學準確的仿真結(jié)果.

在本文的實驗案例中,內(nèi)積倒數(shù)映射的作用還沒有得到完全地展現(xiàn).在之后的工作中,可以針對具體的布料材質(zhì)可以由力學實驗測得其相應的極限應變,例如布料橫紋(即布料經(jīng)向)、直紋(布料緯向)與斜紋(布料經(jīng)緯45度方向)的拉伸剛度就具有很大的差異,并通過材料力學的相關理論確定其不同方向的應變限制閾值,從而取代人為設定的閾值.再結(jié)合本文所提出的基于內(nèi)積倒數(shù)映射的應變限制,可以精確地仿真出各向異性的布料材質(zhì),從而進一步地提高布料仿真實驗的準確性與真實感.

1Hutchinson D,Preston M,Hewitt T.Adaptive refinement for mass/spring simulations.Boulic R,Hégron G.Computer Animation and Simulation’96.Vienna,Austria:Springer,1996.31-45.

2Villard J,Borouchaki H.Adaptive meshing for cloth animation.Engineering with Computers,2005,20(4):333-341.[doi:10.1007/s00366-005-0302-1]

3Simnett TJR,Laycock SD,Day AM.An edge-based approach to adaptively refining a mesh for cloth deformation.Proceedings of 2009 EG UK Theory and Practice of Computer Graphics.Cardiff University,United Kingdom.2009.77-84.

4Lee Y,Yoon SE,Oh S,et al.Multi-resolution cloth simulation.Computer Graphics Forum,2010,29(7):2225-2232.[doi:10.1111/cgf.2010.29.issue-7]

5Narain R,Samii A,O’Brien JF.Adaptive anisotropic remeshing for cloth simulation.ACM Transactions on Graphics (TOG),2012,31(6):152.

6Provot X.Deformation constraints in a mass-spring model to describe rigid cloth behaviour.Graphics interface.Canadian Information Processing Society.Toronto,Canada.1995.147-147.

7Hauth M,Etzmu? O,Stra?er W.Analysis of numerical methods for the simulation of deformable models.The Visual Computer,2003,19(7-8):581-600.[doi:10.1007/s00371-003-0206-2]

8Thomaszewski B,Pabst S,Stra?er W.Continuum-based strain limiting.Computer Graphics Forum,2009,28(2):569-576.[doi:10.1111/cgf.2009.28.issue-2]

9Ma GH,Ye JT,Li JT,et al.Anisotropic strain limiting for quadrilateral and triangular cloth meshes.Computer Graphics Forum,2015,35(1):89-99.

10Baraff D,Witkin A.Large steps in cloth simulation.Proceedings of the 25th Annual Conference on Computer Graphics and Interactive Techniques.New York,NY,USA.1998.43-54.

11Papadopoulo T,Lourakis MIA.Estimating the Jacobian of the singular value decomposition:Theory and applications.Proceedings of European Conference on Computer Vision.Berlin Heidelberg,Germany.2000.554-570.