先放上效果
驚現塞拉醬
算法是Weta Digital根據siggraph2003的論文加以改進，改進之前使用的是Kajiya and Kay’s 模型，它能量不守恆，也就是說不是基於物理的，不准確
電鏡下真實頭發絲纖維的照片，我們發現上面有很多重疊的角質層叫做毛小皮也叫毛鱗片，他們相對根部的傾斜角度大約為3°，近似模型如下圖
頭發纖維的模型
R為反射（reflection），T為穿透（transmission），也就是折射
這里假設光有三種傳播方式R， TT， TRT
R是直接反射，
TT是經過兩次折射
TRT是穿透（折射）進頭發內部，在頭發內部進行反射，在穿透（折射）出，詳細如圖
人頭發纖維主要由兩部分組成：表皮與皮質。表皮作為外鞘保護內部皮質層，在光散射方面表皮很重要，因為他是折射的邊緣，
頭發核心部有顏色的東西叫做髓質，在皮質與髓質中的色素決定了頭發的顏色，
我們研究頭發的光學屬性有兩種要素：1.折射參數η（大約為1.55）2.橫截面的吸收率σa
幾何圖如上圖，
頭 發方向從根部到尖部，我們假定u是頭發的tan值，向量v和w是一個右手標准正交基，光線入射方向為ωi，散射出的光的方向為ωr，角度分別為 θi，θr（0°垂直於頭發，u為90°）。圍繞着頭發的方位角表示為φi和φr（v為0°w為+90°），差角(θi-θr)/2表示為θd相關方位角 為（φr-φi）/2表示為φd。φr-φi表示為φ。平均值θh = (θi+θr)/2為半角。
頭發的吸光性主要取決於真黑色素（eumelanin）和偽黑色素（pheomelanin），他們的濃度我們分別定義為ρe和ρp，他們的吸收率分別為σa，e和σa，p，光譜吸收率為μa = ρeσa，e + ρpσa，p。
光線離開頭發前變化（折射，反射）次數為p。p ∈ {R = 0, TT = 1, TRT = 2, TRRT = 3, ...}
總反射函數S是所有的縱向散射（longitudinal scattering）函數Mp 和方位角散射函數Np之積的和
我們的最終目標就是把0-3的p所有的S相加，也就是R, TT, TRT , TRRT所有光線之和。
縱向散射（longitudinal scattering）
接下來是高斯Mp縱向散射（longitudinal scattering）函數的求法。
在頭發是完美圓弧的光滑表面的情況下θr = −θi，如圖所示
在光滑的表面反射，形成一個單一的錐形（左圖）。因為表面粗糙，Mp模擬了一個沒有高光的表面，這種像錐形一樣的散射會被縮放，更偏向於法平面（右圖）。光散射范圍更寬，能量也隨之減少。
因為表面粗糙而不能產生完美反射的高光，Marschner使用了半角的高斯函數來求Mp
inline fixed G(fixed beta, fixed theta)
{
return pow(E, -(theta*theta) / (2 * beta * beta)) / (sqrt(2 * PIE)* beta);
}
G函數是一個以縱向傾角θ為參數的標准高斯函數，β是粗糙程度（在縱向的光滑錐形上的標准偏差角），αp是皮質層傾角的一個簡單函數。但是求得的Mp能量不守恆，有如下原因
1. g函數中的θ ∈ {−∞,∞}，但是θh ∈{−π/2,π/2}。使用θh代替θ會使反射能量增倍。
2. 當θi = −θr時，從錐形高光偏光−θi到θr涉及到縮放錐形。在Marschner中只有大約1/cos ²θd
3. 偏轉光照角度移動了相當大的能量導致成角大於θ ∈{−π/2,π/2}（角過大，永遠不會接收到，導致能量損失）
我們使用一種能量守恆的Mp
這種方法在球形重新分配反射角度，使用球形高斯卷積（spherical- Gaussian convolution）得到
G經過一系列的變形化簡得到最終的Mp：
Mp = (1 / ((pow(E, V) - pow(E, -V))*V))*pow(E, (sin(-Theta_i)*sin(Theta_r) / V))*I0(cos(-Theta_i)*cos(Theta_r) / V);
在這里v = β²是粗糙度的平方
I0函數為第一類貝塞爾函數
得到的能量守恆的縱向散射曲線如下
錐形高光是不對稱的，off-specular peak近似於二維BRDF所有粗糙度是能量守恆的
方位角散射（Azimuthal Scattering）
接下來計算方位角散射（Azimuthal Scattering）Np
Marschner利用在發絲是光滑圓情況下的Bravais特性，研究在發絲纖維內部的方位角散射。
修正折射指數η’
方位角改變方向Φ與在發絲中的偏移h ∈ {−1,1}和p值有關，反射模式為：
γi = arcsin(h) ， γt = arcsin(h/η’ )
由於h不可求，所以我們采用另一種可以模擬h的方法
Marschner使用立方近似來模擬h，但是尚未驗證其准確性，當p值較大或反射指數較低時（頭發在水中）
粗糙度方位角散射（Roughened Azimuthal Scattering）函數Np
我們用高斯分布模擬一個變形法平面表面粗糙度的效果。在發絲纖維中的每一個偏移h 產生了一個方位角的連續的分布Dp(φ−Φ(p,h))，這個分布是由φ中的高斯和光滑頭發中的離散方位角Φ(p,h)產生的。如下圖所示：
粗糙離散角散射（TT）：在光滑發絲上的偏移h能產生確定的方位角Φ(p,h)，可以求出方位角的高斯分布Dp(φ−Φ(p,h))。
全部的出射光線集成在一起：
A(p,h)是光線吸收和菲涅爾產生的衰減。我們使用一個新的標准高速分布函數Dp稱之為Gaussian detector
for (fixed i = -1; i < 1; i += 0.1)
{
Mu_a_d = Mu_a / cos(i / Eta_d);
ref_1 = acos(cos(Theta_d) * cos(asin(i)));
f = -_SC + (1 + _SC)*pow(2, -10 * dot(ref_1, _Eta));
T = exp(-2 * Mu_a_d*(1 + cos(2 * asin(i / Eta_d))));
A = pow((1 - f), 2) * pow(f, p - 1) * T;
// Np += A * Dp(Phi - Phi_p_h(p, i, Eta_d));
Np += A * Dp(Phi - Phi_p_h(i, Eta_d));
}
Np *= 0.5;
到此我們所有公式全部完事，我們可以用for循環累加不同p值的S得到最終結果
經過論文中的測量σa,e = {0.419,0.697,1.37}σa,p ={0.187,0.4,1.05}得到結果相對真實。
上圖論文中的實現結果
博主目前能力不夠不能在Unity中弄出這樣的頭發（感覺他這個頭發像是曲面細分生成的），只能用個好一點的模型來湊數了：
本文為部分翻譯的加上博主自身的理解與實現，如有不正確的地方歡迎批評指正
參考：
1. An Energy-Conserving Hair Reflectance Model
2. Light Scattering from Human Hair Fibers
-----by wolf96