【PBR系列六】基于物理的环境光照(上)：漫反射辐照度(Diffuse irradiance)

本文核心知识主要参照learnopengl-cn文章总结归并，并根据个人学习方向进行了筛选摘抄，如有错误或不完整之处，可参照原文阅读。

基于图像的光照（IBL）是对光源物体的技巧集合，与直接光照不同，它将周围环境当成一个大光源。IBL通常结合cubemap环境贴图，cubemap通常采集自真实的照片或从3D场景生成，这样可以将其用于光照方程：将cubemap的每个像素当成一个光源。这样可以更有效地捕获全局光照和常规感观，使得被渲染的物体更好地融入所处的环境中。

当基于图像的光照算法获得一些（全局的）环境光照时，它的输入被当成更加精密形式的环境光照，甚至是一种粗糙的全局光照的模拟。这使得IBL有助于PBR的渲染，使得物体渲染效果更真实。

在介绍IBL结合PBR之前，先回顾一下反射方程：

如之前所述，我们的主目标是解决所有入射光wi通过半球Ω的积分∫。与直接光照不同的是，在IBL中，每一个来自周围环境的入射光ωi都可能存在辐射，这些辐射对解决积分有着重要的作用。为解决积分有两个要求：

需要用某种方法获得给定任意方向向量ωi的场景辐射。
解决积分需尽可能快并实时。

对第一个要求，相对简单，采用环境cubemap。给定一个cubemap，可以假设它的每个像素是一个单独的发光光源。通过任意方向向量ωi采样cubemap，可以获得场景在这个方向的辐射。

获取任意方向向量ωi的场景辐射很简单，如下：

vec3 radiance = texture(_cubemapEnvironment, w_i).rgb;

对要求二，解决积分能只考虑一个方向的辐射，要考虑环境贴图的半球Ω的所有可能的方向ωi，但常规积分方法在片元着色器中开销非常大。为了有效解决积分问题，可采用预计算或预处理的方法。因此，需要深究一下反射方程：

可将上述的kd和ks项拆分：

拆分后，可分开处理漫反射和镜面反射的积分。先从漫反射积分开始。

一、漫反射辐照度(Diffuse irradiance)

仔细分析上面方程的漫反射积分部分，发现Lambert漫反射是个常量项（颜色c，折射因子kd和π）并且不依赖积分变量。因此，可见常量部分移出漫反射积分：

因此，积分只依赖ωi（假设p在环境贴图的中心）。据此，可以计算或预计算出一个新的cubemap，这个cubemap存储了用卷积（convolution）计算出的每个采样方向（或像素）ωo的漫反射积分结果。

卷积（convolution）是对数据集的每个入口应用一些计算，假设其它所有的入口都在这个数据集里。此处的数据集就是场景辐射或环境图。因此，对cubemap的每个采样方向，我们可以顾及在半球Ω的其它所有的采样方向。

为了卷积环境图，我们要解决每个输出ωo采样方向的积分，通过离散地采样大量的在半球Ω的方向ωi并取它们辐射的平均值。采样方向ωi的半球是以点p为中心以ωo为法平面的。

这个预计算的为每个采样方向ωo存储了积分结果的cubemap，可被当成是预计算的在场景中所有的击中平行于ωo表面的非直接漫反射的光照之和。这种cubemap被称为 辐照度图（Irradiance map）。

辐射方程依赖于位置p，假设它在辐照度图的中心。这意味着所有非直接漫反射光需来自于同一个环境图，它可能打破真实的幻觉（特别是室内）。渲染引擎用放置遍布场景的反射探头（reflection probe）来解决，每个反射探头计算其所处环境的独自的辐照度图。这样，点p的辐射率（和辐射）是与其最近的反射探头的辐照度插值。这里我们假设总是在环境图的中心采样。反射探头将在其它章节探讨。

下面是cubemap环境图(下图左）和对应的辐照度图（下图右）：

通过存储每个cubemap像素卷积的结果，辐照度图有点像环境的平均颜色或光照显示。从这个环境图采样任意方向，可获得这个方向的场景辐照度。

1.1 球体图（Equirectangular map）

球体图（Equirectangular map）有些文献翻译成全景图，它与cubemap不一样的是：cubemap需要6张图，而球体图只需要一张，并且存储的贴图有一定形变：

cubemap是可以通过一定算法转成球体图的，详见Converting a Cubemap into Equirectangular Panorama。

1.2 从球体图到立方体图

直接从球体图采样出环境光照信息是可能的，但它的开销远大于直接采样立方体图（cubemap）。因此，需要将球体图先转成立方体图，以便更好地实现后面的逻辑。当然，这里也会阐述如何从作为3D环境图的球体图采样，以便大家有更多的选择权。

为了将球体图映射到立方体图，首先需要构建一个立方体模型，渲染这个立方体模型的顶点着色器如下：

#version 450layout (location = 0) in vec3 inPos;
layout (location = 1) in vec3 inNormal;
layout (location = 2) in vec2 inUV;layout (binding = 0) uniform UBO
{mat4 projection;mat4 model;
} ubo;layout (location = 0) out vec3 outUVW;out gl_PerVertex
{vec4 gl_Position;
};void main()
{outUVW = inPos;gl_Position = ubo.projection * ubo.model * vec4(inPos.xyz, 1.0);
}

在像素着色器中，将会对变形的球体图的每个部位映射到立方体的每一边，具体实现如下：

#version 450layout (binding = 2) uniform samplerCube samplerEnv;
layout (binding = 3) uniform samplerCube irradianceEnv;//辐照度贴图测试用layout (location = 0) in vec3 inUVW;layout (location = 0) out vec4 outColor;layout (binding = 1) uniform UBOParams {vec4 lights[4];float exposure;float gamma;
} uboParams;void main()
{vec3 color = texture(samplerEnv, inUVW).rgb;//vec3 color = texture(irradianceEnv, inUVW).rgb;// Exposure mappingcolor = vec3(1.0) - exp(-color * 2.0);// Gamma correctioncolor = pow(color, vec3(1.0f / uboParams.gamma));outColor = vec4(color, 1.0);
}

现在，在之前渲染的球体上渲染环境贴图，效果应该如下图：

1.3 立方体贴图的卷积（辐射度图）

1.3.1 基础原理

辐射度图提供了漫反射部分的积分，该积分表示来自非直接的所有方向的环境光辐射之和。由于辐射度图被当成是无方向性的光源，所以可以将漫反射镜面反射合成环境光。

如本节教程开头所述，我们的主要目标是计算所有间接漫反射光的积分，其中光照的辐照度以环境立方体贴图的形式给出。我们已经知道，在方向 wi 上采样 HDR 环境贴图，可以获得场景在此方向上的辐射度 L(p,wi) 。虽然如此，要解决积分，我们仍然不能仅从一个方向对环境贴图采样，而要从半球 Ω 上所有可能的方向进行采样，这对于片段着色器而言还是过于昂贵。

然而，计算上又不可能从 Ω 的每个可能的方向采样环境光照，理论上可能的方向数量是无限的。不过我们可以对有限数量的方向采样以近似求解，在半球内均匀间隔或随机取方向可以获得一个相当精确的辐照度近似值，从而离散地计算积分 ∫ 。

然而，对于每个片段实时执行此操作仍然太昂贵，因为仍然需要非常大的样本数量才能获得不错的结果，因此我们希望可以预计算。既然半球的朝向决定了我们捕捉辐照度的位置，我们可以预先计算每个可能的半球朝向的辐照度，这些半球朝向涵盖了所有可能的出射方向 wo ：

给定任何方向向量 wi ，我们可以对预计算的辐照度图采样以获取方向 wi 的总漫反射辐照度。为了确定片段上间接漫反射光的数量（辐照度），我们获取以表面法线为中心的半球的总辐照度。那么在最终的PBR着色器中获取场景辐照度的方法就简化为：

vec3 irradiance = texture(irradianceMap, N);

1.3.2 实现步骤

现在，为了生成辐照度贴图，我们需要将环境光照求卷积，转换为立方体贴图。假设对于每个片段，表面的半球朝向法向量 N ，对立方体贴图进行卷积等于计算朝向 N 的半球 Ω 中每个方向 wi 的总平均辐射率。
samplerEnv是从天空盒中采样而来的 HDR 立方体贴图。有很多方法可以对环境贴图进行卷积，但是对于本教程，我们的方法是：对于立方体贴图的每个纹素，在纹素所代表的方向的半球 Ω 内生成固定数量的采样向量，并对采样结果取平均值。数量固定的采样向量将均匀地分布在半球内部。注意，积分是连续函数，在采样向量数量固定的情况下离散地采样只是一种近似计算方法，我们采样的向量越多，就越接近正确的结果。反射方程的积分 ∫ 是围绕立体角 dw 旋转，而这个立体角相当难以处理。为了避免对难处理的立体角求积分，我们使用球坐标 θ 和 ϕ 来代替立体角。

对于围绕半球大圆的航向角 ϕ ，我们在 0 到 2π 内采样，而从半球顶点出发的倾斜角 θ ，采样范围是 0 到 12π 。于是我们更新一下反射积分方程：

求解积分需要我们在半球 Ω 内采集固定数量的离散样本并对其结果求平均值。分别给每个球坐标轴指定离散样本数量 n1 和 n2 以求其黎曼和，积分式会转换为以下离散版本：

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

给定每个片段的积分球坐标，对半球进行离散采样，辐射辐照度贴图的顶点着色器过程代码如下：

#version 450layout (location = 0) in vec3 inPos;layout(push_constant) uniform PushConsts {layout (offset = 0) mat4 mvp;
} pushConsts;layout (location = 0) out vec3 outUVW;out gl_PerVertex {vec4 gl_Position;
};void main()
{outUVW = inPos;gl_Position = pushConsts.mvp * vec4(inPos.xyz, 1.0);
}

辐射辐照度贴图的片元着色器顶点着色器过程代码如下：

//使用卷积从环境地图生成辐照度立方体#version 450layout (location = 0) in vec3 inPos;
layout (location = 0) out vec4 irradiance ;
layout (binding = 0) uniform samplerCube samplerEnv;layout(push_constant) uniform PushConsts {layout (offset = 64) float deltaPhi;layout (offset = 68) float deltaTheta;
} consts;#define PI 3.1415926535897932384626433832795void main()
{vec3 N = normalize(inPos);vec3 up = vec3(0.0, 1.0, 0.0);vec3 right = normalize(cross(up, N));up = cross(N, right);const float TWO_PI = PI * 2.0;const float HALF_PI = PI * 0.5;vec3 color = vec3(0.0);uint sampleCount = 0u;for (float phi = 0.0; phi < TWO_PI; phi += consts.deltaPhi) {for (float theta = 0.0; theta < HALF_PI; theta += consts.deltaTheta) {// 球面到笛卡尔(切空间)vec3 tempVec = cos(phi) * right + sin(phi) * up;// 相切空间到世界空间vec3 sampleVector = cos(theta) * N + sin(theta) * tempVec;color += texture(samplerEnv, sampleVector).rgb * cos(theta) * sin(theta);sampleCount++;}}irradiance  = vec4(PI * color / float(sampleCount), 1.0);
}

其中推入常量结构体为：

     // Pipeline layoutstruct PushBlock {glm::mat4 mvp;// Sampling deltasfloat deltaPhi = (2.0f * float(M_PI)) / 180.0f;float deltaTheta = (0.5f * float(M_PI)) / 64.0f;} pushBlock;

我们以一个固定的推入管线中的deltaPhi 和deltaTheta 增量值遍历半球，减小（或增加）这个增量将会增加（或减少）精确度。

在两层循环内，我们获取一个球面坐标并将它们转换为 3D 直角坐标向量，将向量从切线空间转换为世界空间，并使用此向量直接采样 HDR 环境贴图。我们将每个采样结果加到 irradiance，最后除以采样的总数，得到平均采样辐照度。请注意，我们将采样的颜色值乘以系数 cos(θ) ，因为较大角度的光较弱，而系数 sin(θ) 则用于权衡较高半球区域的较小采样区域的贡献度。

现在，完成这个流程之后，我们应该得到了一个预计算好的辐照度图，可以直接将其用于IBL 计算。为了查看我们是否成功地对环境贴图进行了卷积，让我们将天空盒的环境采样贴图替换为辐照度贴图：

1.4 PBR和非直接辐射度光照（indirect irradiance lighting）

首先，得声明预计算出的辐射度图的立方体采样器：

layout (binding = 2) uniform samplerCube samplerIrradiance;

通过表面的法线，获得环境光可以简化成下面的代码：

vec3 ambient= texture(samplerIrradiance, N).rgb;

尽管如此，在之前所述的反射方程中，非直接光依旧包含了漫反射和镜面反射两个部分，所以我们需要加个权重给漫反射。下面采用了菲涅尔方程来计算漫反射因子：

vec3 F = F_Schlick(max(dot(N, V), 0.0), F0);
vec3 kD = 1.0 - F;
vec3 irradiance = texture(samplerIrradiance, N).rgb;
vec3 diffuse = irradiance * ALBEDO;

其中F_Schlick为：

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

由于环境光来自在半球内所有围绕着法线N的方向，没有单一的半向量去决定菲涅尔因子。为了仍然能模拟菲涅尔，这里采用了法线和视线的夹角。之前的算法采用了受表面粗糙度影响的微平面半向量，作为菲涅尔方程的输入。这里，我们并不考虑粗糙度，表面的反射因子被视作相当大。

非直接光照将沿用直接光照的相同的属性，所以，期望越粗糙的表面镜面反射越少。由于不考虑表面粗糙度，非直接光照的菲涅尔方程强度被视作粗糙的非金属表面（下图）（为了演示目的略微夸大）：

为了缓解这个问题，可在Fresnel-Schlick方程注入粗糙度项：

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

考虑了表面粗糙度后，菲涅尔相关计算最终如下：

vec3 F = F_SchlickR(max(dot(N, V), 0.0), F0, roughness);
vec3 kD = 1.0 - F;
vec3 irradiance = texture(samplerIrradiance, N).rgb;
vec3 diffuse = irradiance * ALBEDO;

如上所述，实际上，基于图片的光照计算非常简单，只需要单一的cubemap纹理采样。大多数的工作在于预计算或卷积环境图到辐射度图。

加入了IBL的漫反射辐照度渲染效果如下：

二、镜面反射 IBL(Specular IBL)

详见本系列下一节