基于物理的毛发渲染-Marschner模型
前言：
本文是我在学习Marschner毛发渲染理论时做的一点梳理，同时呢，也有一些对学习过程中遇到的有意思的公式背后的物理原理的心得理解，如果你和我一样出于兴趣或是别的什么原因，刚开始接触毛发渲染理论，那么不妨看看本文，应该会对你有所帮助。PS. 又是忙碌的一个月，各种挤时间，好在利用封控在家省下的通勤时间，于月底前赶完进度，提交了本文，也算落下了一件心事 :)
回顾KajiyaKay模型：
Kajiya Kay发丝模型是一种不透明的单圆柱，在垂直于发丝的朝向（fiber direction）上观察高光走向。该模型提出由平行入射光束所激发的散射中（scattering），高光反射部分（specular）遵循圆锥（cone）分布，圆锥的朝向就是光线沿发丝轴向的朝向，圆锥的中心位于入射光束与发丝表面的接触位置；而漫反射（diffuse）部分则正比于观察角度与发丝表面夹角的余弦。具体可参考图。
图1
下面来说说Kajiya Kay模型的缺点：其一是能量不守恒，也就是说渲染方程并非是基于物理的；其二是无法利用该渲染方程复现出发丝的许多其他光学特性，这主要是因为理论模型本身的缺陷（既不认为发丝是一种透明的电介质，同时也没有考虑到发丝表面的特殊几何构造）
Marschner模型基础：
Marschner等人在KajiyaKay的基础上，结合自己和同行对人体发丝的进一步观察实验，提出了更为全面和准确的物理与数学模型。
图2
简便起见，我们先从解剖学角度切入。人体的毛发可以大体上分为3个主要成分，它们分别是表面的角质层（cuticle），填充内部大量空间的皮质层（cortex）以及最内层的髓质（medulla）。首先说说角质，它是包裹在发丝最外层的是一种组织，很薄，形如瓦片般层层堆叠，互相依附从而构成了发丝的外层几何形态，上图是毛发的经向切面，注意两侧外沿的锯齿形太，它们正是角质层在2D平面上投影的几何形状。相对于发丝的轴向（root -> tip）方向，这种阶梯状的构造会在局部形成很小的角度偏折α，使局部微表面的法线沿着逆发丝朝向（root到tip为正）同样偏移α度。在宏观上观察的结果就是，头发的高光主部（primary specular peak）会比原先预计出现的位置更加靠向发丝根部一些，偏移的程度按照反射定律推算应该在2α度左右。
实验显示，在入射光的另一侧存在第二波瓣（secondary lobe），属于透过发丝上下两个表面后再次进入空气的光的分量。引入第二波瓣可谓革新了KajiyaKay模型对发丝提出的不透明的假设。实际上充斥在角质层内部的大量皮质组织是具有透光性的（特别是对棕色或金色等浅色发丝而言），但是皮质层内部也沉淀有少许色素（pigment），会对某些波长μ的光波产生吸收作用，从而表现出对出射光的染色。说到染色，不得不提在实验中观察到的第二高光主部（secondary specular peak），它们一般位于入射光的同一侧，相较第一高光主部来说，会更靠近发梢一处（靠近tip），宏观上看起来就像一条染了色的次级高光带。论文认为这个第二高光来自于经历了内壁反射后再次折射出的光线分量，具体光路是这样的：首次折射进入皮质的光线在发丝另一个表面上发生了反射（角质与皮质，角质与空气之间），而后再次穿透发丝外壁，折射进入空气。由于多次在发丝内部传播，第二高光的光强会有所削弱，同时被皮质染色。此外论文中还提到黑色头发这条第二高光带不明显（可能是黑色素对光波吸收能力太强所致吧）。
最后一种头发的成分组织叫髓质（medulla），它是位于发丝中心，且更加致密，相比皮质含有更多的色素。一般认为人类发丝的髓质占比较低，因此很多论文在讨论毛发的光照特性时会忽略这部分组织，另一方面，近年来的研究发现，对于自然界中的大部分动物来说，髓质是其毛发的主要组成部分，也是光照计算不可忽略的重要成分，这个后续有空再说。
回到图所示的毛发横剖面，Marschner依据之前的观察与分析，将毛发的散射光归类成三个最主要的成分，分别命名为R项, TT项以及TRT项。其中字母R代表反射（Reflection），字母T则代表折射（Transmission），字母的排列则代表了光路，非常简明形象。我们知道单独一个R代表了入射既反射的第一高光主部；TT则代表发丝反面的透射波瓣分部，因为它来自入射光经过两次折射后形成的散射；最后TRT表示的是第二高光部，暗示其光路轨迹由一次透射（进入皮质层）+ 一次反射（在发丝对侧的角质层前）+ 一次折射（重新进入空气）构成。特别注意的是，由于角质层微表面固有的倾斜角α，会让第一高光和第二高光沿着相反的方向偏移，从而在宏观视觉层面显示出明显分离开来的两种颜色的高光带，这在图中也非常明显的表示了出来。
基础数学建模：
图3
在我们真正深入理解Marschner提出毛发渲染理论前，最好先对毛发与光交互的物理过程有一个比较清晰的数学模型，定义其中的重要符号与变量。如图所示，我们看到右上方有发丝的3D模型及其横截面，一张巨大的矩形平面（横截面或法平面）位于图中间，此外还有许多分布于其上的矢量和角度，具体而言它们是：
u （轴向矢量）：表示沿着发根朝向发梢的方向，对于弯曲的头发，它代表在当前研究的截面所处位置的切线方向。v-w 平面 （法平面）：该平面和方向矢量u合并构成右手坐标系，平面中v轴方向沿着发丝横截面中的长轴方向，w方向沿发丝截面的短轴方向。（发丝截面被认为是椭圆）ωi（光源方向）：定义了计算过程中使用的入射光方向，ωi和u共同决定了入射光平面（图中包括了矢量ωi的矩形边框）θi（入射角，或入射天顶角）：定义了入射光与法平面之间的夹角，该夹角在入射光平面上，从该平面与法平面交线开始，转到光源方向为止。φi（入射方位角）：定义为从法平面v轴转向入射平面与法平面交线为止。ωr （散射光方向）：通常这是我们想要求取的目标矢量，我们对其方向和强度都是感兴趣，不同情况下，它可以用来表示R,TT和TRT中的任意一种出射光矢量。θr （散射角，或散射天顶角）：类似θi，从散射光平面与法平面交线开始，转到散射光矢量为止。φr （散射方位角）：类似φi，从法平面v轴转到两平面交线为止。
如何在宏观上描述这种既包含了透射又存在反射的物理过程呢？通常人们使用BRDF和BTDF来描述在反射与透射光照计算时，出射光亮度与入射光辉度之比，这两个分量合并在一起也有个名字，叫BSDF(双向散射分布函数)，不妨定义为S：
1
其中分子中的L表示沿出射方向单位立体角，单位弧长，单位时间内发散出的光能量（通常也叫光强度）；而分母E则是单位弧长段占有的面积，单位时间内接受到的光能量（也叫辉度）。对于分母E我们还能将其转换为来自各个方向入射光亮度的积分，因此有 dE(ωi) = DL(ωi)cos(θi)dω 的关系，其中D表示发丝直径。综合起来，可以得到如下渲染方程的定义：
2
注意方程的积分区间是全空域（包含了正面的反射+散射，以及背面的透射+散射）。由于毛发直径一般是个常量，该公式中的D一般可以忽略。但是也可以看出，更加粗的毛发会带来更强的反射光泽。
构建散射方程BSDF：
公式（1）和（2）中的散射方程S(ωi，ωr)是一个4维的函数，在不考虑强度的前提下，任何一个方向矢量ω都需要天顶角分量θ和方位角分量φ来联合表示，毕竟我们的要探讨的问题和现象发生在三维空间中的发丝身上，当然期望散射方程S能够在三维空间中正确描述光线在R，TT，TRT分量上的的散布性质。为了降低函数S的复杂度，Marschner一文将三维光场的计算拆分到了2个2D平面上：
3
如式所示，等号右边的第一项是狄拉克函数，只有当θi = -θr时才返回1，其余情况为0。这是基于前人对匀质圆柱体色散研究形成的认知：既光束以入射角度θi照射圆柱体表面时所形成的高光散射只会分布在以-θi为夹角，中心在入射点的圆锥面上（cone surface），可参考图。由此可知，在计算有价值的反射光线时，只需要关注-θi附近的反射分布即可，这消除了维度θr的影响。上式中间项N表示的是发丝垂截面上不同方位角散射光线形成的概率分布，梳理这项稍微复杂点，首先是符号η'，它代表修正后的折射系数（IOR或者Index of Refraction），是关于投影前折射系数η与投影角θ的一个简单函数。因为Bravais' Low告诉我们，对电介质表面来说，如果入射和折射向量被投影到了一个包含法线的平面上（此处为垂截面），其投影向量仍然满足Snell定律，只不过要使用依据投影角修正过的折射系数η'(θ)替换原有的系数η。N项中的另外2个参数分别是前面提到的入射方位角和出射方位角，也都是在垂截面上的度量值，如果假设发丝截面是正圆形，那么对于φi来说，也有相应的公式可以推导出φr的方位，这部分是后话了。最后我们看到等式右边还除以了两个cos(θ)，这是针对对立体角投影的修正，每个cos(θ)各自作用于一个方向上的投影。
4
上式是更加一般的形式，Marschner用M项代替了原先的狄拉克函数，用以专职描述发丝横切面（轴向切面，参考图）上的散射分布；另一方面N项维持形式不变，专职描述垂截面（法平面）上光线散射分布的情况，此处不再赘述。
有意思的地方来了，为了找出N项的正确形式，Marschner一文引出了如下图所示的法平面光路分析：
图4
这张图有几层含义，我们从最简明直接的光路传播开始。首先注意到光线从点L出发的，入射到发丝截面上的一点A后发生了折射和反射。其中的反射矢量沿着AR0返回空气介质，折射矢量则沿着AB射向发丝内壁的另一侧，通过折射形成TT分量，沿BR1返回空气，最后在内壁形成的反射矢量BC则于发丝另一表面C点处形成TRT分量，完成出射任务。我们定义入射角为γi，折射角为γt，那么有以下若干关键角度关系：内壁反射角 ∠ABC = 2γt，形成TT分量的折射出射角为γi（沿矢量BR1方向出射），同样形成TRT分量的入射和出射角亦是γt和γi。
我们定义光源矢量为AL（既从A指向L），顺时针方向为方位角的负方向，那么有：
R：AL - 2γiTT：AL + π + 2*(γt - γi)TRT：AL + π + 2*(γt - γi) + (π + 2γt)
简单解释一下，首次的反射方向AR0，也就是R分量的方向非常容易理解，它等于在光源矢量AL的基础上绕着顺时针方向旋转两倍的γi角度，因为定义了顺时针方向为角度负方向，所以我们得到AR0等于在AL方向上作用一个 -2γi角度的选择。
当发生一次折射时，从图中可以看出光线会从折射前的状态（入射光线LA）沿着顺时针方向（角度逆方向）偏折一个 -(γi - γt) 度，或者记成(γt - γi)，在形成TT前，光线必然经理了2回类似的折射，所以最终TT分量的出射角包含有2*(γt - γi) ，至于上面为何还有一个落单的π？那是因为AL方向是光源方向，要得到入射方向，必须对AL矢量180°掉转朝向。
最后的TRT分量类似上面提及的TT，也经历了两次折射，同时由于包含了一次反射，需要使用一个(π + 2γt)转向操作来描述，其意义在于先对入射光线AB掉头成为BA，然后让BA沿着逆时针（正角度）方向旋转2γt度，怎么样，结合图示很容易理解吧！
Marschner使用了一个非常聪明的公式统一了上式光线路径的转变：
5
它的参数p代表迭代次数，比如对于R来说p = 0，TT时 p = 1，TRT则意味着p = 2，Marschner在论文里认为当 p > 2 时产生的影响可以忽略不计。至于上式中的变量h，其定义是入射光线到圆心的最短距离，可以参考图示中的线段h。由于预设了发丝半径为单位1（unit length)，由三角函数我们可以很容易得到 sin(γi) = h，同时通过Snell定律也能得到 η' * sin(γt) = 1 * sin(γi) = h。由此可知，引入变量h的意义在于将彼此独立的变量γt和γi统一到一个变量名下，以便后续计算。
反过来，如果我们已知入射和出射方位角φi和φr，那么角度变化量 φ = φr - φi也是已知的了，通过解算方程 Φ(p, h) - φ = 0，可以得到h的具体数值。论文中提及，当p = 0 or 1 时，这样的方程求解只有一个根，但是当 p = 2 时，方程会有1个或3个根，而更多得根代表了更多的有效光路（ray path），光路互相叠加会带来能量的不连续，Marschner认为这些多光路的解析解在物理上代表了焦散（caustic）现象，同时也是我们在宏观上感到发丝细节处带有闪烁（glint）的源头。
图中的第二层含义是能量分布的守恒性，既对于一块平坦的微表面而言，所有入射光能等于所有出射光能（不考虑吸收的前提下）。参考如下公式：
6
若定义单位弧长的辉度为 ，在采用了半径为单位1的发丝模型中，单位长度的辐照度E可以写作：E(h) = /2，其中 2 = 2 * 1 unit，既两倍的圆的半径。那么上面双等式中间项表示的自然是对投影到发丝截面“边条”上（此时为1维）的辐射照度E在发丝直径线段上的微分，它等于左侧项所表示的：在以h确立的出射立体角方向上，出射光亮度对方位角的微分。多少有点拗口，但是本质意思就是利用了能量守恒定理来把弧度相关的出射光亮度与弧度相关的入射光辉度联系在了一起（目前还没考虑吸收，当然电介质一般认为吸收不多啦）。
有了上述联系，再加上吸收系数T以及菲涅尔项F就能得到比较完整的N项了。增加后大致上上式可转换为如下形式：
7
当p = 0时，A项的形式可以方便获得：
8
此时没有任何光路经过发丝内部的皮质层，故忽略吸收分量，只余下比较传统的Fresnel项，它是个关于表面介质系数以及入射角度的函数，记作F。
当p > 0时，Marschner给出了如下近似：
9
该式模拟了p次菲涅尔反射的耗损，同时最后一项T对进入发丝内部的光路长度做了估算，并据此对光能进行必要的衰减（模拟吸收），特别得，依据三角函数可近似获取内部光路长度的计算公式，它等于发丝的半径乘以 2 + cos(2γt)，本文随后会详细阐明这个公式的物理意义，此处先略过。
综合起来，完整的垂截面（法平面）上的散射方程子项N可以做如下表示：
10
上式第一项告诉我们，N项将首先按照 p（折射）次数划分为Nr，Ntt以及Ntrt三项。上式中第二项的求和区间为所有可能的光路，此处不用h的原因是因为给定一个h可能会产生多个解（root），而每个解都将对应一条合法的光路，这点前文已经提及过。
解释了略显麻烦的垂切面分量N之后，再来看轴切面分量M会发现非常简单。前文已经说过，基于圆柱的观察结果是所有的散射光线分布在镜面反射方向的一个圆椎上，考虑到角质层特有的几何构造会产生一定量的有向偏移，Marschner等人做了两件简单暴力且有效的事情，第一就是通过实验观察记录轴向散射的具体分布，从而解算出若干具有指标意义的常量；第二，采用归一化高斯分布函数作为基础，采用第一步获取到的参数进行拟合，从而模拟出预期散射分布表达式。
11
具体参考上式，其中函数g为刚刚提及的归一化高斯函数，θh是入射和出射天顶角的半角值（这再次说明圆柱散射只在镜像反射角附近有意义，也就是只需要知道入射或反射角中的一个量即可，此处作者选用了它们二者和的一半来表征），αr表示轴向单次反射波瓣的偏移角，一般为-5°~-10°之间，αtt表示TT的出射波瓣偏移，大约为-αr/2，而最后的αtrt自然是TRT的出射波瓣偏移，被作者定义为大约-3αr/2左右的度数。最后βr，βtt和βtrt表示的是相关标准差（standard deviation）具体物理量亦可通过原文查询，此处不再赘述。
我们来总结一下到目前为止梳理的部分重要内容：Marschner理论认为一束打到发丝上的光线会被分散到以入射点为中心，沿着发丝轴向的一个圆锥面上，而在面上的散射光线可以被归纳为三个主要的分量，它们是来自一次反射的R分量，经过两次透射后从另一侧出射的TT分量，以及透射并在内部上形成反射而后再次出射的TRT分量。对于每一个分量，Marschner将复杂的空间关系拆分到两个2D平面上，它们分别是垂直于发丝的法平面投影以及平行与发丝轴向的横切面投影，各自以N和M表示，控制参数分别是天顶角θ和方位角φ。于是我们有了如下形式的BSDF方程：
12
实时渲染之路
从理论来到实践还有不少路要走，特别是在实时渲染领域，我们没有那么多资源去逐根逐点计算光路，也不会有很多算力可以分给复杂的采样积分运算，所有就必须灵活取舍，找到可接受的近似来加速渲染。
多重散射（Multiple scattering）
这个概念是之前梳理Marshner论文时不曾提及的部分，然而对于非光照模型来说它又至关重要。多重散射简言之就是那些来自发丝与发丝之间的多次弹射并最终落入摄像机的那部分光能。头发并不是一根根孤立的个体，而是一簇簇的发丝聚合，多重散射的本质就是以发丝聚合体为单位所发出的一种漫反射（diffuse），理论上它没有特别得传播方向，理想状况下就如同普通散射一般，多重散射会与发丝簇的宏观法线以及入射光线方向夹角的余弦相关，此外还可以纳入阴影强度 (shadow)，并考虑上头发的基础颜色。业界相对流行的模拟方案是UE4提出的基于Kajiya Kay模型改进来的计算公式：
13
简单罗列下几个关键量的含义：
albedo：头发基础色。n：被UE4定义为fake normal，可用公式 n = normalize(v - N * (N dot v)) 计算得到，该公式中的矢量N指代发丝宏观法线。shadow：来自采样获得的光照衰减系数，是个0~1之间的数值。Luma：从albedo中提取的lum值，可参考公式：dot(albedo, (.3,.59,.11))计算获得。
参考代码实现：
//c(olor) or albedo, m(etallic), l(ight), v(iew), n(ormal), s(hadow)float3 KajiyaDiffuse(float3 C, float m, float3 l, float3 v, float3 n, float s){ float KajiyaDiffuse = 1 - abs(dot(n, l)); float3 fakeNormal = normalize(v - n * dot(v, n)); n = fakeNormal; float Wrap = 1; float NoL = saturate((dot(n ,l) + Wrap) / pow2(1 + Wrap)); //漫反射项的另一种近似，考虑了表面金属度 float diffuseScatter = (1 / UNITY_PI) * lerp(NoL, KajiyaDiffuse, 0.33) * m; float3 Luminance = (0.3f, 0.59f, 0.11f); float luma = dot(C, Luminance); float3 ScatterTint = pow(abs(C / luma), 1 - s); return sqrt(abs(C)) * diffuseScatter * ScatterTint;}轴向散射函数M
文献里一般会叫 Longitudinal scattering function，在Marschner的原文里给出的建议是使用高斯函数G作为描述散射在圆锥面上光线的分布。UE4在2016的文章里也提到了他们的方案，对于第一次反射分量Mr，他们采用了效果更好但是更加昂贵的Weta模型，同时参考如下：
14
简单说明下参数：
v = λ^2I (cosθicosθr/v) 是第一类修正贝塞尔函数（Modified Bessel function of first kind）λ 一般认为是与roughness有关的参数，也可以认为是角质层倾斜角宽度（longitudinal width）
UE4之所以不惜代价在计算Mr时采用Weta的方案，主要在于这个模型是能量守恒的，对于基于物理的环境光，能给出相比高斯G更加真实的辐射亮度估计。
对于Mtt和Mtrt这两项，UE4与Unity的HDRP都采用了基于高斯模型，具体参考如下公式：
15
同样说明一下参数：
λ ： 与roughness相关的参数α ：轴向的偏移角，一般可以暴露给艺术家调节
基本上就是利用了高斯函数的一些特性，例如：
连续性，光线分布不间断。期望(既偏移角α)处达到极致，用以模拟高光在角质层瓦片结构的驱使下的角度偏移。用关联粗糙度的参数控制方差，使得相对粗糙的发质能够得到更加均匀的散射，而光滑的发质高光集中的特性。
参考代码实现：
Mp = Hair_G(roughness, ThetaH - Alpha);//Gaussian Distribution for M term//B->λ， theta->αinline float Hair_G(float B, float Theta){ return exp(-0.5 * square(Theta) / (B*B)) / (SQRT2PI * B);}垂面（方位角）散射函数N
文献里这部分叫做Azimuthal scattering function，这个函数与出射，入射光线在垂截面上的投影矢量有关，如果涉及到透射过程，那么还会和发丝的折射率相关（不考虑吸收）。和M项类似，求取N项前会拆分成Nr，Ntt和Ntrt三项分别进行。
Nr
Nr在三个垂面散射函数中，属于相对简单同时很独特的一种，值得单独拿出来说说。以UE4为例：Nr(φ) ∝ cos(φ/2)，既入射光线经表面的直接反射(R)分量在法平面上的投影所形成的概率分布，将正比于方位角φ一半的余弦值，具体有如下公式：
16
参数：
φ = φr - φi， 既入射和出射光线在法平面上投影线段的夹角
在讨论上述公式的几何意义前首先需要明确一个之前没有清晰阐明的认知：我们所考察的每一条入射光束（如下图矢量LO）在到达发丝表面后都将会以入射点（o）为中心，向四周散射开去，形成一个圆锥，如下图所示：
图5
首先我们引出三条典型的反射光线oa，ob和oc，同时定义入射光线Lo与发丝轴线u这两个向量构成平面A。假设中间的反射光线oa也位于平面A内，说明向量Lo, oa和轴线u共面，又因为轴线u垂直于法平面，说明平面A也垂直于法平面，从而我们可以判断：入射光线Lo与出射光线oa在垂截面（法平面）上投影矢量形成的方位角将会完全相等，也就是：φr - φi = 0. 依据Nr的定义（cos(0) = 1），反射光线oa在它的众多兄弟姐妹中（ob，oc等）是亮度最高的那一个（取到极大值）。至于图中向左右侧偏移的反射光线ob和oa，它们必然不在平面A上，且随着向左右偏移得越远，反射矢量在法平面上的投影向量所形成的方位角也会慢慢远离它们的兄弟oa所形成的方位角，结果就是φ = φr - φi将会逐渐增大，直到发丝的另一侧，此时φ = π，对应cos(φ/2) = cos(π/2) = 0，从而形成最弱的反射强度（取到极小值）。这就是为什么Nr使用cos(φ/2)近似是可行的，它模拟了光线随方位角变化而衰减的物理现象，虽然不一定精确。
Ntt和Ntrt
一旦折射光进入发丝内部后，光路解算起来就会变得相对复杂许多，相信大家在读Marschner原文有关内容时就有所体会。总体来说提升的复杂度落在两处，其一是衰减项（Attenuation term），其二是方位角分布函数Dp：
17
上式积分号中的A就是衰减项，它通常可以拆分为两部分，一部分因折射定律引起，通过折射率η和入射角计算光线折射或反射的占比，可由菲涅尔项F定义；另一部分作为吸收项（absorption term），用于估算光能在发丝内部传播时被色素吸收从而衰减的程度，它是一个由光线波长μ与长度h（入射光所在直线到圆心的距离）定义的函数T，后续会展开。
至于积分式中的方位角分布函数Dp，我们可以认为它是一个均值为0的高斯分布，当D的入参为0，既 φ - Phi(p, h) = 0 时取到极大值。这里被减数φ表示待测的方位角变化量（φr - φi），减数是一个函数，它通过确定内部反射次数（由p确定）以及原始入射角度（由h确定），可以计算出当前光路的实际出射方向。于是在一次积分片段中，当给定一个h，我们能通过Dp估算出该光路对我们的目标方位角φr的贡献度，将全部的光路加总起来（通过对h在-1到1区间积分来实现）后，就能获得目标出射方向的最终光场分布了。以上就是d‘Eon在中通过Weta模型求解发丝内部光路所得的积分表达式，已经非常的简洁了，只是对于实时渲染的需求来说任然太过复杂和昂贵，我们需要更简洁近似。
菲涅尔衰减F
先说菲涅尔衰减F，这里主要参考UE4的经验，在对比了两种相近的模型：d'Eon对Cook Torrence标准模型适用于圆柱形坐标系的修改版，和Marschner利用Bravais指数进行推导的模型，最终选择了标准模型：
18
f作为经典菲涅尔计算公式，具有如下入参形式：
19
其中cos(arcsin(h))是为了计算h边对角的余弦值，即为求直角三角形另一条直角边的值，因此计算公式也可以写成这样：
20
当然在实际计算中，一般选择Schlick近似来模拟菲涅尔公式f：
21
其中η为头发的IOR，具有测定值1.55，故F0大约为0.0465。
为了能求取h，我们需要先定义和推导以下公式：
22
符号η'参考上文梳理Marschner投影折射率的段落，它是在法平面投影上的折射率。上式第一项α是给出的一个定义，后续会用。而中间项η'的定义则是UE4提出的改良版折射率计算公式，从中我们可以看出，η'是有关原始折射系数η意思入射出射角关系θd的函数。鉴于η = 1.55，η'可以近似成(c/X) + c'X 这样两项之和，X~cosθd。
float eta_prime = 1.19 / CosThetaD + 0.36 * CosThetaD;
UE4依据Weta模型用于计算htt的表达式：
23
应用了一些我们熟悉的三角函数技巧（主要是二倍角公式），将之变化为如下形式，注意此时等式两倍皆为正值：
24
然而这个形式UE4任然不满足，认为其既复杂又昂贵，于是一顿拟合和近似后给出了如下的最终形式：
25
这个公式中的所有参数我们或者已知，或者在上面已经有了近似的获取方法，所以此时我们完成了在p = 1情况下h求解的梳理，虽然这个近似htt无法在负数域正确表示（近似过程中舍弃了符号sign），好在后续对h的应用中并没有出现对符号的任何需求，相安无事。
float a = rcp(eta_prime);float Htt = CosHalfPhi * (1 + a * (0.6 - 0.8 * CosPhi));
至于p = 2时，也就是Htrt，UE4直接给了个常量 Htrt = sqrt(3) /2 来参与Fresnel衰减计算。
总结一下，各部分菲涅尔项计算代码可参考如下：
//R//补充说下，R项不涉及折射，所以无需考虑投影折射率η'和间距h的影响，只需回归纯粹的BRDF公式中F项在三维空间中的定义即可//入参sqrt(0.5f + 0.5f * LdotV)计算所得是VdotH的值，既cos(L与V夹角的一半)，符合Schlick Fresnel关于cos(angle)参数的定义 float FRValue = Fresnel(sqrt(saturate(0.5f + 0.5f * LdotV)));//TT//通过计算获得的Htt和cosθd构成投影方向cos(angle)float FTTValue = FAttenTT(Fresnel(cosThetaD * sqrt(saturate(1 - Htt * Htt))));//TRT//可见这种情况下Htrt被常数替代了float FTRTValue = FAttenTRT(Fresnel(cosThetaD * 0.5f));//Fp in Ap//p = R, FR = fresnel. p = TT, FTT = fattentt. p = TRT, FTRT = fattentrtfloat Fresnel(float cosAngle, float ior = 1.55){ float F0 = pow2((1 - ior) / (1 + ior)); return F0 + (1 - F0) * pow5(1 - cosAngle);}float FAttenTT(float Fresnel){ return pow2(1 - Fresnel);}float FAttenTRT(float Fresnel){ return pow2(1 - Fresnel) * Fresnel;}体积吸收项T
这项的最主要作用是为头发提供染色。我们回顾下电介质固有色的物理意义：当光线入射到电介质表面，少部分光能被直接反射，大量的光线会通过折射进入物质内部，其中一部分波长的光子容易被散布在物体内部的金属离子或其他组织（色素）吸收从而完成染色，最后在若干次反射后再次折射出物体表面，形成散射光被人眼或摄像机捕捉到。由此可见，吸收项T一定会带有光波长μ相关的参数以便对输出的颜色进行控制，此外T项还需要考虑在物质内传播的路径长度，这和入射角相关（θ和h），也和反射次数有关（TT或TRT）。最后非常明确的一点，R项的体积吸收项 T = 1，没有衰减！
以下是Marschner在论文中提出的吸收项T的形式：
26
结合图并利用余弦函数我们可以很轻易知道在垂截面上，单支投影光路的长度为 2cos(γt)，单次入射的最大光路（只考虑到p = 2）为 2 * 2cos(γt) = 4cos(γt)。这个函数随着γ变化有正有负会比较麻烦，所以Marschner使用 2 + 2cos(2γt)做了近似，确保γ无论处于什么角度，光路的长度总是处在0~4倍的发丝半斤之间。这部分投影光路长度的估计对应上式幂指数的分子上：1+cos(2γt)部分，注意尝试2已经被提出。分子上的cos(θt)也很好理解，因为刚才的长度只是法平面投影长度，考虑到光路在轴向上的倾角必然造成实际光路长度的放大，我们自然需要此弥补上这些轴向长度损失。
上式中的σ被Marschner定义为体积吸收因子，也就是路径传播越长，对颜色中某些分量的衰减就越大的意思。
至于EXP这种指数形式，可以参考下图对T(σ，h)的模拟结果，横轴为γt，纵轴为输出的衰减系数：
图6
当折射角γt = 0时，表示光线垂直入射，光线经过圆形到另一侧内部后原路返回，整个过程走了4 * radius的长度，认为吸收幅度最大，所有输出的衰减系数最低。当入射角γt = π/2时，代表入射光线是以掠射角度照射头发的，此时可以认为折射不发生，或者只在很小的局部发生，故没有颜色衰减，返回值如图所示为1。
OK，明白了原理之后，我们再看看UE4对此项做出的改进，首先他们做优化的起点是皮克斯（Pixar）提出的简化版T项：
27
其中ζ(C)是输入颜色的整体波长。
UE4认为折射角γt是可以进一步求解和替换的，参考如下公式使其变成包含已知量h和η'的方程：
28
然后使用头发的基础色（baseColor）C替换ζ(C)，并且带入上式后得到了一个新的T项（只针对tt分量）：
29
带入γt有些技巧，需要构建一个对边是h，临边是η'，斜边为sqrt(h^2 + η'^2)的三角形，而且要带入之前α = 1/η'的定义，此处不再赘述。
对于p = 2，既TRT分量来说，UE4直接舍去了幂指数分子的部分，改用常数项模拟：
30
最终计算所得的T项是一个三维的衰减系数，结合式所示，它与菲涅尔衰减组成了我们在式中的A项。
工程中相关shader代码如下：
float3 Ttt(float h, float3 thetaD, float cosPhi, float cosHalfPhi, float3 C){ float power = sqrt(1 - pow2(hTT(cosPhi, cosHalfPhi, thetaD)) * pow2(1 /ModifiedIor(thetaD))) / (2 * cos(thetaD)); return pow(C, power);}float3 Ttrt(float3 thetaD, float3 C){ float power = 0.8 / cos(thetaD); return pow(C, power);}方位角分布函数Dp
有关方位角分布函数的具体物理意义前文已经详细描述过，此处不再赘述，对其的实现有很多种：基于正态分布的，基于Cauchy分布的，也有基于Logistic分布的，UE4在2016年的文章中提出了自己的分布公式，据说灵感来自于观察皮克斯先前提出的基于Logistic分布的分布函数，通过产生各种高斯函数的逼近，最终确定下了如下以cosφ为变量的近似高斯分布表达式，分别是当p = 1和p = 2时的形式：
31
而下图是原文给出的Ntt输出曲线：
图7
不难发现函数大约在φ = π时取到极大值，也就是当入射方位角φi与出射方位角φr相差了一个π时，可以想象，此情况属于光线直接对穿发丝，最终能够透射出光能确实是最大的。
Ntt和Ntrt这两个分布函数代码实现如下：
//Ntt & Ntrt//UE4 approxfloat NTT(float cosPhi){ return exp(-3.65 * cosPhi - 3.98);}float NTRT(float cosPhi){ return exp(17 * cosPhi - 16.78);}总结
本文梳理了基于Marschner模型的理论依据和部分工程实现，但对毛发渲染来说仍然是管中窥豹，后续有时间我会跟进毛发渲染的一些其他领域，如毛发几何生成(卡带,线段等)，毛发顶点/骨骼动画，毛发面片之间AlphaBlend，以及一些其他比较流行的渲染理论（如基于KajiyaKay模型对Marschner的逼近等）。总之一家之谈，有错是必然的，欢迎各位读者指正错误。
marchner
Reference
An Energy-Conserving Hair Reflectance Model, d’Eon(Weta) 2011
Light Scattering from Human Hair Fibers, Marschner 2003
Rendering Fur With Three Dimensional Texture, Kajiya, Kay 1989
Physically Based Hair Shading in Unreal, Karis 2016
Marschner-Hair-Unity
Marschner-Hair-UE4
页:
[1]