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

解析全局光照Global Illumination Explained

前言：Global Illumination全局光照技术是实时渲染的必然发展方向。我参考了一些研究成果，琢磨了一下，让更多的人可以理解这项“古老”的技术。

Front Line

虽然说，如今的GPU速度已经非常之快，并行浮点计算能力已经相当的强大，但是由于很多人依旧使用传统的Phong局部光照模型对场景进行Shading处理，所以成像效果上还远远不能够与离线渲染的图像相比。如让我们目瞪口呆的Unreal Engine3，技术规格依旧是48bit Color Phong HDR。虽然人们已经开始普遍的使用Shading Language，但是光栅化、对屏幕空间差值的操作依旧没有变化。这一点不同于Pixar RenderMan。重要的原因就是，无法精确的从物理上模拟光线在场景中的交互过程。对于PC机平台来说，由于硬件性能的限制，实时光线跟踪还是不现实的事情，但是我相信未来人们一定可以在PC机平台上置身于可以和电影画面媲美的场景中交互、游戏。

经典全局光照模型的公式如下：

 

       场景的每一个物体，既有可能是发光体（萤火虫、火把），也有可能是接受光线的物体（举火把的人）。意义等同于，物体的每个点都可能会互相影响。让我们想象，有一个模型S由许多个三角形组成，以每个顶点Vertex为中心，以这点的向量Normal为法向量，构造一个半径为1的半球体Normalized Hemisphere。这样的做法类似于Ambient Occlusion的步骤。当我们从y发出一条光线以ωi的方向入射的时候，函数h返回光线与模型的第一个交点。这里可以进行并行优化。θi是指ωi和Normal的夹角，cosθi将射入的辐射光线投影到Normalized Hemisphere上。模型上每一个点y的辐射度L等于，这个点本身的辐射度Le与反射的辐射度之和。可以使用在半球上取一些点采样，进行测试，或者使用更精确不过更加耗时的Monte Carlo采样法。我们也可以简单的把上面的积分式写成这样的简单式子：

 

       fr是BRDF函数，表现了模型表面的颜色光泽等特性。一般来说这个函数的计算依赖于辐射的入射方向ωi，辐射的反射方向ωr，以及y。为了简化BRDF的表达，我们将公式简化如下，这样表面将只存在各向同性漫反射：

       

       pd(y)是相应位置的漫反射贴图数值。我个人觉得大家如果不知道如何实现BRDF，可以参考GPU Gems1上的那篇《Spatial BRDFs》一文。

        现在我们有了4个量，（S, fr, Le,Ψ），最后一个是什么呢？这是处理微平面反射向量集合的函数，可以从理论上分析得到。此时，再试图分解这个公式已经没有多少意义，不如直接用一个近似模型简化：

 

       我们直接使用这个公式得到，穿过像素Pmn光的平均辐射度值 —— 我们知道“摄象机”上有一个个像素，把yf想象为一个针孔，ωyf := P – yf作为方向，XPmn是像素特征函数，δ是Kronecker Delta函数。（是什么东西……）

        Alexander Keller在他的《Instant Radiosity》中展现了一种算法，目的是为了快速计算光照特性实现实时渲染。想象场景由有限多的虚拟光源组成。对于像素Pmn，我们计算通过这个像素的平均辐射度：

 

       最后一项TmnL决定了最少有一条辐射光线通过像素Pmn。当场景中有M个光源的时候，我们可以简单的进行累加：

 

       对于最右边的求和式，

　　Li是辐射度，Pi是第i个光源的位置。Tmn的数值可以用普通图形卡计算得到，直接计算当场景中只有一个光源存在时每个像素的亮度数值。这个算法的最大好处就是在性能和质量上可以寻求平衡，通过控制采样点数和光源数目，也就是控制多通道渲染的次数。

Behind the Enemy

       Monte Carlo方法的大名想必大家都听说过。它需要一个随机数发生器，而且这个发生器对统计结果得质量来说非常重要。这里省略了详细的推导过程，大家可以看原文[1]。

 

       首先想象从光源发出N个粒子，从光源开始，赋予初始的辐射强度Le，射向场景，p-bar≈||Tfd||约等于场景的平均反射率。这些虚拟的光子不断地反射折射，不断地点亮，直到自己强度为0。伪代码如下：

 

       其中，入射方向可以使用下列公式计算得出：

 

       这种方法的思路被用于目前大部分的Monte Carlo光线跟踪渲染引擎中，比如赫赫有名的Mental Ray，而且可以同时使用GPU与CPU进行加速渲染，缺点是多次累计填充率极高。

Continue Fighting

       可以看我的那篇斯图加特大学光线投射引擎的分析文章。

 Conclusion

       事实上，正如[2]中所说，实时光线跟踪还需要克服下面几点难点：

       Parallesim并行：无论是何种方法，为了在有限计算能力下得到最大的性能提升，并行处理都是必须的。这一点在GPU上得到充分表现，而CPU上的并行计算依赖于有经验的设计和出色的算法。

       Efficiency效率：我们需要在画质和速度中追求一个平衡点，找到最适合的采样数目，避免额外的而对画质的提升没有帮助的计算。

       Independence From Geometry独立于几何体：事实上这才是真正我们需要关注的，毫无疑问，体渲染光线投射在实现上更加的有效，可是它只能处理Volumn Texture，体积大处理速度慢来源少；传统的模型体积小来源广泛，优化分割算法成熟，但是却更适合光栅化处理，速度更快，效果也不差。而且真实的场景中，Raytracer应该无论是Volumn Textures和Indexed Models都可以一视同仁的处理。

       我相信会有越来越多的程序依赖或者逐渐采用全局光照技术，从GPU Gems1/2中的静态/动态遮挡算法的进化，到NVIDIA OpenGL SDK 10里Render To 3D Texture的展示，还有采用Volumn Rendering技术的新一代游戏的出现，全局光照技术将应用的越拉越广泛。

相关技术文档下载

[1]Alexander Keller Instant Radiosity

[2]Wald et al. Realtime Ray Tracing and its use for Interactive Global Illumination 

以后我会详细的收集相关技术文章给大家打包下载

有几本书我很感兴趣，不知道什么地方有的下载，麻烦直到的朋友告诉我一声

Programming Mental Ray和Realtime Rendering

一、什么是全局光照

全局光照（Global Illumination,简称 GI）， 作为图形学中比较酷的概念之一，是指既考虑场景中来自光源的直接光照，又考虑经过场景中其他物体反射后的间接光照的一种渲染技术。

即可以理解为：全局光照 = 直接光照(Direct Light) + 间接光照(Indirect Light)

图1 Direct illumination （原图连接）

图2 Global illumination = Direct illumination +Indirect illumination

虽说实际应用中只有漫反射全局照明的模拟算法被称为全局照明算法，但其实理论上说反射、折射、阴影都属于全局光照的范畴，因为模拟它们的时候不仅仅要考虑光源对物体的直接作用还要考虑物体与物体之间的相互作用。也是因为，镜面反射、折射、阴影一般不需要进行复杂的光照方程求解，也不需要进行迭代的计算。因此，这些部分的算法已经十分高效，甚至可以做到实时。不同于镜面反射，光的漫反射表面反弹时的方向是近似“随机”，因此不能用简单的光线跟踪得到反射的结果，往往需要利用多种方法进行多次迭代，直到光能分布达到一个基本平衡的状态。
实际应用中只有漫反射全局照明的模拟算法被称为全局照明算法
PBR中的BRDF主要模拟的是镜面高光，漫反射通常使用Lambert着色方程模拟

二、全局光照的主要算法流派

经过几十年的发展，全局光照现今已有多种实现方向，常见的全局光照主要流派列举如下：

• Ray tracing 光线追踪 [1968]
• Path tracing 路径追踪 [1986]
• Photon mapping 光子映射
• Point Based Global Illumination 基于点的全局光照
• Radiosity 辐射度
• Metropolis light transport 梅特波利斯光照传输
• Spherical harmonic lighting 球谐光照
• Ambient occlusion 环境光遮蔽
• Voxel-based Global Illumination 基于体素的全局光照
• Light Propagation Volumes Global Illumination
• Deferred Radiance Transfer Global Illumination
• Deep G-Buffer based Global Illumination
• 等。

三、渲染方程 The Rendering Equation [1986]

在前人的研究基础上，Kajiya于1986年进一步建立了渲染方程的理论，并使用它来解释光能传输的产生的各种现象。这一方程描述了场景中光能传输达到稳定状态以后，物体表面某个点在某个方向上的辐射率（Radiance）与入射辐射亮度等的关系。

可以将渲染方程理解为全局光照算法的基础，Kajiya在1986年第一次将渲染方程引入图形学后，随后出现的很多全局光照的算法，都是以渲染方程为基础，对其进行简化的求解，以达到优化性能的目的。渲染方程根据光的物理学原理，以及能量守恒定律，完美地描述了光能在场景中的传播。很多真实感渲染技术都是对它的一个近似。渲染方程在数学上的表示如下：


图 渲染方程描述了从x点沿某一方向看的光放射的总额。

四、光线投射 Ray Casting [1968]

光线投射（Ray Casting），作为光线追踪算法中的第一步，其理念起源于1968年，由Arthur Appel在一篇名为《 Some techniques for shading machine rendering of solids》的文章中提出。其具体思路是从每一个像素射出一条射线，然后找到最接近的物体挡住射线的路径，而视平面上每个像素的颜色取决于从可见光表面产生的亮度。

五、光线追踪 Ray Tracing [1979]

1979年，Turner Whitted在光线投射的基础上，加入光与物体表面的交互，让光线在物体表面沿着反射，折射以及散射方式上继续传播，直到与光源相交。这一方法后来也被称为经典光线跟踪方法、递归式光线追踪（Recursive Ray Tracing）方法，或 Whitted-style 光线跟踪方法。

光线追踪方法主要思想是从视点向成像平面上的像素发射光线，找到与该光线相交的最近物体的交点，如果该点处的表面是散射面，则计算光源直接照射该点产生的颜色；如果该点处表面是镜面或折射面，则继续向反射或折射方向跟踪另一条光线，如此递归下去，直到光线逃逸出场景或达到设定的最大递归深度。

一般的光线追踪算法有两大缺点，一是表面属性比较单一，很难去丰富光接触到物体表面以后发生的多种光学效应，二是忽略了漫反射。我们可以通过模型修正的方式来缓解，首先，物体的表面属性是混合的，比如它有一部分是反射，一部分折射，一部分漫反射。我们这么理解各个部分：当光线接触该表面以后，有一定概率发生反射，折射和漫反射。发生这些光学效应的概率我们可以事先定义好。然后我们多次计算光线跟踪，每次按照概率决定光线的反射属性，就能把各种光学效应都考虑进去。

以下这张图示可以很好的说明光线追踪方法的思路：

图 Ray Tracing Illustration First Bounce

图 典型效果图

光线追踪的伪代码：

for each pixel of the screen
{
	Final color = 0;
      	Ray = { starting point, direction };
      	Repeat
	{
		for each object in the scene
          	{
                 	determine closest ray object/intersection;
          	}
     	if intersection exists
      	{
             	for each light inthe scene
             	{
                    	if the light is not in shadow of anotherobject
                    	{
                           	addthis light contribution to computed color;
                    	}
		}
	}
      	Final color = Final color + computed color * previous reflectionfactor;
      	reflection factor = reflection factor * surface reflectionproperty;
      	increment depth;
      } until reflection factor is 0 or maximumdepth is reached
}

六、路径追踪 Path Tracing [1986]

Kajiya也于1986年提出了路径追踪算法的理念，开创了基于蒙特卡洛的全局光照这一领域。根据渲染方程， Kajiya 提出的路径追踪方法是第一个无偏（Unbiased）的渲染方法。路径追踪的基本思想是从视点发出一条光线，光线与物体表面相交时根据表面的材质属性继续采样一个方向，发出另一条光线，如此迭代，直到光线打到光源上（或逃逸出场景），然后用蒙特卡洛的方法，计算其贡献，作为像素的颜色值。

路径追踪 = 光线追踪+ 蒙特卡洛方法

具体算法流程：(1)从视点出发，到成像平面上的每一个像素向场景发出一条虚拟的光线。 (2)当光线与物体相交时按照概率确定光学效应，即是折射，反射还是散射（漫反射）。 (3)根据不同的光学效应继续跟踪和计算，直到得到最终结果。如果光学效应为漫反射，随机选择一个反射方向进行跟踪。 (4)重复前面的过程，把每次渲染出来的图像像素叠加混合，直到渲染出的结果达到满意程度。

【转】MonteCarloPathTracing的参考代码链接 https://github.com/ybbbbt/MonteCarloPathTracing
很香，真实可用，环境也配好了

图 蒙特卡洛法多次迭代的渲染图

图 基于路径追踪渲染的效果图

六、光线投射（Ray Casting ），光线追踪（Ray Tracing），路径追踪（Path Tracing）三者的的区别：

• Ray Tracing：这其实是个框架，而不是个方法。符合这个框架的都叫raytracing。这个框架就是从视点发射ray，与物体相交就根据规则反射、折射或吸收。遇到光源或者走太远就停住。一般来说运算量不小。
• Ray Casting：其实这个和volumetric可以脱钩。它就是ray tracing的第一步，发射光线，与物体相交。这个可以做的很快，在Doom 1里用它来做遮挡。
• Path Tracing：是ray tracing + 蒙特卡洛法。在相交后会选一个随机方向继续跟踪，并根据BRDF计算颜色。运算量也不小。还有一些小分类，比如Bidirectional path tracing。

参考文献：

1、【《Real-Time Rendering 3rd》 提炼总结】(八) 第九章 · 全局光照:光线追踪、路径追踪与GI技术进化编年史
2、光线投射,光线追踪与路径追踪的概念与区别
3、蒙特卡洛光线追踪
4、蒙特卡洛光线追踪 计算机图形学冬季作业