【翻译】Coding Labs :: Physically Based Rendering - Cook-Torrance

介绍

在本篇中，我们将探索和理解Cook-Torrance BRDF在PBR中的应用。我假设你熟悉上一篇文章的内容，如果没有，可以考虑从那里开始，然后再来这个更高级的话题。

Cook-Torrance BRDF 可以插入渲染方程作为其中的 $f_r$ 。使用这个BRDF，我们以两种不同的方式模拟光的行为，区分为 “漫反射” 和 “镜面反射”。它的想法是：我们所模拟的材质将在所有方向上反射一定量的光(Lambert)，并以镜面反射的方式反射另一部分的光(如镜子)。因此，Cook-Torrance BRDF 并没有完全取代我们之前的模型，而是说，对于我们新的模型，我们能够根据正在模拟的材质来指定有多少的辐射属于漫反射而有多少属于镜面反射。BRDF函数看起来会是这样：
$f_r=k_df_{lambert}+k_sf_{cook-torrance}$

其中： $f_{lambert}$ 描述了漫反射部分， $f_{cook-torrance}$ 描述了镜面反射部分， $k_d$ 是入射光被漫反射的量， $k_s$ 是被镜面反射的量。如果材质表现出很强的漫反射行为，那么我们的 $k_d$ 就很高；而如果它表现得更像一面镜子，那么 $k_s$ 就高。
在这里插入图片描述
需要注意的是：现实中的材质反射的能量都不会超过它们所接收到的。这意味着整个半球的所有 BRDF 权重之和必须小于或等于 1，我们将其写为： $\forall \omega_i \int_\Omega f_r(n\cdot \omega_o)d\omega_o \leq1$ 。这个公式可以解读为：对于所有入射光线，所有出射的BRDF的权重乘以角度余弦的总和必须小于等于1。如果结果大于1，那么就表示表面会发射能量（就像光源一样），这对于非自发光的材质来说是我们不想要的。同样，为了确保能量守恒，我们也会强制 $k_d+k_s \leq 1$ 。

Cook-Torrance BRDF

在上一篇中我们知道了 $f_{lambert} = \frac{c}{\pi}$ ，所以接下来为了完成 $f_r$ 我们需要定义 $f_{cook-torrance}$ ：
$f_{cook-torrance} = \frac{DFG}{4(\omega_o\cdot n)(\omega_i\cdot n)}$

其中： $D$ 是分布函数， $F$ 是菲涅尔函数， $G$ 是几何函数。这三个函数是Cook-Torrance的支柱，每个函数都试图模拟了表面材料的一个特定方面的行为。Cook 和 Torrance 使用的表面近似模型属于 微表面模型(Microfacet Model) 的范畴，该模型的核心想法是粗糙的表面可以看作是一组很小的微表面。而这些微表面被看作是非常小的完美反射器，它们的行为被统计学模型所描述，在我们的例子中， $D$ 和 $G$ 都是基于此想法的。

分布函数 D

$D$ ：分布函数，它用于描述给定点上微表面的统计学上的朝向。例如：假设有20%的面是朝着“向量m”的，那么将“m”填入这个分布函数中将会得到 0.2 。在文献中有些对其他分布函数的描述如 Phong 和 Beckmann，但我们所使用的分布函数叫做 GGX 【1】，定义如下：
$\begin{aligned} D(m,n,\alpha) & = \frac{\alpha^2\chi(m\cdot n)}{\pi((m\cdot n)^2(\alpha^2+tan^2(\theta_m)))^2} \\ & = \frac{\alpha^2\chi(m\cdot n)}{\pi((m\cdot n)^2(\alpha^2+(\frac{1-(m\cdot n)^2}{(m\cdot n)^2})))^2} \\ & = \frac{\alpha^2\chi(m\cdot n)}{\pi((m\cdot n)^2\alpha^2+{1-(m\cdot n)^2})^2} \\ \end{aligned}$

其中：

$\chi$ 是“指示函数”，即 $\chi(x)=\begin{cases} 1 & \text{if }x > 0 \\ 0 & \text{if } x \leq0 \end{cases}$
$\alpha$ 是表面的粗糙度（范围[0,1]）。
$m$ 是我们想要检查的方向。

译者补充：
关于 $n$ 是什么作者并没有明确的解释，我认为这里是表面的法线方向，即宏观上面的方向。
$\theta_m$ 是 $m$ 与 $n$ 的夹角，关于此处“正切表示”转为“点乘表示”的细节如下：
$tan^2\theta_m=\frac{sin^2\theta_m}{cos^2\theta_m}=\frac{1-cos^2\theta_m}{cos^2\theta_m}=\frac{1-(m\cdot n)^2}{(m\cdot n)^2}$

如果你将 “半向量” 带入 $m$ （译者注：所谓“半向量”指的是光线方向与视线方向中间的向量，即 $H=\frac{L+V}{\lVert L+V \lVert}$ ），那么不同粗糙度的结果看起来如下：（从左到右粗糙度依次是 0.1，0.3，0.6，0.8，1.0）
在这里插入图片描述
你可以看到，当粗糙度较低时（最左边的图像），只有较小范围内的点有可能朝向半向量，但是这些范围内有很高的比例的微平面确实朝向半向量，这让其输出了白色。而当粗糙度变高时，越来越大的范围的点都有概率拥有朝向半向量的微平面，但只有少数的微平面真的朝向了半向量，这让其输出了灰色。

HLSL代码如下：

float chiGGX(float v)
{
    return v > 0 ? 1 : 0;
}

float GGX_Distribution(float3 n, float3 h, float alpha)
{
    float NoH = dot(n,h);
    float alpha2 = alpha * alpha;
    float NoH2 = NoH * NoH;
    float den = NoH2 * alpha2 + (1 - NoH2);
    return (chiGGX(NoH) * alpha2) / ( PI * den * den );
}

几何函数 G

$G$ ：几何函数，它用于描述由于微面相互遮蔽而导致的光的衰减。这又是一个统计学上的近似，它模拟了 “在给定点微平面被彼此遮挡” 或 “光线在多个微平面上反弹并在到达人眼之前失去能量” 这样的概率。几何衰减来自分布函数，所以我们的 $G$ 与 GGX 的 $D$ 相配：
$G(\omega_i,\omega_o,m,n,\alpha)=G_p(\omega_o,m,n,\alpha)G_p(\omega_i,m,n,\alpha)\\ G_p(\omega,m,n,\alpha)=\chi(\frac{\omega\cdot m}{\omega\cdot n})\frac{2}{1+\sqrt{1+\alpha^2\frac{1-(\omega\cdot n)^2}{(\omega\cdot n)^2}}}$
GGX的几何函数：（从左到右粗糙度依次是：0.0，0.2，0.5，0.8，1.0）
在这里插入图片描述
$G_p$ 的HLSL代码如下：

float chiGGX(float v)
{
    return v > 0 ? 1 : 0;
}

float GGX_PartialGeometryTerm(float3 v, float3 n, float3 h, float alpha)
{
    float VoH2 = saturate(dot(v,h));
    float chi = chiGGX( VoH2 / saturate(dot(v,n)) );
    VoH2 = VoH2 * VoH2;
    float tan2 = ( 1 - VoH2 ) / VoH2;
    return (chi * 2) / ( 1 + sqrt( 1 + alpha * alpha * tan2 ) );
}

菲涅尔函数 F

$F$ ：菲涅尔函数，它用于模拟光以不同角度与表面相互作用的方式。该函数接受输入的入射光和一些表面的属性。实际的公式相当复杂，而且对于导体和绝缘体是不同的。但 Shlick C. (1994) 提供了一个很好的近似版本可以非常快速地计算，我们将用它来作为我们的模型（对于绝缘体）：
$F=F_0+(1-F_0)(1-cos\theta)^5\\ F_0=(\frac{\eta_1-\eta_2}{\eta_1+\eta_2})^2$

其中： $\theta$ 是视线方向与半向量的夹角， $\eta_1$ 和 $\eta_2$ 是两种介质的 IOR（index of refraction，折射率）。这个公式是完整版本的一个近似，描述了光在不同折射率的绝缘体中的传输。如果想模拟导体材质，就要使用不同的公式，它的一个很好的近似是【2】：
$F=\frac{(\eta-1)^2+4\eta(1-cos\theta)^5+k^2}{(\eta+1)^2+k^2}$

其中： $\eta$ 是导体的折射率， $k$ 是导体的吸收系数。该公式假定光线来自的介质是空气或真空。

绝缘体材质的菲涅尔：（从左到右折射率依次是 1.2，1.5，1.8，2.4）
在这里插入图片描述
导体材质的菲涅尔：（吸收系数 k=2，从左到右折射率依次是 1.2，1.5，1.8，2.4）

导体材质的菲涅尔：（吸收系数 k=4，从左到右折射率依次是 1.2，1.5，1.8，2.4）

如你所见，这两个函数给出了非常不同的结果，因此我们必须根据所模拟的材质是否为绝缘体来从两个公式中选择其中的一个。不过目前由Unreal等一流引擎所引领的趋势是使用 Cook-Torrance 论文中引入的概念进行进一步的近似：“法向入射反射率”。为了让数学简单，其想法是预先计算材质在法向入射时的表现，然后基于视角对这个值进行插值。Schlick的近似就是基于此的：

float3 Fresnel_Schlick(float cosT, float3 F0)
{
  return F0 + (1-F0) * pow( 1 - cosT, 5);
}

公式中的 $F_0$ 就是“法向入射反射率”，我们用以下方式计算 $F_0$ ：

// Calculate colour at normal incidence
float3 F0 = abs ((1.0 - ior) / (1.0 + ior));
F0 = F0 * F0;
F0 = lerp(F0, materialColour.rgb, metallic);

代码的最后一行耍了一个小技巧，从而让它可以同时处理导体与绝缘体。其想法是：如果材质是绝缘体，那么结果仍旧是使用IOR计算出的 $F_0$ ；但如果材质是导体，则其被替换为导体的albedo颜色。这会有一点令人困惑，尤其是在创作资产时。因此在这个实现里记住这一点很重要：如果金属性是1则将忽略IOR的值而使用albedo颜色

Cook-Torrance

现在，我们已经有了公式中所有的碎片，可以将它们拼在一起完成我们的BRDF了：
$f_r = k_d \frac{c}{\pi} + k_s \frac{DFG}{4(\omega_o \cdot n)(\omega_i \cdot n)}$

正如上一篇所说，插入进渲染方程中的BRDF是由蒙特卡洛方法近似计算的。在之前我们是预计算了积分的结果然后将其存到了一张环境CubeMap中的，但不幸的是我们在 Cook-Torrance’s BRDF 中不能这么做——因为这个公式依赖了两个向量而非先前的一个。这意味着我们没有办法预计算然后将结果存到CubeMap中了，我们必须在每一帧都计算。
$\begin{aligned} L_o(p,\omega_o) & = \int\limits_{\Omega} (k_d \frac{c}{\pi} + k_s \frac{DFG}{4(\omega_o \cdot n)(\omega_i \cdot n)}) L_i(p,\omega_i) (\mathbf n \cdot \omega_i) d\omega_i \\ & = k_d \frac{c}{\pi} \int\limits_{\Omega} L_i(p,\omega_i) (\mathbf n \cdot \omega_i) d\omega_i + k_s \int\limits_{\Omega} \frac{DFG}{4(\omega_o \cdot n)(\omega_i \cdot n)}) L_i(p,\omega_i) (\mathbf n \cdot \omega_i) d\omega_i \end{aligned}$

这个式子从加号划分为了两部分：左侧的部分就像之前那样可以使用预计算；而右侧的部分是新的，且无法预计算。因此我们可以对左侧部分使用老方法，但我们还需要一种方法来计算右侧的积分：
$k_s \int\limits_{\Omega} \frac{DFG}{4(\omega_o \cdot n)(\omega_i \cdot n)}) L_i(p,\omega_i) (\mathbf n \cdot \omega_i) d\omega_i$

对右侧的积分求值会有很高的开销，因为它要求在每一帧都运行一次蒙特卡罗估算。但幸运的是，有两种主要的方法可以优化它，使其可以在实时帧率中运行：

One is to save inside the envmap part of the formula convolved with the actual radiance map; this will also take advantage of the lower mips of the envmap.
第二种是使用 “重要性采样” 而不是蒙特卡罗方法。

要想真正解释我们如何使用这些方法以及为什么，每项技术都可以轻易展开成另一篇文章。因此我这里只是简单地呈现出最终结果，以及过程中几个关键点。Also, we will limit ourselves to use the mip levels without any prefiltering，这有点棘手，但会减少我们在重要性采样过程中需要的样本数量。

重要性采样

本质上，重要性采样是通过引入 “引导式采样” 来改进蒙特卡罗算法。它的想法是：我们可以定义一个概率分布函数描述了我们想要在何处采样更多以及在何处采样更少。这会产生一些影响稍微改变了蒙特卡罗公式，但不幸的是我们没有时间看了，因为它太长了。我们就直接给出使用重要性采样的估算器吧：
$\hat{f}(x)=\frac{1}{N} \sum_{i=1}^{N} \frac{f(x)}{p(x)}$
对于GGX，我们使用这个函数作为概率分布函数： $p(h,n,\alpha)=D(h,n,\alpha)|h\cdot n|$ 。另外，我们要对估算器运行 N 次，但是我们要给它什么样的输入呢？之前我们是使用离散的方式扫描了整个半球，但现在我们只采样 “重要” 的区域，那么我们该怎样生成正确的向量呢？生成采样向量的方法的完整推导相当长（不过值得学习，特别是如果你想使用自己的概率分布函数），所以在这一篇中我就直接使用下面的结果而不证明了：
$\theta_s = arctan(\frac{\alpha\sqrt{\xi_1}}{\sqrt{1-\xi_1}})\\ \phi_s = 2\pi\xi_2$
其中 $\theta_s$ 和 $\phi_s$ 是采样向量在球形坐标的极角和方位角。为了得到这两个角，我们使用 [0,1] 范围内的两个均匀随机生成的变量 $\xi_1$ 和 $\xi_2$ 。

最后的步骤

终于，我们可以把所有东西放在一起组合为一个实际可用的公式了：
$k_s \int\limits_{\Omega} \frac{DFG}{4(\omega_o \cdot n)(\omega_i \cdot n)}) L_i(p,\omega_i) (\mathbf n \cdot \omega_i) d\omega_i \approx k_s \frac{1}{N} \sum_{i=1}^{N} \frac{L_i(\omega_i,p)F(\eta_i,\eta_t,h,\omega_o)G( \omega_i,\omega_o,h,n)sin(\theta_i)}{4|\omega_o \cdot n||h \cdot n|}$
我们翻译的代码如下：

float3 GGX_Specular( TextureCube SpecularEnvmap, float3 normal, float3 viewVector, float roughness, float3 F0, out float3 kS )
{
    float3 reflectionVector = reflect(-viewVector, normal);
    float3x3 worldFrame = GenerateFrame(reflectionVector);
    float3 radiance = 0;
    float  NoV = saturate(dot(normal, viewVector));

    for(int i = 0; i < SamplesCount; ++i)
    {
        // Generate a sample vector in some local space
        float3 sampleVector = GenerateGGXsampleVector(i, SamplesCount, roughness);
        // Convert the vector in world space
        sampleVector = normalize( mul( sampleVector, worldFrame ) );

        // Calculate the half vector
        float3 halfVector = normalize(sampleVector + viewVector);
        float cosT = saturatedDot( sampleVector, normal );
        float sinT = sqrt( 1 - cosT * cosT);

        // Calculate fresnel
        float3 fresnel = Fresnel_Schlick( saturate(dot( halfVector, viewVector )), F0 );
        // Geometry term
        float geometry = GGX_PartialGeometryTerm(viewVector, normal, halfVector, roughness) * GGX_PartialGeometryTerm(sampleVector, normal, halfVector, roughness);
        // Calculate the Cook-Torrance denominator
        float denominator = saturate( 4 * (NoV * saturate(dot(halfVector, normal)) + 0.05) );
        kS += fresnel;
        // Accumulate the radiance
        radiance += SpecularEnvmap.SampleLevel( trilinearSampler, sampleVector, ( roughness * mipsCount ) ).rgb * geometry * fresnel * sinT / denominator;
    }

    // Scale back for the samples count
    kS = saturate( kS / SamplesCount );
    return radiance / SamplesCount;        
}

通过蒙特卡罗估算器，我们可以得出表面反射了多少光，然后我们可以把它加到已使用Lambert计算好的漫反射分量上。我们要做的最后一件事就是计算 $k_d$ 和 $k_s$ 。因为这些权重代表了有多少比例的光被反射，所以很自然会想到使用菲涅尔。事实上，我们就是使用了菲涅尔作为 $k_s$ （译者注：此处逻辑见上面的代码），然后根据能量守恒定律可得 $k_d=1-k_s$ 。注意：由于在GGX的计算中我们已经为辐射度乘算菲涅尔了，所以在最终的公式中我们不需要再乘一遍了。

float4 PixelShaderFunction(VertexShaderOutput input) : SV_TARGET0
{
    ...

    float3 surface = tex2D(surfaceMap_Sampler, input.Texcoord).rgb;   
    float ior = 1 + surface.r;
    float roughness = saturate(surface.g - EPSILON) + EPSILON;
    float metallic = surface.b;

    // Calculate colour at normal incidence
    float3 F0 = abs ((1.0 - ior) / (1.0 + ior));
    F0 = F0 * F0;
    F0 = lerp(F0, materialColour.rgb, metallic);
        
    // Calculate the specular contribution
    float3 ks = 0;
    float3 specular = GGX_Specular(specularCubemap, normal, viewVector, roughness, F0, ks );
    float3 kd = (1 - ks) * (1 - metallic);
    // Calculate the diffuse contribution
    float3 irradiance = texCUBE(diffuseCubemap_Sampler, normal ).rgb;
    float3 diffuse = materialColour * irradiance;

    return float4( kd * diffuse + /*ks **/ specular, 1);     
}