浙江工業(yè)職業(yè)技術(shù)學(xué)院 李 偉

引言

長(zhǎng)期以來(lái)，研究人員一直尋求布料仿真的高效、穩(wěn)定、實(shí)時(shí)的算法[1,2,3]，目前已提出了多種模型來(lái)解決這一問(wèn)題。這些研究成果廣泛應(yīng)用于3D電影、網(wǎng)絡(luò)游戲虛擬現(xiàn)實(shí)等領(lǐng)域。人物角色一直是網(wǎng)絡(luò)游戲和虛擬現(xiàn)實(shí)中最重要的對(duì)象之一，而服裝又覆蓋人體9O%區(qū)域以上，因此視覺(jué)上真實(shí)穩(wěn)定的動(dòng)態(tài)布料仿真對(duì)網(wǎng)絡(luò)游戲和虛擬現(xiàn)實(shí)的重要性不言而喻，但具有高度真實(shí)感的布料仿真往往計(jì)算成本很高，實(shí)時(shí)性欠缺。

本文采用典型的質(zhì)點(diǎn)-彈簧模型，利用自適應(yīng)混合積分方法來(lái)對(duì)微分方程求解，巧妙運(yùn)用狀態(tài)方程中的局部結(jié)構(gòu)參數(shù)，計(jì)算得出自適應(yīng)時(shí)間步長(zhǎng)。仿真實(shí)驗(yàn)表明，該方法在保證系統(tǒng)穩(wěn)定性同時(shí)加速了求解過(guò)程。

1.改進(jìn)的質(zhì)點(diǎn)-彈簧模型

目前研究人員提出了各種基于物理性能建模的圖形表達(dá)。代表性模型分為兩大類，即粒子系統(tǒng)模型和連續(xù)介質(zhì)模型。本文中采用的布料模型與choi[4]的模型相類似，是一種粒子系統(tǒng)模型。這種模型的突出優(yōu)點(diǎn)是對(duì)彎曲力和壓縮力統(tǒng)一處理，因而得以避免了“后屈曲失穩(wěn)”現(xiàn)象的發(fā)生，從而提供了很有吸引力的穩(wěn)定結(jié)果。

本文將布料抽象為一個(gè)由m×n個(gè)虛擬質(zhì)點(diǎn)組成的網(wǎng)格，質(zhì)點(diǎn)之間用無(wú)質(zhì)量的、自然長(zhǎng)度為零的彈簧相連接。其連接關(guān)系有三種：“結(jié)構(gòu)彈簧”，“剪切彈簧”，“彎曲彈簧”。

網(wǎng)格中每一個(gè)質(zhì)點(diǎn)與其周圍相鄰的四個(gè)質(zhì)點(diǎn)通過(guò)剛性的結(jié)構(gòu)彈簧（Stretch spring）相連。與其對(duì)角線上的質(zhì)點(diǎn)之間則通過(guò)剛性較小的剪切彈簧（Shear Spring）相連。最后，質(zhì)點(diǎn)與其空間上間隔一點(diǎn)相鄰的質(zhì)點(diǎn)間則通過(guò)非線性的彎曲彈簧（bend spring）相連。

布料的物理模型是基于質(zhì)點(diǎn)一彈簧模型的粒子系統(tǒng)，粒子間的連接是通過(guò)彈簧力。布料的類型和運(yùn)動(dòng)取決于粒子間的彈簧力和布料的拓?fù)浣Y(jié)構(gòu)。作用在粒子上的力可分為內(nèi)力和外力，內(nèi)力包括彈簧力，而外力則包括重力、風(fēng)力等各種阻尼力。

1.1 內(nèi)力分析

用fi,y表示質(zhì)點(diǎn)i和質(zhì)點(diǎn)j之間的彈簧作用在質(zhì)點(diǎn)i上的力，這里彈簧遵循h(huán)ooke定律，則有

這里xi,y=xi-xj，L是彈簧原長(zhǎng)，ki,y是彈簧的彈性系數(shù)

該模型計(jì)算質(zhì)點(diǎn)的受力非常簡(jiǎn)單直觀，并且用三種彈簧模擬了三種力：拉伸力、剪切力和彎曲抗力。但該模型的缺陷在于彎曲彈簧的處理，使得布料很難生成褶皺。Choi采用了力矩平衡方程來(lái)求得彎曲抗力，成功解決了壓桿穩(wěn)定性問(wèn)題，同時(shí)給布料帶來(lái)的褶皺效果要比前者好的多，不過(guò)其計(jì)算相當(dāng)復(fù)雜，這樣在要求實(shí)時(shí)計(jì)算的條件下是不可取的，因此采用了更為簡(jiǎn)便的計(jì)算方式對(duì)彎曲抗力進(jìn)行求解。彎曲抗力的計(jì)算公式如下：

注意，上述彎曲力計(jì)算公式表明彎曲彈簧僅僅在彈簧被壓縮（即彎曲發(fā)生）的時(shí)候才起作用。這里的ki,y要比結(jié)構(gòu)彈簧的剛度小的多，因此這個(gè)力對(duì)結(jié)構(gòu)力的影響并不大，只用來(lái)產(chǎn)生布料的有效彎曲，即褶皺。

另外，考慮到結(jié)構(gòu)彈簧和剪切彈簧的一致性，對(duì)其合并處理。最終采取的模型表示方法，其彈簧連接關(guān)系與Choi一致，結(jié)構(gòu)力和剪切力計(jì)算方法類似方程（1）。

這個(gè)模型結(jié)合了provot和Choi兩種代表性質(zhì)點(diǎn)彈簧模型的優(yōu)點(diǎn)，一方面它保持了provot模型的計(jì)算簡(jiǎn)單的特點(diǎn)，另一方面它采用了choi模型的方法，大幅提高了布料的穩(wěn)定特性，可有效生成各種褶皺。

1.2 外力分析

以上給出了布料模型和內(nèi)力的計(jì)算方法。但這是一個(gè)理想化的模型，它忽略了布料運(yùn)動(dòng)過(guò)程中能量的耗散?，F(xiàn)實(shí)的布料運(yùn)動(dòng)會(huì)受到一些外力的影響，比如空氣阻力、摩擦力等等。為了使布料仿真具有足夠真實(shí)感的動(dòng)態(tài)表現(xiàn)，布料物理模型必須考慮到這些外力的影響。為了讓布料運(yùn)動(dòng)更具真實(shí)性，在布料的運(yùn)動(dòng)模擬中添加一個(gè)阻尼模型是有必要的。

阻尼力是用來(lái)抵抗粒子之間剪切彈簧引起的平面內(nèi)運(yùn)動(dòng)，該力取決與相互作用粒子間的相對(duì)運(yùn)動(dòng)速度，我們采用的阻尼模型如下：

這里fd是連接質(zhì)點(diǎn)i,y之間的阻尼力，Cd是阻尼系數(shù)，vi,vy是兩質(zhì)點(diǎn)的運(yùn)動(dòng)速度。

這個(gè)阻尼模型在布料的模擬中廣為采用，在應(yīng)用上也取得了良好的效果，有效地提高了算法的穩(wěn)定性，增強(qiáng)了布料模擬的真實(shí)感。

2.數(shù)值積分

布料受到的作用力與其運(yùn)動(dòng)狀態(tài)的關(guān)系可以由動(dòng)力學(xué)方程來(lái)描述。一般來(lái)說(shuō)，布料粒子系統(tǒng)的二階動(dòng)力學(xué)方程可以表示為：

這里，M是彈簧一質(zhì)點(diǎn)系統(tǒng)的質(zhì)量矩陣，向量x包括了系統(tǒng)中所有質(zhì)點(diǎn)的位置。F代表了所有的非保守力，比如摩擦力、約束力等外力。表示了所有的保守力，比如重力、結(jié)構(gòu)力、剪切力等，每種保守力都關(guān)聯(lián)著一種勢(shì)能，比如重力能、彈性能等。

基于物理的布料變形動(dòng)畫(huà)可以歸結(jié)為沿時(shí)間軸求解布料動(dòng)力學(xué)方程的初值問(wèn)題。動(dòng)力學(xué)方程中位置和速度都是時(shí)間的未知函數(shù)，然而根據(jù)位置、速度和力三者之間的物理關(guān)系，我們可以將上述具有位置二階偏導(dǎo)的動(dòng)力學(xué)方程分解為兩個(gè)具有一階偏導(dǎo)的運(yùn)動(dòng)方程，即:

這里f代表合力，為了獲得布料粒子系統(tǒng)在每一幀的位置和速度，我們需要利用數(shù)值積分方法對(duì)上式進(jìn)行求解。一般來(lái)說(shuō)，只要已知當(dāng)前t時(shí)刻的布料系統(tǒng)中所有粒子的運(yùn)動(dòng)狀態(tài)(xt,vt)，并由方程構(gòu)造顯式、隱式或半隱式積分式，那么就可以計(jì)算出一個(gè)時(shí)間步長(zhǎng)h后布料系統(tǒng)的運(yùn)動(dòng)狀態(tài)(xt+h,vt+h)，在每個(gè)時(shí)間步長(zhǎng)內(nèi)，顯式方法直接計(jì)算，而隱式方法需要迭代求解。由于布料物理模型的數(shù)值積分求解是影響布料仿真系統(tǒng)性能的一個(gè)重要因素，因此對(duì)數(shù)值積分求解器進(jìn)行改進(jìn)和優(yōu)化始終是快速布料動(dòng)畫(huà)研究的一個(gè)重要內(nèi)容。

2.1 混合積分

目前已有布料仿真的實(shí)踐中，顯式方法存在通常穩(wěn)定性問(wèn)題，所有提到的隱式求解方法只有一階，因此精度較差，容易導(dǎo)致運(yùn)動(dòng)緩慢并可能丟失褶皺等細(xì)節(jié)。其他方法甚至不對(duì)系統(tǒng)求解，只是給出一些合理的近似。

當(dāng)然，我們不是只有顯式和隱式兩種積分方法可以選擇，我們也可以把兩種方法結(jié)合到一起同時(shí)運(yùn)用，這就是混合積分法。混合積分法的核心思想就是將常微分方程系統(tǒng)解耦為剛性部分和非剛性部分，對(duì)于剛性部分使用隱式的方法，對(duì)于非剛性部分則使用顯式的方法，從而把隱式法的求解穩(wěn)定性和顯式法的計(jì)算簡(jiǎn)單性很好的結(jié)合起來(lái)。我們將方程（1）所示系統(tǒng)的力拆分為線性部分和非線性部分，考慮到布料在空間運(yùn)動(dòng)的剛性特征，線性力同布料的剛性特征密切相關(guān)，以隱式方法求解；非線性力部分則直接采用顯式解法。

針對(duì)布料的質(zhì)點(diǎn)—彈簧模型，我們對(duì)結(jié)構(gòu)彈簧進(jìn)行隱式求解，而剪切或彎曲彈簧則顯式求解。不妨將系統(tǒng)合力f解耦為f=gI+gH，gI代表剛性部分，gH代表非剛性部分，則有：

由上述方程知系統(tǒng)的剛性部分運(yùn)用后向歐拉法計(jì)算而非剛性部分則采用前后向歐拉法計(jì)算。

半隱式積分法中每一時(shí)間步長(zhǎng)都需要進(jìn)行系統(tǒng)的單步牛頓迭代計(jì)算，因而顯著降低了每一時(shí)間步長(zhǎng)的計(jì)算效率。混合積分法中不再需要顯式的計(jì)算相鄰質(zhì)點(diǎn)間的系數(shù)矩陣Jacobian矩陣。更重要的是，Jacobian矩陣是巨型稀疏矩陣，導(dǎo)致矩陣的向量積計(jì)算相對(duì)于共軛梯度法的巨大計(jì)算成本降低了很多。

2.2 自適應(yīng)混合積分法

針對(duì)目前布料仿真領(lǐng)域廣泛使用的半隱式計(jì)算框架，在混合積分的基礎(chǔ)上，我們對(duì)其做進(jìn)一步改進(jìn)，提出了一種新的數(shù)值求解方法，即自適應(yīng)混合積分法。我們不再費(fèi)力分析推測(cè)對(duì)系統(tǒng)的哪一部分應(yīng)用何種積分求解，相反我們依賴當(dāng)前模擬參數(shù)和一定的穩(wěn)定性判斷準(zhǔn)則迅速自主靈活決定積分方法的運(yùn)用。另外，在空間上參數(shù)局部變化的區(qū)域，則依賴彈簧的連接類型來(lái)確定。與其他方法不同，我們的方法沒(méi)有犧牲準(zhǔn)確性。與混合積分方法相比，該法計(jì)算簡(jiǎn)單，降低了計(jì)算成本，改善了系統(tǒng)的稀疏性。

一般說(shuō)來(lái)我們對(duì)彎曲彈簧顯式求解，但是彎曲彈簧剛度變化范圍幅度通常很大。當(dāng)模擬布料對(duì)象具有較小的彎曲剛度時(shí)，這種顯式求解方法的結(jié)果很不錯(cuò)。但是當(dāng)模擬的布料對(duì)象具有大剛度特征時(shí)，這種處理方法就顯得不合時(shí)宜了，此時(shí)隱式求解就顯得更加合理。

設(shè)定一剪切或彎曲彈簧，定義其剛度為ks、阻尼系數(shù)為kd，原始長(zhǎng)度為L(zhǎng)。理想狀態(tài)下我們期待建立一個(gè)穩(wěn)定性準(zhǔn)則使我們能夠在模擬過(guò)程中自動(dòng)確定合適的時(shí)間步長(zhǎng)如何對(duì)系統(tǒng)進(jìn)行解耦。另外，針對(duì)自適應(yīng)模擬技術(shù)，即每一空間網(wǎng)格間距h是變化的，或者說(shuō)局部網(wǎng)格參數(shù)m,L,ks,kd均是可變的，我們需要一個(gè)在時(shí)間、空間上局部均可用的判斷準(zhǔn)則。提出的一個(gè)自適應(yīng)的時(shí)間步長(zhǎng)計(jì)算方法如下：

圖1 彈簧間的連接關(guān)系模型

圖2 不同數(shù)值積分方法計(jì)算速度比較

圖3 風(fēng)力作用下的布料

在任一時(shí)間步長(zhǎng)，利用當(dāng)前空間網(wǎng)格間距h和m,ks,kd取值，我們計(jì)算某一類型彈簧連接的一對(duì)粒子作用在其上的力時(shí)，我們利用上式進(jìn)行判斷，如果判斷為真，則跳過(guò)相關(guān)Jacobian矩陣的計(jì)算，否則就如往常一樣，計(jì)算相應(yīng)的Jacobian矩陣。這樣我們就可以主動(dòng)控制系統(tǒng)的解耦計(jì)算而不是被動(dòng)決定。

此外，該準(zhǔn)則僅對(duì)線性化結(jié)構(gòu)或剪切彈簧有效，對(duì)小剛度的彎曲彈簧，我們一般用顯式方法處理，不會(huì)遇到穩(wěn)定性問(wèn)題。

實(shí)驗(yàn)結(jié)果表明，我們的自適應(yīng)混合積分法是穩(wěn)定的，這一點(diǎn)和標(biāo)準(zhǔn)的半隱式積分法幾乎不分上下。比較這兩種算法的計(jì)算效率，自適應(yīng)混合積分法大約降低了2O%左右的計(jì)算時(shí)間，且隨著網(wǎng)格質(zhì)點(diǎn)數(shù)增多計(jì)算成本節(jié)約更多。

3.實(shí)驗(yàn)結(jié)果和分析

基于上述方法，在Windows XP sp2操作系統(tǒng)和Intel pentium 4 2.4G處理器的軟硬件環(huán)境下，在Visual C++ 6.O平臺(tái)并結(jié)合Open GL圖形庫(kù)，對(duì)布料進(jìn)行模擬。

為了測(cè)試本文局部自適應(yīng)混合積分方法和其他常用數(shù)值積分方法對(duì)計(jì)算速度的提升差異，我們利用一1O×1O布料網(wǎng)格做單擺運(yùn)動(dòng)仿真實(shí)驗(yàn)，實(shí)驗(yàn)中分別應(yīng)用顯式龍格-庫(kù)塔法（Runge-Kutta法）、半隱式歐拉法結(jié)合共軛梯度法（CG法）以及本文方法各自模擬1OO幀，記錄總的計(jì)算時(shí)間，從而求出生成每幀的平均計(jì)算時(shí)間，利用這個(gè)參數(shù)來(lái)衡量不同數(shù)值方法的計(jì)算速度。另外根據(jù)仿真中總步長(zhǎng)數(shù)還可以求解得到平均單步計(jì)算時(shí)間。實(shí)驗(yàn)結(jié)果如圖（2）所示。由實(shí)驗(yàn)結(jié)果可知，盡管顯式方法單步計(jì)算時(shí)間比本文方法要小，但考慮到本文方法在時(shí)間步長(zhǎng)上要遠(yuǎn)大于它，因此本文方法的速度還是會(huì)達(dá)到超過(guò)顯示方法，布料規(guī)模越大這種實(shí)時(shí)性能提升的效果越明顯。

圖3給出了16×12個(gè)質(zhì)點(diǎn)，1O14根彈簧構(gòu)成的布料旗幟在風(fēng)力作用下效果。

實(shí)驗(yàn)結(jié)果表明，自適應(yīng)混合積分法保持了很好的數(shù)值穩(wěn)定性，仿真系統(tǒng)始終保持穩(wěn)定，無(wú)散亂失真現(xiàn)象。通過(guò)仿真過(guò)程中自適應(yīng)混合積分方法的使用，能夠有效而逼真地動(dòng)態(tài)仿真織物，不僅保證了精度，同時(shí)提高了計(jì)算效率。

[1]Ari Stern and Mathieu Desbrun.Discrete geometric mechanics for varia tional integrators(a primer for cs).ACM SIGGRAPH ’06 Course Notes on Discrete Di erential Geometry,2006.

[2]R.Bridson,S.Marino,and R.Fedkiw.Simulation of clothing with folds and wrinkles. In ACM SIGGRAPH/Eurographics Symposium Computer Animation,pages 28-36.ACM Press,2003.

[3]E.Boxerman and U. Ascher,“Decomposing cloth,”in Proc.ACM SIGGRAPH/Eurographics Symp.on Comput.Anim.,pp.153-161,2004.

[4]K.Choi and H.Ko.Stable but responsive cloth.In Proceedings of the 29thannual conference on Computer graphics and interactive techniques,pages 604-611.ACM Press,2002.

[5]D.Baraf,A.Witkin,and M.Kass.Untangling cloth.In ACM Trans.Graphics,pages 862{870.ACM Press,2003.

[6]J.Teran,S.Salinas-Blemker,V.Ng,and R.Fedkiw.Finite volume methods for the simulation of muscle tissue.Eurographics/ACM Symp.on Comp.Animation,2003.

[7]P.Volino and N.Magnenat-Thalmann,“Implicit midpoint integration and adaptive damping for eff i cient cloth simulation,”Computer Animation and Virtual Worlds,vo l.16,pp.163-175,2005.