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

发表于 更新于 分类于 阅读次数： 本文字数： 9.7k 阅读时长 ≈ 9 分钟

前言

书接前文，让我们来继续实现specular IBL。上回书说到我们可以使用暴力的手段，在片段着色器里实时计算积分来求得最终的光强，只不过采样的时候通过一些数学方法(重要性采样)来提高效率。但是Epic提供了一个更高效的方案，虽然损失了一些精度，但是图形学的本质就是欺骗嘛。

此处列一下后面需要引用到的以前的文章里介绍过的方程。

\[ L_{specular}(p, \vec{v}) = \int_{\vec{l}_i \in \Omega} \color{red}{f_s}\color{black}{} \ast L(p, \vec{l}_i) \ast (\vec{n} \cdot \vec{l}_{i}) \ast d\vec{l}_i \tag{0.1} \]

重要性采样近似方程。

\[ L_{specular}(p, \vec{v}) \approx \frac{1}{N} \sum_{i = 0}^{i \lt N} \frac { f_s \ast L(p, \vec{l_i}) \ast (\vec{n} \cdot \vec{l_i}) } { p_{s}(p, \vec{l_i}, \vec{v}) } \tag{0.14} \]

近似

观察公式\((0.1)\)，现在我们来思考一下，为什么非要在片段着色器里实时地去算积分呢？

笔者认为主要是因为每个片段对应的顶点信息不同，主要是其法线方向不同，而且我们的相机在每一时刻的位置也不同，观察方向也不同，所以就会导致对于光源环境立方体贴图，每次采样点位置不一样，采集的光源拿过来积分，得到的积分结果也不一样。

但如果不考虑光源的话，我们似乎可以求出任意法线方向和观察方向上brdf的积分值。所以我们先把公式\((0.1)\)拆分成光源部分和brdf部分，这一步叫做Split Sum Approximation。

而拆分方法就是，把光源项\(L(p, \vec{l_i})\)去掉，剩下的作为brdf部分，并且设另外一部分为\(L_d(p, \vec{v})\)，则有

\[ L_{specular}(p, \vec{v}) \approx L_c(p, \vec{v}) \ast ( \frac{1}{N} \sum_{i = 0}^{i \lt N} \frac { f_s \ast (\vec{n} \cdot \vec{l_i}) } { p_{s}(p, \vec{l_i}, \vec{v}) } ) \tag{1.1} \]

那么

\[ L_c(p, \vec{v}) \approx \frac { \frac{1}{N} \sum_{i = 0}^{i \lt N} \frac { f_s \ast L(p, \vec{l_i}) \ast (\vec{n} \cdot \vec{l_i}) } { p_{s}(p, \vec{l_i}, \vec{v}) } } { \frac{1}{N} \sum_{i = 0}^{i \lt N} \frac { f_s \ast (\vec{n} \cdot \vec{l_i}) } { p_{s}(p, \vec{l_i}, \vec{v}) } } \tag{1.2} \]

展开brdf函数\(f_s\)和pdf函数\(p_s\)得

\[ L_c(p, \vec{v}) \approx \frac { \sum_{i = 0}^{i \lt N} \frac { F \ast G \ast L(p, \vec{l_i}) \ast (\vec{v} \cdot \vec{h_i}) } { (\vec{n} \cdot \vec{v}) \ast (\vec{n} \cdot \vec{h_i}) } } { \sum_{i = 0}^{i \lt N} \frac { F \ast G \ast (\vec{v} \cdot \vec{h_i}) } { (\vec{n} \cdot \vec{v}) \ast (\vec{n} \cdot \vec{h_i}) } } \tag{1.3} \]

有先知做过实现，如果将菲涅尔系数F作为常量约去，与不约去直接计算，获得的结果从视觉上看不到差别，所以这里的分子和分母里的F项也可以约去。

注意到分子和分母有相同的部分，其实分子只比分母多了一个\(L(p, \vec{l_i})\)，所以我们记共同部分为\(W\)，即

\[ W =\frac { G \ast (\vec{v} \cdot \vec{h_i}) } { (\vec{n} \cdot \vec{v}) \ast (\vec{n} \cdot \vec{h_i}) } \tag{1.4} \]

那么\((1.3)\)变为

\[ L_c(p, \vec{v}) \approx \frac { \sum_{i = 0}^{i \lt N} L(p, l_i) \ast W } { \sum_{i = 0}^{i \lt N} W } \tag{1.5} \]

ok，现在，又有先知做了个实验，发现当\(\vec{n} = \vec{v}\)时，有\(W = G\)(分子分母都消去了)，此时无肉眼可见的差别，嗯，就这么定了。

ok，现在，双有先知来做实验了，发现\(W = \vec{n} \cdot \vec{l_i}\)时，与\(W = G\)时相比，也没有啥由此可见的差别，嗯嗯，就这么定了。。

(笔者并不是很能理解这部分，所以先这么定Orz。)

那么\(L_c\)就变成了如下模样。

\[ L_c(p, \vec{v}) \approx \frac { \sum_{i = 0}^{i \lt N} L(p, l_i) \ast (\vec{n} \cdot \vec{l_i}) } { \sum_{i = 0}^{i \lt N} (\vec{n} \cdot \vec{l_i}) } \tag{1.6} \]

代入\((1.1)\)，得

\[ L_{specular}(p, \vec{v}) \approx ( \underbrace { \frac { \sum_{i = 0}^{i \lt N} L(p, l_i) \ast (\vec{n} \cdot \vec{l_i}) } { \sum_{i = 0}^{i \lt N} (\vec{n} \cdot \vec{l_i}) } }_{LD} ) \ast ( \underbrace { \frac{1}{N} \sum_{i = 0}^{i \lt N} \frac { f_s \ast (\vec{n} \cdot \vec{l_i}) } { p_{s}(p, \vec{l_i}, \vec{v}) } }_{FGD} ) \tag{1.7} \]

这便是Epic对反射率方程进行近似的大致过程了，省略了一些笔者也不明白的细节Orz。

其中左边括号内的部分称为LD项，笔者称之为光源积分贴图，而右边括号内的部分称为FGD项，官方叫法是brdf积分贴图。

显然这两项都是可以预处理出来的，下面分别来看。

LD项

\[ LD(p, \vec{n}) = \frac { \sum_{i = 0}^{i \lt N} L(p, l_i) \ast (\vec{n} \cdot \vec{l_i}) } { \sum_{i = 0}^{i \lt N} (\vec{n} \cdot \vec{l_i}) } \tag{1.8} \]

从公式中也可以看出，LD项还是使用重要性采样来从光源立方体贴图上采样N个点。 所以我们的原料是一张环境立方体贴图。

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

拥有了这张表，我们就可以将任何一个确定的方向，作为这张立方体贴图的纹理坐标，然后在这个坐标点上进行采样，对应的就是这个方向上的LD项的值。

所以这整个预处理过程，无非就是在绘制一个skybox，输入其顶点和纹理坐标，传到片段着色器，在片段着色器里采样立方体贴图源材料作为光源。

而对于每个像素，输出的便是使用\((1.8)\)计算得到的积分值。

回想一下我们在计算diffuse IBL的时候，也是类似的过程。

还有一个问题需要考虑，我们知道重要性采样最终得到的坐标，跟粗糙度\(\alpha = roughness^2\)有关，所以对于每一个特定的方向，根据粗糙度的不同，自然采样点也不同，计算出的LD值也不同，那么保存所有情况的结果就不可能仅使用一个立方体贴图就能存的下的。

这个时候我们利用mipmap来做粗糙度的映射，即将不同的粗糙度采样的结果，分别存在到不同的mipmap层级里去，设置当前mipmap层级为CurMipLevel，最大mipmap层级为MaxMipLevel，则有

\[ roughness = \frac{CurMipLevel}{MaxMipLevel} \tag{1.9} \]

其中MaxMipLevel可以根据当前贴图的像素分辨率w来确定。

\[ MaxMipLevel = \lceil \log_2(w) \rceil \tag{1.10} \]

将使用以上公式计算好的roughness传到相应的shader之后，就可以在片段着色器里进行积分计算了。代码如下。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
// pbr_ibl_ld_bake.frag

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

// 原材料 -- 光源立方体贴图
uniform samplerCube     u_texture;
// 粗糙度 -- 根据mipmap层级在cpu计算好传过来的
uniform float           u_roughness;

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

// 积分
vec3 ibl_ld_integral(vec3 n, float roughness) {
    // 总光强值积分
    vec3 res = vec3(0.0);

    // 假设 view = normal 和 reflect = normal
    vec3 v = n;
    vec3 r = n;

    // (1.8)分子
    float tot_weight = 0.0;

    for(uint i = 0u; i < sampler_num; ++i) {
        // hammersley 采样
        vec2 uv = hammersley(i, sampler_num);
        // 重要性采样
        vec3 h = importance_sampling_ggx(uv, roughness, n);
        // 从 v 和 h 逆算 l -- 光源方向
        vec3 l = normalize(2.0 * dot(v, h) * h - v);

        // 常用cos值
        float n_o_v = max(dot(n, v), 0.0);
        float n_o_l = max(dot(n, l), 0.0);
        float n_o_h = max(dot(n, h), 0.0);
        float h_o_v = max(dot(h, v), 0.0);

        if(n_o_l > 0.0) {
            // 法线分布函数
            float d = brdf_d_tr_ggx(n_o_h, roughness);

            // 计算 pdf
            float pdf = d * n_o_h / (4.0 * h_o_v);
            // 根据pdf 获得 lod -- 降噪
            float lod = pbr_ibl_lod_get(pdf, 512, 512, sampler_num);

            // 从光源采样，然后代入(1.8)的分子进行积分
            res += textureLod(u_texture, l, lod).rgb * n_o_l;
            // res += texture(u_texture, l).rgb * n_o_l;

            // (1.8)分子项
            tot_weight += n_o_l;
        }

    }

    // 公式(1.8)
    res /= tot_weight;

    return res;
}

最终我们将渲染到RT上的图像保存到本地，可以看到类似下图(不记得是在哪盗的图了，自己的结果当时没有保存Orz)。

LD项积分贴图

上图虽然看起来和原图像差别不大，但其实质已经改变了，当我们将其绘制到一个skybox上的时候，如果从某个纹理坐标处采样，得到的颜色值就不是那个方向的光源了，而是\(\vec{n}\)或者\(\vec{v}\)或才\(\vec{R}\)(反射向量)为这个方向时的总的光源积分项(即LD项)了。

总结一下过程。

输入：一个立方体贴图，用作光源

输出：一个有多级mipmap的立方体贴图

过程：绘制立方体，在像素着色器里根据\((1.8)\)进行积分计算。

该过程在实时运行前烘焙一下，生成立方体贴图以做后续计算备用即可。

FGD项

\[ FGD(\vec{v}) = \frac{1}{N} \sum_{i = 0}^{i \lt N} \frac { f_s \ast (\vec{n} \cdot \vec{l_i}) } { p_{s}(p, \vec{l_i}, \vec{v}) } \tag{1.11} \]

展开brdf函数\(f_s\)和pdf函数\(p_s\)得

\[ FGD(\vec{v}) = \frac{1}{N} \sum_{i = 0}^{i \lt N} \frac { F \ast G \ast (\vec{v} \cdot \vec{h_i}) } { (\vec{n} \cdot \vec{v}) \ast (\vec{n} \cdot \vec{h_i}) } \tag{1.12} \]

我们之前有介绍过，菲涅尔方程F为

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

使用这个继续展开\((1.12)\)，得

\[ \begin{array}{l} FGD(\vec{v}) &=& \frac{1}{N} \sum\limits_{i = 0}^{i \lt N} (F_0 + (1 - F_0) * (1 - (\vec{n} \cdot \vec{v}))^5) \ast \frac { G \ast (\vec{v} \cdot \vec{h_i}) } { (\vec{n} \cdot \vec{v}) \ast (\vec{n} \cdot \vec{h_i}) } \\ &=& \frac{1}{N} \sum\limits_{i = 0}^{i \lt N} ( F_0 * [1 - (1-(\vec{v}\cdot\vec{h_i}))^5] + [(1-(\vec{v}\cdot\vec{h_i}))^5] ) \ast \frac { G \ast (\vec{v} \cdot \vec{h_i}) } { (\vec{n} \cdot \vec{v}) \ast (\vec{n} \cdot \vec{h_i}) } \\ &=& F_0 \ast \underbrace { \{ \frac{1}{N} \sum\limits_{i = 0}^{i \lt N} ( [1 - (1-(\vec{v}\cdot\vec{h_i}))^5]) \ast \frac { G \ast (\vec{v} \cdot \vec{h_i}) } { (\vec{n} \cdot \vec{v}) \ast (\vec{n} \cdot \vec{h_i}) } \}}_{scale} + \underbrace { \{ \frac{1}{N} \sum\limits_{i = 0}^{i \lt N} (1-(\vec{v}\cdot\vec{h_i}))^5 \ast \frac { G \ast (\vec{v} \cdot \vec{h_i}) } { (\vec{n} \cdot \vec{v}) \ast (\vec{n} \cdot \vec{h_i}) } \}}_{bias} \end{array} \tag{1.13} \]

则又得到两个可以预计算的求和项，设左边为scale，右边为bias，即

\[ scale = \frac{1}{N} \sum_{i = 0}^{i \lt N} ( [1 - (1-(\vec{v}\cdot\vec{h_i}))^5]) \ast \frac { G \ast (\vec{v} \cdot \vec{h_i}) } { (\vec{n} \cdot \vec{v}) \ast (\vec{n} \cdot \vec{h_i}) } \tag{1.14} \]

\[ bias = \frac{1}{N} \sum_{i = 0}^{i \lt N} (1-(\vec{v}\cdot\vec{h_i}))^5 \ast \frac { G \ast (\vec{v} \cdot \vec{h_i}) } { (\vec{n} \cdot \vec{v}) \ast (\vec{n} \cdot \vec{h_i}) } \tag{1.15} \]

那么，\((1.13)\)可记为

\[ FGD(\vec{v}) = F_0 \ast scale + bias \tag{1.16} \]

显然，我们只需要在切线空间来计算这个结果即可，因为其值与光源无关，所以不管法线朝向哪里，其最终预处理出来的能量值都是一样的，即这个积分结果是不受法线影响的。所以我们假设\(\vec{n} = (0.0, 0.0, 1.0)\)，把法线固定成朝上的单位法线即可。

那么剩下影响这个积分值的变量，就只有粗糙度roughness(影响重要性采样以及G函数)和观察方向\(\vec{v}\)(影响\(1.13\)的计算和G函数)了。

对于观察方向\(\vec{v}\)，这里有一个技巧是，因为\(\vec{n}\)已经固定下来了，所以如果我们知道\((\vec{n} \cdot \vec{v})\)的话，那么\(v\)是可以反推出来的。所以我们不妨将\((\vec{n} \cdot \vec{v})\)作为输入，这样就只需要传递一个值，而非一个向量了。

所以，对于FGD项，我们传入两个值roughness和\((\vec{n} \cdot \vec{v})\)，就可以得到scale和bias了。

2个变量，得到2个值，似乎一张2D贴图就可以保存所有的数据了。我们使用一张2D贴图来保存这些信息。

x坐标表示\((\vec{n} \cdot \vec{v})\)。

y坐标表示roughness。

对于(x, y)处的颜色值(r, g, b)。

红色通道r保存scale。

绿色通道g保存bias。

代码如下。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
// pbr_ibl_fgd_bake.frag

vec2 ibl_dgf_integral(float n_o_v, float roughness) {
    float scale = 0.0;
    float bias  = 0.0;

    // 固定法线方向
    vec3 n = vec3(0.0, 0.0, 1.0);
    // 逆推观察向量
    vec3 v = vec3(
        sqrt(1.0 - n_o_v * n_o_v),
        0.0,
        n_o_v
    );

    for(uint i = 0u; i < sampler_num; ++i) {
        // hammersley 采样
        vec2 uv = hammersley(i, sampler_num);
        // 重要性采样
        vec3 h = importance_sampling_ggx(uv, roughness, n);
        // 从 v 和 h 逆算 l -- 光源方向
        vec3 l = normalize(2.0 * dot(v, h) * h - v);

        // 常用cos值
        float n_o_l = max(dot(n, l), 0.0);
        float n_o_h = max(dot(n, h), 0.0);
        float h_o_v = max(dot(h, v), 0.0);

        if(n_o_l > 0.0) {
            // 根据直接光源/IBL 选择相应k值
            float k = brdf_g_k_ibl(roughness);

            // 几何函数
            float g = brdf_g_smith(n_o_v, n_o_l, k);

            // 开始代入公式(1.13)计算
            float f = pow(1.0 - h_o_v, 5.0);
            float g_vis = g * h_o_v / (n_o_h * n_o_v);

            scale += (1.0 - f)  * g_vis;
            bias  += f          * g_vis;

        }

    }

    scale /= float(sampler_num);
    bias  /= float(sampler_num);

    return vec2(scale, bias);
}

最终得到的2D贴图如下所示。

FDG项积分贴图

从左往右为\((\vec{n} \cdot \vec{v})\)，从上往下为roughness。此图官方叫做lookup texture, LUT。

根据预处理结果计算

通过上一步近似过程，我们可以把反射率方程\((1.7)\)近似为以下形式。

\[ L_{specular}(p, \vec{v}) \approx LD \ast (F_0 \ast scale + bias) \tag{1.17} \]

所以我们在运行时计算specular IBL时，只需要将LD项和FGD项的贴图作为输入，查询对应的值，并使用\((1.17)\)进行计算即可。便无需实时进行积分计算了。代码如下。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
// pbr_ibl_epic.frag

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

// LD项
uniform samplerCube u_ibl_ld_map;
// FGD项
uniform sampler2D   u_ibl_dgf_map;

// ld项cubemap的最大LOD数 -- 通过roughness确定mipmap层级用
uniform int         u_max_reflection_lod;

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

vec3 pbr_ibl_specular_epic(vec3 n, vec3 v, vec3 albedo, vec3 f0, float roughness, float metallic) {
    vec3 res = vec3(0.0, 0.0, 0.0);

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

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

    // 反射方向
    vec3 R = reflect(-v, n);

    // 采样LD项
    vec3 ld = textureLod(u_ibl_ld_map, R, roughness * u_max_reflection_lod).rgb;
    //vec3 ld = texture(u_ibl_ld_map, R, 0).rgb;

    // 采样FGD项
    vec2 dgf = texture(u_ibl_dgf_map, vec2(n_o_v, roughness)).rg;

    // 代入公式(1.17)
    res = ld *(f * dgf.x + dgf.y);

    return res;
}

需要注意的是，代码中使用反射方向\(\vec{R}\)来对LD项进行采样，这是因为虽然LD项在生成是假设了\(\vec{n} = \vec{v} = \vec{R}\)，但是这里把反射方向当作是整体的光源入射方向来采样，即\(\vec{l_i} = \vec{R}\)，这大概是因为在反射方向上brdf的波瓣最大，且对不同的入射光线\(l_i\)来说变化不大(这里不是很懂Orz)。

IBL_epic效果图

最终运行得到的效果如图所示。从后向前，粗糙度roughness依次减少；从左向右，金属度metallic依次增加。

结语

以上及以前的几篇便是读者所有的对pbr粗浅的认识了，终于完结了Orz。能力有限，很多地方未能求得甚解，只能等以后有机会再深入学习和补充了Orz。

\(^o^)/

参考

Alternative Take on the Split Sum Approximation for Cubemap Pre-filtering

Real Shading in Unreal Engine 4