【cg】【pbr】基于物理的渲染实现篇之间接光照(上)

发表于 2019-11-27 更新于 2022-07-29 分类于 game_graphics 阅读次数：本文字数： 8.5k 阅读时长 ≈ 8 分钟

前言

上篇使用Cook-Torrance BRDF模型实现了精准光源直接光照的计算。本篇继续基于PBR的理论基础，通过对微平面进行半球领域近似积分来计算周围环境作用于物体上的间接光源。

本人能力有限，对其中的某些步骤尚有疑惑之处，希望以后通过进一步的深入学习可以真正地理解每一个过程，共勉。

下面将需要在本篇引用到的理论篇推导过的公式列举一下，方便后面引用。

\[ L(p, \vec{v}) = \int_{\vec{l}_i \in \Omega} f(p, \vec{l}_i, \vec{v}) \ast L(p, \vec{l}_i) \ast (\vec{n} \cdot \vec{l}_{i}) \ast d\vec{l}_i \tag{0.19} \]

\[ f(p, \vec{l}, \vec{v}) = f_{lambert} (p, \vec{l}, \vec{v}) + f_{cook\_torrance} (p, \vec{l}, \vec{v}) \tag{0.21} \]

\[ f_{lambert} (p, \vec{l}, \vec{v}) = k_d \ast \frac {C_{diffuse}} {\pi} \tag{0.22} \]

IBL简介

基于图像的光照(Image Based Lighting)是计算PBR间接光照的一种方法，用于计算除了精准光源之外的其他环境光照。

该方法将物体周围的整体环境作为一个大的光源作用于物体。所以我们需要对物体的每一个平面的半球领域进行积分才能获得最终的光照强度，即对每一个微平面来说，它会接收任意方向上的光源。

我们同样将间接光照分为漫反射项和镜面反射项两部分进行计算。

对于漫反射部分，由于其低频域的特性(可以理解为是很模糊的)，我们将预烘焙确定的积分项到一个较小分辨率的立方体贴图上，然后在片段着色器里对此立方体再采样进行计算。

对于镜面反射部分，分别使用了两种不同的方案实现。

第一种方案是使用重要性采样加速的蒙特卡洛积分对反射率方程进行近似，通过较少的采样数实现其实时性，即可以在每个片段着色器里采样并进行积分计算，其中主要涉及到重要性采样的原理和实现，以及一些概率论相关的采样分布映射的计算。

第二种方案是Epic公司在ue4里采用的实现方案，它将前一种方案里的蒙特卡洛积分进一步近似地划分成两个不同的积分项，即LD项和FGD项，这两项都可以通过预先计算得出结果。我们将LD项预烘焙到一张Cubemap上，将FGD项预烘焙到一张2D纹理上，然后再在片段着色器里直接采样进行计算，不需要实时地在片段着色器里求积分，性能很高，这部分内容限于篇幅可能会放到下篇来推导和实现。

由Cook-Torrance BRDF公式\((0.21)\)可将反射率方程\((0.19)\)改写为如下格式

\[ L(p, \vec{v}) = \int_{\vec{l}_i \in \Omega} f_d \ast L(p, \vec{l}_i) \ast (\vec{n} \cdot \vec{l}_{i}) \ast d\vec{l}_i + \int_{\vec{l}_i \in \Omega} f_s \ast L(p, \vec{l}_i) \ast (\vec{n} \cdot \vec{l}_{i}) \ast d\vec{l}_i \tag{1.1} \]

其中简记 \(f_d = f_{lambert} (p, \vec{l}, \vec{v})\), \(f_s = f_{cook\\_torrance} (p, \vec{l}, \vec{v})\)。

此格式将反射率积分分为左边的漫反射部分和右边的镜面反射部分。下面分别来求。

Diffuse IBL

本小节所要实现的便是漫反射部分，即

\[ L_{diffuse} (p, \vec{v}) = \int_{\vec{l}_i \in \Omega} f_d \ast L(p, \vec{l}_i) \ast (\vec{n} \cdot \vec{l}_{i}) \ast d\vec{l}_i \tag{1.2} \]

又由公式\((0.22)\)，并提取积分项里的常量，可得

\[ L_{diffuse} (p, \vec{v}) = k_d \ast \frac {C_{diffuse}} {\pi} \int_{\vec{l}_i \in \Omega } L(p, \vec{l}_i) \ast (\vec{n} \cdot \vec{l}_{i}) \ast d\vec{l}_i \tag{1.3} \]

如果我们直接用公式\((1.3)\)进行积分的话，会发现代码很难写，因为使用向量作为微元进行积分的话，是很难说对向量如何进行步进采样的。

所以我们进一步将该公式转换为极角坐标系的形式，然后进行以极角坐标系为基础的逐角度积分，就变得很惬(qia)意了。

设入射方向\(\vec{l}_i\)与法线\(\vec{n}\)的夹角为\(\theta\)，设入射方向偏离\(x\)轴的夹角为\(\phi\)。则有

\[ \vec{n} \cdot \vec{l}_i = \cos{\theta} \tag{1.4} \]

\[ d\vec{l}_i = \sin{\theta} d\phi \ast d\theta \tag{1.5} \]

那么，公式\((1.3)\)可转换为以极角坐标系为基础的逐角度积分，即

\[ L_{diffuse} (p, \vec{v}) = k_d \ast \frac {C_{diffuse}} {\pi} \iint_{\theta = 0, \phi = 0} ^{\theta <\frac{\pi}{2}, \phi < 2\pi} L(p, \theta_i, \phi_i) \ast \cos{\theta_i} \ast \sin{\theta_i} d\phi_i \ast d\theta_i \tag{1.6} \]

这是一个在半球领域的连续型二重积分，我们只有转换为离散的积分才可以使用代码近似地计算出结果。

这里使用黎曼和(Riemann Sum)进行离散化近似。对于二重积分而言，黎曼和使用一个步进值，步进\(\theta\)和\(\phi\)，然后求其平均值，以获得和原来二重积分对应的面积相近似的面积值。

假设在\(\theta\)方向上步进了\(N_\theta\)步，则\(\theta\)每步的步进值为

\[ d\theta_i = \frac{\frac{\pi}{2}}{N_\theta} = \frac{\pi}{2N_\theta} \tag{1.7} \]

同理，设在\(\phi\)方向上步进了\(N_\phi\)步，则\(\phi\)每步的步进值为

\[ d\phi_i = \frac{2\pi}{N_\phi} \tag{1.8} \]

将\((1.7)\)、\((1.8)\)代入\((1.6)\)得使用黎曼和进行计算时的公式为

\[ \require{cancel} \\ L_{diffuse} (p, \vec{v}) \approx k_d \frac {C_{diffuse}} {\cancel{\pi}} \ast \frac{\cancel{\pi}}{\cancel{2}N_\theta} \frac{\cancel{2}\pi}{N_\phi} \sum_{j = 0}^{N_\phi} \sum_{i = 0}^{N_\theta} L(p, \theta_i, \phi_j) \ast \cos{\theta_i} \ast \sin{\theta_i} \tag{1.9} \]

即

\[ L_{diffuse} (p, \vec{v}) \approx k_d C_{diffuse} \frac{\pi}{N_\theta N_\phi} \sum_{j = 0}^{N_\phi} \sum_{i = 0}^{N_\theta} L(p, \theta_i, \phi_j) \ast \cos{\theta_i} \ast \sin{\theta_i} \tag{1.10} \]

其中除了\(k_d\)、\(C_{diffuse}\)项为外部变量外，剩下的都可以进行预处理积分。

实现

有了上面的理论基础，下面就可以开始实现了，但是这个实现需要建立在一个很好的理解之上的，不然可能得不到想要的结果。

为了更好地理解这个过程，我问了自己两个问题，首先，我们现在拥有什么？其次，我们想要输出什么？

好在我确实可以回答出来这两个问题，所以我实现了，虽然现在看来显得有些轻松愉快，但在第一次coding的时候还是踉跄了许久。

首先，我们现在拥有的，无非就是一张立方体贴图，一张代表着周围环境的光源信息，这是我们珍贵的源材料。

其次，我们想要输出的，无非也是一张立方体贴图，一张记录着每个方向上最终的光强积分的表。

拥有了这张表，我们就可以将任何一个确定的方向，作为这张立方体贴图的纹理坐标，然后在这个坐标点上进行采样，对应的就是这个方向上的间接漫反射光强的积分(即公式\((1.10)\)的预处理部分)。

所以这整个实现过程，无非就是在绘制一个skybox，输入其顶点和纹理坐标，传到片段着色器，在片段着色器里采样立方体贴图源材料作为光源。只不过每个片段输出的颜色要通过\((1.10)\)进行积分。

下面来看下我的片段着色器中的实现。

/* pbr_ibl_diffuse_bake.frag */

#version 330 core

const float pi = acos(-1.0);

// ================================================================================
// in & out

/* 最终输出的颜色 -- 积分结果 */
out vec4 r_color;

/* 输入 -- 立方体贴图的纹理坐标 */
in O_VS {
	vec3 tex_coord;
} i_fs;

// ================================================================================
// uniform

/* 输入 -- 立方体贴图 -- 光源材料 */
uniform samplerCube	u_texture;

// ================================================================================
// pre dec

/* TBN矩阵的各向量 -- 切线空间转世界空间用 */
vec3 t_normal;
vec3 t_up;
vec3 t_right;

/* 极坐标转轴坐标 */
vec3 spherical_to_cartesian(float phi, float theta);

/* 积分 -- delta为每次步进值 */
vec3 ibl_diffuse_integral(float delta);

// ================================================================================

void main() {

    /* 切线空间向量基 -- 计算TBN矩阵 */
    t_normal = normalize(i_fs.tex_coord);           // N
    t_up = vec3(0.0, 1.0, 0.0);
    t_right = normalize(cross(t_up, t_normal));     // T
    t_up = normalize(cross(t_normal, t_right));     // B

    /* 黎曼和积分 */
    vec3 t_color = ibl_diffuse_integral(0.05);

    r_color = vec4(t_color, 1.0);

}

// ================================================================================

vec3 spherical_to_cartesian(float phi, float theta) {
    return vec3(
        sin(theta) * cos(phi),
        sin(theta) * sin(phi),
        cos(theta)
    );
}

// ================================================================================

vec3 ibl_diffuse_integral(float delta) {
    vec3 res = vec3(0.0);
    int cnt = 0;

    // 对半球领域积分 -- 黎曼和近似
    for(float phi = 0.0; phi < 2.0 * pi; phi += delta) {
        for(float theta = 0.0; theta < 0.5 * pi; theta += delta) {
            // 获得轴坐标系下的方向向量 -- 但是是在切线空间的 -- 以当前片段的切线为中心轴的半球领域
            vec3 t_uv_tg = spherical_to_cartesian(phi, theta);
            // 要转化成世界空间的方向，才能输入的光源立方体贴图的正确位置采样 -- 乘以TBN矩阵
            vec3 t_uv = t_uv_tg.x * t_right + t_uv_tg.y * t_up + t_uv_tg.z * t_normal;

            // 这一步就是对光源进行采样  -- 正确的方向
            vec4 t_Li = texture(u_texture, normalize(t_uv));

            // 积分
            res += t_Li.rgb * cos(theta) * sin(theta);

            // cnt实际上是 N_theta * N_phi
            ++cnt;
        }
    }

    // 最终结果
    res = pi * (1.0 / float(cnt)) * res;
    return res;
}

一切都很顺理成章，只要保证方向正确！

注意实现中有一步切换空间转世界空间的，如果读者实现过或者熟悉法线贴图的话就无需解释了。不熟悉也没关系，笔者确有一篇推导和实现法线贴图的文章，还没整理出来，先在这里打个广告。

但是这里还是需要解释下，之所以要有这一步转换，是因为我们在黎曼和近似的时候所拿到的\(\theta\)和\(\phi\)，不管是哪个片段(像素)，其值都是一样的，都是[0.0,\(\frac{\pi}{2}\)]和[0.0,\(\pi\)]。这肯定是不对的，这样一来每一个片段所得到的积分结果其实都是法线为\(\vec{n} = (0.0, 0.0, 1.0)\)方向上的积分结果。即这里步进时取到的\(\theta\)和\(\phi\)是将当前片段的法线做为主朝向得到的相对值，但是我们对光源立方体贴图的采样却必须要用绝对的方向才可以拿到正确的值，所以把这个方向通过TBN矩阵转换成世界空间下的方向就可以了。

最终由一个光源立方体贴图(左图)得到的积分值存到RT上或者保存到硬件上之后，呈现出来的便像右图一样模(mǒ)糊(hū)。这也说明了间接漫反射光是一种低频域的信号，模(mǒ)糊(hū)就对了。

使用

将上一步得到的结果放到片段着色器里采样，并乘以\(k_d \ast C_{diffuse}\)，便可得到最终的辐射率。

但是需要注意的是，由于之前的\(k_d\)项是依赖于菲涅尔方程的(见[直接光照篇]公式\((1.11)\))，而菲涅尔方程则依赖于微观法线\(\vec{h}\)(见[直接光照篇]公式\((1.9)\))。

但是微观法线只在预计算时有不同的光源方向才可确定，但是这里并没有一个具体的光源方向而是一个半球领域，所以这里使用宏观法线\(\vec{n}\)代替微观法线\(\vec{h}\)。

但是这样一来，菲涅尔系数就失去了粗糙度对微观法线数量的影响，所以我们使用加入粗糙度作为输入参数的菲涅尔方程，如下所示。

\[ t = max⁡(1-roughness, F_0) \]

\[ F_{schlick\_roughness} (p, \vec{l}, \vec{v}) = F_0 + (t - F_0) * (1 - (\vec{n} \cdot \vec{v}))^5 \tag{1.11} \]

最终使用ibl_diffuse预处理积分计算间接漫反射光照的代码如下(省略了不重要的细节)。

// ================================================================================
// in & out

out vec4 r_color;

// ================================================================================
// uniform

/* 采样上述步骤得到的ibl_diffuse积分图 -- 一个立方体贴图 */
uniform samplerCube u_ibl_diffuse_map;

// ================================================================================

/* 带roughness的F函数 */
vec3 brdf_f_fresnel_schlick_roughness(float n_o_v, vec3 f0, float roughness) {
    vec3 t_s = max(vec3(1.0 - roughness), f0);
    return f0 + (t_s - f0) * pow(1.0 - n_o_v, 5.0);
}

/* ibl_diffuse */
vec3 pbr_ibl_diffuse(vec3 n, vec3 v,
                     vec3 albedo, vec3 f0,
                     float roughness, float metallic,
                     vec3 irradiance) {

    float n_o_v = max(dot(n, v), 0.0);

    vec3 f90 = brdf_f_f0(f0, albedo, metallic);                         // 获取真正的 f0 -- 有金属性影响
    vec3 f = brdf_f_fresnel_schlick_roughness(n_o_v, f90, roughness);   // 菲涅尔方程 -- 带roughness的

    // get k_s and k_d
    vec3  k_s = f;                                                      // 镜面反射系数 -- 等于菲涅尔方程的值
    vec3  k_d = (vec3(1.0) - k_s) * (1.0 - metallic);                   // 漫反射系数 -- 考虑金属度

    return k_d * albedo * irradiance;
}

// ================================================================================

void main()
{

    // 获得光源方向 -- 根据观察方向 和 法线方向
    vec3 t_light_dir = reflect(-t_view_dir, t_normal);

    // 从光源方向采样ibl_diffuse积分
    vec3 t_ibl_diffuse_irradiance = texture(u_ibl_diffuse_map, t_light_dir).rgb;

    // 计算间接漫反射光照
    vec3 t_ibl_diffuse = pbr_ibl_diffuse(
        t_normal, t_view_dir,
        t_c_diffuse, t_f0,
        t_roughnes, t_metallic,
        t_ibl_diffuse_irradiance
    );

    r_color = vec4(t_ibl_diffuse, 1.0);
}

最终添加了间接漫反射光照的效果如下图所示。

可以拿来跟直接光照篇作一下对比，似乎差距不大，但是可以看到直接光照篇全黑的地方，在这里似乎有了一些亮堂的地方，谁又知道这么细微的差别，背后却多了这么多次积分，这也许就是图形学的魅力所在吧。

结语

本文中提到的有关法线贴图的部分，仅涉及到一些小小的矩阵运算，近期如果需求不多的话会尽快整理出来。

本来打算把Specular IBL也给一口气写了，但是为了尽可能清楚详细地把问题表述清楚，加上自己文笔不是很好，所以可能需要更多的篇幅来填充内容。所以后面关于间接镜面反射光照的重头戏，打算另启一到两篇来表述。

\(^o^)/