基于Cosserat理論的頭發動態模擬

摘 要:針對基于彈性桿理論的超螺旋模型中對動力學方程的求解影響頭發模擬實時性的問題，采用Cosserat彈性桿理論對頭發建模，引入角速度變量，并加入頭發運動平衡的固有約束，得到改進的拉格朗日動力學運動方程;然后，將頭發單體離散化，用角速度和四元數作為狀態變量簡化動力學方程，用半顯式的歐拉方程加速方程的執行，大大降低了系統的運行時間，提高了模擬的速度，在滿足實時性要求的同時提高了頭發模擬的真實度。

關鍵詞:Cosserat理論; 頭發建模; 離散; 歐拉方程

中圖法分類號:TP391.9

文獻標志碼:A

文章編號:1001-3695(2010)02-0757-03

doi:10.3969/j.issn.1001-3695.2010.02.0099

Simulation of dynamic hair based on theory of Cosserat

TANG Yonga， GAO Ying-huia， FENG Li-yingb， LV Meng-yaa， LI Xiao-yana

(a.School of Information Science Engineering， b.Library， Yanshan University， Qinhuangdao Hebei 066004， China)

Abstract:The real-time simulation of hair in the super-helix model based-on the elastic rod is a big problem because of complicatedly solving dynamical equation.Applied Cosserat rod theory to model hair. Obtained the equation of motion by combining the angular velocity with the holonomic constraints in the dynamic equilibrium configuration for an elastic rod. Successfully improved the Lagrangian dynamic equation. Further more，discretized the rod into finite element.Brought the angular velocity together with the quaternions to be the state variables. Decoupled the resulting system of equations. It had gotten an efficient numerical solution of Lagrangian.Used semi-explicit Euler equations to speed up the Lagrangian dynamic equation. Finally， increased the velocity of the simulation，improved the reality of the simulation， satisfied the requirement of fast and realistic simulation of hair.

Key words:Cosserat theory; hair model; discretization; Euler equation

0 引言

頭發的模擬在虛擬現實、電子游戲、廣告設計、電影電視制作等領域都有著極其廣泛的應用，但如何恰當地對頭發建模并真實快速地模擬頭發的動態行為一直是計算機圖形學研究的難點之一。

20世紀80年代以來，學者們在頭發的建模和模擬方面做了很多工作。2001年，Hadap等人[1]在流體力學模型的基礎上對頭發建模，通過相鄰段之間位置的變化來模擬頭發運動，由于流體與頭發屬性的固有差異，模擬效果不夠真實。2003年，Bando等人[2]將頭發構造成寬松連接的粒子系統，這種近似于粒子系統的模擬方式結構簡單、計算方便，使繪制速度大幅度提高，但只能實現對直發的模擬。2005年Choe等人[3]擴展了質點彈簧系統，用連接兩個剛體段中點的扭轉彈簧實現頭發的扭轉形變，但復雜控制矩陣的求解降低了模擬的實時性。2006年Volino等人[4]用基于空間自由形態變形的晶格機械模型實時模擬了各種復雜的發型，但沒有對單根頭發建模，真實度不夠。2008年Selle等人[5]采用在四面體中構造高度彈簧的方法，避免了模擬時物體的倒塌，并實現了頭發的扭轉，但由于直接對所有發絲建模，超長的處理時間使得該模型不適合交互模擬。另外，我國的一些學者也對頭發的研究做出了自己的貢獻，但這些方法大都沒有考慮頭發的物理屬性，難以實現真實的頭發運動。

2006年，Bertails等人[6]提出了基于經典基爾霍夫彈性桿理論的超螺旋模型，將頭發絲構造成具有有限自由度的分段螺旋體，并通過拉格朗日力學實現對頭發絲運動的計算。在經典力學基礎上建立的超螺旋模型，由于利用了彈性桿與頭發極其相似的特點而成為迄今為止最精確的頭發建模模型。但是該模型沒有實現頭發的拉伸特性和固有約束，更重要的是，采用半顯式的牛頓方程對動力學方程求解嚴重降低了模擬的實時性;為提高運行效率，選取較少的離散段也影響了模擬的真實性。在前人的工作基礎上，以真實快速模擬頭發為目標，用Cosserat彈性桿理論[7～9]對頭發建模，將頭發的線速度和角速度分開考慮，改進拉格朗日運動方程，用半顯式歐拉方程控制時間步進行，以期獲得真實快速的頭發模擬和簡單有效的求解算法。

1 構造頭發動力學模型

1.1 建立頭發單體模型

根據Cosserat理論，將頭發定義為又長又細的變形體，其長度遠遠大于橫截面半徑。用中心線 r(σ)=(rx(σ)，ry(σ)，rz(σ))T(σ∈[0，1])特征化頭發單體的連續結構，表示出頭發在空間中的位置，圍繞頭發單體建立標準正交坐標系，如圖1所示。圖中，n3(σ)代表曲線的切線方向r′(σ);n1(α)和n2(σ)分別是曲線的法線和副法線。空間導數‖r′(σ)‖的長度表示中心線r(σ)在任一點σ處的拉伸度。引入曲線的Darboux矢量u，滿足n′k=u×nk(k=1，2，3)，則曲線的彎曲度和扭轉度可表示為uk=u#8226;nk(k=1，2，3)。其中，u1和u2是曲線的彎曲度，u3是曲線的扭轉度。為確定uk，用四元數q=(q1，q2，q3，q4)T(qi∈R，‖q‖=1)表示nk，則曲線在局部坐標系中的彎曲和扭轉為

uk=2‖q‖2Bkq#8226;q′;k=1，2，3(1)

1.2 引入曲線的角速度

用(σ)表示質點在任意位置σ處的線速度，那么質點在σ處的角速度ω滿足:k=ω×nk(k=1，2，3) ，則頭發單體橫截面圍繞nk旋轉的角速度為ωk=ω#8226;nk。容易得出角速度在局部坐標系中的表達式為

ωk=2‖q‖2Bkq#8226;;k=1，2，3(2)

角速度的引入為動力學方程的簡化奠定了基礎。

1.3 加入約束的拉格朗日運動方程

拉格朗日運動方程中頭發的動態平衡結構受到兩個完整性約束的限制，即

Cp≡rL-R=0(3)

Cq≡q21+q22+q23+q24-1=0(4)

式(3)是對一端固定、另一端在指定位置的彈性桿的端部約束;式(4)要求四元數具有單位長度，能夠呈現真實的旋轉。通過引入拉格朗日乘子λ和μ，得到加入約束后的拉格朗日運動方程的表達式為

ddtTi-Tgi+Vi+Di+λ#8226;Cpgi+uCqgi=∫10Fedσ(5)

其中:gi∈{rx，ry，rz，q1，q2，q3，q4}是頭發單體在標準正交坐標系中的坐標;Fe是額外力;T、V和D分別為頭發的連續動能、勢能和散逸能。將曲線的線速度和角速度分開考慮后，頭發所具有的各種能量均由平移能和旋轉能兩部分組成，用公式表示為

T=Tt+Tr=12∫10pπr2#8226;rdσ+12∫10∑3k=1Ikkω2kdσ(6)

V=Vs+Vb=12∫10Ks(‖r′‖-1)2dσ+12∫10∑3k=1Kkk(uk-k)2dσ(7)

D=Dt+Dr=12∫10γtv(rel)#8226;v(rel)dσ+12∫10γrω0′#8226;ω0′dσ(8)

其中:p為頭發密度;r為頭發的橫截面半徑;I為慣性張量;Es和E分別為拉伸和彎曲楊氏模量;G為扭轉剪切模量;γt和γr分別為平移和旋轉內摩擦系數。I11=I22=pπr2/4，I33=pπr2/2，Ks=Esπr2，K11=K22=Eπr2/4，K33=Gπr2/4。

2 動力學方程的簡化

對式(5)中頭發平衡狀態的數值計算要求先對該方程進行數值積分，因此，應用有限元方法將頭發離散為單位元，頭發的離散為用角速度和四元數作為狀態變量解決動力學方程奠定了基礎。

2.1 離散頭發單體

將頭發的中心線r(σ)離散為具有N個節點的鏈，每段ri+1-ri的長度可以不等。中心線元素ri+1-ri的方向由四元數q(j∈[1，N-1])j表示，方向元素j起點在中心線元素i的中點，終點在元素i+1的中點(圖2)，則方向的離散空間導數表示為q′j≈1lj(qj+1-qj)，lj=12(‖r0i+2-r0i+1‖+‖r0i+1-r0i‖)。

假設頭發的拉伸剛度很大，那么離散后中心線的空間導數可以近似為r′i≈1li(ri+1-ri)，li=‖r0i+1-r0i‖為中心線元素i的靜止長度。其中:r0i為質點的初始位置。

2.2 簡化動力學方程

拉格朗日運動方程式(5)描述了頭發單體的動態平衡，用離散后得到的單位元素能量取代式(5)中的連續能量，得到離散拉格朗日運動方程，用來描述單位元素的動態平衡。對應的動力學方程體系如下:

M +f(，g)=(Fe

τe)T(9)

矩陣M(，g)包含了={，}和g={r，q}中的非線性項;f(，g)對應于積分勢能和散逸能得到的非線性剛性函數;Fe和τe分別代表頭發所受外力和扭轉。但該方程的求解是相當費時的，難以滿足實時模擬的需要，因此需對該方程體系進行簡化。

式(9)中，M具有分塊對角形式M=Mr00Mq。其中，Mr是質量矩陣，控制質點ri的狀態方程;Mq是四元數qj的質量矩陣。假設質量居中在ri上，則Mr是對角矩陣，可以得到質點的運動方程為

i=Fi/mii=vi(10)

其中:Fi是質點所受內力和外力的合力;vi是質點的速度;mi是質點的質量。四元數qj的質量矩陣Mq不是對角化的，因此用角速度ωj和四元數qj作為狀態變量簡化運動方程:

j=I-1(τj-ωj×Iωj)j=12Qj0ωj(11)

其中:τj是內扭轉和外扭轉的和;I是慣性張量;Qj是qj的四元數矩陣。

通過式(10)和(11)的定義，整個系統的處理過程為:在預處理階段積分勢能、約束能和散逸能得到控制內力和扭轉的剛性函數f(，r，，q)，然后執行運動方程式(10)和(11)，得到曲線在下一個時間步的位置rt+h和方向qt+h。

用Ft表示當前形態下頭發所受內力與外力的和，用τt表示內扭轉，忽略外扭轉，采用半顯式的歐拉方程積分運動方程式(10)，推動質點的平移運動:

vt+hi=vti+hmiFt(12)

rt+hi=rti+hvt+hi(13)

這種設計為系統的執行提供了高效的穩定性。

對旋轉運動方程式(11)的數值積分也用半顯式的歐拉方程解決:

ωt+hj=I-1(tj-ωtj×(Iωtj))h+ωtj(14)

t+hj=12Qtj0ωt+hjh+qtj(15)

qt+hj=t+hj/‖t+hj‖(16)

式(16)對四元數進行顯式規則化，實現約束式(4)。

以上設計加快了方程的執行速率，大大降低了系統的運行時間。

3 碰撞處理

由于頭發數量龐大，碰撞處理成了模擬過程中最耗時的部分。這里只考慮頭發與頭發之間的碰撞檢測和響應，并且允許發生碰撞的頭發段產生一定的重疊。在允許的圓柱體模型中，定義段i與j間允許的沖突深度δij為它們切線方向的函數:

δij=(1+g|ti#8226;tj|)λ(17)

其中:ti和tj分別是段i和j的切線方向;g和λ是常量系數。

當有S個段同時與段i重疊時，采用增加外力的方法來解決碰撞。設δij為段i和j的當前滲透深度，依據δij的值將S分為兩個集合S1和S2。其中，S1由δij≤δij的線段組成，即在允許的極限范圍內的束段;S2由剩下的線段組成，即發生碰撞，必須立即解決的束段。對S1和S2采用不同的方法計算外力。

對S1中的每個段j，計算外力Fwwij，這個力將加到段i所受力上。對于Fwwij，先計算因為碰撞產生的速度變化，然后從速度變化中得到Fwwij。法向方向上:

v1i-v1j=-ε(v0i-v0j)v1i+v1j=v0i+v0j(18)

由式(18)可以得到段i的速度變化量Δvwwij:

Δvwwij=v1i-v0i=12(1+ε)(v0j-voi)(19)

其中:ε為經驗恢復系數ε=δij/δij-1，它隨當前的滲透深度變化。當δij接近δij時，負ε逐漸減少滲透段的速度。

通過Δvwwij可以得到Fwwij:Fwwij=(mi/h)Δvwwij#8226;nij。

對于S2中的段，滲透深度大于極限深度δij，它們必須立即解決。需要增加的外力為Fwwij=mi2δij-δijh2nij。因此，應用到段i的所有外力由重疊段提供的所有力的合力計算，合力表示為

Fwwi=∑j∈SFwwij(20)

最后將式(20)加到式(12)的Ft中，得到頭發因碰撞產生的動態變化。

4 實驗與分析

首先用Cosserat彈性桿理論構造頭發絲單體結構，然后用簡化后的動力學方程計算頭發的運動狀態。實驗的微機配置為AMD Athlon(tm) 64×2 Dual2.11 GHz CPU，1GB內存，GeForce7300LE顯卡; 軟件開發平臺為VC++與開放性圖形庫OpenGL相結合。其中，采用B樣條插值算法實現輔助發絲的生成，用Zinke等人[10]的散射光照模型渲染發絲。在同等規模下，用本文方法和文獻[6]方法分別進行頭發模擬，實驗數據如表1所示。

由表1可知，繪制速度主要受頭發單體數量以及離散段數的影響。與文獻[6]中用半顯式的牛頓方程解決動力學方程方法相比，用角速度和四元數作為狀態變量，用半顯式的歐拉方程執行動力學方程，有效減少了系統運行時間，加速了頭發的模擬，達到快速模擬頭發的要求。

圖3為采用本文方法模擬的頭發運動效果。其中頭發單體數量為50 000根，離散段數為50，獲得的運行時間為0.542fps。實驗中各個參數取值如表2所示。

5 結束語

基于Cosserat彈性桿理論對頭發建模，從頭發的物理屬性上考慮頭發在外力作用下的彎曲和扭轉形變。引入了角速度變量，加入了頭發的固有約束，改進了拉格朗日運動方程，最后用角速度和四元數表示狀態方程，大大簡化了動力學方程的執行。實驗結果表明，動力學方程的有效簡化，使得系統能夠保證在將頭發單體離散為近百個頭發段后仍能滿足快速模擬的需要，提高了動態頭發模擬的實時性和真實性。頭發的運動屬性是極其復雜和無規則的，本文僅利用頭發與彈性桿物理屬性的相似點對頭發進行建模，對于頭發的特有屬性如易產生靜電力等還有必要作進一步研究。

參考文獻:

[1]HADAP S，MAGNENAT-THALMANN N. Modeling dynamic hair as a continuum[J].Computer Graphics Forum， 2001，20(3):329-338.

[2]BANDO Y， CHEN Yu-bing，NISHITA T.Animating hair with loosely connected particles[J].Computer Graphics Forum， 2003， 22(3):411-418.

[3]CHOE B， CHOI M， KO H S， et al. Simulating complex hair with robust collision handling[C]//Proc of ACM SIGGRAPH Symposium on Computer Animation. New York:ACM Press， 2005:153-160.

[4]VOLINO P， MAGNENAT-THALMANN N. Real-time animation of complex hairstyles[J] . IEEE Trans on Visualization and Computer Graphics， 2006，12(2):212-223.

[5]SELLE A， LENTINE M， FEDKIWR， et al. A mass spring model for hair simulation[J]. ACM Trans on Graphics，2008，27(3):64.1-64.11.

[6]BERTAILS F， AUDOLY B， CANI M P， et al. Super-helices for predicting the dynamics of natural hair[C]// Proc of ACM SIGGRAPH Conference on ACM Computer Graphics and Interactive Techniques.New York:ACM Press， 2006:1180-1187.

[7]DICHMANN D J， MADDOCKS J H， PEGO R L. Hamiltonian dynamics of an elastica and the stability of solitary waves[J].Archive for Rational Mechanics and Analysis，1996，135(4):357-396.

[8]SPILLMANN J， TESCHNER M. Cossetat nets[J].IEEE Trans on Visualization and Computer Graphics. 2009， 15(2):325-338.

[9]劉延柱.彈性細桿的非線性力學[M].北京:清華大學出版社，2006:137-213.

[10]ZINKE A， YUKSEL C， WEBER A， et al.Dual scattering approximation for fast multiple scattering in hair[J]. ACM Trans on Grap-hics，2008，27(3):1-32.