完整推导PBR IBL

Posted on 2025-06-06 Waline:

https://google.github.io/filament/Filament.html#endnote-ibltypes1

本文将完整推导 PBR IBL (image-based lighting)，将数学公式、近似和蒙特卡洛积分讲清楚。

PBR

我们使用的 PBR 在 diffuse BRDF 使用 Lambert 模型，在 specular BRDF 使用微表面模型（GGX）。即

\[ L_o( \boldsymbol \omega_r) = \int_{\Omega} L_i( \boldsymbol \omega_i) f_r( \boldsymbol \omega_r, \boldsymbol \omega_i) (\boldsymbol \omega_i \cdot \boldsymbol n)_{+} \, \mathrm d \boldsymbol \omega_i. \]

\[ f_r(\boldsymbol \omega_i, \boldsymbol \omega_r) = (1-m)(1 - F_{\text{avg}}) af_d(\boldsymbol \omega_i, \boldsymbol \omega_r) + f_s(\boldsymbol \omega_i, \boldsymbol \omega_r) \]

\[ f_d(\boldsymbol \omega_i, \boldsymbol \omega_r) = \frac{1}{\pi}, \]

\[ f_s(\boldsymbol \omega_i, \boldsymbol \omega_r) = \frac{D(\boldsymbol \omega_h, \boldsymbol n, \alpha) F(\boldsymbol \omega_r, \boldsymbol \omega_h) G(\boldsymbol \omega_i, \boldsymbol \omega_r, \alpha)}{4 (\boldsymbol n \cdot \boldsymbol \omega_i)_+ ( \boldsymbol n \cdot \boldsymbol \omega_r)_{+}}, \]

\[ F(\boldsymbol \omega_r, \boldsymbol \omega_h) = F_0 + (F_{90} - F_0) (1 - (\boldsymbol \omega_r \cdot \boldsymbol \omega_h)_{+})^5. \]

记 $k_d = (1-m)(1 - F_{\text{avg}})a, k_s = 1$，那么 $f_r = k_df_d + k_s f_s$。

我们需要计算的 radiance 结果为

\[ \begin{aligned} L_o( \boldsymbol \omega_r) &= \int_{\Omega} L_i( \boldsymbol \omega_i) f_r( \boldsymbol \omega_r, \boldsymbol \omega_i) (\boldsymbol \omega_i \cdot \boldsymbol n)_{+} \, \mathrm d \boldsymbol \omega_i. \\ &= k_d \int_{\Omega} L_i( \boldsymbol \omega_i) f_d(\boldsymbol \omega_i, \boldsymbol \omega_r) (\boldsymbol \omega_i \cdot \boldsymbol n)_{+} \, \mathrm d \boldsymbol \omega_i + k_s \int_{\Omega} L_i( \boldsymbol \omega_i) f_s(\boldsymbol \omega_i, \boldsymbol \omega_r) (\boldsymbol \omega_i \cdot \boldsymbol n)_{+} \, \mathrm d \boldsymbol \omega_i \\ &= k_d L_{\text{diffuse}} + k_s L_{\text{specular}}. \end{aligned} \]

Diffuse

Irradiance map

Diffuse 部分为

\[ \begin{aligned} L_{\text{diffuse}} &= \int_{\Omega} L_i(\boldsymbol \omega_i) \cdot \frac{1}{\pi} (\boldsymbol \omega_i \cdot \boldsymbol n)_+ \, \mathrm d \boldsymbol \omega_i \\ &= \frac{1}{\pi} \int_{\Omega} L_i(\boldsymbol \omega_i) (\boldsymbol \omega_i \cdot \boldsymbol n)_+ \, \mathrm d \boldsymbol \omega_i \\ &= \frac{1}{\pi} E(\boldsymbol n) \end{aligned} \]

其中 $E(\boldsymbol n)$ 为可预计算的 irradiance 项。我们可以将 irradiance 存储成 cubemap，里面包含每个法线方向的 irradiance 信息：

预计算时，我们对 cubemap 的每个法线方向卷积计算该积分
运行时，我们直接根据物体表面的法线方向，采样该 cubemap，得到具体的 irradiance 值。

其他存储预计算 irradiance 项 $E(\boldsymbol n)$ 的方式包括 sphere map 和球谐函数，这里不再赘述。

蒙特卡洛积分

考虑蒙特卡洛积分

\[ \begin{aligned} E(\boldsymbol n) &= \int_{\Omega} L_i(\boldsymbol \omega_i) (\boldsymbol \omega_i \cdot \boldsymbol n)_+ \, \mathrm d \boldsymbol \omega_i \\ &\approx \frac{1}{N} \sum_{i=1}^N \frac{L_i(\boldsymbol \omega_i) (\boldsymbol \omega_i \cdot \boldsymbol n)_{+}}{p(\boldsymbol \omega_i)}. \end{aligned} \]

朴素的均匀半球采样可设

\[ p(\boldsymbol \omega_i) = \frac{1}{ 2 \pi} \]

并且在代码中实现时，均匀地采样每个半球方向，但这样收敛较慢。

一个可能更好的的方法是，让 PDF 正比于 $(\boldsymbol \omega_i \cdot \boldsymbol n)_{+}$，并且根据该值偏向采样方向，即余弦加权的半球重要性采样。不妨设

\[ p(\boldsymbol \omega_i) = (\boldsymbol \omega_i \cdot \boldsymbol n)_{+}, \]

那么

\[ E(\boldsymbol n) \approx \frac{1}{N} \sum_{i=1}^N L_i(\boldsymbol \omega_i). \]

算法

对 irradiance cubemap 的每个像素（可并行）：
- 确定该像素对应的世界空间法线方向 $\boldsymbol n$
- 构造该位置的局部 TBN 空间
- $c=0$
- 在该局部空间生成 $N$ 个余弦加权采样的方向 $\boldsymbol \omega_1, \dots, \boldsymbol \omega_N$
- 对每个方向 $\boldsymbol \omega_i$：
  - 从环境贴图中采样 $L_i(\boldsymbol \omega_i)$，并累加到 $c$
- 得到该法线方向的预计算 irrdiance $L_i = c / N$。

不过，由于 irradiance 预计算本身的计算量就不大，且 diffuse 的计算结果本来也比较均匀，因此使用均匀采样并不会明显影响性能或结果。

Specular

Split sum

使用 split sum（这里的 split sum 是不约去 $(\boldsymbol \omega_i \cdot \boldsymbol n)_+$ 的版本，更准确）近似得到

\[ \begin{aligned} L_{\text{specular}} &= \int_{\Omega} L_i(\boldsymbol \omega_i) f_s(\boldsymbol \omega_i, \boldsymbol \omega_r) (\boldsymbol \omega_i \cdot \boldsymbol n)_+ \, \mathrm d \boldsymbol \omega_i \\ &= \left (\frac{\int_{\Omega} L_i(\boldsymbol \omega_i) f_s(\boldsymbol \omega_i, \boldsymbol \omega_r) (\boldsymbol \omega_i \cdot \boldsymbol n)_+ \, \mathrm d \boldsymbol \omega_i}{\int_{\Omega} f_s(\boldsymbol \omega_i, \boldsymbol \omega_r) (\boldsymbol \omega_i \cdot \boldsymbol n)_+ \, \mathrm d \boldsymbol \omega_i} \right) \left (\int_{\Omega} f_s(\boldsymbol \omega_i, \boldsymbol \omega_r) (\boldsymbol \omega_i \cdot \boldsymbol n)_+ \, \mathrm d \boldsymbol \omega_i \right) \\ & \approx \left (\frac{\int_{\Omega} L_i(\boldsymbol \omega_i) (\boldsymbol \omega_i \cdot \boldsymbol n)_+ \, \mathrm d \boldsymbol \omega_i}{\int_{\Omega} (\boldsymbol \omega_i \cdot \boldsymbol n)_+ \, \mathrm d \boldsymbol \omega_i} \right) \left (\int_{\Omega} f_s(\boldsymbol \omega_i, \boldsymbol \omega_r) (\boldsymbol \omega_i \cdot \boldsymbol n)_+ \, \mathrm d \boldsymbol \omega_i \right) \\ & = L_c \cdot k_{\text{BRDF}} \end{aligned} \]

其中

\[ L_c = \frac{\int_{\Omega} L_i(\boldsymbol \omega_i) (\boldsymbol \omega_i \cdot \boldsymbol n)_+ \, \mathrm d \boldsymbol \omega_i}{\int_{\Omega} (\boldsymbol \omega_i \cdot \boldsymbol n)_+ \, \mathrm d \boldsymbol \omega_i}. \\ \]

PFEM (Pre-filtered environment map)

注意到，在 split sum 近似后：

$L_c$ 中不再含有 $\boldsymbol \omega_r$。常用的工程近似手段是，令 $\boldsymbol \omega_r = \boldsymbol n = \boldsymbol R$，其中 $\boldsymbol R$ 是新定义的变量名，即假设我们的镜面反射角度总是为表面法线。这会引入误差，最明显的现象就是在掠射角时无法得到拖长的反射。
$L_c$ 中不再含有 $\alpha$。常用的工程近似手段是，预先定义一个离散的 $\alpha$ 等级列表，为每个粗糙度等级生成一个独立的 $L_c$ mipmap level，后续在高粗糙度下采样 mipmap level 高的 $L_c$，而在低粗糙度下采样 mipmap level 低的 $L_c$，以模拟不同粗糙度下粗糙/光滑镜面反射的情况。
- 在采样时，我们可以根据 PDF 的值采样不同 mipmap 的环境贴图。具体而言，当 PDF 较大时，此时表面较光滑，样本较集中，因此采样更清晰的环境贴图（mipmap level 低）；当 PDF 较小时，采样 mipmap level 高的环境贴图。PDF 的值会与 $\alpha$ 有关系，但环境贴图的 mipmap 与 $L_c$ 的 mipmap 没有直接联系。

计算 $L_c$ 项，使用离散和代替积分得到

\[ \begin{aligned} L_c &= \frac{\int_{\Omega} L_i(\boldsymbol \omega_i) (\boldsymbol \omega_i \cdot \boldsymbol R)_+ \, \mathrm d \boldsymbol \omega_i}{\int_{\Omega} (\boldsymbol \omega_i \cdot \boldsymbol R)_+ \, \mathrm d \boldsymbol \omega_i}. \\ &\approx \frac{\frac{1}{N} \sum_{i=1}^N \frac{L_i(\boldsymbol \omega_i) (\boldsymbol \omega_i \cdot \boldsymbol R)_{+}}{p(\boldsymbol \omega_i)}}{\frac{1}{N} \sum_{i=1}^N \frac{ (\boldsymbol \omega_i \cdot \boldsymbol R)_{+}}{p(\boldsymbol \omega_i)}} \\ &\approx \frac{\sum_{i=1}^N L_i(\boldsymbol \omega_i) (\boldsymbol \omega_i \cdot \boldsymbol R)_+ }{\sum_{i=1}^N (\boldsymbol \omega_i \cdot \boldsymbol R)_+} \\ \end{aligned} \]

其中 $\boldsymbol \omega_i$ 的采样不应使用朴素的均匀半球采样，因为这是 specular 项，均匀半球采样难以收敛。

具体地，应该使用 GGX 加权重要性采样（如果你使用的 specular BRDF 的 $D, G$ 不是 GGX，则需要使用对应的重要性采样方法），让 PDF 正比于法线分布函数 $D$。

算法

对每个离散的粗糙度 $\alpha$ 等级（可并行）：
- 对 pre-filtered environment map 的每个像素（可并行）：
  - 确定该像素对应的世界空间法线方向 $\boldsymbol n$
  - 令 $\boldsymbol \omega_r = \boldsymbol n = \boldsymbol R$
  - 构造该位置的局部 TBN 空间
  - $c=0, w = 0$
  - 在该局部空间生成 $N$ 个 GGX 加权采样的方向 $\boldsymbol \omega_1, \dots, \boldsymbol \omega_N$
    - 具体地，可以先生成 GGX 加权采样的微表面法线 $\boldsymbol \omega_h$，再使用反射公式 $\boldsymbol \omega_i = \text{reflect}(\boldsymbol \omega_r, \boldsymbol \omega_h) = 2(\boldsymbol \omega_r \cdot \boldsymbol \omega_h) \boldsymbol \omega_h - \boldsymbol \omega_r$ 生成采样方向
  - 对每个方向 $\boldsymbol \omega_i$：
    - 从当前的采样 PDF 确定要采样的环境贴图的 mipmap 等级
    - 从环境贴图中采样 $L_i(\boldsymbol \omega_i)$
    - 将 $L_i(\boldsymbol \omega_i)(\boldsymbol \omega_i \cdot \boldsymbol R)_+$ 累加到 $c$
    - 将 $(\boldsymbol \omega_i \cdot \boldsymbol R)_+$ 累加到 $w$
  - 得到该法线方向的预计算 $L_c(\boldsymbol R, \alpha) = c / w$。

分离菲涅耳

接下来我们计算 BRDF 项。记

\[ F_c = (1 - (\boldsymbol \omega_r \cdot \boldsymbol \omega_h)_{+})^5. \]

我们表演一个菲涅耳项的花活：

\[ \begin{aligned} k_\text{BRDF} &= \int_{\Omega} f_s(\boldsymbol \omega_i, \boldsymbol \omega_r) (\boldsymbol \omega_i \cdot \boldsymbol n)_+ \, \mathrm d \boldsymbol \omega_i \\ &= \int_{\Omega} f_s(\boldsymbol \omega_i, \boldsymbol \omega_r) \frac{F(\boldsymbol \omega_r, \boldsymbol \omega_h)}{F(\boldsymbol \omega_r, \boldsymbol \omega_h)} (\boldsymbol \omega_i \cdot \boldsymbol n)_+ \, \mathrm d \boldsymbol \omega_i \\ &= \int_{\Omega} \frac{f_s(\boldsymbol \omega_i, \boldsymbol \omega_r)}{F(\boldsymbol \omega_r, \boldsymbol \omega_h)} (F_0 + (F_{90} - F_0) F_c) (\boldsymbol \omega_i \cdot \boldsymbol n)_+ \, \mathrm d \boldsymbol \omega_i \\ &= \int_{\Omega} \frac{f_s(\boldsymbol \omega_i, \boldsymbol \omega_r)}{F(\boldsymbol \omega_r, \boldsymbol \omega_h)} (F_0(1 - F_c )+ F_{90} F_c) (\boldsymbol \omega_i \cdot \boldsymbol n)_+ \, \mathrm d \boldsymbol \omega_i \\ &= F_0 \int_{\Omega} \frac{f_s(\boldsymbol \omega_i, \boldsymbol \omega_r)}{F(\boldsymbol \omega_r, \boldsymbol \omega_h)} (1 - F_c) (\boldsymbol \omega_i \cdot \boldsymbol n)_+ \, \mathrm d \boldsymbol \omega_i + F_{90} \int_{\Omega} L_i(\boldsymbol \omega_i) \frac{f_s(\boldsymbol \omega_i, \boldsymbol \omega_r)}{F(\boldsymbol \omega_r, \boldsymbol \omega_h)} F_c (\boldsymbol \omega_i \cdot \boldsymbol n)_+ \, \mathrm d \boldsymbol \omega_i \\ &= F_0 \cdot \text{DFG}_1(( \boldsymbol n \cdot \boldsymbol \omega_r)_{+}, \alpha) + F_{90} \cdot \text{DFG}_2(( \boldsymbol n \cdot \boldsymbol \omega_r)_{+}, \alpha) \end{aligned} \]

其中

\[ \text{DFG}_1(( \boldsymbol n \cdot \boldsymbol \omega_r)_{+}, \alpha) = \int_{\Omega} \frac{DG}{4 ( \boldsymbol n \cdot \boldsymbol \omega_r)_{+}} (1 - F_c) \, \mathrm d \boldsymbol \omega_i, \]

\[ \text{DFG}_2(( \boldsymbol n \cdot \boldsymbol \omega_r)_{+}, \alpha) = \int_{\Omega} \frac{DG}{4 ( \boldsymbol n \cdot \boldsymbol \omega_r)_{+}} F_c \, \mathrm d \boldsymbol \omega_i, \]

是两个只与 $\boldsymbol n, \boldsymbol \omega_r$ 和 $\alpha$ 有关的函数。由于我们使用的 GGX 是各向同性的，因此 $\boldsymbol \omega_r$ 的取值在其与 $\boldsymbol n$ 的夹角一定时具有任意性，所以，我们只需用两个变量 $(\boldsymbol n \cdot \boldsymbol \omega_r)_{+}$ 和 $\alpha$ 作为 DFG 的输入即可。

BRDF LUT

我们将使用蒙特卡洛积分计算 $\text{DFG}_1$ 和 $\text{DFG}_2$。使用 GGX 加权重要性采样，使用的 PDF 为

\[ p(\boldsymbol \omega_i) = D \cdot \frac{(\boldsymbol n \cdot \boldsymbol \omega_h )_+}{4 (\boldsymbol \omega_r \cdot \boldsymbol \omega_h)_+}. \]

那么

\[ \begin{aligned} \text{DFG}_1 &= \int_{\Omega} \frac{DG}{4 ( \boldsymbol n \cdot \boldsymbol \omega_r)_{+}} (1 - F_c) \, \mathrm d \boldsymbol \omega_i, \\ & \approx \frac{1}{N} \sum_{i=1}^N \frac{\frac{DG}{4 ( \boldsymbol n \cdot \boldsymbol \omega_r)_{+}} (1 - F_c)}{p(\boldsymbol \omega_i)} \\ &= \frac{1}{N} \sum_{i=1}^N \frac{G \cdot (\boldsymbol \omega_r \cdot \boldsymbol \omega_h)_+ (1 - F_c)}{(\boldsymbol n \cdot \boldsymbol \omega_r)_+ (\boldsymbol n \cdot \boldsymbol \omega_h)_+}. \end{aligned} \]

同理

\[ \text{DFG}_2 = \frac{1}{N} \sum_{i=1}^N \frac{G \cdot (\boldsymbol \omega_r \cdot \boldsymbol \omega_h)_+ F_c}{(\boldsymbol n \cdot \boldsymbol \omega_r)_+ (\boldsymbol n \cdot \boldsymbol \omega_h)_+}. \]

算法

对 BRDF 2D LUT 的每个像素 $(u, v) \in [0,1]^2$（可并行）：
- $(\boldsymbol n \cdot \boldsymbol \omega_r)_{+} = u, \alpha = v$
- 在局部 TBN 空间中有 $\boldsymbol n = (0, 0, 1)$，故不失一般性，可取 $\boldsymbol \omega_r = (\sqrt{1 - u^2}, 0, u)$
- $l_1 = 0, l_2 = 0$
- 在该局部空间生成 $N$ 个 GGX 加权采样的方向 $\boldsymbol \omega_1, \dots, \boldsymbol \omega_N$
- 对每个方向 $\boldsymbol \omega_i$：
  - $F_c = (1 - (\boldsymbol \omega_r \cdot \boldsymbol \omega_h)_{+})^5.$
  - $G_{\text{vis}} = \frac{G \cdot (\boldsymbol \omega_r \cdot \boldsymbol \omega_h)_+}{(\boldsymbol n \cdot \boldsymbol \omega_r)_+ (\boldsymbol n \cdot \boldsymbol \omega_h)_+}$
  - 将 $G_{} (1-F_c) $ 累加到 $l_1$
  - 将 $G_{} F_c $ 累加到 $l_2$
- 得到该 LUT 在该 uv 处的预计算结果 $(l_1 / N, l_2 / N)$。

若要可视化该 LUT，可使用红绿颜色对应显示每个像素保存的两个分量，这样 LUT 就是一个双通道红绿图了。

总结

我们使用三个预计算数据结构

1 个 cubemap $E(\boldsymbol n)$，
1 个有 mipmap 的 cubemap $L_c(\boldsymbol R, \alpha)$，
1 个 2D LUT $\text{DFG}(( \boldsymbol n \cdot \boldsymbol \omega_r)_{+}, \alpha)$

实现了 IBL：

一般 $F_{90}$ 可用 $1$ 近似。

在 shading 采样 $L_c$ 时要注意一个细节，我们一般使用

\[ \boldsymbol R = \text{reflect}(\boldsymbol \omega_r, \boldsymbol n) \]

而非 $\boldsymbol n$ 采样 $L_c$。这是因为，我们只关心当前视角下 specular BRDF lobe 中能量集中的部分，这部分的中心是观察方向 $\boldsymbol \omega_r$ 的镜面反射方向。