实时渲染Real-time Rendering：Image Based Lighting

本文在前一篇PBR渲染器的基础上，进一步深入理解并实现了Image Based Lighting的全局光照。Image Based Lighting预先烘培辐照度贴图、构建BRDF积分查找表，然后在实时渲染中直接查询查找表快速计算渲染方程，实现的效果高效且真实。这种技术已经集成在虚幻4引擎中。

• 渲染方程的求解
• hdr文件转成cubemap
• 预计算漫反射积分
• 预计算镜面反射积分
• 计算渲染方程
• 实现效果
• 参考资料

Image Based Lighting

  在前面一篇基于物理的渲染中我们仅仅考虑了直接光照部分，仅仅考虑直接光照时求解渲染方程非常简单，只需需根据给定的光源直接计算并叠加，而不需要在法线轴向的半球方向做积分。Image Based Lighting（以下简称IBL）翻译过来的意思就是基于图像的光照，这种技术专门计算来自于周围环境的间接光照，周围环境的光影信息存储在一张环境贴图中，这个环境贴图可以是预先生成好的，也可以是运行时动态生成好的。IBL技术使得物体的光照效果与周围的环境更加融合，大大增强沉浸感。

一、渲染方程的求解

  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_id\omega_i \tag {1}$$

  我们实现基于物理的光照效果的本质就是要计算上面的渲染方程，它是对半球方向$\Omega$的所有$\omega_i$积分。在前面一偏PBR中我们求解这个积分方程非常容易，因为仅考虑直接光照且光源为非面积光源时，上述的被积函数部分是一个狄拉克函数，即除了在光源方向上函数值不为$0$，在半球内的其余定义域该被积函数均为$0$，因此没有必要求半球积分，或者说半球积分的值就等于在光源方向上的被积函数值。但是在IBL中，就没有这么简单了。我们把环境贴图中的每一个像素都当作一个光源，这使得我们必须要求解半球方向的积分，即在整个半球方向采样光照信息。如果我们直接暴力求解渲染方程$(1)$的话，在每一个片元着色器上我们都要计算渲染方程的半球积分，耗费极大的性能，这对于实时性应用来说几乎不可能。

  因此，为了高效地求解渲染积分方程，IBL方法将渲染方程中的积分项预先计算出来并存储到指定的纹理图片上，在后面进行光照计算的时候根据相关的信息索引预先计算好的积分值，从而快速高效地求解积分。这其实就是以空间换时间的思想。为了方便预先计算积分值，我们把公式$(1)$写成如下的形式：

$$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_id\omega_i\\ =\int_{\Omega}k_d\frac{c}{\pi}L_i(p,\omega_i)n\cdot \omega_id\omega_i +\int_{\Omega}(\frac{DFG}{4(\omega_o\cdot n)(\omega_i\cdot n)})L_i(p,\omega_i)n\cdot \omega_id\omega_i \tag {2}$$

  上面公式$(2)$的两个项分别对应光照的漫反射部分和镜面反射部分，我们把这两部分分开来，然后分别单独计算漫反射的积分值、镜面反射的积分值，这个其实相当于对立方体环境贴图做一个特殊的卷积操作，计算得到的结果分别存储到各自指定查找表中供后续的渲染使用。这个就是IBL算法的核心思想，接下来我们就围绕这个核心思想展开相应的实现细节。

二、hdr文件转成cubemap

  在实现积分预计算之前我们还有一件事要完成，就是获取周围环境贴图的辐射率值，因为我们把环境贴图的每一个像素都当作一个光源。在进行积分计算时，我们要根据给定的入射方向$\omega_i$去获取该方向上辐射过来的辐射率值。如果环境贴图是一个立方体贴图的话，这很容易实现，因为立方体贴图就是根据三维向量进行索引的。但是在IBL中我们的环境贴图并不是一个普通的环境贴图，普通的环境贴图的像素值在低动态范围（Low Dynamic Range），即每个像素的分量不超过1.0。IBL同样是PBR渲染方法，为了能够捕捉真实的物理光影和细节，我们的环境贴图应该是高动态范围的（High Dynamic Range）。

  目前有一种以.hdr为后缀的图像格式保存了高动态范围的像素值，这种格式的文件采用一种巧妙的方法将超出$0.0$到$1.0$的正确的颜色值保存下来。高动态范围的环境贴图就采用了这种格式保存，stb_image支持.hdr文件加载。这里收集了一些hdr格式的环境贴图，stb_image加载.hdr文件也非常简单，如下所示：

float *data = stbi_loadf(path.c_str(), &m_width, &m_height, &m_channel, 0);

  加载后得到的图片如下图1所示，这是一个二维的图片，类似于广角镜头拍摄得到的画面。.hdr格式保存的不是六张独立的图片，它将6张图片从球面投影到一个二维的平面上，从而构成下面这张略带扭曲的矩形纹理图。为了后续方便获取指定像素的纹理值，我们需要把这张矩形的hdr图转换成立方体贴图。

图1 .hdr贴图

  从hdr转换到cubemap的过程其实就是利用球面投影到二维平面，可以理解为球面的纹理映射为二维的纹理（或者说球面上的点映射到二维平面上的点）。这个映射过程就是利用了球坐标和笛卡尔坐标的关系，对于一个球心在原点的单位球体上的点$(x,y,z)$，我们也可以用采用天顶角$\theta$（与$xz$平面的夹角）和方位角$\phi$（在$xz$平面上与$x$轴的夹角）来唯一地表示，他们的关系如下：

$$x=cos(\theta)cos(\phi)\\ y=sin(\theta) \\ z=cos(\theta)sin(\phi) \tag {3}$$

  其中$\theta$的取值范围为$[-\pi/2,+\pi/2]$，$\phi$的取值范围为$[-\pi,+\pi]$，我们很容易地可以分别将其映射到二维纹理坐标uv的$[0,1]$，这样我们就将一个三维立方体贴图纹理坐标映射到了二维的纹理坐标，这个其实就是构建.hdr时的映射过程。现在我们要再现这个过程，获取三维纹理坐标对应的二维纹理坐标，然后索引上面的hdr贴图，从而得到三维纹理坐标对应的像素值，并将其保存到cubemap当中。我们采用的方法就是绘制6次立方体，每次保存1个立方体面的像素值：

void IBLAuxiliary::convertToCubemap(int width, int height,
	unsigned int hdrTexIndex, unsigned int cuebmapTexIndex)
{
	// manager.
	TextureMgr::ptr texMgr = TextureMgr::getSingleton();
	ShaderMgr::ptr shaderMgr = ShaderMgr::getSingleton();
	// load shader.
	unsigned int shaderIndex = shaderMgr->loadShader("convertToCubemap",
		"./glsl/convertToCubemap.vert", "./glsl/convertToCubemap.frag");
	// load cube mesh.
	Mesh::ptr cubeMesh = std::shared_ptr<Mesh>(new Cube(1.0f, 1.0f, 1.0f));
	// load framebuffer.
	FrameBuffer::ptr framebuffer = std::shared_ptr<FrameBuffer>(new FrameBuffer(width, height,
		"convertDepth", {}, true));
	// projection matrix and view matrix.
	glm::mat4 captureProjectMatrix = glm::perspective(glm::radians(90.0f), 1.0f, 0.1f, 10.0f);
	glm::mat4 captureViewMatrix[] =
	{
		glm::lookAt(glm::vec3(0.0f), glm::vec3(+1.0f, 0.0f, 0.0f), glm::vec3(0.0f, -1.0f, 0.0f)),
		glm::lookAt(glm::vec3(0.0f), glm::vec3(-1.0f, 0.0f, 0.0f), glm::vec3(0.0f, -1.0f, 0.0f)),
		glm::lookAt(glm::vec3(0.0f), glm::vec3(0.0f,+1.0f, 0.0f), glm::vec3(0.0f,  0.0f,+1.0f)),
		glm::lookAt(glm::vec3(0.0f), glm::vec3(0.0f,-1.0f, 0.0f), glm::vec3(0.0f,  0.0f,-1.0f)),
		glm::lookAt(glm::vec3(0.0f), glm::vec3(0.0f, 0.0f,+1.0f), glm::vec3(0.0f, -1.0f, 0.0f)),
		glm::lookAt(glm::vec3(0.0f), glm::vec3(0.0f, 0.0f,-1.0f), glm::vec3(0.0f, -1.0f, 0.0f)),
	};

	// convert.
	framebuffer->bind();
	glDisable(GL_BLEND);
	glDisable(GL_CULL_FACE);
	glEnable(GL_DEPTH_TEST);
	glDepthFunc(GL_LEQUAL);
	glClearColor(0.0f, 1.0f, 0.0f, 1.0f);
	glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
	GLuint cubemapId = texMgr->getTexture(cuebmapTexIndex)->getTextureId();
	Shader::ptr shader = shaderMgr->getShader(shaderIndex);
	shader->bind();
	shader->setInt("hdrMap", 0);
	shader->setMat4("projectMatrix", captureProjectMatrix);
	texMgr->bindTexture(hdrTexIndex, 0);
	for (unsigned int i = 0; i < 6; ++i)
	{
		shader->setMat4("viewMatrix", captureViewMatrix[i]);
		glFramebufferTexture2D(GL_FRAMEBUFFER, GL_COLOR_ATTACHMENT0,
			GL_TEXTURE_CUBE_MAP_POSITIVE_X + i, cubemapId, 0);
		cubeMesh->draw(false, 0);
	}
	shader->unBind();
	texMgr->unBindTexture(hdrTexIndex);
	framebuffer->unBind();
}

  相应的，在着色器中我们将三维纹理坐标投影到二维纹理坐标，索引hdr值，输出为片元着色器的值。

=======================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
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，如下图2所示，方便我们后续的使用。有了这个高动态范围的环境贴图，接下来我们就展开漫反射积分和镜面反射积分的预计算。

图2 hdr立方体贴图

三、预计算漫反射积分

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

$$\int_{\Omega}k_d\frac{c}{\pi}L_i(p,\omega_i)n\cdot \omega_id\omega_i =k_d\frac{c}{\pi}\int_{\Omega}L_i(p,\omega_i)n\cdot \omega_id\omega_i \tag {4}$$

  上面的积分公式中，积分变量为$\omega_i$，和$\frac{c}{\pi}$项都与$\omega_i$无关，故可以提出积分外。我们实际上要求积分的被积函数是$L_i(p,\omega_i)n\cdot \omega_i$。假设点$p$在原点，则对于每一个$n$确定的半球方向，积分值仅仅取决于$\omega_i$。因此，我们遍历所有的$n$，然后预先计算其积分值$\int_{\Omega}L_i(p,\omega_i)n\cdot \omega_id\omega_i$到一张cubemap当中，最后渲染时直接用$n$去采样这个cubemap得到预先计算好的积分值。这个就是预计算漫反射积分的思路。

  $n$就是给定法线向量，它的取值范围就是所有方向，相当于一个球心在原点的单位球体上的所有点。为了遍历所有的$n$，我们送入一个立方体进行预计算的渲染，不需要球体，因为立方体的顶点归一化normalize之后也就是对应的球面上的点。然后在片元着色器根据这个$n$进行半球方向的积分，最后存储到cubemap当中。同样的，我们需要绘制6次，每次绘制保存到cubemap的六个面中的一面。

  然后就是关于公式$(4)$当中的积分数值计算。首先需要注意的是公式$(4)$当中的积分变量$\omega_i$是立体角，为了方便计算，我们需要把他转成以天顶角$\theta$和方位角$\phi$为变量的表示形式。前面在光线追踪器Ray Tracer：进阶篇已经提到过微分立体角$d\omega_i$如何转成用$d\theta$和$d\phi$表示，即$d\omega_i=sin\theta d\theta d\phi$。因此，公式$(4)$中的积分项转成如下的二重积分：

$$k_d\frac{c}{\pi}\int_{\Omega}L_i(p,\omega_i)n\cdot \omega_id\omega_i = k_d\frac{c}{\pi}\int_{\phi=0}^{2\pi}\int_{\theta=0}^{\frac12\pi}L_i(p,\phi_i,\theta_i)cos\theta sin \theta d\theta d\phi \tag {5}$$

  这里要提一下，公式$(4)$中的$n\cdot \omega_i$就是$cos\theta$，半球方向内$\theta$取值$[0,\frac12\pi]$，$\phi$取值$[0,2\pi]$。公式$(5)$等号两边完全是等价的。然后我们需要采用数值方法近似计算公式$(5)$的积分值，在半球方向做均匀的离散采样$(\phi,\theta)$，计算采样的方向对应的被积函数值，最后叠加起来，叠加的结果就是黎曼积分和。公式$(5)$的黎曼积分近似公式为：

$$k_d\frac{c}{\pi}\int_{\phi=0}^{2\pi}\int_{\theta=0}^{\frac12\pi}L_i(p,\phi_i,\theta_i)cos\theta sin \theta d\theta d\phi \approx k_d\frac{c}{\pi}\Sigma_0^{2\pi}\Sigma_0^{\frac12\pi}L_i(p, \phi_i, \theta_i)cos\theta sin\theta d\theta d\phi \tag{6}$$

  对于公式$(6)$中的黎曼积分，我们需要确定$d\theta$和$d\phi$，即数值积分步长。这个步长越小，则计算结果越接近于理论值，但耗费的时间也越多。此外还需要提的是，我们是在切线空间的半球方向进行采样的，在切线空间获取采样的方向向量之后，我们需要把它转换到世界空间，这里我们没有构建变换矩阵，直接计算切线空间的三个基向量。计算近似的黎曼积分的着色器代码如下所示：

============================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
in vec3 worldPos;
out vec4 fragColor;

uniform samplerCube environmentMap;

const float PI = 3.14159265359;

void main()
{
	vec3 normal = normalize(worldPos);
	vec3 irradiance = vec3(0.0f);
	
	// tangent space.
	vec3 up = vec3(0.0f, 1.0f, 0.0f);
	vec3 right = cross(up, normal);
	up = cross(normal, right);
	
	float sampleDelta = 0.025;
	float sampleDeltaSquared = sampleDelta * sampleDelta;
	for(float phi = 0.0f; phi < 2.0 * PI;phi += sampleDelta)
	{
		for(float theta = 0.0f; theta < 0.5 * PI;theta += sampleDelta)
		{
			vec3 sampleDir = vec3(sin(theta) * cos(phi), sin(theta) * sin(phi), cos(theta));
			// tangent space to world space.
			sampleDir = sampleDir.x * right + sampleDir.y * up + sampleDir.z * normal;
			
			irradiance += texture(environmentMap, sampleDir).rgb * cos(theta) * sin(theta) * sampleDeltaSquared;
			++ nSamples;
		}
	}
	irradiance = (1.0f / PI) * irradiance ;
	fragColor = vec4(irradiance, 1.0f);
}

  在片元着色器中，我们取采样步长$d\phi$和$d\theta$均为$0.025$，然后是两重循环，获取采样方向向量的$(\phi,\theta)$后我们需要把它转换成笛卡尔坐标下的方向向量，用这个方向相应去获取该方向对应的hdr立方体纹理。最后结果除以一个$\pi$，对应公式$(6)$中的$\pi$，而剩下的漫反射因子$k_d$和反照率$c$由渲染时确定。上面的过程相当于对hdr立方体贴图做一个卷积操作，存储结果的纹理我们称之为辐照度纹理（Irradiance Map）。由于辐照度纹理没有高频的细节，因此我们不需要设置太大分辨率，在cpu端创建一个$64\times 64$大小的cubemap，然后调用着色器绘制6次写入这个cubemap当中，具体如下：

void IBLAuxiliary::convoluteDiffuseIntegral(int width, int height, unsigned int cubemapTexIndex, unsigned int irradianceTexIndex)
{
	// manager.
	TextureMgr::ptr texMgr = TextureMgr::getSingleton();
	ShaderMgr::ptr shaderMgr = ShaderMgr::getSingleton();
	// load shader.
	unsigned int shaderIndex = shaderMgr->loadShader("diffuseIntegral",
		"./glsl/diffuseIntegral.vert", "./glsl/diffuseIntegral.frag");
	// load cube mesh.
	Mesh::ptr cubeMesh = std::shared_ptr<Mesh>(new Cube(1.0f, 1.0f, 1.0f));
	// load framebuffer.
	FrameBuffer::ptr framebuffer = std::shared_ptr<FrameBuffer>(
		new FrameBuffer(width, height, "irradianceDepth", {}, true));
	// projection matrix and view matrix.
	glm::mat4 captureProjectMatrix = glm::perspective(glm::radians(90.0f), 1.0f, 0.1f, 10.0f);
	glm::mat4 captureViewMatrix[] =
	{
		glm::lookAt(glm::vec3(0.0f), glm::vec3(+1.0f, 0.0f, 0.0f), glm::vec3(0.0f, -1.0f, 0.0f)),
		glm::lookAt(glm::vec3(0.0f), glm::vec3(-1.0f, 0.0f, 0.0f), glm::vec3(0.0f, -1.0f, 0.0f)),
		glm::lookAt(glm::vec3(0.0f), glm::vec3(0.0f,+1.0f, 0.0f), glm::vec3(0.0f,  0.0f,+1.0f)),
		glm::lookAt(glm::vec3(0.0f), glm::vec3(0.0f,-1.0f, 0.0f), glm::vec3(0.0f,  0.0f,-1.0f)),
		glm::lookAt(glm::vec3(0.0f), glm::vec3(0.0f, 0.0f,+1.0f), glm::vec3(0.0f, -1.0f, 0.0f)),
		glm::lookAt(glm::vec3(0.0f), glm::vec3(0.0f, 0.0f,-1.0f), glm::vec3(0.0f, -1.0f, 0.0f)),
	};

	// begin to convolute.
	framebuffer->bind();
	glDisable(GL_BLEND);
	glDisable(GL_CULL_FACE);
	glEnable(GL_DEPTH_TEST);
	glDepthFunc(GL_LEQUAL);
	glClearColor(0.0f, 1.0f, 0.0f, 1.0f);
	glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
	GLuint irradianceTexId = texMgr->getTexture(irradianceTexIndex)->getTextureId();
	Shader::ptr shader = shaderMgr->getShader(shaderIndex);
	shader->bind();
	shader->setInt("environmentMap", 0);
	shader->setMat4("projectMatrix", captureProjectMatrix);
	texMgr->bindTexture(cubemapTexIndex, 0);
	for (unsigned int i = 0; i < 6; ++i)
	{
		shader->setMat4("viewMatrix", captureViewMatrix[i]);
		glFramebufferTexture2D(GL_FRAMEBUFFER, GL_COLOR_ATTACHMENT0,
			GL_TEXTURE_CUBE_MAP_POSITIVE_X + i, irradianceTexId, 0);
		cubeMesh->draw(false, 0);
	}
	shader->unBind();
	texMgr->unBindTexture(cubemapTexIndex);
	framebuffer->unBind();
}

  以上的预计算或者说预烘培步骤只需执行一次，然后我们就得到了一张辐照度立方体贴图，如下图3所示，可以看到类似对原来的立方体贴图做了一个模糊处理，但这并不是模糊处理。

图3 预烘培得到的辐照度纹理

四、预计算镜面反射积分

  接下来我们就把目标转到镜面反射的积分预计算上面。我们要预计算的就是渲染方程$(2)$中的镜面反射部分，如下所示：

$$\int_{\Omega}(\frac{DFG}{4(\omega_o\cdot n)(\omega_i\cdot n)})L_i(p,\omega_i)n\cdot \omega_id\omega_i =\int_{\Omega}f_r(p,\omega_i,\omega_o)L_i(p,\omega_i)n\cdot \omega_id\omega_i \tag {7}$$

  被积函数中的$f_r(p,\omega_i,\omega_o)$是之前提到的brdf函数，它与积分变量$\omega_i$有关，因此它不是一个常量，我们不能像前面的漫反射积分那样把brdf函数当作常数项并提出积分外。除此之外，我们还注意到最终的积分值还取决于出射方向或者说观察方向$\omega_o$，因此公式$(7)$的积分值取决于$n$和$\omega_o$这两个三维的向量。前面的漫反射积分仅取决于$n$，所以我们可以很容易地遍历所有的$n$，然后将积分值存储到立方体贴图中。但是在镜面反射积分这里，积分值取决于$n$和$\omega_o$，这使得问题变得复杂起来，这是因为我们不能同时用$n$和$\omega_o$去索引立方体贴图。Epic Games公司提出了提出了一种近似的方案，这种方法就是将公式$(7)$中的积分分成两部分：

$$\int_{\Omega}f_r(p,\omega_i,\omega_o)L_i(p,\omega_i)n\cdot \omega_id\omega_i \approx \int_{\Omega}L_i(p,\omega_i)d\omega_i * \int_{\Omega}f_r(p,\omega_i,\omega_o)n\cdot \omega_i d\omega_i \tag {8}$$

  其中$\int_{\Omega}L_i(p,\omega_i)d\omega_i$是对半球方向的辐射率进行积分，$\int_{\Omega}f_r(p,\omega_i,\omega_o)n\cdot
\omega_id\omega_i$则对半球方向的brdf函数进行积分，我们将原来的积分分成了这两部分，预计算然后保存到两张纹理当中。可以看到前面一部分的积分仅仅取决于法线向量$n$，后面一部分表面上取决于$n$和$\omega_o$，后面深入了解之后就会发现不用$n$和$\omega_o$去索引brdf函数的积分值。

1、Pre-Filtered Environment Map

  首先我们来看前面一部分的积分值，即$\int_{\Omega}L_i(p,\omega_i)d\omega_i$，这部分预计算得到的纹理为预过滤的环境贴图，它返回的是反射的像素值。我们知道，给定一个物体表面，若其表面的粗糙程度越高，则反射的内容越模糊。因此，为了把物体的粗糙度考虑进去，我们将考虑物体的粗糙度划分成几个等级，每个等级预计算一遍积分，并存储到cubemap的一个mipmap当中。我们为生成的预过滤环境贴图构建mipmap，越粗糙则存储到mipmap的等级越高，正如如下图4所示：

图4 不同粗糙度对应的mipmap

  表面越粗糙则反射得越模糊得现象其实涉及到一个反射波瓣的问题，所以的反射波瓣就是反射光线的分布范围。如下图5所示，对于一个光滑的完美镜面，其反射波瓣就是反射向量，越粗糙则反射波瓣越大，且基本上都是以反射向量为中心。

图5 不同粗糙度下的反射波瓣

  因此，当我们进行积分采样时，我们没有必要再在半球方向做一个均匀的采样，因为超出反射波瓣的采样是无效的。我们应该偏向反射波瓣采样，这里就涉及到了一个重要性采样的问题，重要性采样的内容在前面的光线追踪器Ray Tracer：进阶篇已经详细地展开过，这里就不再赘述。采用了重要性采样方法，我们求解积分方程的数值方法不再是黎曼积分法，而是蒙特卡洛积分法，蒙特卡洛积分同样在光线追踪渲染器中提及过。

  为了加速蒙特卡洛积分方法的收敛速度，Epic Games公司提出使用超均匀分布序列（Low-discrepancy Sequence）——Hammersley序列。相对于普通的伪随机数，Hammersley序列的随机数分布更加均匀，如下图6所示，将其应用到蒙特卡洛采样能够提升收敛速度。

图6 伪随机数vs超均匀分布随机数

  我们将使用的是二维的Hammersley序列。一个二维的Hammersley序列$H_N=\{x_1,…,x_N\}, N\geq 1$是散落分布在单位正方形内的点集，其数学定义为：

$$H_N=\{x_i= (\begin{matrix} i/N\\ \Phi_2(i) \end{matrix} ),\ \ for\ \ i=0,...,N-1 \} \tag {9}$$

  其中$\Phi_2(i)$是Van der Corput序列，它输入$i$，然后将$i$的二进制编码以小数点为对称做一个镜像操作，返回$[0,1)$的浮点数，其数学定义为：

$$\Phi_2(i)=\frac{a_0}{2}+\frac{a_1}{2^2}+...+\frac{a_r}{2^{r+1}} \tag {10}$$

  其中的$a_0a_1…a_n$是$i$的二进制编码每一位二进制位，即$i=a_0+a_1\cdot 2+a_2\cdot 2^2+a_3\cdot 2^3+…+a_r\cdot 2^r$。然后我们需要将这个二维的序列转换成我们对半球方向的三维采样，同样我们利用球面坐标和笛卡尔坐标之间的关系，首先将Hammersley序列$x_i=(u,v)^T\in H_N$映射到$(\phi,\theta)$，然后转换成笛卡尔坐标下的向量形式。一个均匀映射和一个余弦映射公式如下：

$$Uniform\ \ mapping= \{\begin{matrix} \theta=cos^{-1}(1-u)\\ \phi=2\pi v \end{matrix} \\ Cosinus\ \ mapping= \{\begin{matrix} \theta=cos^{-1}(\sqrt{(1-u)})\\ \phi=2\pi v \end{matrix} \tag {11}$$

  均匀映射就是将序列映射到一个均匀的分布，余弦映射则将序列映射到一个更偏向于半球中心轴上的分布。两者的区别如下图7所示。Hammersley序列可以通过位移操作快速地实现，具体代码见下面。

图7 均匀映射vs余弦映射

float radicalInverseVdc(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), radicalInverseVdc(i));
}

  紧接着我们要用Hammersley序列实现我们的重要性采样。前面已经提到过，我们将考虑物体的粗糙度，因为不同粗糙度下的反射波瓣大小不同。我们将粗糙度换分成5个等级，每个等级根据当前的粗糙度进行重要性采样。因此做重要性采样时我们需要根据粗糙度确定当前的反射波瓣大小，反射波瓣越大则采样范围越大。我们将结合之前PBR提到的法线分布函数，法线分布函数给定一个法线向量，它返回微平面法线与给定法线朝向一致的分布概率。Trowbridge-Reitz GGX法线分布函数的数学定义为：

$$NDF_{GGXTR}(n,h,\alpha)=\frac{\alpha^2}{\pi((n\cdot h)^2(\alpha^2-1)+1)^2} \tag {12}$$

  我们将法线分布函数与余弦映射结合起来做重要性采样：

$$Important\ \ sampling= \{\begin{matrix} \theta=cos^{-1}(\sqrt{\frac{1-u}{u(\alpha^2-1)+1}})\\ \phi=2\pi v \end{matrix} \tag {13}$$

  同样这也是Epic Games公司提出的映射方法，注意到与公式$(11)$中的余弦映射相比，公式$(13)$多了一个分母$u(\alpha^2-1)+1$取自公式$(12)$中的法线分布函数。当粗糙度$\alpha$增大时，余弦值减小，$\theta$取值范围越大，反射波瓣也就越大，这个就是公式$(13)$的大体思路。对于给定的Hammersley序列、法线向量N、粗糙度roughness，一个重要性采样代码如下所示：

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);
	vec3 H;
	H.x = cos(phi) * sinTheta;
	H.y = sin(phi) * sinTheta;
	H.z = cosTheta;
	
	// from tangent space to world space.
	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 = H.x * tangent + H.y * bitangent + H.z * N;
	return normalize(sampleVec);
}

  最后我们就利用重要性采样进行数值积分的计算，如下所示：

$$\int_{\Omega}L_i(p,\omega_i)d\omega_i \approx \frac1N\Sigma_{k=1}^NL_i(l_k) \tag {14}$$

  然后在实现时，Ephic Games公司发现再乘上一个权重$cos\theta_{l_k}$效果更加，因而实现的积分计算代码如下：

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

uniform float roughness;
uniform samplerCube environmentMap;

const float PI = 3.14159265359;

float radicalInverseVdc(uint bits);
vec2 hammersley(uint i, uint N);
vec3 importanceSampleGGX(vec2 Xi, vec3 N, float roughness);

void main()
{
	vec3 N = normalize(worldPos);
	vec3 V = N;

	const uint sampleCount = 1024u;
	float totalWeight = 0.0f;
	vec3 prefilteredColor = vec3(0.0f);
	for(uint i = 0u;i < sampleCount;++ i)
	{
		vec2 Xi = hammersley(i, sampleCount);
		// sample halfway vector.
		vec3 H = importanceSampleGGX(Xi, N, roughness);
		// reflect vector.
		vec3 L = normalize(2.0 * dot(V, H) * H - V);
		
		float NdotL = max(dot(N, L), 0.0);
		if(NdotL > 0.0f);
		{
			prefilteredColor += texture(environmentMap, L).rgb * NdotL;
			totalWeight += NdotL;
		}
	}
	prefilteredColor = prefilteredColor / totalWeight;
	
	fragColor = vec4(prefilteredColor, 1.0f);
}

  在cpu端生成多个mipmap层次的cubemap，对于每一个粗糙度等级，我们将预计算的结果存储到对应mipmap等级的cubemap纹理当中，最后我们就得到如图4所示的多个mipmap等级的Pre-Filtered Environment Map。

void IBLAuxiliary::convoluteSpecularIntegral(int width, int height, unsigned int cubemapTexIndex,
	unsigned int prefilteredTexIndex)
{
	// manager.
	TextureMgr::ptr texMgr = TextureMgr::getSingleton();
	ShaderMgr::ptr shaderMgr = ShaderMgr::getSingleton();
	// load shader.
	unsigned int shaderIndex = shaderMgr->loadShader("prefilterEnvMap",
		"./glsl/prefilterEnvMap.vert", "./glsl/prefilterEnvMap.frag");
	// load cube mesh.
	Mesh::ptr cubeMesh = std::shared_ptr<Mesh>(new Cube(1.0f, 1.0f, 1.0f));
	// projection matrix and view matrix.
	glm::mat4 captureProjectMatrix = glm::perspective(glm::radians(90.0f), 1.0f, 0.1f, 10.0f);
	glm::mat4 captureViewMatrix[] =
	{
		glm::lookAt(glm::vec3(0.0f), glm::vec3(+1.0f, 0.0f, 0.0f), glm::vec3(0.0f, -1.0f, 0.0f)),
		glm::lookAt(glm::vec3(0.0f), glm::vec3(-1.0f, 0.0f, 0.0f), glm::vec3(0.0f, -1.0f, 0.0f)),
		glm::lookAt(glm::vec3(0.0f), glm::vec3(0.0f,+1.0f, 0.0f), glm::vec3(0.0f,  0.0f,+1.0f)),
		glm::lookAt(glm::vec3(0.0f), glm::vec3(0.0f,-1.0f, 0.0f), glm::vec3(0.0f,  0.0f,-1.0f)),
		glm::lookAt(glm::vec3(0.0f), glm::vec3(0.0f, 0.0f,+1.0f), glm::vec3(0.0f, -1.0f, 0.0f)),
		glm::lookAt(glm::vec3(0.0f), glm::vec3(0.0f, 0.0f,-1.0f), glm::vec3(0.0f, -1.0f, 0.0f)),
	};

	// begin to filter.
	GLuint prefilteredTexId = texMgr->getTexture(prefilteredTexIndex)->getTextureId();
	Shader::ptr shader = shaderMgr->getShader(shaderIndex);
	shader->bind();
	shader->setInt("environmentMap", 0);
	shader->setMat4("projectMatrix", captureProjectMatrix);
	texMgr->bindTexture(cubemapTexIndex, 0);
	unsigned int maxMipLevels = 5;
	for (unsigned int mip = 0; mip < maxMipLevels; ++mip)
	{
		unsigned int mipWidth = width * std::pow(0.5, mip);
		unsigned int mipHeight = height * std::pow(0.5, mip);
		std::stringstream ss;
		ss << mip;
		FrameBuffer::ptr framebuffer = std::shared_ptr<FrameBuffer>(
			new FrameBuffer(mipWidth, mipHeight, "prefilteredDepth" + ss.str(), {}, true));
		framebuffer->bind();
		glDisable(GL_BLEND);
		glDisable(GL_CULL_FACE);
		glEnable(GL_DEPTH_TEST);
		glEnable(GL_TEXTURE_CUBE_MAP_SEAMLESS);
		glDepthFunc(GL_LEQUAL);
		glClearColor(0.0f, 1.0f, 0.0f, 1.0f);

		float roughness = (float)mip / (float)(maxMipLevels - 1);
		shader->setFloat("roughness", roughness);
		for (unsigned int i = 0; i < 6; ++i)
		{
			shader->setMat4("viewMatrix", captureViewMatrix[i]);
			glFramebufferTexture2D(GL_FRAMEBUFFER, GL_COLOR_ATTACHMENT0,
				GL_TEXTURE_CUBE_MAP_POSITIVE_X + i, prefilteredTexId, mip);
			glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
			cubeMesh->draw(false, 0);
		}
		framebuffer->unBind();
	}
	shader->unBind();
	texMgr->unBindTexture(cubemapTexIndex);
}

2、Pre-computing the BRDF

  接下来我们把目标方法第二部分的积分计算，即brdf函数积分：

$$\int_{\Omega}f_r(p,\omega_i,\omega_o)n\cdot \omega_id\omega_i \tag {15}$$

  公式$(15)$的积分值取决于三个变量，即$n\cdot \omega_o$、表面粗糙度以及$F_0$，$F_0$是菲涅尔方程的一个输入值。三个变量太多了，为了简化且方便预计算，我们设法将一些变量提出积分符号外：

$$\int_{\Omega}f_r(p,\omega_i,\omega_o)n\cdot \omega_id\omega_i =\int_{\Omega}f_r(p,\omega_i,\omega_o)\frac{F(\omega_o,h)}{F(\omega_o,h)}n\cdot \omega_id\omega_i\\ =\int_{\Omega}\frac{f_r(p,\omega_i,\omega_o)}{F(\omega_o,h)}F(\omega_o,h)n\cdot \omega_id\omega_i \tag{16}$$

  将菲涅尔项$F(\omega_o,h)=(F_0+(1-F_0)(1-\omega_o\cdot h)^5)$代入上式，得：

$$\int_{\Omega}\frac{f_r(p,\omega_i,\omega_o)}{F(\omega_o,h)}(F_0+(1-F_0)(1-\omega_o\cdot h)^5)n\cdot \omega_id\omega_i \\ = \int_{\Omega}\frac{f_r(p,\omega_i,\omega_o)}{F(\omega_o,h)}(F_0(1-(1-\omega_o\cdot h)^5)+(1-\omega_o\cdot h)^5)n\cdot \omega_id\omega_i \tag {17}$$

​   然后将公式$(17)$得到的结果分成两部分：

$$\int_{\Omega}\frac{f_r(p,\omega_i,\omega_o)}{F(\omega_o,h)}(F_0(1-(1-\omega_o\cdot h)^5)+(1-\omega_o\cdot h)^5)n\cdot \omega_id\omega_i\\ =F_o\int_{\Omega}f_r(p,\omega_i,\omega_o)(1-(1-\omega_o\cdot h)^5)n\cdot \omega_id\omega_i +\int_{\Omega}f_r(p,\omega_i,\omega_o)(1-\omega_o\cdot h)^5 n\cdot \omega_i d\omega_i \tag {18}$$

  注意，公式$(18)$中的是去掉了菲涅尔项的brdf函数，因为菲涅尔项被消去了。现在我们把$F_0$提出积分外了，公式$(18)$中的两项积分取决于$n\cdot \omega_o$和粗糙度roughness。我们构建这样的一个二维查找表，它的横轴坐标取值为$n\cdot \omega_o$，纵轴坐标取值为粗糙度roughness，我们渲染屏幕空间大小的四边形，遍历$n\cdot \omega_o$和粗糙度的所有取值，计算其对应的公式$(18)$中的两项积分的结果，存储为纹理的像素值。最后渲染时使用纹理坐标$(n\cdot \omega_o, roughness)$去索引像素值。

float geometrySchlickGGX(float NdotV, float roughness)
{
	float a = roughness;
	float k = (a * a) / 2.0f;
	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;
}

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

	float A = 0.0;
	float B = 0.0;
	
	vec3 N = vec3(0.0, 0.0, 1.0);
	
	const uint sampleCount = 1024u;
	for(uint i = 0u;i < sampleCount;++i)
	{
		vec2 Xi = hammersley(i, sampleCount);
		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);
			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(sampleCount);
	B /= float(sampleCount);
	return vec2(A, B);
}
void main()
{
	vec2 value = integrateBRDF(Texcoord.x, Texcoord.y);
	fragColor = vec4(value, 0.0f, 1.0f);
}

  在cpu端创建一个二维纹理，并送入一个屏幕大小的四边形进行预计算的渲染。

void IBLAuxiliary::convoluteSpecularBRDFIntegral(int width, int height, unsigned int brdfLutTexIndex)
{
	// manager.
	TextureMgr::ptr texMgr = TextureMgr::getSingleton();
	ShaderMgr::ptr shaderMgr = ShaderMgr::getSingleton();
	// load shader.
	unsigned int shaderIndex = shaderMgr->loadShader("genBrdfLUT",
		"./glsl/genBrdfLUT.vert", "./glsl/genBrdfLUT.frag");
	// load quad mesh.
	Mesh::ptr quadMesh = std::shared_ptr<Mesh>(new ScreenQuad());
	FrameBuffer::ptr framebuffer = std::shared_ptr<FrameBuffer>(
		new FrameBuffer(width, height, "brdfDepth", {}, true));
	Shader::ptr shader = shaderMgr->getShader(shaderIndex);

	framebuffer->bind();
	glDisable(GL_BLEND);
	glDisable(GL_CULL_FACE);
	glDisable(GL_DEPTH_TEST);
	glDepthFunc(GL_LEQUAL);
	glClearColor(0.0f, 1.0f, 0.0f, 1.0f);
	glClear(GL_COLOR_BUFFER_BIT);
	GLuint brdfLutTexId = texMgr->getTexture(brdfLutTexIndex)->getTextureId();
	glFramebufferTexture2D(GL_FRAMEBUFFER, GL_COLOR_ATTACHMENT0, GL_TEXTURE_2D, brdfLutTexId, 0);

	shader->bind();
	quadMesh->draw(false, 0);
	shader->unBind();
	framebuffer->unBind();
}

  然后我们就可以得到下面这张纹理。

图8 brdf积分查找表

五、计算渲染方程

  在前面的步骤中我们预计算获取了Irradiance Map、Pre-Filtered Environment Map以及BRDF Lookup Texture，最后我们就直接查找这些纹理，用以我们的光照计算。光照计算部分直接就是一开始提到的渲染方程，渲染方程中的积分项从纹理中直接获取，不再需要实时计算。

// ambient lighting.
vec3 ambientS = fresnelSchlickRoughness(max(dot(normal, viewDir), 0.0f), F0, roughness);
vec3 ambientD = vec3(1.0f) - ambientS;
ambientD *= (1.0 - metallic);
vec3 irradiance = texture(irradianceMap, normal).rgb;

vec3 R = normalize(reflect(-viewDir, normal));
const float MAX_REFLECTION_LOD = 4.0;
vec3 prefilteredColor = texture(prefilteredMap, R, roughness * MAX_REFLECTION_LOD).rgb;

vec2 envBrdf = texture(brdfLutMap, vec2(max(dot(normal, viewDir), 0.0f), roughness)).rg;
vec3 envSpecular = prefilteredColor * (ambientS * envBrdf.x + envBrdf.y);

vec3 ambient = (albedo * irradiance * ambientD + envSpecular) * ao;

fragColor.xyz = ambient + fragColor.xyz * shadow + pointLightRadiance;

六、实现效果

参考资料：

$[1]$ Hammersley Points on the Hemisphere

$[2]$ Real Shading in Unreal Engine 4, by Brian Karis, Epic Games

$[3]$ https://learnopengl.com/PBR/IBL/Diffuse-irradiance

$[4]$ https://learnopengl.com/PBR/IBL/Specular-IBL