全局光照算法：IBL

1. 介绍

基于图像的光照（Image based lighting, IBL）是一类光照技术的集合。其光源不是可分解的直接光源，而是将周围环境整体视为一个大光源。现代渲染引擎中使用的IBL有四种常见类型:

• 远程光探头，用于捕捉"无限远"处的光照信息，可以忽略视差。远程探头通常包括天空， 远处的景观特征或建筑物等。它们可以由渲染引擎捕捉, 也可以高动态范围图像的形式从相机获得.
• 局部光探头，用于从特定角度捕捉世界的某个区域。捕捉会投影到立方体或球体上, 具体取决于周围的几何体。局部探头比远程探头更精确，在为材质添加局部反射时特别有用.
• 平面反射，用于通过渲染镜像场景来捕捉反射。此技术只适用于平面，如建筑地板，道路和水。
• 屏幕空间反射，基于在深度缓冲区使用光线行进方法渲染的场景，来捕捉反射。SSR效果很好，但可能非常昂贵。

IBL 通常使用（取自现实世界或从3D场景生成的）环境立方体贴图 (Cubemap) ，我们可以将立方体贴图的每个像素视为光源，在渲染方程中直接使用它。这种方式可以有效地捕捉环境的全局光照，使物体更好地融入环境。

全景图也有其它形式的记录方法，例如Spherical Map（我们常见的 HDRI Map多用 Spherical Map 表示），其对于地平线的分辨率要高于正上方的天空，比较适合室外这种天空啥都没有的环境。

理论上，物体的每一个细分的表面都应该对应着自己独有的一个半球体光照环境，而不是整个物体共享一个 Cubemap，但是这样的话对于每个细分表面都算Cubemap，性能开销还不如光线追踪——所以使用此技术的一个前提是：周围的物体足够远。

1.1 积分方程的分解

先快速回顾一下基于Cook-Torrance BRDF的反射积分方程：
$$L_o(p,w_o)={\int }_{\Omega }{k}_d\frac{c}{\pi }+k_s\frac{DFG}{4(w_0\cdot n)(w_i\cdot n)}{L}_i(p,w_i)n\cdot w_{i}d{w}_i$$

以上公式应该都是很熟了，每个符号就无需解释了

来自周围环境的入射光都可能具有一些辐射度，使得解决积分变得不那么简单。这为解决积分提出了两个要求：

• 给定任何方向向量 $w_i$ ，我们需要一些方法来获取这个方向上场景的辐射度。
• 求解积分需要快速且实时

但实际上，第一点和第二点似乎是冲突的——因为如果真的在渲染时，考虑每个方向的辐射度，那就很难做到快速实时。为了以更有效的方式解决积分，我们需要对其大部分结果进行预计算。仔细研究反射方程，我们发现BRDF 的漫反射 $k_d$ 和镜面 $k_s$ 项是相互独立的——可以将积分分成两部分：
$$L_o(p,w_o)={\int }_{\Omega }{k}_d\frac{c}{\pi }{L}_i(p,w_i)n\cdot w_{i}d{w}_i+{\int }_{\Omega }{k}_s\frac{DFG}{4(w_0\cdot n)(w_i\cdot n)}{L}_i(p,w_i)n\cdot w_{i}d{w}_i$$

这样，就可以分开研究漫反射和镜面反射。

2. Diffuse IBL

仔细观察漫反射积分，我们发现漫反射兰伯特项是一个常数项，不依赖于任何积分变量。基于此，我们可以将常数项移出漫反射积分：
$$L_o(p,w_o)=k_d\frac{c}{\pi }{\int }_{\Omega }{L}_i(p,w_i)n\cdot w_{i}d{w}_i$$

因为有上面的“周围物体足够远”的假设，这里p也和积分变量没什么关系了。但如果是室内场景，可以通过在场景中放置多个反射探针，来解决此问题——每个反射探针单独预计算其周围环境的辐照度图。这样，位置 p 处的辐照度（以及辐射度）是：最近的几个反射探针之间的辐照度插值。

这时候积分只和入射方向有关 了，而Cubemap本身就可以记录整个球面的入射角度，所以 Cubemap在这里很好的派上了用场。为了避免运行时做“积分”，我们可以在运行前预处理这张Cubemap。

这个预计算的立方体贴图，在每个采样方向 $w_o$ 上存储其积分结果，可以理解为：场景中所有能够击中面向 $w_o$ 的表面的间接漫反射光的预计算总和。这样的立方体贴图被称为辐照度图 Irradiance Map。

Irradiance Map：使用表面法线n作为索引，存储的是辐照度值 irradiance，要求的分辨率极低。可以使用预过滤的高光环境贴图的最低等级Mip来存储。
Irradiance Map直接和光照挂钩，所以这张图存的时候一定要存线性空间，并且支持分量大于1的像素值。.hdr 或者 .exr 格式的图像（radiance HDR）就可以做到这点。

2.1 立方体贴图的卷积

如果只是使用上诉的cubeMap，那我们为了得到尽量正确的积分结果，还是要对半球的各个方向进行采样。为了避免这一点，我们需要对cubemap进行预过滤，让我们可以在实时运行时，通过一次采样，得到所需的辐照度。

既然半球的朝向决定了我们捕捉辐照度的位置，我们可以预先计算每个可能的半球朝向的辐照度，这些半球朝向涵盖了所有可能的出射方向 $w_o$。

真的对每一个可能方向去采样，理论上根本做不完，也没什么必要，所以这里肯定是用一种离散积分的方法，比如最简单的黎曼积分。为了避免对难处理的立体角求积分，我们使用球坐标 $\theta$ 和 $\varphi$ 来代替立体角 $w$——这时候积分方程就变成了：

$$L_o(p,{\varphi }_o,{\theta }_o)=k_d\frac{c}{\pi }{\int }_{\varphi =0}^{2\pi }{\int }_{\theta =0}^{\frac{1}{2}\pi }{L}_i(p,{\varphi }_i,{\theta }_i)cos(\theta )sin(\theta )d\varphi d\theta$$

当我们离散地对两个球坐标轴进行采样时，每个采样近似代表了半球上的一小块区域，如上图所示。注意，由于球的一般性质，当采样区域朝向中心顶部会聚时，天顶角 θ 变高，半球的离散采样区域变小。为了平衡较小的区域贡献度，我们使用 sinθ 来权衡区域贡献度，这就是多出来的 sin 的作用。

将上诉公式转换为黎曼和形式的离散版本：

一般情况下，连续形式可按照如下方法，转换成离散形式：

所以，参考离散版本，实际代码如下（这里存取的是辐照度E，注意！）：

vec3 irradiance = vec3(0.0);  
vec3 up    = vec3(0.0, 1.0, 0.0); 
vec3 right = normalize(cross(up, normal)); 
up         = normalize(cross(normal, right));
float sampleDelta = 0.025; 
float nrSamples = 0.0; 
for (float phi = 0.0; phi < 2.0 * PI; phi += sampleDelta) 
{    
	for(float theta = 0.0; theta < 0.5 * PI; theta += sampleDelta)    
	{        
		// spherical to cartesian (in tangent space)        
		vec3 tangentSample = vec3(sin(theta) * cos(phi),  sin(theta) * sin(phi), cos(theta));        
		// tangent space to world space (for cubemap to use)        
		vec3 sampleVec = tangentSample.x * right + tangentSample.y * up + tangentSample.z * N; 
	    irradiance += texture(environmentMap, sampleVec).rgb * cos(theta) * sin(theta);        
	    nrSamples++;    
    } 
} 
irradiance = PI * irradiance * (1.0 / float(nrSamples));

这段程序可以在运行图形应用程序之前就执行，然后存在硬盘上，程序执行直接读；也可以动态生成，然后存在显存上（和Mipmap 的生成差不多）。

2.2 应用

实际使用很简单：

vec3 kS = fresnelSchlickRoughness(max(dot(N, V), 0.0), F0, roughness); 
vec3 kD = 1.0 - kS;
vec3 irradiance = texture(irradianceMap, N).rgb;
vec3 diffuse    = irradiance * albedo;
// 间接光照的漫反射部分
vec3 ambient    = (kD * diffuse) * ao; 

为什么上诉代码中，是乘以albedo，而不是$albedo/\pi$，原因见上诉：离散版本的推导。

2.3 关于菲涅尔项

由于环境光来自半球内围绕法线 N 的所有方向，因此没有一个确定的半向量，来计算菲涅耳效应。

为了模拟菲涅耳效应，我们用法线和视线之间的夹角计算菲涅耳系数。然而，之前我们是以受粗糙度影响的微表面半向量作为菲涅耳公式的输入，但我们目前没有考虑任何粗糙度，表面的反射率总是会相对较高。间接光和直射光遵循相同的属性，因此我们期望较粗糙的表面在边缘反射较弱。由于我们没有考虑表面的粗糙度，间接菲涅耳反射在粗糙非金属表面上看起来有点过强。

可以通过在 Sébastien Lagarde 提出的 Fresnel-Schlick 方程中加入粗糙度项来缓解这个问题：

vec3 fresnelSchlickRoughness(float cosTheta, vec3 F0, float roughness)
{
    return F0 + (max(vec3(1.0 - roughness), F0) - F0) * pow(1.0 - cosTheta, 5.0);
}   

此时的结果如下：

自己也写过，照着learnOpenGL这个教程，但是找不到了，就直接照搬吧。

2.4 后续研究

考虑到运行时性能和效果，还可以做几点增强：

• 用 Irradiance Map算出 SH，提升性能。全局可以使用一个SH 模拟天光照明（例如UE4中的 SkyLight），也可以组建成网格，实现动态物体的 GI（Unity中的 Light Probe Group，UE 中的 ILC）。An Efficient Representation for Irradiance Environment Maps
• 使用 Hammersley 随机采样（球面上的均匀采样）实现 Importance Sampling，提升离线资产处理的性能。Hammersley Points on the Hemisphere

2.5 Irradiance Environment Maps

Irradiance Map 和 Light Map

当初第一次看RTR4时，看到十一章的第五节——Diffuse Global Illumination时，就特别疑惑，搞不清楚light map和后文诸如H基技术的区别。特别是关于法线的描述。

这里给出一个不知道对不对的理解：

• light map : 对于每一个物体的表面，以较低的分辨率预计算它的每个像素的光照情况，然后存储在纹理中。这个时候如果把墙旋转一个角度，那么如下所示的两张图都失效了。所以，light map更多用于
  

• Irradiance Map：则是只预计算周围环境的光照情况。

SH基础

球谐函数在球域$S$上定义了一个正交基。使用如下参数化：⬇️
$$s=(x,y,z)=(\sin \theta \cos (\varphi ),\sin \theta \cos \varphi ,\cos \varphi )$$
基函数定义为：⬇️
$$Y_l^m(\theta ,\varphi )=K_l^m{e}^{im\varphi }{P}_l^{|m|}(\cos \theta ),l\in N，-l\le m\le l$$

其中，$P_l^m$是相关的勒让德多项式，$K_l^m$是归一化常数：⬇️
$$K_l^m=\sqrt{\frac{(2l+1)}{4\pi }\frac{(l-|m|)!}{(l+|m|)!}}$$
上述定义形成了一个复基complex basis，实值基是由简单的变换给出的：⬇️

$l$的低值（称为频带指数band index）表示球上的低频基函数，频带$l$的基函数在$x,y,z$上化为$l$阶多项式。可以用简单的递推公式进行计算。因为SH基是正交的，所以定义在S上的标量函数$f$可以通过如下积分⬇️，投影Projection得到系数：⬇️
$$f_l^m=\int f(s)y_l^m(s)\mathrm{d}s$$
这些系数提供了n阶重建Reconstruction函数：⬇️
$$\overline{f}(s)=\sum _{l=0}^{n-1}\sum _{m=-l}^l{f}_l^m{y}_l^m(s)$$
Properties： SH投影的一个重要性质是它的旋转不变性，例如：给定$g(s)=f(Q(s))$，$Q$是S上的一个任意的旋转函数，则：⬇️（这类似于一维傅里叶变换的位移不变性）
$$\overline{g}(s)=\overline{f}(Q(s))$$
这个不变性意味着，当在一组旋转的样本处采样f时，SH投影不会产生任何锯齿。

SH基的正交性提供了一个有用的性质，即给定任意两个球上的函数$a$和$b$，它们的投影满足：⭐️

换句话说，对带限函数的乘积进行积分，可将其简化为投影系数的点积和。

使用SH进行投影

如果我们使用3个频带的SH系数，也就是说，我们只需要计算得到9个$L_{lm}$，而不是对每个像素进行积分。这些系数的计算方法如下：
$$L_{lm}={\int }_{\theta =0}^{\pi }{\int }_{\varphi =0}^{2\pi }L(\theta ,\varphi )Y_{lm}(\theta ,\varphi )sin\theta d\theta d\varphi$$

它的离散形式如下：
$$L_{lm}=\frac{\pi }{N1}\cdot \frac{2\pi }{N2}\sum _{pixels(\theta ,\varphi )}envmap[pixel]\cdot Y_{lm}(\theta ,\varphi )$$
实际上，更加精确的公式，应该是实际考虑每个cubemap上的像素所代表的矩形区域投影到单位球面的面积：
$$L_{lm}=\sum _{pixels(\theta ,\varphi )}envmap[pixel]\cdot Y_{lm}(\theta ,\varphi )\Delta w_i$$

最终我们得到9个SH系数（27个float）。

这里可以简单讲下为什么我们只需要这么点数据（27个float）就可以记录光照情况：最大的原因就是我们的假设——光照环境无限远，每个着色点 $P_i$ 对于光照贴图来说，都是一样的，这个时候，$L(p,\theta ,\varphi )$ 就变成了 $L(\theta ,\varphi )$，所以我们只需要在这个cube map（或者探针）的中心 $P_0$，算一次投影结果，就可以推广到其他着色点。

重建

我们重建的光照度量是辐照度E，而对于E和L的SH系数有如下关系：
$$E_{lm}=A_l{L}_{lm}$$
公式中的$A_l$的定义如下：

这个 $A_l$ 的物理意义是什么？考虑diffuse项的计算，光照积分里面是 $L\cdot cos$，我们上面计算的只是 $L$，而没有余弦项。所以 $A_l$ 应该是这个余弦项的SH系数 （根据SH基的正交性可以推断出）。目前还有一个问题是，我们是要投影 $\int (n\cdot w)dw$，而这个积分哪怕基于环境光无限远的假设，除了$(\theta ,\varphi )$，我们还要考虑 $n$，这意味着我们需要$width\ast height\ast NU{M}_{SH}$个系数（其实就是三张纹理），这根本不节省带宽。考虑球谐函数的旋转不变性（法线的各不相同，本质就是旋转），我们可以对其进行简化，最终得到 $A_l$。具体推导可以见论文复现：A Non-Photorealistic Lighting Model For Automatic Technical Illustration

所以重建公式变成了：
$$E(\theta ,\varphi )=\sum _{l,m}{A}_l{L}_{lm}{Y}_{lm}(\theta ,\varphi )$$
对于实际渲染，我们可以使用下列公式来计算E：

由于我们只考虑$l\le 2$，辐照度就是一个（归一化）表面法线坐标的二次多项式。因此，对于$n^t=(x,y,z,1)$，我们可以有：
$$E(n)=n^{t}Mn$$
M是一个对称的4x4矩阵。下面的方程对渲染特别有用，因为我们只需要一个矩阵-向量乘法和一个点乘法来计算E：

总结

我们在Diffuse IBL的假设以及思想上，结合球谐函数，将环境贴图的低频光照信息投影到SH基上，这样就可以极大节省带宽，因为我们不需要存储一张预过滤图了，而是存储几十个vector就行了（以$l=2$为例，我们只需要存储9个SH vector系数即可）

实时运行过程中，只需要以法线为索引，就可以快速重建辐照度E，下面是UE4的源码：

// filament根据预缩放的SH重建辐照度的GLSL代码
vec3 irradianceSH(vec3 n) {
    // uniform vec3 sphericalHarmonics[9]
    // 我们只使用前两个波段以获得更好的性能
    return
    //另外, 由于使用 Kml 进行了预缩放, SH系数可视为颜色, 
    //特别地sphericalHarmonics[0]直接就是平均辐照度.
          sphericalHarmonics[0]
        + sphericalHarmonics[1] * (n.y)
        + sphericalHarmonics[2] * (n.z)
        + sphericalHarmonics[3] * (n.x)
        + sphericalHarmonics[4] * (n.y * n.x)
        + sphericalHarmonics[5] * (n.y * n.z)
        + sphericalHarmonics[6] * (3.0 * n.z * n.z - 1.0)
        + sphericalHarmonics[7] * (n.z * n.x)
        + sphericalHarmonics[8] * (n.x * n.x - n.y * n.y);
}

3. Specular IBL

将重点关注反射方程的镜面部分：
$$L_o(p,w_o)={\int }_{\Omega }{k}_s\frac{DFG}{4(w_0\cdot n)(w_i\cdot n)}{L}_i(p,w_i)n\cdot w_{i}d{w}_i$$
很明显，镜面部分要复杂的多，不仅受入射光方向影响，还受视角影响。如果试图解算所有入射光方向加所有可能的视角方向的积分，二者组合数会极其庞大，实时计算太昂贵。

进行预计算？但是这里的积分依赖于$w_i$和$w_o$，我们无法用两个方向向量采样预计算的立方体图。Epic Games提出了一个解决方案，他们预计算镜面部分的卷积，为实时计算作了一些妥协，这种方案被称为分割求和近似法（split sum approximation）——将预计算分成两个单独的部分求解，再将两部分组合起来，得到预计算结果。分割求和近似法将镜面反射积分拆成两个独立的积分：
$$L_o(p_{,}{w}_o)={\int }_{\Omega }{L}_i(p,w_i)d{w}_i\ast {\int }_{\Omega }{f}_r(p,w_i,w_o)n\cdot w_{i}d{w}_i$$

3.1 第一部分：光照部分

$$L_o(p_{,}{w}_o)={\int }_{{\Omega }_{f_r}}{L}_i(p,w_i)d{w}_i$$

卷积的第一部分被称为预滤波环境贴图，它类似于辐照度图，是预先计算的环境卷积贴图，但这次考虑了粗糙度。这部分看起来和上面 Diffuse IBL非常接近，唯一不一样的是：它的积分域从整个半球，变为了BRDF的覆盖范围，也就是Specular Lobe / BRDF Lobe。于是积分域就和 Lobe “撑起来的胖瘦程度”有关了。而Lobe和BRDF项的 Roughness有直接关系——越粗糙，高光越分散（极端情况就是diffuse了）。Roughness 是变量，因此需要得到一系列不同Roughness 所对应的Cubemap。

diffuse IBL中对粗糙度的考虑是菲涅尔项，但在求解积分的时候，并没有考虑，和这里是不同的。

这里轮到Mipmapping来救场了：用不同的 mipmaps离散的表示不同的 Roughness，借助着 Trilinear Filtering 三线性纹理过滤，来插值得到真正 Roughness所对应的光照强度。在实时渲染中，可以预处理原始环境贴图，得到的Mipmap 过的环境贴图被称为Pre-filtered Environment Map（预处理环境贴图），如下图所示（来自 LearnOpenGL）：

虽然积分中没有$w_o$（视线向量）的身影，但采样的球面积分域和出射角有关，我们还没有考虑——一个lobe除了胖瘦程度，还有朝向！但正如之前所说，我们已经和Diffuse IBL一样有了法线N作为索引，来采样这个预滤波环境贴图，不能在考虑第二个向量了，因此Epic Games假设视角方向V——也就是镜面反射方向R——总是等于输出采样方向N，以作进一步近似。翻译成代码如下：

vec3 N = normalize(w_o);
vec3 R = N;
vec3 V = R;

显然，这种近似会导致：在视线几乎垂直于法线的掠射方向上，会无法获得很好的掠射镜面反射

3.2 光照部分的实现

一些基础技术，这里直接给出实现，具体原理请百度。

低差异序列：Hammersley 序列

float RadicalInverse_VdC(uint bits) 
{
    bits = (bits << 16u) | (bits >> 16u);
    bits = ((bits & 0x55555555u) << 1u) | ((bits & 0xAAAAAAAAu) >> 1u);
    bits = ((bits & 0x33333333u) << 2u) | ((bits & 0xCCCCCCCCu) >> 2u);
    bits = ((bits & 0x0F0F0F0Fu) << 4u) | ((bits & 0xF0F0F0F0u) >> 4u);
    bits = ((bits & 0x00FF00FFu) << 8u) | ((bits & 0xFF00FF00u) >> 8u);
    return float(bits) * 2.3283064365386963e-10; // / 0x100000000
}
// ----------------------------------------------------------------------------
vec2 Hammersley(uint i, uint N)
{
    return vec2(float(i)/float(N), RadicalInverse_VdC(i));
}  

GGX重要性采样

有别于均匀或纯随机地（比如蒙特卡洛）在积分半球 Ω 产生采样向量，我们的采样会根据粗糙度，偏向微表面的半向量的宏观反射方向。采样过程将与我们之前看到的过程相似：

1. 开始一个大循环，生成一个随机（低差异）序列值，用该序列值在切线空间中生成样本向量，
2. 将样本向量变换到世界空间，并对场景的辐射度采样。

const uint SAMPLE_COUNT = 4096u;
for(uint i = 0u; i < SAMPLE_COUNT; ++i)
{
    vec2 Xi = Hammersley(i, SAMPLE_COUNT);   
}

此外，要构建采样向量，我们需要一些方法定向和偏移采样向量，以使其朝向特定粗糙度的镜面波瓣方向。我们可以如理论教程中所述使用 NDF，并将GGX NDF结合到 Epic Games 所述的球形采样向量的处理中：

vec3 ImportanceSampleGGX(vec2 Xi, vec3 N, float roughness)
{
    float a = roughness*roughness;

    float phi = 2.0 * PI * Xi.x;
    float cosTheta = sqrt((1.0 - Xi.y) / (1.0 + (a*a - 1.0) * Xi.y));
    float sinTheta = sqrt(1.0 - cosTheta*cosTheta);

    // from spherical coordinates to cartesian coordinates
    vec3 H;
    H.x = cos(phi) * sinTheta;
    H.y = sin(phi) * sinTheta;
    H.z = cosTheta;

    // from tangent-space vector to world-space sample vector
    vec3 up        = abs(N.z) < 0.999 ? vec3(0.0, 0.0, 1.0) : vec3(1.0, 0.0, 0.0);
    vec3 tangent   = normalize(cross(up, N));
    vec3 bitangent = cross(N, tangent);

    vec3 sampleVec = tangent * H.x + bitangent * H.y + N * H.z;
    return normalize(sampleVec);
}

着色器

#version 330 core
out vec4 FragColor;
in vec3 localPos;

uniform samplerCube environmentMap;
uniform float roughness;

const float PI = 3.14159265359;

float RadicalInverse_VdC(uint bits);
vec2 Hammersley(uint i, uint N);
vec3 ImportanceSampleGGX(vec2 Xi, vec3 N, float roughness);

void main()
{       
    vec3 N = normalize(localPos);    
    vec3 R = N;
    vec3 V = R;

    const uint SAMPLE_COUNT = 1024u;
    float totalWeight = 0.0;   
    vec3 prefilteredColor = vec3(0.0);     
    for(uint i = 0u; i < SAMPLE_COUNT; ++i)
    {
        vec2 Xi = Hammersley(i, SAMPLE_COUNT);
        vec3 H  = ImportanceSampleGGX(Xi, N, roughness);
        vec3 L  = normalize(2.0 * dot(V, H) * H - V);

        float NdotL = max(dot(N, L), 0.0);
        if(NdotL > 0.0)
        {
            prefilteredColor += texture(environmentMap, L).rgb * NdotL;
            totalWeight      += NdotL;
        }
    }
    prefilteredColor = prefilteredColor / totalWeight;

    FragColor = vec4(prefilteredColor, 1.0);
}  

当然，也可以使用Diffuse IBL中均匀采样的方法，但这样的效率太低。

3.3 第二部分：BRDF 部分

推导

$$L_o(p_{,}{w}_o)={\int }_{\Omega }{f}_r(p,w_i,w_o)n\cdot w_{i}d{w}_i$$
这个方程要求我们在$n\cdot {\omega }_o$ 、表面粗糙度、菲涅尔系数 $F_0$ 上计算BRDF方程的卷积。这等同于在纯白的环境光或者辐射度恒定为1.0的设置下，对镜面BRDF求积分。对3个变量做卷积有点复杂，不过我们可以把$F_0$移出镜面BRDF方程：

F为菲涅耳方程。将菲涅耳分母移到 BRDF 下面可以得到如下等式：

用 Fresnel-Schlick近似公式替换右边的F 可以得到：

让我们用$\alpha =(1-w_o\cdot h{)}^5$，以便更轻松地求解$F_0$：

然后我们将菲涅耳函数F分拆到两个积分里：

接下来，我们将$\alpha$替换回其原始形式，从而得到最终分割求和的BRDF方程：

公式中的两个积分分别表示$F_0$的比例和偏差 。注意，这里的$f_r$中不计算F项。积分式子里面留下来了夹角（$n$ 和 $w_o$）和粗糙度。我们将卷积后的结果存储在2D查找纹理（Look Up Texture, LUT）中，这张纹理被称为 BRDF 积分贴图。

着色器

BRDF卷积着色器在2D 平面上执行计算，直接使用其2D纹理坐标作为卷积输入（NdotV 和 roughness）。代码与预滤波器的卷积代码大体相似，不同之处在于，它现在根据 BRDF的几何函数和 Fresnel-Schlick近似来处理采样向量：

vec2 IntegrateBRDF(float NdotV, float roughness)
{
    vec3 V;
    V.x = sqrt(1.0 - NdotV*NdotV);
    V.y = 0.0;
    V.z = NdotV;

    float A = 0.0;
    float B = 0.0;

    vec3 N = vec3(0.0, 0.0, 1.0);

    const uint SAMPLE_COUNT = 1024u;
    for(uint i = 0u; i < SAMPLE_COUNT; ++i)
    {
        vec2 Xi = Hammersley(i, SAMPLE_COUNT);
        vec3 H  = ImportanceSampleGGX(Xi, N, roughness);
        vec3 L  = normalize(2.0 * dot(V, H) * H - V);

        float NdotL = max(L.z, 0.0);
        float NdotH = max(H.z, 0.0);
        float VdotH = max(dot(V, H), 0.0);

        if(NdotL > 0.0)
        {
            float G = GeometrySmith(N, V, L, roughness);
            // 我们就是基于NDF进行重要性采样的
            // 所以这里除了F，实际上也不需要计算D项。
            // 所以fr只剩下了分母，和几何项G。
            float G_Vis = (G * VdotH) / (NdotH * NdotV);
            float Fc = pow(1.0 - VdotH, 5.0);

            A += (1.0 - Fc) * G_Vis;
            B += Fc * G_Vis;
        }
    }
    A /= float(SAMPLE_COUNT);
    B /= float(SAMPLE_COUNT);
    return vec2(A, B);
}
// ----------------------------------------------------------------------------
void main() 
{
    vec2 integratedBRDF = IntegrateBRDF(TexCoords.x, TexCoords.y);
    FragColor = integratedBRDF;
}

如你所见，BRDF卷积部分是从数学到代码的直接转换。我们将角度$\theta$和粗糙度作为输入，以重要性采样产生采样向量，在整个几何体上结合BRDF的菲涅耳项对向量进行处理，然后输出每个样本上$F_0$的系数和偏差，最后取平均值。

关于几何项

与IBL 一起使用时，BRDF的几何项略有不同，因为k变量的含义稍有不同：

由于BRDF卷积是镜面IBL积分的一部分，因此我们要在 Schlick-GGX几何函数中使用$k_{IBL}$：

float GeometrySchlickGGX(float NdotV, float roughness)
{
    float a = roughness;
    float k = (a * a) / 2.0;

    float nom   = NdotV;
    float denom = NdotV * (1.0 - k) + k;

    return nom / denom;
}
// ----------------------------------------------------------------------------
float GeometrySmith(vec3 N, vec3 V, vec3 L, float roughness)
{
    float NdotV = max(dot(N, V), 0.0);
    float NdotL = max(dot(N, L), 0.0);
    float ggx2 = GeometrySchlickGGX(NdotV, roughness);
    float ggx1 = GeometrySchlickGGX(NdotL, roughness);

    return ggx1 * ggx2;
}  

请注意，虽然 k 还是从 a 计算出来的，但这里的 a 不是 roughness 的平方——如同最初对 a 的其他解释那样——在这里我们假装平方过了。我不确定这样处理是否与 Epic Games 或迪士尼原始论文不一致，但是直接将 roughness 赋给 a 得到的 BRDF 积分贴图与 Epic Games 的版本完全一致。

结果

3.4 实时运行阶段

也是公式到代码的直接复刻：

vec3 F = FresnelSchlickRoughness(max(dot(N, V), 0.0), F0, roughness);

vec3 kS = F;
vec3 kD = 1.0 - kS;
kD *= 1.0 - metallic;     

vec3 irradiance = texture(irradianceMap, N).rgb;
vec3 diffuse    = irradiance * albedo;

const float MAX_REFLECTION_LOD = 4.0;
vec3 prefilteredColor = textureLod(prefilterMap, R,  roughness * MAX_REFLECTION_LOD).rgb;   
vec2 envBRDF  = texture(brdfLUT, vec2(max(dot(N, V), 0.0), roughness)).rg;
vec3 specular = prefilteredColor * (F * envBRDF.x + envBRDF.y);

vec3 ambient = (kD * diffuse + specular) * ao; 

请注意，specular没有乘以$k_s$，因为已经乘过了菲涅耳系数。 现在，在一系列粗糙度和金属度各异的球上运行此代码：

参考

[1] LearnOpenGLCN
[2] Real Time Rendering 4th.
[3] Filament白皮书
[4] 学姐的笔记
[5] Games202