[OpenGL]使用OpenGL实现基于物理的渲染模型PBR（下）

一、简介

在[OpenGL]使用OpenGL实现基于物理的渲染模型PBR（中）中介绍了基于物理的渲染(Physically Based Rendering, PBR) 的基本概念，并且实现了基于图像光源（IBL）的PBR中的漫反射部分。本文中将继续介绍IBL中的镜面反射部分。
按照本文代码实现后，可以实现以下效果：

二、基于IBL的PBR

1. 什么是IBL

IBL（Image-Based Lighting，基于图像的光照） 是一种使用环境贴图（Environment Map）提供间接光照的渲染技术，广泛应用于物理渲染（PBR）、电影特效 和 实时渲染（如游戏、VR等）。它通过 预计算的光照信息 让物体在复杂环境中表现出更加真实的全局光照（Global Illumination, GI） 效果。
简单来讲，IBL是一种存储环境光照信息的方法。在 IBL 中默认待渲染的模型接受来自环境中四面八方的光照，环境光使用一个 环境贴图 将各个方向的 irrdiance 存储到一张 texture 中（类似于 skybox），那么在渲染时就可以根据该 环境贴图 得到不同方向环境光对目标着色点的光照贡献。
根据渲染方程：
$$L_o\left(p,wo\right)={\int }_{{\mathrm{\Omega }}_{\mathrm{np}}}fr\left(p,wi,wo\right)\ast Li\left(p,wi\right)\ast n\cdot widwi$$
在渲染时，对于目标点 $p$、当前的视线向量$wo$，我们可以根据 $wi$ 采样环境贴图对应方向的值，而环境贴图 $wi$方向的采样值即为 $Li(p,wi)$。因此我们可以使用数值积分方法求解上述渲染方程，即可得到目标值 $L_o$。

2. IBL中的漫反射部分

如[OpenGL]使用OpenGL实现基于物理的渲染模型PBR（上）中我们讲的那样，在 pbr 中$fr(p,wi,wo)$可以分为 漫反射 和 镜面反射 两部分，如以下公式所示：
$$fr\left(p,wi,wo\right)=kd\ast f_{lambert}+fs\ast f_{Cook-Torrance}$$
因此渲染方程可以写为：
$$L_o\left(p,wo\right)={\int }_{\mathrm{\Omega }}(kd\ast f_{lambert}+fs\ast f_{Cook-Torrance})\ast Li\left(p,wi\right)\ast n\cdot widwi$$
其中
$$f_{lambert}=\frac{c}{\pi }$$
如果我们只考虑 漫反射 部分，那么可以得到：
$$\begin{gathered} L_{o,lambert}\left(p,wo\right)={\int }_{{\mathrm{\Omega }}_{\mathrm{np}}}(kd\ast f_{lambert})\ast Li\left(p,wi\right)\ast n\cdot widwi \\ ={\int }_{{\mathrm{\Omega }}_{\mathrm{np}}}(kd\ast \frac{c}{\pi })\ast Li\left(p,wi\right)\ast n\cdot widwi \\ =kd\ast \frac{c}{\pi }{\int }_{{\mathrm{\Omega }}_{\mathrm{np}}}Li\left(p,wi\right)\ast n\cdot widwi \end{gathered}$$
因此，假如只考虑漫反射部分，任意的目标着色点 $p$ 对应的积分值只与点 $p$ 法向 $n_p$ 所在的半球 ${\Omega }_{np}$ 有关，跟当前的视线向量 $wo$ 、点 $p$ 的材质属性都没关系。
那么，我们可以根据环境贴图预先计算得到法向为 $n_p$ 对应的半球 ${\Omega }_{np}$ 范围内的积分值：
$$\frac{1}{\pi }{\int }_{{\mathrm{\Omega }}_{\mathrm{np}}}Li\left(p,wi\right)\ast n\cdot widwi$$
即，预计算得到一个查找表（look up table），该查找表的仅仅包含一个维度的自变量，即目标着色点的法向$n_p$，因变量为目标积分值（离线预计算得到）$\frac{1}{\pi }{\int }_{{\mathrm{\Omega }}_{\mathrm{np}}}Li\left(p,wi\right)\ast n\cdot widwi$。
如下所示：
$$L_{o,lambert}(n_p)=lookuptable(n_p)=\frac{1}{\pi }{\int }_{{\mathrm{\Omega }}_{\mathrm{np}}}Li\left(wi\right)\ast n\cdot widwi$$
然后根据目标着色点的 漫反射系数$kd$ 和 颜色$c$，即可得到 $L_{o,lambert}=kd\ast c\ast L_{o,lambert}(n_p)$。

3. IBL中的镜面反射部分

对于渲染方程中的 镜面反射 部分：
$$\begin{gathered} L_{o,Cook-Torrance}\left(p,wo\right)={\int }_{{\mathrm{\Omega }}_{\mathrm{np}}}(fs\ast f_{Cook-Torrance})\ast Li\left(p,wi\right)\ast n\cdot widwi \\ =ks\ast {\int }_{{\mathrm{\Omega }}_{\mathrm{np}}}\frac{D\ast F\ast G}{4\ast \left(wo\cdot n\right)\left(wi\cdot n\right)}Li\left(p,wi\right)\ast n\cdot widwi \end{gathered}$$
可以看到，对于镜面部分任意点 $p$ 受到环境光的影响不仅仅跟 $p$ 对应的半球 ${\Omega }_{np}$ 相关，还受到视角方向 $wo$、$p$ 的材质属性的影响，如果如处理漫反射一样构建一个查找表，其输入维度不仅需要包含法向$n_p$，还需要包含 $wo$，着色点材质属性 $\alpha$（粗糙度）等，维度太大。因此不能简单地如处理漫反射一样构建一个look up table 表示任意点 $p$ 在任意视角方向 $wo$ 下的环境光镜面反射。因此我们需要对渲染方程进行简化近似，以达到可以实时计算，或者可以预计算一个look up table 的程度。

3.1 镜面反射渲染方程

渲染方程中镜面反射的部分公式如下：
$$\begin{gathered} L_{o,Cook-Torrance}\left(p,wo\right)={\int }_{{\mathrm{\Omega }}_{\mathrm{np}}}(fs\ast f_{Cook-Torrance})\ast Li\left(p,wi\right)\ast n\cdot widwi \\ =ks\ast {\int }_{{\mathrm{\Omega }}_{\mathrm{np}}}\frac{D\ast F\ast G}{4\ast \left(wo\cdot n\right)\left(wi\cdot n\right)}Li\left(p,wi\right)\ast n\cdot widwi \\ ={\int }_{{\mathrm{\Omega }}_{\mathrm{np}}}\frac{D\ast F\ast G}{4\ast \left(wo\cdot n\right)\left(wi\cdot n\right)}Li\left(p,wi\right)\ast n\cdot widwi\text{\hspace{1em}(1)} \end{gathered}$$
其中 $ks$ 即 $F$，因此在公式中省略了积分前的 $ks$。镜面反射项主要 分为 D，F，G三项。其中：

• D 为法向分布函数，Normal Distribution Function，微表面中法线朝向半程向量 ℎ 的表面片段所占的比例。使用以下公式计算：
  $$D=ND{F}_{GGX-TR}\left(n,h,\alpha \right)=\frac{{\alpha }^2}{\pi {\left({\left(n\cdot h\right)}^2\left({\alpha }^2-1\right)+1\right)}^2}$$
  其中 $\alpha$为着色点粗糙度，$n$为着色点法向，$h$为半程向量。
• F 为菲涅尔项，Fresnel Rquation，描述了光在表面上的反射率随视角变化的情况，即上面的$ks$。可以使用以下公式近似表示：
  $$F=F_{Schlick}\left(h,v,F_0\right)=F_0+\left(1-F_0\right){\left(1-\left(h\cdot v\right)\right)}^5$$
  其中 $F_0$ 为着色点的基础反射率，$h$为半程向量，$v$为视线向量。
• G 为几何函数，Geometry Function，计算微表面之间的遮挡效应。G 项比价复杂，计算公式如下：
  $$G=G_{Schlick-GXX}\left(n,v,l,k\right)=G_{sub}\left(n,v,k\right)\ast G_{sub}\left(n,l,k\right)$$
  其中：
  $$\begin{gathered} k=(\alpha +1{)}^2/8，针对直接光照 \\ k={\alpha }^2/2，针对IBL光照 \end{gathered}$$
  $$\begin{gathered} G_{sub}\left(n,v,k\right)=G_{Schlick-GGX}\left(n,v,k\right)=\frac{n\cdot v}{\left(n\ast v\right)\left(1-k\right)+k} \\ {G}_{sub}\left(n,l,k\right)=G_{Schlick-GGX}\left(n,l,k\right)=\frac{n\cdot l}{\left(n\ast l\right)\left(1-k\right)+k} \end{gathered}$$
  $\alpha$为着色点粗糙度，$n$为着色点法向，$l$为光线向量，$v$为视线向量。

3.2 渲染方程简化

首先明确我们的目标，即如何快速的计算公式：
$$L_o\left(n_p,wo,\alpha \right)={\int }_{{\mathrm{\Omega }}_{\mathrm{np}}}fr(wi,n,wo,\alpha )Li\left(p,wi\right)\ast n\cdot widwi\text{\hspace{1em}(2)}$$
其中 $n_p$ 为目标着色点的法向，$\alpha$为目标着色点的粗糙度，$wo$ 为视线向量。

我们将渲染场景分为 待渲染的模型 和 光照贴图 两部分。
使用公式 (3) 近似表示 目标公式(2)：
$$L_o\left(n_p,wo,\alpha \right)=pre\_filtered(n_p,wo,\alpha )\ast BRDF(n_p,wo,\alpha )\text{\hspace{1em}(3)}$$
其中 pre_filtered项（又称LD项） 是一个跟光照贴图相关的函数， BRDF 项（又称 DFG项）为只与 待渲染场景相关的函数。
接下来我们将会详细介绍如何定义、求解 BRDF 项 和 pre_filtered 项。

3.3 BRDF 项

3.3.1 BRDF 项的定义

我们假设渲染场景中存在一个待渲染的模型，并且光照贴图各个位置的值都为 1，即 $Li(wi)\equiv 1$。
那么此时的渲染结果为：
$$\begin{gathered} \overline{L_o\left(n_p,wo,\alpha \right)}={\int }_{{\mathrm{\Omega }}_{\mathrm{np}}}fr(wi,n,wo,\alpha )\ast Li\left(p,wi\right)\ast n\cdot widwi \\ ={\int }_{{\mathrm{\Omega }}_{\mathrm{np}}}fr(wi,n,wo,\alpha )\ast 1\ast n\cdot widwi\text{\hspace{1em}(4)} \end{gathered}$$
可以看出 公式(4) 是一个只与待渲染模型自身属性相关的函数，与外界的光照贴图无关。因此可以将公式(4) 视作 BRDF 项。那么就有：
$$BRDF(n_p,wo,\alpha )={\int }_{{\mathrm{\Omega }}_{\mathrm{np}}}fr(wi,n,wo,\alpha )\ast n\cdot widwi\text{\hspace{1em}(5)}$$

3.3.2 BRDF 项的数值求解

根据 基于重要性采样的蒙特卡洛数值积分 可以得到：
$$\begin{gathered} BRDF(n_p,wo,\alpha )={\int }_{{\mathrm{\Omega }}_{\mathrm{np}}}fr(wi,n,wo,\alpha )\ast n\cdot widwi \\ =\frac{1}{N}\sum _{j=1}^{j=N}\frac{fr(w{i}_j,n,wo,\alpha )\ast (n\cdot w{i}_j)}{pdf(w{i}_j)}\text{\hspace{1em}(6)} \end{gathered}$$
其中 $w{i}_j$ 为我们根据概率密度分布函数 $pdf(w{i}_j)$ 采样得到的 $wi$ 向量。

我们将 $fr(wi,n,wo,\alpha )$ 函数展开可以得到：
$$\begin{gathered} BRDF(n_p,wo,\alpha )=\frac{1}{N}\sum _{j=1}^{j=N}\frac{fr(w{i}_j,n,wo,\alpha )\ast (n\cdot w{i}_j)}{pdf(w{i}_j)} \\ =\frac{1}{N}\sum _{j=1}^{j=N}\frac{D(n_p,w{h}_j,\alpha )\ast F(w{h}_j,wo,F_0)\ast G(n_p,wo,w{i}_j,\alpha )\ast (n\cdot w{i}_j)}{4\ast (wo\cdot n_p)(w{i}_j\cdot n)\ast pdf(w{i}_j)}\text{\hspace{1em}(7)} \end{gathered}$$
其中$w{h}_j$为半程向量，可以根据$w{i}_j$和$wo$计算得到：
$$w{h}_j=\frac{(w{i}_j+wo)}{|w{i}_j+wo|}\text{\hspace{1em}(8)}$$
由于 公式(7) 中的 微表面法向分布函数 $D(n_p,w{h}_j,\alpha )$ 是以$w{h}_j$为自变量，因此我们需要根据$wh$的概率密度分布函数$pdf(wh)$采样$w{h}_j$，而不是根据$pdf(wi)$。
根据公式
$$\begin{gathered} pdf(wh)=D(n_p,wh,\alpha )\ast (n\ast wh) \\ pdf(wi)=pdf(wh)\ast |\frac{dwh}{dwi}|=pdf(wh)\ast \frac{1}{4(wh\cdot wo)}=D(n_p,wh,\alpha )\ast (n\ast wh)\ast \frac{1}{4(wh\cdot wo)} \end{gathered}$$
可以得到：
$$pdf(wi)=D(n_p,wh,\alpha )\ast (n\cdot wh)\ast \frac{1}{4(wh\cdot wo)}\text{\hspace{1em}(9)}$$
将公式(9)代入公式(7) 可以得到：
$$\begin{gathered} BRDF(n_p,wo,\alpha )=\frac{1}{N}\sum _{j=1}^{j=N}\frac{D(n_p,w{h}_j,\alpha )\ast F(w{h}_j,wo,F_0)\ast G(n_p,wo,w{i}_j,\alpha )\ast (n\cdot w{i}_j)}{4\ast (wo\cdot n_p)(w{i}_j\cdot n)\ast pdf(w{i}_j)} \\ =\frac{1}{N}\sum _{j=1}^{j=N}\frac{\sout{D(n_p,w{h}_j,\alpha )}\ast F(w{h}_j,wo,F_0)\ast G(n_p,wo,w{i}_j,\alpha )\ast \sout{(n_p\cdot w{i}_j)}}{\bcancel{4}\ast (wo\cdot n_p)\sout{(w{i}_j\cdot n)\ast D(n_p,w{h}_j,\alpha )}\ast (n_p\cdot w{h}_j)\ast \frac{1}{\bcancel{4}\ast (w{h}_j\cdot wo)}} \\ =\frac{1}{N}\sum _{j=1}^{j=N}\frac{F(w{h}_j,wo,F_0)\ast G(n_p,wo,w{i}_j,\alpha )\ast (w{h}_j\cdot wo)}{(wo\cdot n_p)\ast (n_p\cdot w{h}_j)}\text{\hspace{1em}(10)} \end{gathered}$$

3.3.3 F(wh,wo,F0) 函数的简化

因为
$$F(wh,wo,F_0)=F_0+\left(1-F_0\right){\left(1-\left(wh\ast wo\right)\right)}^5$$
整理得到：
$$\begin{gathered} F(wh,wo,F_0)=F_0+\left(1-F_0\right){\left(1-\left(wh\ast wo\right)\right)}^5 \\ =F_0\ast (1-(1-(wh\cdot wo){)}^5)+(1-(wh\cdot wo){)}^5\text{\hspace{1em}(11)} \end{gathered}$$
我们仔细观察下公式(11)，其自变量为 $wh,wo$ 和 $F_0$，其中 $F_0$ 为一个系数，我们可以暂时不管它。主要关注于自变量 $wh$ 和 $wo$ 。

从公式(11) 可以看出，函数 $F()$ 的自变量有三个分别是 $wh,wo$ 和 $F_0$。但是根据公式(10)，我们仅仅需要根据 $pdf(wh)$ 采样得到 $N$ 个 $w{h}_j$，然后计算 $N$ 个 $w{h}_j$ 对应的 $\frac{1}{N}{\sum }_{j=1}^{j=N}\cdots \ast (1-(w{h}_j\ast wo){)}^5\ast \dots$ ，而在采样 $wh$ 时，$wh$ 在径向满足概率分布函数 $pdf(wh)=D(wh)\ast (wh\cdot n)$，但是在周向是对称的，如下图所示：

那么说明，只要保证 $n\cdot wo$（即下图中的$\theta$） 为定值，不管我们选择下图中$wo1$还是$wo2$进行计算，采样$N$个$w{h}_j$得到的期望 $\frac{1}{N}{\sum }_{j=1}^{j=N}\cdots \ast (1-(w{h}_j\cdot wo){)}^5\ast \dots$ 就相等。

这说明，对于函数$F(wh,wo,\alpha )$，我们无需知道 $wo$ 的真实值是多少，只需要知道$n\cdot wo$，然后根据$n$（在采样$wh$时$n$为$(0,1,0)$或者$(0,0,1)$等单位向量）任取一个 $w{o}^{\prime }$，使得 $n\cdot w{o}^{\prime }=n\cdot wo$，那么根据$w{o}^{\prime }$计算得到的 $BRDF(n_p,w{o}^{\prime },\alpha )$与目标值 $BRDF(n_p,w{o}^{\prime },\alpha )$ 就相等。

因此，我们可以将 函数$F(wh,wo,\alpha )$，写成
$$F(wh,wo,F0)=F(wh,w{o}^{\prime },F0),其中w{o}^{\prime }满足n\cdot w{o}^{\prime }=n\cdot wo\text{\hspace{1em}(12)}$$
那么根据公式(12)，函数 $F(wh,wo,F0)$ 的自变量就变为 $wh$，$n\cdot wo$（因为可以根据$n\cdot wo$得到$w{o}^{\prime }$） 和 $F_0$ 了。

我们令：
$$F_c(wh,wo)=(1-(wh\cdot wo){)}^5$$
那么：
$$F(wh,n\cdot wo,F_0)=F(wh,n\cdot wo,F_0)=F_0\ast (1.0-Fc(wh,n\cdot wo))+F_c(wh,n\cdot wo)$$

3.3.4 BRDF 项的简化

将 $F(wh,n\cdot wo,F_0)$ 带入公式(10)中得到：
$$\begin{gathered} BRDF(n_p,wo,\alpha )=\frac{1}{N}\sum _{j=1}^{j=N}\frac{\left(F_0\ast (1.0-Fc(wh,n\cdot wo))+F_c(w{h}_j,n\cdot wo)\right)\ast G(n_p,wo,w{i}_j,\alpha )\ast (w{h}_j\cdot wo)}{(wo\cdot n_p)\ast (n_p\cdot w{h}_j)} \\ =F_0\ast \frac{1}{N}\ast \sum _{j=1}^{j=N}\frac{(1.0-Fc(wh,n\cdot wo))\ast G(n_p,wo,w{i}_j,\alpha )\ast (w{h}_j\cdot wo)}{(wo\cdot n_p)\ast (n_p\cdot w{h}_j)}+\frac{1}{N}\ast \sum _{j=1}^{j=N}\frac{F_c(w{h}_j,n\cdot wo)\ast G(n_p,wo,w{i}_j,\alpha )\ast (w{h}_j\cdot wo)}{(wo\cdot n_p)\ast (n_p\cdot w{h}_j)}\text{\hspace{1em}(13)} \end{gathered}$$
同样的根据前面的分析，由于我们是根据$pdf(wh)$采样得到大量的 $w{h}_j$，来计算上式，而且 $pdf(wh)$ 函数在周向上是对称的，因此我们只要根据 $n\cdot wo$ 确定任意一个 $w{o}^{\prime }$，使得 $n\cdot w{o}^{\prime }=n\cdot wo$ 即可，无需使用 $wo$ 计算上式。同时由于函数$G_{sub}\left(n,wo,k\right)=\frac{n\ast wo}{\left(n\ast wo\right)\left(1-k\right)+k}$
因此函数$G(n,wo,k)=G(n\cdot wo,\alpha )$，
那么公式(13) 可以写成：
$$BRDF(n_p,wo,\alpha )=BRDF(n_p\cdot wo,\alpha )=F_0\ast \underset{scale}{\underbrace{\frac{1}{N}\ast \sum _{j=1}^{j=N}\frac{F_{scale}(w{h}_j,n\cdot wo)\ast G(n\cdot wo,w{i}_j,\alpha )\ast (w{h}_j\cdot w{o}^{\prime })}{(wo\cdot n_p)\ast (n_p\cdot w{h}_j)}}}+\underset{bias}{\underbrace{\frac{1}{N}\ast \sum _{j=1}^{j=N}\frac{F_{bias}(w{h}_j,n\cdot wo)\ast G(n_p\cdot wo,w{i}_j,\alpha )\ast (w{h}_j\cdot w{o}^{\prime })}{(wo\cdot n_p)\ast (n_p\cdot w{h}_j)}}},其中w{o}^{\prime }可以根据n_p\cdot wo确定\text{\hspace{1em}(14)}$$
从公式(13) 可以看出，目标函数 BRDF 项是可以看作：
$$BRDF(n_p\cdot wo,\alpha )=F_0\ast scale(n_p\cdot wo,\alpha )+bais(n_p\cdot wo,\alpha )\text{\hspace{1em}(15)}$$
我们只要预计算得到 $scale(n_p\cdot wo,\alpha )$ 和 $bais(n_p\cdot wo,\alpha )$ 两个二维的查找表，然后根据 目标着色点的 $F_0$ 即可得到 BRDF项。
其中：
$$\begin{gathered} scale(n_p\cdot wo,\alpha )=\frac{1}{N}\ast \sum _{j=1}^{j=N}\frac{(1.0-Fc(wh,n\cdot wo))\ast G(n\cdot wo,w{i}_j,\alpha )\ast (w{h}_j\cdot w{o}^{\prime })}{(wo\cdot n_p)\ast (n_p\cdot w{h}_j)} \\ bais(n_p\cdot wo,\alpha )=\frac{1}{N}\ast \sum _{j=1}^{j=N}\frac{F_c(w{h}_j,n\cdot wo)\ast G(n_p\cdot wo,w{i}_j,\alpha )\ast (w{h}_j\cdot w{o}^{\prime })}{(wo\cdot n_p)\ast (n_p\cdot w{h}_j)}\text{\hspace{1em}(16)} \end{gathered}$$

3.3.5 预计算 BRDF 项的代码实现

根据公式(16)，计算 BRDF 项的代码如下，在代码中 $Fc(wh,n\cdot wo)\to Fc$，$n_p\to N$，$w{o}^{\prime }\to V$，$wh\to H$，$wi\to L$：

// ----------------------------------------------------------------------------
/**
input: n dot wo (NdotV), alpha (roughness)
output: scale(n dot wo,alpha), bias(n dot wo, roughness)
*/
vec2 IntegrateBRDF(float NdotV, float roughness)
{// 假设 N = (0,0,1)// 根据 N dot V，选择 V' = (sqrt(1.0 - NdotV*NdotV), 0, NdotV),// 可以保证 N dot V = N dot V'vec3 N = vec3(0.0, 0.0, 1.0);vec3 V;V.x = sqrt(1.0 - NdotV*NdotV);V.y = 0.0;V.z = NdotV;float scale = 0.0;float bias = 0.0; const uint SAMPLE_COUNT = 1024u;for(uint i = 0u; i < SAMPLE_COUNT; ++i){// generates a sample vector that's biased towards the// preferred alignment direction (importance sampling).vec2 Xi = Hammersley(i, SAMPLE_COUNT); // 生成一个 [0.0,1.0] 范围内的随机数vec3 H = ImportanceSampleGGX(Xi, N, roughness); // 根据重要性采样 pdf(wh) = D(wh)*(wh dot n) 采样 whvec3 L = normalize(2.0 * dot(V, H) * H - V); // 采样得到的 wh 对应的 wifloat NdotL = max(L.z, 0.0); // n dot wifloat NdotH = max(H.z, 0.0); // n dot whfloat VdotH = max(dot(V, H), 0.0); // wo dot whif(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);scale += (1.0 - Fc) * G_Vis;bias += Fc * G_Vis;}}scale /= float(SAMPLE_COUNT);bias /= float(SAMPLE_COUNT);return vec2(scale, bias);
}

最终我们可以预计算不同$(n_p\cdot wo,\alpha )$输入时对应的 $scale(n_p\cdot wo,\alpha )$ 和 $bais(n_p\cdot wo,\alpha )$，得到这两个 二维查找表后 根据公式(15) 得到 BRDF 项对应的值。
在代码实现时，往往将 $scale$ 和 $bias$ 两个分量作为 $(r,g)$ 分量，存储到一个 二维纹理中，得到的最终纹理结果如下所示：

3.4 pre_filtered项

3.4.1 pre_fileterd 项的定义

接下来介绍如何快速计算 公式 (3) 中 pre_filtered 项。根据公式(3)和公式(5)我们可以得到：

$$pre\_filtered(n_p,wo,\alpha )=\frac{L_o\left(n_p,wo,\alpha \right)}{BRDF(n_p,wo,\alpha )}=\frac{{\int }_{{\mathrm{\Omega }}_{\mathrm{np}}}fr(wi,n,wo,\alpha )Li\left(p,wi\right)\ast n\cdot widwi}{{\int }_{{\mathrm{\Omega }}_{\mathrm{np}}}fr(wi,n,wo,\alpha )\ast n\cdot widwi}\text{\hspace{1em}(17)}$$

3.4.2 pre_fileterd 项的数值求解

根据 基于重要性采样的蒙特卡洛数值积分 可以得到：
$$\begin{gathered} pre\_filtered(n_p,wo,\alpha )=\frac{{\int }_{{\mathrm{\Omega }}_{\mathrm{np}}}fr(wi,n,wo,\alpha )\ast Li\left(p,wi\right)\ast n\cdot widwi}{{\int }_{{\mathrm{\Omega }}_{\mathrm{np}}}fr(wi,n,wo,\alpha )\ast n\cdot widwi} \\ =\left(\frac{1}{N}\sum _{j=1}^{j=N}\frac{fr(w{i}_j,n,wo,\alpha )\ast Li\left(p,wi\right)\ast (n\cdot w{i}_j)}{pdf(w{i}_j)}\right)/\left(\frac{1}{N}\sum _{j=1}^{j=N}\frac{fr(w{i}_j,n,wo,\alpha )\ast (n\cdot w{i}_j)}{pdf(w{i}_j)}\right) \\ =\left(\sum _{j=1}^{j=N}\frac{fr(w{i}_j,n,wo,\alpha )\ast Li\left(p,wi\right)\ast (n\cdot w{i}_j)}{pdf(w{i}_j)}\right)/\left(\sum _{j=1}^{j=N}\frac{fr(w{i}_j,n,wo,\alpha )\ast (n\cdot w{i}_j)}{pdf(w{i}_j)}\right)\text{\hspace{1em}(17)} \end{gathered}$$
其中 $w{i}_j$ 为我们根据概率密度分布函数 $pdf(w{i}_j)$ 采样得到的 $wi$ 向量。与 求解 BRDF 项时类似，我们将
$$pdf(wi)=pdf(wh)\ast \frac{1}{4(wh\cdot wo)}=D(n_p,wh,\alpha )\ast (n\ast wh)\ast \frac{1}{4(wh\cdot wo)}$$
代入公式(17)可以得到：
$$\begin{gathered} pre\_filtered(n_p,wo,\alpha )=\left(\sum _{j=1}^{j=N}\frac{fr(w{i}_j,n,wo,\alpha )\ast Li\left(wi\right)\ast (n\cdot w{i}_j)}{pdf(w{i}_j)}\right)/\left(\sum _{j=1}^{j=N}\frac{fr(w{i}_j,n,wo,\alpha )\ast (n\cdot w{i}_j)}{pdf(w{i}_j)}\right) \\ =\left(\sum _{j=1}^{j=N}\frac{\sout{D(n_p,w{h}_j,\alpha )}\ast F(w{h}_j,wo,F_0)\ast G(n_p,wo,w{i}_j,\alpha )\ast Li\left(wi\right)\ast \sout{(n\cdot w{i}_j)}}{\bcancel{4}\ast (wo\cdot n_p)\sout{(w{i}_j\cdot n)\ast D(n_p,w{h}_j,\alpha )}\ast (n_p\cdot w{h}_j)\ast \frac{1}{\bcancel{4}(w{h}_j\cdot wo)}}\right)/\left(\sum _{j=1}^{j=N}\frac{\sout{D(n_p,w{h}_j,\alpha )}\ast F(w{h}_j,wo,F_0)\ast G(n_p,wo,w{i}_j,\alpha )\ast \sout{(n\cdot w{i}_j)}}{\bcancel{4}\ast (wo\cdot n_p)\sout{(w{i}_j\cdot n)\ast D(n_p,w{h}_j,\alpha )}\ast (n_p\cdot w{h}_j)\ast \frac{1}{\bcancel{4}(w{h}_j\cdot wo)}}\right) \\ =\left(\sum _{j=1}^{j=N}\frac{F(w{h}_j,wo,F_0)\ast G(n_p,wo,w{i}_j,\alpha )\ast Li\left(wi\right)\ast (w{h}_j\cdot wo)}{(wo\cdot n_p)\ast (n_p\cdot w{h}_j)}\right)/\left(\sum _{j=1}^{j=N}\frac{F(w{h}_j,wo,F_0)\ast G(n_p,wo,w{i}_j,\alpha )\ast (w{h}_j\cdot wo)}{(wo\cdot n_p)\ast (n_p\cdot w{h}_j)}\right)\text{\hspace{1em}(18)} \end{gathered}$$
公式(18)依旧是一个比较复杂的公式，并且其输入自变量包含$n_p,wo,\alpha$ 和 $F0$，依旧不能满足建立一个预计算表达到快速计算的目的。

3.4.3 pre_filtered 项的简化

接下来将介绍一系列假设，以简化公式(18)。
首先，公式(18)中，除了$Li(wi)$中的$wi$根据实际的$n_p$和采样得到的$wh$计算得到。假设其他部分的 $n_p$ 与 $wo$ 相同，注意，这里是假设$n_p$等于wo（改变$n_p$），而不是假设wo等于n_p。
因为我们使用真实的$n_p$和$wh$计算$wi$，那么在pre_filetred 公式中 $n_p$ 的变化只影响函数 $fr(wi,wo,n_p,\alpha )$。假设$n_p=wo$，即假设不同 $n$ 时，着色点表面的 lobe形状不会发生变化（只是发生旋转）。（然而实际上 $n_p$ 的变化会影响 lobe 的形状）。如下图所示：

其中 $\overline{lobe}$ 可以理解为使用 $Li(wi)=1$ 时计算得到的绝对值，$lobe$ 为使用 IBL $Li(wi)$ 计算得到的绝对值。
在假设 $wo$ 与 $n_p$ 相同的条件下，公式(18) 可以写成：
$$\begin{gathered} pre\_filtered(n_p,wo,\alpha )==\left(\sum _{j=1}^{j=N}\frac{F(w{h}_j,wo,F_0)\ast G(wo,wo,w{i}_j,\alpha )\ast Li\left(w{i}_j\right)\ast (w{h}_j\cdot wo)}{(wo\cdot wo)\ast (wo\ast w{h}_j)}\right)/\left(\sum _{j=1}^{j=N}\frac{F(w{h}_j,wo,F_0)\ast G(wo,wo,w{i}_j,\alpha )\ast (w{h}_j\cdot wo)}{(wo\cdot wo)\ast (wo\ast w{h}_j)}\right) \\ =\frac{{\sum }_{j=1}^{j=N}F(w{h}_j,wo,F_0)\ast G(wo,wo,w{i}_j,\alpha )\ast Li\left(w{i}_j\right)}{{\sum }_{j=1}^{j=N}F(w{h}_j,wo,F_0)\ast G(wo,wo,w{i}_j,\alpha )}\text{\hspace{1em}(19)} \end{gathered}$$
然后，由于$F()$ 函数对结果的影响不大（此处只给出结论，不给出推导证明，读者可以参考深入理解 PBR/基于图像照明 (IBL)-3.2 speular。
可以将$F()$函数省略，那么公式(19)简化为：
$$pre\_filtered(n_p,wo,\alpha )=\frac{{\sum }_{j=1}^{j=N}G(wo,wo,w{i}_j,\alpha )\ast Li\left(w{i}_j\right)}{{\sum }_{j=1}^{j=N}G(wo,wo,w{i}_j,\alpha )}\text{\hspace{1em}(19)}$$
最后，使用 $(n_p\cdot wi)$ 近似表示函数 $G(wo,wo,w{i}_j,\alpha )$，就得到的 pre_filtered 项的最终计算公式：
$$pre\_filtered(n_p,wo,\alpha )=\frac{{\sum }_{j=1}^{j=N}(n_p\cdot w{i}_j)\ast Li\left(w{i}_j\right)}{{\sum }_{j=1}^{j=N}(n_p\cdot w{i}_j)}\text{\hspace{1em}(20)}$$
需要注意的是，上式中的 $wi$ 是我们根据$n_p$、$wo$和使用重要性采样中的 $pdf(wh)$ 采样得到的 $wh$，计算得到的 $wi$，也就是说 { n p , w o , w h j } → { w i j } \{n_p, wo, wh_j\} \rightarrow{\{wi_j\}} {np​,wo,whj​}→{wij​}。而$wo$和$w{h}_j$是相对于$n_p=(0,1,0)$（或者$n_p=(0,0,1)$）的局部坐标系的。
从下图中可以看出，只要$wr=reflect(n_p,wo)$ 相同，那么在全局坐标系中根据$pdf(wi)$采样得到的 $wi$ 就相同。

上图中，黑色虚线的向量为采样得到的 $w{h}_j$，绿色虚线的向量为对应的 $w{i}_j$。

因此公式(20)可以写为：
$$pre\_filtered(wr,\alpha )=\frac{{\sum }_{j=1}^{j=N}(n_p\cdot w{i}_j)\ast Li\left(w{i}_j\right)}{{\sum }_{j=1}^{j=N}(n_p\cdot w{i}_j)}\text{\hspace{1em}(21)}$$
其中 $w{i}_j$ 为将 $wr$ 视作 $n_p$，基于$pdf(wh)$采样得到$w{h}_j$，然后计算得到的$w{i}_j$。

我们可以预计算不同$(wr,\alpha )$输入时对应的 $pre\_filtered(wr,\alpha )$，得到这二维查找表后 根据公式(21) 得到 pre_filtered 项对应的值。

3.4.4 预计算 pre_filtered 项的代码实现

根据公式(21)，计算 pre_filtered 项的代码如下，在代码中 $n_p\to N$，$wr\to R$，$wo\to V$，$wh\to H$，$wi\to L$：

/**
input: wr (WorldPos), alpha (roughness)
output: prefiltered(wr, roughness)
*/
void main()
{vec3 N = normalize(WorldPos); // 预计算时, wr = n// make the simplifying assumption that V equals R equals the normal vec3 R = N; // wr = nvec3 V = R; // wo = n = wrconst uint SAMPLE_COUNT = 1024u;vec3 prefilteredColor = vec3(0.0);float totalWeight = 0.0;for(uint i = 0u; i < SAMPLE_COUNT; ++i){// generates a sample vector that's biased towards the preferred alignment direction (importance sampling).vec2 Xi = Hammersley(i, SAMPLE_COUNT); // Xi 是一个 [0.0,1.0] 范围内部的随机数vec3 H = ImportanceSampleGGX(Xi, N, roughness); // H 是根据 pdf(wh) 采样得到的 wh_jvec3 L  = normalize(2.0 * dot(V, H) * H - V); // 采样得到的 wh_j 对应的 wi_jfloat NdotL = max(dot(N, L), 0.0); // n_p * wi_jif(NdotL > 0.0){// sample from the environment's mip level based on roughness/pdf// 根据 roughness 从 environmentMap 不同的 mipmap 中采样，仅仅用于加速计算// 即使 environmentMap 不使用 mipmap, 如果采样数目足够大，也可以保证最终的 prefilteredColor 值收敛float D   = DistributionGGX(N, H, roughness);float NdotH = max(dot(N, H), 0.0);float HdotV = max(dot(H, V), 0.0);float pdf = D * NdotH / (4.0 * HdotV) + 0.0001;  // pdf(wi_j)float resolution = 512.0; // resolution of source cubemap (per face)float saTexel  = 4.0 * PI / (6.0 * resolution * resolution);float saSample = 1.0 / (float(SAMPLE_COUNT) * pdf + 0.0001);float mipLevel = roughness == 0.0 ? 0.0 : 0.5 * log2(saSample / saTexel); prefilteredColor += textureLod(environmentMap, L, mipLevel).rgb * NdotL;totalWeight      += NdotL;}}prefilteredColor = prefilteredColor / totalWeight;FragColor = vec4(prefilteredColor, 1.0);
}

3.5 总结

根据上面的分析，我们可以预计算得到两个查找表 $pre\_filtered(wr,\alpha )$, $BRDF(n\cdot wo,\alpha )$。前者只与 渲染场景中的光照贴图有关，后者只与待渲染的模型有关。在得到这两个查找表后，在渲染时对于目标点$n_p$的最终着色，可以根据公式(3)计算得到：
$$L_o\left(n_p,wo,\alpha \right)=pre\_filtered(wo,\alpha )\ast BRDF(n_p\cdot wo,\alpha )\text{\hspace{1em}(3)}$$

三、全部代码及模型文件

使用OpenGL实现IBL漫反射+镜面反射 PBR的全部代码以及模型文件可以在OpenGL使用OpenGL实现基于物理的渲染模型PBR（下） 中下载。
下载源代码后使用以下命令编译运行：

mkdir build
cd build
cmake ..
make
./OpenGL_PBR

渲染结果如下：

四、参考

[1].LearnOpenGL-PBR-IBL-漫反射辐照
[2].深入理解 PBR/基于图像照明 (IBL)
[3].漫谈split sum approximate
[4].Alternative Take on the Split Sum Approximation for Cubemap Pre-filtering