0.04毫米的真實

01-23

最近花了一段時間仔細研究了毛髮模擬和渲染的相關技術,這裡會寫一系列文章分析相關演算法，作為個人的技術筆記，同時簡單科普一下這個看起來比較神秘的小領域。

這裡先羅列一下近十年這個領域最重要的幾篇論文。

首先是毛髮模擬的基本模型，Discrete Elastic Rods，以下簡稱DER。這是目前毛髮動力學中使用的基本模型，用微分幾何的一些方法描述了如何計算離散曲率和扭矩等毛髮物理學特性。包括了兩篇siggraph論文和Miklos Bergou的博士論文。

Discrete Elastic Rods - Columbia University

Discrete Viscous Threads - Columbia University

Phd Thesis

接著是毛髮動力模擬中最困難的部分，滑動摩擦解算。這部分是表現大規模真實毛髮運動的關鍵特性，也就是頭髮自然糾纏的來源。由法國Inria研究所的Gilles Daviet發表，他同時是目前最好的開源滑動摩擦解算器的作者。

Inverse Dynamic Hair Modeling with Frictional Contact

A Hybrid Iterative Solver for Robustly Capturing Coulomb Friction in Hair Dynamics

接下來是毛髮渲染的部分。今年比較重要的一個發表就是Disney對毛髮渲染演算法的改進，直接把毛髮寬度積分幹掉，換成在spp時的隨機積分，渲染速度提高了10倍。這就是瘋狂動物城中大量毛髮渲染的基礎。

An Energy-Conserving Hair Reflectance Model

Importance Sampling for Physically-Based Hair Fiber Models

A Practical and Controllable Hair and Fur Model for Production Path Tracing

Disney在前兩年還發布了一個大規模毛髮模擬的adonis演算法，既保存了毛髮運動的複雜特性，也將模擬速度和穩定度提高了不少。

Adaptive Nonlinearity for Collisions in Complex Rod Assemblies

最後還有一篇關於結合流體粒子與毛髮動力的論文

Coupling Hair with Smoothed Particle Hydrodynamics Fluids

第一部分基本毛髮模型

下面先主要介紹一下DER演算法。作為基本演算法，DER是開創性的，從整體上描述了彈性桿（毛髮）的數學模型和這個模型中曲率、扭矩等物理特性的計算方法。

DER的很多概念來自於微分幾何，本人微分幾何的專業知識有限，只是按照自己的理解做簡單說明，很多地方並不嚴謹（當然不會影響到實現代碼）。

幾何描述

首先說明兩個微分幾何的概念，活動架標和平行移動。

如圖所示，在一條空間曲線上的每一個點，可以建立一組正交的坐標。圖中橙色和藍色的箭頭就表示了曲線上每一個點上的正交坐標的兩個軸，第三個軸就是那一點的曲線切線。這樣隨著曲線變動的坐標稱為活動架標。

想像一下，如果有一個向量（如圖中的綠色箭頭）是用這種坐標描述的，那這個向量順著曲線移動時也按照活動架標的變動而變動。這樣這個向量在局部坐標里其實沒有變動，但因為活動架標在曲線上變動了，所以向量在曲線上的移動也變動了，這就是所謂的平行移動。

平行移動和活動架標將不同點的曲率平移在一起，從而可以在曲線上作向量微分。

當然這裡只是直觀地說明一下,實際上活動架標是通過平行移動計算出來的,而不是相反,下面會講到計算的方法。

這裡要注意一下，上面描述的活動架標是有扭矩的，也就是會順著曲線旋轉。實際上初始狀態的活動架標扭矩為0，如圖中所示的黑白箭頭，這個作為參考架標。實際架標比參考架標多了扭矩，用橙藍箭頭表示，和標轉架標之間有個旋轉角度 $heta$ ，作為彈性桿的第四個自由度。

用計算機來計算參數曲線的積分會比較麻煩，渲染也不能用傳統的多邊形渲染，所以需要把參數曲線分解成離散的點桿結構。

如圖所示，點桿結構的參考架標是在每條連接桿上離散的。計算的方法是選一頭開始，先做任意垂直於桿的箭頭T1，相對應的黑色箭頭T2是T1和桿的切線的叉積。接著將第一組坐標順著桿做平行移動，依次得到剩下的坐標。

這種離散平行移動在計算上其實是一個三位空間旋轉，這個旋轉滿足兩個條件。如圖所示，平行移動的兩個關聯桿為E1和E2，條件1是旋轉作用在E1上能得到E2，條件2是旋轉作用在cross(E1,E2)上得到cross(E1,E2)本身。這個旋轉可以用一個四元數乘法來表示。這樣初始坐標經過每一個節點時做一次旋轉得到下一個坐標，直到末端，就得到了一組參考架標。

動力系統

接下來有了這套幾何工具，我們可以開始描述彈性桿的動力系統了。

上面說到這是一系列連接著的點桿結構。點可以在三位空間變動，所以有xyz三個自由度。桿有一個可以旋轉的旋轉角度 $heta$ 自由度。綜合起來點桿構成了四個自由度的廣義坐標 $q(x^{1},x^{2},x^{3}, heta )$ 。用向量表示整個彈性桿的坐標就是[點，桿，點，桿。。點]。

$[x_{1}, heta_{1},x_{2}, heta_{2}..x_{n}, heta_{n} ]$ 。如果彈性桿有10節，那坐標向量就有11*3+10=43個分量。

整個彈性桿的拉格朗日動力系統可以用如下公式描述

$M(q,t)frac{partial dot q}{partial dot v} =F(q,t)$

$dot q$ 是 $q$ 的時間導數，就是廣義速度。M是廣義質量矩陣。這其實就是拉格朗日版的F=ma

作歐拉積分的時候寫成這樣的形式

$q_{k+1} = q_{k}+hdot q_{k+1}$

$Mfrac{dot q_{k+1}-dot q_{k}}{h}=F(q_{k+1,tk+1})$

$sim F(q_{k},t_{k+1})+hfrac{partial F}{partial q} dot q_{k+1}$

看起來有點麻煩，其實是一個複雜版本的牛頓迭代法。還記得嗎，就是用導數不斷迭代逼近求值。

這裡解這個系統就要求兩個重要的量，一個是系統受力，一個是系統受力的導數（實際是梯度）。然後不斷迭代就能得到系統的變換了。

這裡另外說一個東西，在論文里有時用能量計算，有時用力計算，其實是一樣的。因為彈性桿的受力都是保守力。保守力就是能量函數的負數梯度。論文里還提到了能量的海森矩陣，其實就是負力的雅可比矩陣，道理是一樣的。因為海森矩陣是二次偏導數矩陣，雅可比矩陣是一次偏導數矩陣。概括來說就是能量和力差了一個導數的關係。

$F=- abla E$

$H_{E} = - J_{F}$

整個迭代過程如下

抓重點來講，就是每個迭代過程求兩個量J和g。假設廣義坐標q的自由度是n,g是給定廣義速度 $dot q_{k+1}$ 情況下系統的受力集合，是有n個分量的一個向量。J是力的梯度矩陣，是nxn的矩陣。

這個系統比較複雜一點，因為要計算一個彈性桿的全局變化，而不是局部變化。如果用給單個質點加速度的方法肯定是錯誤的，因為單個質點的運動受到整個彈性桿的制約，牽動附近的質點和桿變動，只有用這種非線性的方法才能求出正確的全局變動。

還是提醒一點，別被這個複雜的公式搞暈了，其實就是牛頓迭代法。只是導數換成了矩陣，不斷地乘矩陣的逆，然後迭代得到新值，逼近真實值。

力的計算

上面把各種概念講完了，接著就集中在如何在上面的系統中計算力和力的雅可比矩陣。

這裡先講彈性桿的內部力，碰撞和滑動摩擦是另一個問題。

內部力基本分為三種，扭曲力，彎折力和拉力。論文中都有詳細描述如何計算。這裡只用扭曲力做一個例子，演示一下如何把論文中的公式變成實際的代碼，計算出力向量和矩陣。

對於力的計算論文中給出來是這樣的，看起來有點懵逼。一個能量的梯度怎麼變成鏈式導數了，其實這裡就是力作用在桿和質點上的不同。

$frac{partial E}{partial heta }$ 是對桿旋轉自由度 $heta$ 的力，可以根據論文中給出的公式計算。但 $frac{partial E}{partial x }$ 是對質點的作用力，需要用鏈式法則求導。

畫個圖來說明一下，如上面所說的，系統的力是n個分量的向量，[點，桿，點，桿...點]的結構。X1是我們當前計算的質點，但它的力同時影響了周圍的11個分量。對兩邊連桿E1和E2的作用力可以用 $frac{partial E}{partial heta }$ 表示，不過E1是正的力，E2相對是負的力。而對於質點X2 X1 X3的力需要用鏈式方法求，就是要再乘上 $frac{partial heta}{partial x }$ 。論文中給出了三個點各自 $frac{partial E}{partial heta }$ 的計算公式。