剖析UE4着色模型源码

前面讲过Marschner着色模型使用了一个基于Bidirectional Scattering Distribution Function（BSDF）的方法来计算细长物体的光学性质。这个BSDF可以分解为M项和N项，其中M项表示沿着纤维方向的反射和透射，N项表示垂直于纤维方向的反射和透射。即对光进行了拆解，但是他的计算仍非常复杂应用在实时渲染不现实，因此我们有必要进行近似逼近化简从而降低开销这也是UE4的做法，其做法就是不断地近似调参直至得到和Marschner着色模型大致的效果，在不保证丢失画面质量的情况下大幅度降低计算开销从而能够使得该基于物理的着色模型应用于实时渲染中，我们会简单讲解一下对BSDF中不同光路R、TT、TRT的M项和N项逼近方法再从UE4源码出发学习一下他的实现代码。

M项逼近

具体的推导原理过于复杂这里就给出各个项的逼近计算公式，首先不论是R,TT还是TRT他们的M项都是类似的，可以使用高斯分布来进行替代。本身Mp的计算应如下计算这也是我们最后推导出来的改进版：

\begin{aligned} & M_p\left(\theta_i, \theta_r\right)=\left(\frac{1}{v e^{\frac{2}{v}}-v}\right) e^{\frac{1-\sin \theta_i \sin \theta_r}{v}} I_0\left(\frac{\cos \theta_i \cos \theta_r}{v}\right) \\ & \boldsymbol{v}=\beta_p^2 \end{aligned}

但是这样的话计算开销过大，所以简化用高斯替代：

M_p\left(\theta_i, \theta_r\right) \approx \frac{1}{\boldsymbol{\beta}_n \sqrt{2 \pi}} e^{-\frac{\left(\sin \theta_i+\sin \theta_r-\alpha_p\right)^2}{2 \beta_p{ }^2}}

其中上式有一个alpha值，R,TT和TRT分别用不同的关于shift函数值，具体格式为

float Shift = 0.035;
float Alpha[] =
{
    -Shift * 2,
    Shift,
    Shift * 4,
};

而其中的β会与粗糙度有关，因此最终Mp的计算code为：

//定义不同光的alpha值
float Shift = 0.035;
float Alpha[] =
{
    -Shift * 2,
    Shift,
    Shift * 4,
};
// 定义不同光的β值
float B[] =
    {
        Area + Pow2(ClampedRoughness),
        Area + Pow2(ClampedRoughness) / 2,
        Area + Pow2(ClampedRoughness) * 2,
    };
// 不同光的计算公式形式类似
float Hair_g(float B, float Theta)
{
    return exp(-0.5 * Pow2(Theta) / (B * B)) / (sqrt(2 * PI) * B);
}
//R
float Mp = Hair_g( B[0] , SinThetaL + SinThetaV - Alpha[0]);
// TT
float Mp = Hair_g( B[1], SinThetaL + SinThetaV - Alpha[1] );
//TRT
float Mp = Hair_g( B[2], SinThetaL + SinThetaV - Alpha[2] );

当然为了得到更真实的效果，我们可以将占贡献主体的R的M项使用原公式，只对TT,TRT进行M项逼近。

N项的逼近

M项逼近很容易，但是N项却不同，不同光的N项公式并不相同有的简单有的复杂。这里我们首先看一下思路：

R的N项逼近：当光线与细长物体垂直相交时，即光线方向与纤维方向相垂直时，可以使用R逼近方法计算N项。R逼近基于Fresnel反射定律，计算垂直方向的反射光线强度。由于这种情况下没有透射光线，因此只需要计算反射光线的强度。
TT的N项逼近：当光线与细长物体平行相交时，即光线方向与纤维方向平行时，可以使用TT逼近方法计算N项。TT逼近基于透射光线的衰减模型，计算光线从一个纤维进入另一个纤维的透射光线的强度。由于这种情况下没有反射光线，因此只需要计算透射光线的强度。
TRT的N项逼近：当光线与细长物体斜向相交时，即光线方向既不与纤维方向平行也不与纤维方向垂直时，可以使用TRT逼近方法计算N项。TRT逼近基于透射-反射-透射（transmission-reflection-transmission）的光线路径，计算光线在细长物体中传播的强度。这种情况下既有反射光线又有透射光线，因此需要同时计算反射和透射光线的强度。

R Path的N项逼近

如果直接计算R的N项，那么计算公式为：

N_r(\theta_i,\theta_r,\phi_r) = \frac{A_r(\theta_i,\phi_r)}{4\cos{\theta_i}}\left(\frac{1}{\cos{\theta_i}} + \frac{1}{\cos{\theta_r}}\right)

因为计算过于复杂所以我们可以使用R逼近方法计算N项。R逼近基于Fresnel反射定律，计算垂直方向的反射光线强度。具体来说，R的N项逼近公式为：

N_R\left(\theta_i, \theta_r, \phi\right)=\left(\frac{1}{4} \cos \frac{\phi}{2}\right) A(0, h)

对于R项的N项逼近公式 $N_R\left(\theta_i, \theta_r, \phi\right)=\left(\frac{1}{4} \cos \frac{\phi}{2}\right) A(0, h)$ ，其中， $\theta_i$ 和 $\theta_r$ 分别为入射光和反射光的夹角， $\phi$ 表示入射光和反射光的平面角度差即方位角度差， $h$ 表示视线方向和入射光/反射光的平均方向的中间向量， $A(0,h)$ 表示由中间向量 $h$ 计算出的散射光强度。

公式中的 $\left(\frac{1}{4} \cos \frac{\phi}{2}\right)$ 是一个常数项，表示在入射角和反射角一定的情况下，散射光的强度随着平面角度差的增大而逐渐减小:

\cos \frac{\phi}{2}=\sqrt{\frac{1}{2}+\frac{1}{2} \cos \phi}

$A(0, h)$ 则是中间向量 $h$ 的函数，表示散射光随着中间向量的不同而发生变化。在实际应用中， $A(0, h)$ 通常会采用一些高效的近似算法来进行计算。这里我们采用如下公式进行近似：

A(0, h)=F\left(\eta, \sqrt{\frac{1}{2}+\frac{1}{2}\left(\omega_i \cdot \omega_r\right)}\right) \begin{aligned} & F(\eta, x)=F_0+\left(1-F_0\right)(1-x)^5 \\ & F_0=\frac{(1-\eta)^2}{(1+\eta)^2} \end{aligned}

因此结合R的N项逼近和M项计算记得R的BSDF代码如下：

float F(n, x)
{
    F0=pow((1-n)/1+n),2);
    return F0 + (1 - F0) * pow((1 - x), 5);
}

vec3 wi, wr;
float Mp = Hair_g(B[0] * BScale, SinThetaL + SinThetaV - Alpha[0]);
float CosHalfPhi = sqrt(0.5 + 0.5 * cos(phi));
float A = F(n, sqrt(0.5 + 0.5 * (dotProduct(wi, wr))));
float Np = 0.25 * CosHalfPhi * A;

s+=Mp*Np;

TT Path的N项逼近

我们在前篇文章给出过N项的公式：

\begin{aligned} &A(0, h)=F(η, \frac{1}{2} \arccos (\omega_j \cdot\omega_r))\\ &A(p, h)=(1-f)^2 f^{p-1} T\left(\mu_a{ }^{\prime}, h\right)^p \end{aligned}\\ N_p(\phi)=\frac{1}{2} \int_{-1}^1 A(p, h) D_p(\phi-\Phi(p, h)) d h

为了能够实现实时计算，我们需要进行简化来减少计算开销，因此我们主要是针对h, $\eta$ ，D，A计算公式进行拟合。

首先h的初始原计算公式如下：

\begin{aligned} & a=\frac{1}{\eta^{\prime}} \\ & h_{T T}=\frac{\operatorname{sign}(\phi) \cos \frac{\phi}{2}}{\sqrt{1+a^2-2 a \operatorname{sign}(\phi) \sin \frac{\phi}{2}}} \\ & \end{aligned}

h_{T T}{ }^2=\frac{\frac{1}{2}+\frac{1}{2} \cos \phi}{1+a^2-2 a \sqrt{\frac{1}{2}-\frac{1}{2} \cos \phi}}

但是还是过于复杂所以使用如下函数近似逼近：

h_{T T} \approx(1+a(0.6-0.8 \cos \phi) \cos \frac{\phi}{2}).

同时对于 $\eta'$ 修正折射率是基于原始材质折射率 $\eta$ 计算得出，因此为了简化我们这里规定将 $\eta$ 使用固定的1.55计算，同时使用线性拟合代替元计算方程：

原始计算公式

\eta^{\prime}=\frac{\sqrt{\eta^2-\sin ^2 \theta_d}}{\cos \theta_d}

近似计算公式

\begin{aligned} & \eta=1.55 \\ & \eta^{\prime} \approx \frac{1.19}{\cos \theta_d}+0.36 \cos \theta_d \\ & 误差：\text { error }<0.68 \% \end{aligned}

而对于A的优化我们主要是需要优化A(p,h)中的吸收因子计算公式，其原始公式为：

\begin{aligned} & T(\theta, \phi)=\exp \left(-2 \mu_a \frac{1+\cos \left(2 \gamma_t\right)}{\cos \theta_t}\right) \\ & \gamma_t=\sin ^{-1}\left(\frac{h}{\eta^{\prime}}\right) \end{aligned}

其中， $\mu_a$ 表示头发的吸收系数， $\theta_t$ 表示入射光线与头发法线的夹角， $\gamma_t$ 表示头发内部的光线传播角度。

这里选择使用Pixar’s absorption尝试简化：

T(\theta, \phi)=\exp \left(-\zeta(C) \frac{\cos \gamma_t}{\cos \theta_d}\right)

简化后这里的 $\zeta(C)$ 是一个根据吸收系数 $\mu_a$ 和头发曲率半径 $C$ 计算得到的常数，一般在渲染前已经确定。但是计算公式仍较为复杂，所以在UE4中进一步优化近似如下：

T(\theta, \phi)=C^{\frac{\sqrt{1-h^2 a^2}}{2 \cos \theta_d}}

其中a和C仍是根据 $\mu_a$ 计算得到的常数，在渲染前确定好，自此T不再是一个复杂的指数计算公式同时减少了复杂的三角运算，从而整体上降低了A的计算开销。

最后我们同样优化一下D，实际上D在原论文中应该是接收三个参数的

D(\phi, s, \mu)=\frac{e^{\frac{\phi-\mu}{s}}}{s\left(1+e^{\frac{\phi-\mu}{s}}\right)^2}

$\phi$ ：单位弧度长度上的横截面积分（cross-sectional area per unit length），表示在弧度长度上纤维横截面积的总和。
$s$ ：控制密度分布的参数，被称为“形状参数”（shape parameter），它越小，密度分布越集中，越大，密度分布越均匀。
$\mu$ ：密度分布的平均值，表示纤维横截面积的平均大小。

但是计算开销太大，因此我们这里考虑到有些因素影响较小因此我们把一些参数设置为常用的常量：

\begin{aligned} & D_{T T}(\phi)=D(\phi, 0.35, \pi) \\ & D_{T T}(\phi) \approx e^{-3.65 \cos \phi-3.98} \end{aligned}

自此TT的N项逼近就大致完成了。

TRT Path的N项逼近

对于TRT我们将h设置为了常数且不包含 $\eta$ 了，因此它的吸收项T也可以大幅简化，对于D我们使用高斯分布近似：

h_{T R T}=\frac{\sqrt{3}}{2}

T_{T R T}(\theta, \phi)=C^{\frac{0.8}{\cos \theta_d}}

\begin{aligned} D_{T R T}(\phi) & =\frac{3}{4} D(\phi, 0.15,0) \\ D_{T R T}(\phi) & \approx e^{17 \cos \phi-16.78} \end{aligned}

自此我们就完成了所有光路的M项和N项逼近，实际上就是基于大量经验调整和函数逼近尽可能在降低开销的情况下得到和原公式差不多的结果，并且仍符合基于物理的特性，以下是最终的对比图我们可以看到效果不错：

虽然近似效果已经很好，但是UE4还又在最后引入了一个多散射近似，用来模拟头发散乱、多时对光线造成的多散射情况，可以看到透光性更强烈，如下所示，我们会在代码中看到其实现：

UE4 HairShading代码剖析

首先思路也是从GBuffer中获取头发的一些属性参数，然后使用封装好的BSDF结合渲染方程进行渲染，这里我们重点分析一下BSDF的封装，即如何实现对BSDF的M项和N项拆解以及近似的。

注意，虽然近似思路大致相同，但是具体的代码实现可能略有出入。

// 这是封装的高斯分布函数，考虑真实效果和开销，这里我们对TT和TRT的M项使用如下公式计算
float Hair_g(float B, float Theta)
{
    return exp(-0.5 * Pow2(Theta) / (B * B)) / (sqrt(2 * PI) * B);
}
float GaussianDetector(float Bp, float Phi)
{
    float Dp = 0;
    for (int k = -4; k <= 4; k++)
    {
        // TODO use symmetry and detect for both Phi and -Phi
        Dp += Hair_g(Bp, Phi - (2 * PI) * k);
    }
    return Dp;
}

// Shlick的菲涅尔函数
float Hair_F(float CosTheta)
{
    // 人类发丝测量折射率写死1.55
    const float n = 1.55;
    const float F0 = Pow2((1 - n) / (1 + n));
    return F0 + (1 - F0) * Pow5(1 - CosTheta);
}

// +++++++++++++ 下面两个函数主要是用来实现最终头发散射近似的 +++++++++++++++++
// Kajiya推导的解析解
float3 KajiyaKayDiffuseAttenuation(FGBufferData GBuffer, float3 L, float3 V, half3 N, float Shadow)
{
    // Use soft Kajiya Kay diffuse attenuation
    float KajiyaDiffuse = 1 - abs(dot(N, L));

    float3 FakeNormal = normalize(V - N * dot(V, N));
    //N = normalize( DiffuseN + FakeNormal * 2 );
    N = FakeNormal;
    
    // Hack approximation for multiple scattering.
    float Wrap = 1;
    float NoL = saturate((dot(N, L) + Wrap) / Square(1 + Wrap));
    float DiffuseScatter = (1 / PI) * lerp(NoL, KajiyaDiffuse, 0.33) * GBuffer.Metallic;
    float Luma = Luminance(GBuffer.BaseColor);
    float3 ScatterTint = pow(GBuffer.BaseColor / Luma, 1 - Shadow);
    return sqrt(GBuffer.BaseColor) * DiffuseScatter * ScatterTint;
}

float3 EvaluateHairMultipleScattering(
    const FHairTransmittanceData TransmittanceData,
    const float Roughness,
    const float3 Fs)
{
    return TransmittanceData.GlobalScattering * (Fs + TransmittanceData.LocalScattering) * TransmittanceData.OpaqueVisibility;
}
// +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
///////////////////////////////////////////////////////////////////////////////////////////////////
// Hair BSDF Reference code
// 首先参考原论文使用未拟合的原公式计算光照着色效果
#define HAIR_REFERENCE 0
#if HAIR_REFERENCE
    struct FHairTemp
    {
        float SinThetaL;
        float SinThetaV;
        float CosThetaD;
        float CosThetaT;
        float CosPhi;
        float CosHalfPhi;
        float VoL;
        float n_prime;
    };

    // 第一类修正贝塞尔函数
    float I0(float x)
    {
        x = abs(x);
        float a;
        if (x < 3.75)
        {
            float t = x / 3.75;
            float t2 = t * t;
            a = +0.0045813;
            a = a * t2 + 0.0360768;
            a = a * t2 + 0.2659732;
            a = a * t2 + 1.2067492;
            a = a * t2 + 3.0899424;
            a = a * t2 + 3.5156229;
            a = a * t2 + 1.0;
        }
        else
        {
            float t = 3.75 / x;
            a = +0.00392377;
            a = a * t - 0.01647633;
            a = a * t + 0.02635537;
            a = a * t - 0.02057706;
            a = a * t + 0.00916281;
            a = a * t - 0.00157565;
            a = a * t + 0.00225319;
            a = a * t + 0.01328592;
            a = a * t + 0.39894228;
            a *= exp(x) * rsqrt(x);
        }
        return a;
    }

    // d'Eon的纵向散射函数：M_p(θ_i,θ_r), v = B*B(即原论文中的v = β^2)
    float LongitudinalScattering(float B, float SinThetaL, float SinThetaV)
    {
        float v = B * B;
        float CosThetaL2 = 1 - SinThetaL * SinThetaL;
        float CosThetaV2 = 1 - SinThetaV * SinThetaV;
        float Mp = 0;
        if (v < 0.1)
        {
            float a = sqrt(CosThetaL2 * CosThetaV2) / v;
            float b = -SinThetaL * SinThetaV / v;
            float logI0a = a > 12 ? a + 0.5 * (-log(2 * PI) + log(1 / a) + 0.125 / a) : log(I0(a));
            Mp = exp(logI0a + b - rcp(v) + 0.6931 + log(0.5 / v));
        }
        else
        {
            Mp = rcp(exp(2 / v) * v - v) * exp((1 - SinThetaL * SinThetaV) / v) * I0(sqrt(CosThetaL2 * CosThetaV2) / v);
        }
        return Mp;
    }

    // 用于拟合版R，TT,TRT路径中M使用的高斯分布
    float GaussianDetector(float Bp, float Phi)
    {
        float Dp = 0;
        for (int k = -4; k <= 4; k++)
        {
            // TODO use symmetry and detect for both Phi and -Phi
            Dp += Hair_g(Bp, Phi - (2 * PI) * k);
        }
        return Dp;
    }

    // d'Eon的吸收项定义为A(η，cosφ/2)
    float3 Attenuation(uint p, float h, float3 Color, FHairTemp HairTemp)
    {
        float3 A;
        if (p == 0)
        {
            //A = F( cos( 0.5 * acos( HairTemp.VoL ) ) );
            A = Hair_F(sqrt(0.5 + 0.5 * HairTemp.VoL));
        }
        else
        {
            // ua is absorption
            // ua = pe*Sigma_ae + pp*Sigma_ap
            float3 Sigma_ae = { 0.419, 0.697, 1.37 };
            float3 Sigma_ap = { 0.187, 0.4, 1.05 };
            //float3 ua = 0.25 * Sigma_ae + 0.25 * Sigma_ap;
            float3 ua = -0.25 * log(Color);
            float3 ua_prime = ua / HairTemp.CosThetaT;
            //float3 ua_prime = ua / sqrt( 1 - Pow2( HairTemp.CosThetaD ) / 2.4 );

            float yi = asin(h);
            float yt = asin(h / HairTemp.n_prime);

            float f = Hair_F(HairTemp.CosThetaD * sqrt(1 - h * h));     // (14)
            //float3 T = exp( -2 * ua_prime * ( 1 + cos(2*yt) ) );
            float3 T = exp(-2 * ua_prime * cos(yt));
            if (p == 1)
                A = Pow2(1 - f) * T;        // (13)
            else
                A = Pow2(1 - f) * f * T * T;    // (13)
        }
        return A;
    }

    // Np = 0.5 * Integral_-1^1( A(p,h) * Dp( Phi - Omega(p,h) ) * dh )
    float Omega(uint p, float h, FHairTemp HairTemp)
    {
        float yi = asin(h);
        float yt = asin(h / HairTemp.n_prime);
        return 2 * p * yt - 2 * yi + p * PI;
    }

    // 方位散射函数
    float3 AzimuthalScattering(uint p, float Bp, float3 Color, FHairTemp HairTemp, uint2 Random)
    {
        float Phi = acos(HairTemp.CosPhi);
        // Np = 0.5 * Integral_-1^1( A(p,h) * Dp( Phi - Omega(p,h) ) * dh )
        float Offset = float(Random.x & 0xffff) / (1 << 16);
        uint Num = 16;
        float3 Np = 0;
        // 不是真的积分，而是使用求和平均近似
        for (uint i = 0; i < Num; i++)
        {
            float h = ((float)(i + Offset) / Num) * 2 - 1;
            Np += Attenuation(p, h, Color, HairTemp) * GaussianDetector(Bp, Phi - Omega(p, h, HairTemp));
        }
        Np *= 2.0 / Num;
        return 0.5 * Np;
    }

    // [d'Eon et al. 2011, "An Energy-Conserving Hair Reflectance Model"]
    // [d'Eon et al. 2014, "A Fiber Scattering Model with Non-Separable Lobes"]
    float3 HairShadingRef(FGBufferData GBuffer, float3 L, float3 V, half3 N, uint2 Random, uint HairComponents)
    {
        // to prevent NaN with decals
        // OR-18489 HERO: IGGY: RMB on E ability causes blinding hair effect
        // OR-17578 HERO: HAMMER: E causes blinding light on heroes with hair
        float ClampedRoughness = clamp(GBuffer.Roughness, 1 / 255.0f, 1.0f);
        float n = 1.55;
        FHairTemp HairTemp;
        
        // N is the vector parallel to hair pointing toward root
        HairTemp.VoL = dot(V, L);
        HairTemp.SinThetaL = dot(N, L);
        HairTemp.SinThetaV = dot(N, V);
        // SinThetaT = 1/n * SinThetaL
        HairTemp.CosThetaT = sqrt(1 - Pow2((1 / n) * HairTemp.SinThetaL));
        HairTemp.CosThetaD = cos(0.5 * abs(asin(HairTemp.SinThetaV) - asin(HairTemp.SinThetaL)));

        float3 Lp = L - HairTemp.SinThetaL * N;
        float3 Vp = V - HairTemp.SinThetaV * N;
        HairTemp.CosPhi = dot(Lp, Vp) * rsqrt(dot(Lp, Lp) * dot(Vp, Vp));
        HairTemp.CosHalfPhi = sqrt(0.5 + 0.5 * HairTemp.CosPhi);

        HairTemp.n_prime = sqrt(n * n - 1 + Pow2(HairTemp.CosThetaD)) / HairTemp.CosThetaD;

        float Shift = 0.035;
        float Alpha[] =
        {
            -Shift * 2,
            Shift,
            Shift * 4,
        };
        float B[] =
        {
            Pow2(ClampedRoughness),
            Pow2(ClampedRoughness) / 2,
            Pow2(ClampedRoughness) * 2,
        };

        float3 S = 0;
        UNROLL for (uint p = 0; p < 3; p++)
        {
            if (p == 0 && (HairComponents & HAIR_COMPONENT_R) == 0) continue;
            if (p == 1 && (HairComponents & HAIR_COMPONENT_TT) == 0) continue;
            if (p == 2 && (HairComponents & HAIR_COMPONENT_TRT) == 0) continue;

            float SinThetaV = HairTemp.SinThetaV;
            float Bp = B[p];
            if (p == 0)
            {
                Bp *= sqrt(2.0) * HairTemp.CosHalfPhi;
                float sa, ca;
                sincos(Alpha[p], sa, ca);
                SinThetaV -= 2 * sa * (HairTemp.CosHalfPhi * ca * sqrt(1 - SinThetaV * SinThetaV) + sa * SinThetaV);
            }
            else
            {
                SinThetaV = sin(asin(SinThetaV) - Alpha[p]);
            }
            float Mp = LongitudinalScattering(Bp, HairTemp.SinThetaL, SinThetaV);
            float3 Np = AzimuthalScattering(p, B[p], GBuffer.BaseColor, HairTemp, Random);

            float3 Sp = Mp * Np;
            S += Sp;
        }
        return S;
    }
#endif

///////////////////////////////////////////////////////////////////////////////////////////////////
// 毛发BSDF
// 依据以下论文中的方案进行曲线拟合
// [Marschner et al. 2003, "Light Scattering from Human Hair Fibers"]
// [Pekelis et al. 2015, "A Data-Driven Light Scattering Model for Hair"]
float3 HairShading( FGBufferData GBuffer, float3 L, float3 V, half3 N, float Shadow, FHairTransmittanceData HairTransmittance, float InBacklit, float Area, uint2 Random )
{
    // to prevent NaN with decals
    // OR-18489 HERO: IGGY: RMB on E ability causes blinding hair effect
    // OR-17578 HERO: HAMMER: E causes blinding light on heroes with hair
    float ClampedRoughness = clamp(GBuffer.Roughness, 1/255.0f, 1.0f);

    //const float3 DiffuseN = OctahedronToUnitVector( GBuffer.CustomData.xy * 2 - 1 );
    const float Backlit = min(InBacklit, HairTransmittance.bUseBacklit ? GBuffer.CustomData.z : 1);

    // 这里的策略是如果HAIR_REFERENCE宏启用就按照d'Eon的方案去做，否则就使用Karis拟合的那套方案
    #if HAIR_REFERENCE
        // todo: ClampedRoughness is missing for this code path
        float3 S = HairShadingRef( GBuffer, L, V, N, Random );
        //float3 S = HairShadingMarschner( GBuffer, L, V, N );
    #else
        const float VoL       = dot(V,L);                                                      
        const float SinThetaL = clamp(dot(N,L), -1.f, 1.f);
        const float SinThetaV = clamp(dot(N,V), -1.f, 1.f);
        float CosThetaD = cos( 0.5 * abs( asinFast( SinThetaV ) - asinFast( SinThetaL ) ) );
        //CosThetaD = abs( CosThetaD ) < 0.01 ? 0.01 : CosThetaD;
        const float3 Lp = L - SinThetaL * N;
        const float3 Vp = V - SinThetaV * N;
        const float CosPhi = dot(Lp,Vp) * rsqrt( dot(Lp,Lp) * dot(Vp,Vp) + 1e-4 );
        const float CosHalfPhi = sqrt( saturate( 0.5 + 0.5 * CosPhi ) );
        //const float Phi = acosFast( CosPhi );
        
        /**
         * η'的拟合
         * 原型：η' = sqrt( η * η - 1 + CosThetaD^2) / CosThetaD;
         * float n_prime = sqrt( n*n - 1 + Pow2( CosThetaD ) ) / CosThetaD;
         * 拟合思路：η即人类发丝折射率写死为1.55, 拟合后的η'如下：
         * η' = 1.19 / CosThetaD + 0.36 * CosThetaD;
        */
        float n = 1.55;
        float n_prime = 1.19 / CosThetaD + 0.36 * CosThetaD;

        float Shift = 0.035;
        float Alpha[] =
        {
            -Shift * 2,
            Shift,
            Shift * 4,
        };  
        float B[] =
        {
            Area + Pow2(ClampedRoughness),
            Area + Pow2(ClampedRoughness) / 2,
            Area + Pow2(ClampedRoughness) * 2,
        };

        /**
         * R路径
         * M：使用高斯分布近似，d'Eon的分布项实时计算运算量太大了：M_p(θ_i, θ_r) ≈ 1 / β * sqrt(2Π) * e^{- (sinθ_i + sinθ_r - α_p)^2 / (2 * β^2)}
         * N：N_R(θ_i, θ_r, φ) = (0.25 * cosφ/2) * A(0，h)
         *                     = (0.25 * cosφ/2) * F(η，sqrt(1/2 + 1/2(w_i . w_r)))
        */
        float3 S = 0;
        if (HairTransmittance.ScatteringComponent & HAIR_COMPONENT_R)
        {
            const float sa = sin(Alpha[0]);
            const float ca = cos(Alpha[0]);
            float Shift = 2 * sa * (ca * CosHalfPhi * sqrt(1 - SinThetaV * SinThetaV) + sa * SinThetaV);
            float BScale = HairTransmittance.bUseSeparableR ? sqrt(2.0) * CosHalfPhi : 1;
            float Mp = Hair_g(B[0] * BScale, SinThetaL + SinThetaV - Shift);
            float Np = 0.25 * CosHalfPhi;
            float Fp = Hair_F(sqrt(saturate(0.5 + 0.5 * VoL)));
            S += Mp * Np * Fp * (GBuffer.Specular * 2) * lerp(1, Backlit, saturate(-VoL));
        }

        /**
         * TT路径
         * M：依旧使用高斯分布
         * N：TT路径的N项比较复杂，基于Weta’s d'Eon与Pixar的工作分别针对"h,η，D,T"进行拟合来简化
        */
        if (HairTransmittance.ScatteringComponent & HAIR_COMPONENT_TT)
        {
            float Mp = Hair_g( B[1], SinThetaL + SinThetaV - Alpha[1] );
            /**
             * Step1: 对h的拟合
             * h的原型计算公式如下：
             * float h = CosHalfPhi * rsqrt( 1 + a*a - 2*a * sqrt( 0.5 - 0.5 * CosPhi ) );
             * float h = CosHalfPhi * ( ( 1 - Pow2( CosHalfPhi ) ) * a + 1 );
             *
             * 最终曲线拟合完的h_tt如下：
            */
            float a = 1 / n_prime;
            float h = CosHalfPhi * ( 1 + a * ( 0.6 - 0.8 * CosPhi ) );

            /**
             * Step2：η'的拟合
             * 原型：η' = sqrt( η * η - 1 + CosThetaD^2) / CosThetaD;
             * 拟合思路：η即人类发丝折射率写死为1.55, 拟合后的η'如下：
             * η' = 1.19 / CosThetaD + 0.36 * CosThetaD;
             * 代码往上翻
            */
        
            float f = Hair_F( CosThetaD * sqrt( saturate( 1 - h*h ) ) );
            float Fp = Pow2(1 - f);

            /**
             * Step3：对于吸收项T的拟合：选择Pixar的方案但没有直接用，还是做了拟合
             *
             * T与γ_t的计算原型如下：
             * T(θ，φ) = e^{-2 * μ_a * (1 + cos(2γ_t)) / (cosθt)}，其中γt = sin^-1(h / η')
             * 代码实现：float yi = asinFast(h); float yt = asinFast(h / n_prime);
             * 
             * 参考Pixar的实现：
             * T(θ，φ) = e^{-epsilo(C) * cosγt / cosθd}
             * 代码实现：float3 Tp = pow( GBuffer.BaseColor, 0.5 * ( 1 + cos(2*yt) ) / CosThetaD );
            */
            float3 Tp = 0;
            if (HairTransmittance.bUseLegacyAbsorption)
            {
                // T(θ，φ) = C^{(sqrt(1 - h^2 * a^2)) / (2 * cosθd)}
                Tp = pow(GBuffer.BaseColor, 0.5 * sqrt(1 - Pow2(h * a)) / CosThetaD);
            }
            else
            {
                // 考虑多散射的情况下
                const float3 AbsorptionColor = HairColorToAbsorption(GBuffer.BaseColor);
                Tp = exp(-AbsorptionColor * 2 * abs(1 - Pow2(h * a) / CosThetaD));
            }

            /**
             * Step4: 对分布项D的拟合
             * 技术原型：Pixar's Logistic Distribution Function
             * D(φ，s, μ) = (e^{(φ - μ) / s}) / (s^{1 + e^{(φ - μ) / s}}^2)
             * 
             * 考虑s_tt实际上贡献很小，因此近似如下：
             * D_TT(φ) = D(φ，0.35，Π) ≈ e^{-3.65cosφ - 3.98}
            */
            float Np = exp( -3.65 * CosPhi - 3.98 );

            /**
             * float t = asin( 1 / n_prime );
             * float d = ( sqrt(2) - t ) / ( 1 - t );
             * float s = -0.5 * PI * (1 - 1 / n_prime) * log( 2*d - 1 - 2 * sqrt( d * (d - 1) ) );
             * float s = 0.35;
             * float Np = exp( (Phi - PI) / s ) / ( s * Pow2( 1 + exp( (Phi - PI) / s ) ) );
             * float Np = 0.71 * exp( -1.65 * Pow2(Phi - PI) );
            */
            S += Mp * Np * Fp * Tp * Backlit;
        }

        /**
         * TRT路径
         * M：高斯分布
         * N: 尽可能拟合为常数，D_trt = e^{17cosφ - 16.78}
        */
        if (HairTransmittance.ScatteringComponent & HAIR_COMPONENT_TRT)
        {
            float Mp = Hair_g( B[2], SinThetaL + SinThetaV - Alpha[2] );
            
            /**
             * Step1 ：对h的拟合
             * h_trt = sqrt(3) / 2
             * float h = 0.75;
            */
            float f = Hair_F( CosThetaD * 0.5 );
            float Fp = Pow2(1 - f) * f;
            
            /**
             * Step2：对于吸收项T的拟合
             * T_TRT(θ，φ) = C^{0.8 / cosθd}
            */
            float3 Tp = pow( GBuffer.BaseColor, 0.8 / CosThetaD );

            /**
             * Step3：对于分布项D的拟合;
             * 技术原型：float Np = 0.75 * exp( Phi / s ) / ( s * Pow2( 1 + exp( Phi / s ) ) );
             * 
             * 拟合 s_TRT写死为0.15，拟合后的计算如下：
             * D_TRT(φ) = 0.75 * D(φ，0.15，0) ≈ e^{17cosφ - 16.78}
            */
            float Np = exp( 17 * CosPhi - 16.78 );

            S += Mp * Np * Fp * Tp;
        }
    #endif

    // 叠加上散射（该项设计的不是很好，叠不叠加基本都看不出来）
    if (HairTransmittance.ScatteringComponent & HAIR_COMPONENT_MULTISCATTER)
    {
        S  = EvaluateHairMultipleScattering(HairTransmittance, ClampedRoughness, S);
        S += KajiyaKayDiffuseAttenuation(GBuffer, L, V, N, Shadow);
    }

    S = -min(-S, 0.0);
    return S;
}