LearnOpenGL学习笔记—PBR：IBL

本节对应官网学习内容：漫反射辐照，镜面IBL

结合英文原站，中文站，以及个人实践进行描述。

本项目地址点击 这里这里这里

文章中的代码基于OpenGL的参考，真正实现时使用了WebGL弄了个网页（为了交作业_(:з」∠)_）
所以文章代码不完全等于实现的代码，但是思路和解说是正确的，可以到git项目里看该项目实现“(终版)”的细致说明

0 引入

在前面一篇PBR的光照中，我们仅仅考虑了直接光照部分，仅仅考虑直接光照时求解渲染方程非常简单，只需需根据给定的光源直接计算并叠加，而不需要在法线轴向的半球方向做积分。

Image Based Lighting（以下简称IBL）是基于图像的光照，这种技术专门计算来自于周围环境的间接光照，其会把周围环境的光影信息存储在一张环境贴图中，

这个环境贴图可以是预先生成好的，也可以是运行时动态生成好的。IBL技术使得物体的光照效果与周围的环境更加融合，大大增强沉浸感。

1 渲染方程的求解

IBL的渲染方程与PBR一样，即特化的渲染方程——反射方程，如下所示：

$L_o(p,{\omega }_o)={\int }_{\Omega }(k_d\frac{c}{\pi }+\frac{DFG}{4({\omega }_o\cdot n)({\omega }_i\cdot n)})L_i(p,{\omega }_i)n\cdot {\omega }_{i}d{\omega }_i$

我们实现基于物理的光照效果，本质就是要计算上面的渲染方程
它是对半球方向$\Omega$的所有${\omega }_i$的积分。

在前面一篇PBR中我们求解这个积分方程非常容易，因为仅考虑直接光照且光源为非面积光源时，上述的被积函数部分是一个狄拉克函数，即除了在光源方向上函数值不为0，在半球内的其余定义域该被积函数均为0，因此没有必要求半球积分。

或者说半球积分的值就等于在光源方向上的被积函数值，所以我们采用之前采用累加的处理方式。

但是在IBL中，就没有这么简单了。我们把环境贴图中的每一个像素都当作一个光源，这使得我们必须要求解半球方向的积分，即在整个半球方向采样光照信息。

如果我们直接暴力求解渲染方程的话，在每一个片元着色器上我们都要计算渲染方程的半球积分，耗费极大的性能，这对于实时性应用来说几乎不可能。

因此，为了高效地求解渲染积分方程，IBL方法将渲染方程中的积分项预先计算出来。

计算出来后并存储到指定的纹理图片上，在后面进行光照计算的时候根据相关的信息索引预先计算好的积分值，从而快速高效地求解积分。

这其实就是以空间换时间的思想。为了方便预先计算积分值，我们把公式写成如下的形式：

$L_o(p,{\omega }_o)={\int }_{\Omega }(k_d\frac{c}{\pi }+\frac{DFG}{4({\omega }_o\cdot n)({\omega }_i\cdot n)})L_i(p,{\omega }_i)n\cdot {\omega }_{i}d{\omega }_i$

$={\int }_{\Omega }{k}_d\frac{c}{\pi }{L}_i(p,{\omega }_i)n\cdot {\omega }_{i}d{\omega }_i+{\int }_{\Omega }(\frac{DFG}{4({\omega }_o\cdot n)({\omega }_i\cdot n)})L_i(p,{\omega }_i)n\cdot {\omega }_{i}d{\omega }_i$

上面公式的两个项，分别对应光照的漫反射部分和镜面反射部分。

我们把这两部分分开来，然后分别单独计算漫反射的积分值、镜面反射的积分值。

这个其实相当于过程是，对立方体环境贴图做一个特殊的卷积操作，计算得到的结果存储到各自指定查找表中，供后续的渲染使用。

这个就是IBL算法的核心思想，接下来我们就围绕这个核心思想展开相应的实现细节。

2 hdr文件转成cubemap

在实现积分预计算之前我们还有一件事要完成，就是获取周围环境贴图的辐射率值。

因为我们把环境贴图的每一个像素都当作一个光源，在进行积分计算时，我们要根据给定的入射方向${\omega }_i$去获取该方向上辐射过来的辐射率值。

如果环境贴图是一个立方体贴图的话，这很容易实现，因为立方体贴图就是根据三维向量进行索引的。

但是在IBL中我们的环境贴图并不是一个普通的环境贴图，普通的环境贴图的像素值在低动态范围（Low Dynamic Range），即每个像素的分量不超过1.0。

IBL同样是PBR渲染方法，为了能够捕捉真实的物理光影和细节，我们的环境贴图应该是高动态范围的（High Dynamic Range）。

辐射度文件的格式（扩展名为 .hdr）存储了一张完整的立方体贴图，所有六个面数据都是浮点数，允许指定 0.0 到 1.0 范围之外的颜色值，以使光线具有正确的颜色强度。

这个文件格式使用了一个聪明的技巧来存储每个浮点值：它并非直接存储每个通道的 32 位数据，而是每个通道存储 8 位，再以 alpha 通道存放指数——虽然确实会导致精度损失，但是非常有效率，不过需要解析程序将每种颜色重新转换为它们的浮点数等效值。

下面是一个示例：

这张环境贴图是从球体投影到平面上，以使我们可以将环境信息存储到一张等距柱状投影图(Equirectangular Map) 中。

类似于广角镜头拍摄得到的画面。.hdr格式保存的不是六张独立的图片，它将6张图片从球面投影到一个二维的平面上，从而构成下面这张略带扭曲的矩形纹理图。

为了后续方便获取指定像素的纹理值，我们需要把这张矩形的hdr图转换成立方体贴图。

从hdr转换到cubemap的过程，其实就是利用球面投影到立方图上的二维平面。

可以理解为球面的纹理映射为立方体上的二维纹理（或者说球面上的点映射到立方体上二维平面上的点）。

这个映射过程就是利用了球坐标和笛卡尔坐标的关系。

对于一个球心在原点的单位球体上的点$(x,y,z)$

我们也可以用采用天顶角$\theta$（与$xz$平面的夹角）和方位角$\varphi$（在$xz$平面上与$x$轴的夹角）来唯一地表示，他们的关系如下：

• $$\begin{gathered} x=cos(\theta )cos(\varphi ) \\ y=sin(\theta ) \\ z=cos(\theta )sin(\varphi ) \end{gathered}$$
  其中$\theta$的取值范围为$[-\pi /2,+\pi /2]$，$\varphi$的取值范围为$[-\pi ,+\pi ]$

我们逆运算很容易地可以分别将其映射到二维纹理坐标uv的$[0,1]$

这样我们就将一个三维立方体贴图纹理坐标映射到了二维的纹理坐标，这个其实就是构建.hdr时的映射过程。

现在我们要再现这个过程，获取三维纹理坐标对应的二维纹理坐标，然后索引上面的hdr贴图，从而得到三维纹理坐标对应的像素值，并将其保存到cubemap当中。

所谓再现这个过程，也就是要将等距柱状投影图转换为立方体贴图，我们需要渲染一个（单位）立方体。

从内部将等距柱状图投影到立方体的每个面，并将立方体的六个面的图像构造成立方体贴图。

此立方体的顶点着色器只是按原样渲染立方体，并将其局部坐标作为 3D 采样向量传递给片段着色器。

代码思路示意

#version 330 core
layout (location = 0) in vec3 position;
layout (location = 1) in vec3 normal;
layout (location = 2) in vec2 texcoord;
layout (location = 3) in vec3 color;

out vec3 worldPos;

uniform mat4 viewMatrix;
uniform mat4 projectMatrix;

void main(){
   
	worldPos = position;
	gl_Position = projectMatrix * viewMatrix * vec4(position,1.0f);
}

而在片段着色器中，我们为立方体的每个部分着色，方法类似于将等距柱状投影图整齐地折叠到立方体的每个面一样。

为了实现这一点，我们先获取片段的采样方向，这个方向是从立方体的局部坐标进行插值得到的。

然后根据之前的投影过程，反过来计算得到$\theta$和$\varphi$，投影到$[0,1]$的范围，对hdr的贴图进行采样到立方体的面上

#version 330 core
in vec3 worldPos;
out vec4 fragColor;

uniform sampler2D hdrMap;

// (1/(pi/2), 1/(pi))
const vec2 invAtan = vec2(0.1591, 0.3183);
vec2 sampleSphericalMap(vec3 v)
{
   
	vec2 uv = vec2(atan(v.z, v.x), asin(v.y));
	// to [0,1].
	uv *= invAtan;
	uv += vec2(0.5f);
	return uv;
}

void main()
{
   
	// map 3d texcoord to 2d texcoord.
	vec2 uv = sampleSphericalMap(normalize(worldPos));
	// sample hdr map.
	vec3 sampler = texture(hdrMap, uv).rgb;
	// save to one face of cubemap.
	fragColor = vec4(sampler, 1.0f);
}

通过上面的步骤，就将一个二维的hdr贴图转换成cubemap，方便我们后续的使用。

有了这个高动态范围的环境贴图，接下来我们就展开漫反射积分和镜面反射积分的预计算。

3 预计算漫反射积分

首先来看渲染方程中的漫反射积分部分，把它单独拎出来：

${\int }_{\Omega }{k}_d\frac{c}{\pi }{L}_i(p,{\omega }_i)n\cdot {\omega }_{i}d{\omega }_i=k_d\frac{c}{\pi }{\int }_{\Omega }{L}_i(p,{\omega }_i)n\cdot {\omega }_{i}d{\omega }_i$

上面的积分公式中，积分变量为${\omega }_i$

$k_d$和$\frac{c}{\pi }$都和积分无关，故可以提出积分外

我们实际上要求积分的被积函数是 $L_i(p,{\omega }_i)n\cdot {\omega }_i$

假设$P$在原点，则对于每一个半球方向，也就是确定的$n$，积分值仅仅取决于
${\omega }_i$。

因此，我们遍历所有的$n$，然后预先计算其积分值${\int }_{\Omega }{L}_i(p,{\omega }_i)n\cdot {\omega }_{i}d{\omega }_i$，存储到到一张cubemap当中.

最后渲染时直接用$n$去采样这个cubemap得到预先计算好的积分值，这个就是预计算漫反射积分的思路。

$n$是给定法线向量，它的取值范围就所有方向，相当于一个球心在原点的单位球体上的所有点。

为了预积分遍历所有的$n$，我们送入一个立方体进行预计算的渲染，不需要球体，因为立方体的顶点归一化normalize之后也就是对应的球面上的点。

然后在片元着色器根据这个$n$进行半球方向的积分，最后存储到cubemap当中。同样的，我们需要绘制6次，每次绘制保存到cubemap的六个面中的一面。

关于公式当中的积分数值计算，首先需要注意的是公式当中的积分变量${\omega }_i$是立体角。

为了方便计算，我们需要把它转成以天顶角$\theta$和方位角$\varphi$为变量的表示形式。

• 在球面坐标中，一个方向向量我们通常采用$(\theta ,\varphi )$来唯一地表示，分别是天顶角和方位角。
  在衡量发光强度和辐射辐射度量学中，立体角有着广泛的应用。
  立体角描述了站在某一点的观察者观测到的物体大小的尺度，它被定义为球表面截取的面积微分与球半径平方之比，单位为球面度。
  显然，立体角是二维圆心角的三维扩展：
  $d\omega =\frac{dA}{r^2}$
• 立体角通常转换为$(\theta ,\varphi )$来表示，在单位球体上，$d\omega =dA$
  我们转换成用$(\theta ,\varphi )$来求$dA$
  微分面积$dA$，可以看成是一个矩形，宽和高分别为对应的弧长
  

因此，公式中的积分项转成如下的二重积分：

$k_d\frac{c}{\pi }{\int }_{\Omega }{L}_i(p,{\omega }_i)n\cdot {\omega }_{i}d{\omega }_i=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\theta d\varphi$

• $n\cdot \omega i$就是$cos\theta$

求解积分需要我们在半球 $\Omega$内采集固定数量的离散样本并对其结果求平均值。

分别给每个球坐标轴指定离散样本数量 $n_1$ 和 $n_2$ 以求其黎曼和，积分式会转换为以下离散版本

对于公式中的黎曼积分，我们需要确定$d\theta$和$d\varphi$，即数值积分步长。

这个步长越小，则计算结果越接近于理论值，但耗费的时间也越多。

此外还需要提的是，我们是在切线空间的半球方向进行采样的，在切线空间获取采样的方向向量之后，我们需要把它转换到世界空间，这里我们没有构建变换矩阵，直接计算切线空间的三个基向量。

计算近似的黎曼积分的着色器代码思路如下

============================Vertex   Shader=======================
#version 330 core
layout (location = 0) in vec3 position;
layout (location = 1) in vec3 normal;
layout (location = 2) in vec2 texcoord;
layout (location = 3) in vec3 color;

out vec3 worldPos;

uniform mat4 viewMatrix;
uniform mat4 projectMatrix;

void main(){
   
	worldPos = position;
	gl_Position = projectMatrix * viewMatrix * vec4(position,1.0f);
}    
============================Fragment Shader=======================
#version 330 core
out vec4 FragColor;
in vec3 worldPos;

uniform samplerCube environmentMap;

const float PI = 3.14159265359;

void main()
{
   		
	// 世界向量充当原点的切线曲面的法线，与WorldPos对齐。
	// 给定此法线，计算环境的所有传入辐射。
    vec3 N = normalize(worldPos);

    vec3 irradiance = vec3(0.0);   
    
    // 计算切线空间
    vec3 up    = vec3(0.0, 1.0, 0.0);
    vec3 right = normalize(cross(up, N));
    up = normalize(cross(N, 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)
        {
   
            // 球面到笛卡尔（在切线空间中）
            vec3 tangentSample = vec3(sin(theta) * cos(phi),  sin(theta) * sin(phi), cos(theta));
            // 切线空间到世界空间
            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));
    
    FragColor = vec4(irradiance, 1.0);
}

• 在片元着色器中，我们取采样步长$d\theta$和$d\varphi$均为0.025，然后是两重循环，获取采样方向向量的 ( θ , ϕ ) (θ,ϕ)