从零开始的OpenGL光栅化渲染器构建6-PBR光照模型

前言

PBR，或者基于物理的渲染(Physically Based Rendering)，它指的是一些在不同程度上都基于与现实世界的物理原理更相符的基本理论所构成的渲染技术的集合。正因为基于物理的渲染目的便是为了使用一种更符合物理学规律的方式来模拟光线，因此这种渲染方式与我们原来的Phong或者Blinn-Phong光照算法相比总体上看起来要更真实一些。除了看起来更好些以外，由于它与物理性质非常接近，因此我们（尤其是美术师们）可以直接以物理参数为依据来编写表面材质，而不必依靠粗劣的修改与调整来让光照效果看上去正常。使用基于物理参数的方法来飙血材质还有一个更大的好处，就是不论光照条件如何，这些材质看上去都会是正确的，而在非PBR的渲染管线当中有些东西就不会那么真实了。

虽然如此，基于物理的渲染依然只是对基于物理原理的现实世界的一种近似，这也就是为什么它被称为基于物理的着色(Physically based Shading)而非物理着色(Physical Shading)的原因。判断一种PBR光照模型是否是基于物理的，必须满足以下三个条件：

基于微表面(Microfacet)的表面模型。
能量守恒。
应用基于物理的BRDF。

理论部分

微平面模型

所有的PBR技术都基于微平面理论。这项理论认为，打到微观尺度之后任何平面都可以用被称为微平面(Microfacet)的细小镜面来进行描绘。根据平面粗糙程度的不同，这些细小镜面的取向排列可以相当不一致：
在这里插入图片描述

产生的效果就是：一个平面越是粗糙，这个平面上的微平面的排列就越混乱。这些微小镜面这样无序取向排列的影响就是，当我们特指镜面光/镜面反射时，入射光线更趋向于向完全不同的方向发散(Scatter)开来，进而产生出分布范围更广泛的镜面反射。而与之相反的是，对于一个光滑的平面，光线大体上会更趋向于向同一个方向反射，造成更小更锐利的反射：

在这里插入图片描述

在微观尺度下，没有任何平面是完全光滑的。然而由于这些微平面已经微小到无法逐像素地继续对其进行区分，因此我们假设一个粗糙度(Roughness)参数，然后用统计学的方法来估计微平面的粗糙程度。我们可以基于一个平面的粗糙度来计算出众多微平面中，朝向方向沿着某个向量 $h$ 方向的比例。这个向量 $h$ 便是位于光线向量 $l$ 和视线向量 $v$ 之间的半程向量(Halfway Vector)。它的计算方法如下：
$\frac{l+v}{||l+v||}$
微平面的朝向方向与半程向量的方向越是一致，镜面反射的效果就越是强烈越是锐利。通过使用一个介于0到1之间的粗糙度参数，我们就能概略地估算微平面的取向情况了：

在这里插入图片描述

能量守恒

微平面近似法使用了这样一种形式的能量守恒(Energy Conservation)：出射光线的能量永远不能超过入射光线的能量（发光面除外）。如上图我们可以看到，随着粗糙度的上升，镜面反射区域会增加，但是镜面反射的亮度却会下降。如果每个像素的镜面反射强度都一样（不管反射轮廓的大小），那么粗糙的平面就会放射出过多的能量，而这样就违背了能量守恒定律。这也就是为什么正如我们看到的一样，光滑平面的镜面反射更强烈而粗糙平面的反射更昏暗。

为了遵守能量守恒定律，我们需要对漫反射光和镜面反射光做出明确的区分。当一束光线碰撞到一个表面的时候，它就会分离成一个折射部分和一个反射部分。反射部分就是会直接反射开而不进入平面的那部分光线，也就是我们所说的镜面光照。而折射部分就是余下的会进入表面并被吸收的那部分光线，也就是我们所说的漫反射光照。

这里还有一些细节需要处理，因为当光线接触到一个表面的时候折射光是不会立即就被吸收的。通过物理学我们可以得知，光线实际上可以被认为是一束没有耗尽就不停向前运动的能量，而光束是通过碰撞的方式来消耗能量。每一种材料都是由无数微小的粒子所组成，这些粒子都能如下图所示一样与光线发生碰撞。这些粒子在每次的碰撞中都可以吸收光线所携带的一部分或者是全部的能量而后转变成为热量。
在这里插入图片描述

一般来说，并非全部能量都会被吸收，而光线也会继续沿着（基本上）随机的方向发散，然后再和其他的粒子碰撞直至能量完全耗尽或者再次离开这个表面。而光线脱离物体表面后将会协同构成该表面的（漫反射）颜色。不过在基于物理的渲染之中我们进行了简化，假设对平面上的每一点所有的折射光都会被完全吸收而不会散开。而有一些被称为次表面散射(Subsurface Scattering)技术的着色器技术将这个问题考虑了进去，它们显著地提升了一些诸如皮肤，大理石或者蜡质这样材质的视觉效果，不过伴随而来的代价是性能的下降。

对于金属(Metallic)表面，当讨论到反射与折射的时候还有一个细节需要注意。金属表面对光的反应与非金属（也被称为介电质(Dielectrics)）表面相比是不同的。它们遵从的反射与折射原理是相同的，但是所有的折射光都会被直接吸收而不会散开，只留下反射光或者说镜面反射光。亦即是说，金属表面只会显示镜面反射颜色，而不会显示出漫反射颜色。由于金属与电介质之间存在这样明显的区别，因此它们两者在PBR渲染管线中被区别处理。

反射光与折射光之间的这个区别使我们得到了另一条关于能量守恒的经验结论：反射光与折射光它们二者之间是互斥的关系。无论何种光线，其被材质表面所反射的能量将无法再被材质吸收。因此，诸如折射光这样的余下的进入表面之中的能量正好就是我们计算完反射之后余下的能量。

我们按照能量守恒的关系，首先计算镜面反射部分，它的值等于入射光线被反射的能量所占的百分比。然后折射光部分就可以直接由镜面反射部分计算得出：

float kS = calculateSpecularComponent(...); // 反射/镜面 部分
float kD = 1.0 - ks;                        // 折射/漫反射 部分

这样我们就能在遵守能量守恒定律的前提下知道入射光线的反射部分与折射部分所占的总量了。按照这种方法折射/漫反射与反射/镜面反射所占的份额都不会超过1.0，如此就能保证它们的能量总和永远不会超过入射光线的能量。

反射率方程

首先给出PBR光照模型的反射率方程：
$L_o(p,w_o) = \int_\Omega f_r(p,w_i,w_o)L_i(p,w_i)n\cdot w_idw_i$

更详细的内容，可以参考games101，第15课，讲解了这个方程是怎么得出的。

如果去除 $f_r(p,\omega_i,\omega_o)$ ，可以得到计算Irradiance的公式：

$\int_{\Omega}L_i(p,\omega_i)n\cdot \omega_i d\omega_i$ ，即点 $p$ 接受到的能量之和。

要正确理解这个方程式，我们需要涉足一些辐射度量学的内容。

辐射通量：辐射通量 $\phi$ 表示的是一个光源所输出的能量，以瓦特为单位。

立体角：立体角用 $\omega$ 表示，它可以为我们描述投射到单位球体上的一个截面的大小或者面积。

在这里插入图片描述

可以把自己想象成为一个站在单位球面的中心的观察者，向着投影的方向看。这个投影轮廓的大小就是立体角。

辐射强度：辐射强度(Radiant Intensity)表示的是在单位球面上，一个光源向每单位立体角所投送的辐射通量。举例来说，假设一个全向光源向所有方向均匀的辐射能量，辐射强度就能帮我们计算出它在一个单位面积（立体角）内的能量大小：

在这里插入图片描述

计算辐射强度的公式如下所示：
$\frac{d\Phi}{dw}$
其中 $I$ 表示辐射通量 $\Phi$ 除以立体角 $\omega$ 。

一个拥有辐射强度 $\Phi$ 的光源在单位面积 $A$ ，单位立体角 $\omega$ 上的辐射出的总能量，即**辐射率(radiance)**的方程如下所示：
$\frac{dI}{dAcos\theta} = \frac{d^2\Phi}{dAd\omega cos\theta}$
在这里插入图片描述

如果我们能够把立体角 $\omega$ 和面积 $A$ 看作是无穷小的，那么我们就能用辐射率来表示单束光线穿过空间中的一个点的通量。这就使我们可以计算得出作用于单个（片段）点上的单束光线的辐射率，我们实际上把立体角 $\omega$ 转变为方向向量 $\omega$ 然后把面 $A$ 转换为点 $p$ ，这样我们就能直接在我们的着色器中使用辐射率来计算单束光线对每个片段的作用了。

事实上，当设计到辐射率时，我们通常关心的是所有投射到点 $p$ 上的光线的总和，而这个和就称为辐射照度或者辐照度(Irradiance)。

再回过头来看反射率方程式：
$L_o(p,w_o) = \int_\Omega f_r(p,w_i,w_o)L_i(p,w_i)n\cdot w_idw_i$
反射率公式计算了点 $p$ 在 $\omega_o$ 方向上被反射出来的辐射率 $L_o(p,w_o)$ 的总和。换句话说： $L_o$ 表示了从 $\omega_o$ 方向上观察，光线投射到点 $p$ 上反射出来的辐照度。

基于反射率基于反射率公式是围绕所有入射辐射率的总和，也就是辐照度来计算的，所以我们需要计算的就不只是是单一的一个方向上的入射光，而是一个以点 $p$ 为求新的半球领域 $\Omega$ 内所有方向上的入射光。一个半球领域(Hemisphere)可以描述为以平面法线 $n$ 为轴所环绕的半个球体：

在这里插入图片描述

反射率方程概括了在半球领域 $\Omega$ 内，碰撞到了点 $p$ 上的所有入射方向 $ω_i$ 上的光线的辐射率，并受到 $f_r$ 的约束，然后返回观察方向上反射光的 $L_o$ 。入射光辐射率可以由光源处获得，此外还可以利用一个环境贴图来测算所有入射方向上的辐射率。

BRDF：最后剩下一个未知符号： $f_r$ ，它被称为BRDF，或双向反射分布函数(Bidirectional Reflective Distribution Function)，它的作用是是基于表面材质属性来对入射辐射率进行缩放或者加权。

BRDF，它接受入射方向 $\omega_i$ ，出射方向 $\omega_o$ ，平面法线 $n$ 以及一个用来表示微平面粗糙程度的参数 $a$ 作为函数的输入参数。BRDF可以近似的求出每一束光线对一个给定了材质属性的平面上最终反射出来的光线所作出的贡献程度。对于一个BRDF，为了实现物理学上的可信度，它必须遵守能量守恒定律，也就是说反射光线的总和永远不能超过入射光线的总量。

目前几乎所有实时渲染管线使用的都是一种被称为Cook-Torrance BRDF的模型。

Cook-Torrance BRDF兼有漫反射和镜面反射两个部分：
$f_r = k_df_{lambert} + k_sf_{cook-torrance}$
这里的 $k_d$ 是入射光线中被折射部分的能量所占的比率，而 $k_s$ 是被反射部分的比率。BRDF的左侧表示的是漫反射部分，这里用 $f_{lambert}$ 来表示。它被称为Lambertian漫反射，这和之前漫反射着色中使用的常数因子类似，用如下的公式来表示：
$f_{lambert} = \frac{c}{\pi}$
$c$ 表示表面颜色，除以 $\pi$ 是为了对漫反射进行标准化，因为前面含有BRDF的积分方程是由 $\pi$ 影响的。

BRDF的右侧表示镜面反射部分，它的形式如下所示：
$f_{cook-torrance} = \frac{DFG}{4(\omega_o \cdot n)(\omega_i \cdot n)}$
Cook-Torrance BRDF的镜面反射部分包含三个函数，此外分母部分还有一个标准化因子。字符D，F与G分别代表着一种类型的函数，各个函数分别用来近似地计算出表面反射特性地一个特定部分。

法线分布函数

估算在受到表面粗糙度的影响下，朝向方向与半程向量一致的微平面的数量。这是用来估算微平面的主要函数。举例来说，假设给定向量 $h$ ，如果我们的微平面中有35%与向量 $h$ 取向一致，则法线分布函数，或者说NDF将会返回0.35。目前有很多种NDF都可以从统计学上来估算微平面的总体取向度，只要给定一些粗糙度的参数。其中Trowbridge-Reitz GGX的分布函数为:
$NDF_{GGXTR}(n,h,a) = \frac{\alpha^2}{\pi((n\cdot h)^2(\alpha^2-1)+1)^2}$
这里 $h$ 表示的是与微平面做比较用的半程向量，而 $\alpha$ 代表表面粗糙度。

如果把 $h$ 当成是不同粗糙度参数下，平面法向量和光线方向向量之间的中间向量的话，我们可以得到如下图示的效果：

在这里插入图片描述

当粗糙度很低（也就是说表面很光滑）的时候，与半程向量取向一致的微平面会高度集中在一个很小的半径范围内。由于这种集中性，NDF最终会生成一个非常明亮的斑点。但是当表面比较粗糙的时候，微平面的取向方向会更加的随机。你将会发现与hℎ向量取向一致的微平面分布在一个大得多的半径范围内，但是同时较低的集中性也会让我们的最终效果显得更加灰暗。

几何函数

描述了微平面自成阴影的属性。当一个平面相对比较粗糙的时候，平面表面上的微平面有可能挡住其他的微平面从而减少表面所反射的光线。几何函数从统计学上近似的求得了微平面间相互遮蔽的比率，这种相互遮蔽会损耗光线的能量。

与NDF类似，几何函数采用一个材料的粗糙度参数作为输入参数，粗糙度较高的表面其微表面间相互遮蔽的概率就越高。我们将要使用的几何函数是GGX与Schlick-Beckmann近似的结合体，因此又称为Schlick-GGX：
$G_{SchlickGGX}(n,v,k) = \frac{n\cdot v}{(n\cdot v)(1-k)+k}$
这里 $k$ 是粗糙度 $\alpha$ 的重映射，取决于我们要用的是针对直接光照还是针对IBL光照的几何函数：
$k_{direct} = \frac{(\alpha+1)^2}{8}\\ k_{IBL} = \frac{\alpha^2}{2}$
为了有效的估算几何部分，需要将观察方向（几何遮蔽，Geometry Obstruction）和光线方向向量（几何阴影，Geometry Shadowing）都考虑进去。我们可以使用史密斯法(Smith’s method)来把两者都纳入其中：
$G(n,v,l,k) = G_{sub}(n,v,k)G_{sub}(n,l,k)$
使用史密斯法与Schlick-GGX作为 $G_{sub}$ 可以得到如下所示不同粗糙度的视觉效果。

在这里插入图片描述

几何函数是一个值域为[0.0, 1.0]的乘数，其中白色或者说1.0表示没有微平面阴影，而黑色或者说0.0则表示微平面彻底被遮蔽。

菲涅尔方程

菲涅尔（发音为Freh-nel）方程描述的是被反射的光线对比光线被折射的部分所占的比率，这个比率会随着我们观察的角度不同而不同。当光线碰撞到一个表面的时候，菲涅尔方程会根据观察角度告诉我们被反射的光线所占的百分比。利用这个反射比率和能量守恒原则，我们可以直接得出光线被折射的部分以及光线剩余的能量。

当垂直观察的时候，任何物体或者材质表面都有一个基础反射率(Base Reflectivity)，但是如果以一定的角度往平面上看的时候所有反光都会变得明显起来。你可以自己尝试一下，用垂直的视角观察你自己的木制/金属桌面，此时一定只有最基本的反射性。但是如果你从近乎90度（译注：应该是指和法线的夹角）的角度观察的话反光就会变得明显的多。如果从理想的90度视角观察，所有的平面理论上来说都能完全的反射光线。这种现象因菲涅尔而闻名，并体现在了菲涅尔方程之中。

菲涅尔方程是一个相当复杂的方程式，不过幸运的是菲涅尔方程可以用Fresnel-Schlick近似法求得近似解：
$F_{Schlick}(h,v,F_0) = F_0 + (1-F_0)(1-(h\cdot v))^5$
$F_0$ 表示平面的基础反射率，它是利用所谓折射指数(Indices of Refraction)或者说IOR计算得出的。然后正如你可以从球体表面看到的那样，我们越是朝球面掠角的方向上看（此时视线和表面法线的夹角接近90度）菲涅尔现象就越明显，反光就越强：
在这里插入图片描述

Cook-Torrance反射率方程

最终，我们可以将上述各项纳入到最终的反射率方程当中去：
$L_o(p,\omega_o) = \int_\Omega(k_d\frac{c}{\pi} + k_s\frac{DFG}{4(\omega_o\cdot n)(\omega_i \cdot n)})L_i(p,\omega_i)n\cdot \omega_i d\omega_i$

实践部分

PBR光照

这一部分介绍如果只计算来自于光源的光线，利用PBR方程如何计算。

假设我们在场景中只有一个光源，位于空间中的某一个点，因而对于一个物体表面点 $p$ 来说，其他可能的入射光线方向上的辐射率为0，示意图如下所示：

在这里插入图片描述

如果从一开始，我们就假设点光源不受光线衰减（光照强度会随着距离变暗）的影响，那么无论我们把光源放在哪，入射光线的辐射率总是一样的（除去入射角 $cos\theta$ 对辐射率的影响之外）。这是因为无论我们从哪个角度观察它，点光源总具有相同的辐射强度，我们可以有效地将其辐射强度建模为其辐射通量: 一个常量向量。

PBR的渲染结果以及加上贴图之后的渲染结果如下：

在这里插入图片描述

漫反射辐照度

仔细研究反射方程，可以发现BRDF的漫反射项和镜面反射项是相互独立的，因此，我们可以将积分分成两个部分：
$L_o(p,\omega_o) = \int_{\Omega}(k_d\frac{c}{\pi})L_i(p,\omega_i)n\cdot \omega_id\omega_i + \int_{\Omega}(k_s\frac{DFG}{4(\omega_o \cdot n)(\omega_i \cdot n)})L_i(p,\omega_i)n\cdot \omega_i d\omega_i$
观察漫反射积分，我们发现漫反射兰伯特项是一个常数项（颜色 $c$ 、折射率 $k_d$ 和 $\pi$ 在整个积分中是常数），不依赖于任何积分变量。基于此，我们可以将常数项移出漫反射积分：
$L_o(p,\omega_o) = k_d\frac{c}{\pi}\int_\Omega L_i(p,\omega_i)n\cdot \omega_id\omega_i$
这给了我们一个只依赖于 $\omega_i$ 的积分（假设 $p$ 位于环境贴图的中心）。有了这些知识，我们就可以计算或预计算一个新的立方体贴图，它在每一个采样方向——也就是纹素）——中存储漫反射积分的结果，这些结果是通过卷积计算出来的。

卷积的特性是，对数据集中的一个条目做一些计算时，要考虑到数据集中的所有其他条目。这里的数据集就是场景的辐射度或环境贴图。因此，要对立方体贴图中的每个采样方向做计算，我们都会考虑半球 $\Omega$ 上的所有其他采样方向。

为了对环境贴图进行卷积，我们通过对半球 $\Omega$ 上的大量方向进行离散采样并对其辐射度取平均值，来计算每个输出采样方向 $\omega_o$ 的积分。用来采样方向 $\omega_i$ 的半球，要面向卷积的输出采样方向 $\omega_o$ 。

在这里插入图片描述

这个预计算的立方体贴图，在每个采样方向 $\omega_o$ 上存储其积分结果，可以理解为场景中所有能够击中面向 $\omega_o$ 的表面的间接漫反射光的预计算总和。这样的立方体贴图被称为辐照度图，因为经过卷积计算的立方体贴图能让我们从任何方向有效地直接采样场景（预计算好的）辐照度。

下面是一个环境立方体贴图及其生成的辐照度图的示例，每个方向 $\omega_o$ 的场景辐射度取平均值。

在这里插入图片描述

首先，读入HDR环境贴图，将它从等距柱状投影图转换为立方体贴图，并将HDR立方体贴图作为天空盒渲染到场景中：

在这里插入图片描述

之后我们需要对立方体贴图进行卷积操作。

立方体贴图的卷积

我们的卷积方法是：对于立方体贴图的每个纹素，在纹素所代表的方向的半球 $\Omega$ 内生成固定数量的采样向量，并对采样结果取平均值。数量固定的采样向量将均匀地分布在半球内部。注意，积分是连续函数，在采样向量数量固定的情况下离散地采样只是一种近似计算方法，我们采样的向量越多，就越接近正确的结果。反射方程的积分 $\int$ 是围绕立体角 $d\omega$ 旋转累积的，而立体角相当难以处理。为了避免对难处理的立体角求积分，我们使用球坐标 $\theta$ 和 $\phi$ 来代替立体角。

在这里插入图片描述

对于立体角和球坐标，我们有这样的关系： $d\omega = sin\theta d\theta d\phi$ 。对于围绕半球大圆的航向角 $\phi$ ，我们在 $0$ 到 $2\pi$ 内采样，而从半球顶点出发的倾斜角 $\theta$ ，采样范围是 $0$ 到 $\frac{1}{2}\pi$ 。于是我们可以更新以下反射积分方程中的漫反射部分：
$L_o(p,\phi_o,\theta_o) = k_d\frac{c}{\pi}\int_{\phi=0}^{2\pi}\int_{\theta=0}^{\frac{1}{2}\pi}L_i(p,\phi_i,\theta_i)cos(\theta_i)sin(\theta_i)d\phi_i d\theta_i$
求解积分需要我们在半球 $\Omega$ 内采集固定数量的离散样本并对其结果求平均值。可以用蒙特卡洛法来对积分进行近似，这里因为对 $\theta$ 和 $\phi$ 都采用了均匀采样，我们令 $\theta$ 的概率密度函数 $pdf(\theta) = \frac{1}{\frac{1}{2}\pi} = \frac{2}{\pi}$ ， $\phi$ 的概率密度函数 $pdf(\phi) = \frac{1}{2\pi}$ 。有下式：
$L_o(p,\phi_o,\theta_o) = k_d\frac{c}{\pi}\frac{1}{n_1}\sum_{j=0}^{n_1}(( \frac{1}{n_2}\sum_{k = 0}^{n_2}\frac{L_i(p,\phi_{ij},\theta_{ik})cos\theta_{ik} sin\theta_{ik}}{pdf(\theta_{ik})} ) / pdf(\phi_{ij}))\\ = k_d\frac{c}{\pi}\frac{1}{n_1}\frac{1}{n_2}\frac{\pi}{2}2\pi \sum_{j=0}^{n_1}\sum_{k=0}^{n_2}L_i(p,\phi_{ij},\theta_{ik})cos\theta_{ik}sin\theta_{ik}\\ = \frac{k_d \cdot c \cdot \pi}{n_1\cdot n_2}\sum_{j=0}^{n_1}\sum_{k=0}^{n_2}L_i(p,\phi_{ij},\theta_{ik})cos\theta_{ik}sin\theta_{ik}$
利用上述公式，我们可以求得辐射度贴图，并利用该贴图进行渲染，结果如下图所示：

在这里插入图片描述

镜面反射辐照度

这一小节，我们只关注镜面反射部分。这一部分的反射方程式如下：
$L_o(p,w_o) = \int_\Omega k_s\frac{DFG}{4(\omega_o \cdot n)(\omega_i \cdot n)}L_i(p,\omega_i)n\cdot \omega_id\omega_i = \int_\Omega f_r(p,\omega_i, \omega_o)L_i(p, \omega_i)n\cdot \omega_id\omega_i$
参考Games202，Lecture3，我们有这样的近似：
$\int_\Omega f(x)g(x)dx \approx \frac{\int_\Omega f(x)dx}{\int_\Omega dx} \cdot \int_\Omega g(x)dx$
这个近似什么时候准确？

Small support（当 $g (x)$ 的积分的范围，积分域比较小的时候，或者当 $f (x)$ 的积分域比较小的时候）
Smooth integrand（ $当 g (x)$ 比较光滑的时候，或者当 $f (x)$ 比较光滑的时候）

因此对于我们上述的反射方程，我们也可以做一个近似：
$L_o(p,\omega_o) \approx \frac{\int_\Omega L_i(p,\omega_i)d\omega_i}{\int_\Omega d\omega_i} * \int_\Omega f_r(p,\omega_i, \omega_o)n\cdot \omega_id\omega_i$
卷积的第一部分被称为预滤波环境贴图，它类似于辐照度图，是预先计算的环境卷积贴图，但这次考虑了粗糙度。因为随着粗糙度的增加，参与环境贴图卷积的采样向量会更分散，导致反射更模糊，所以对于卷积的每个粗糙度级别，我们将按顺序把模糊后的结果存储在预滤波贴图的mipmap中。例如，预过滤的环境贴图在其5个mipmap级别中存储5个不同粗糙度值得预卷积结果，如下图所示：

在这里插入图片描述

我们使用Cook-Torrance BRDF的法线分布函数(NDF)生成采样向量及其散射强度，该函数将法线和视角方向作为输入。由于我们在卷积环境贴图时事先不知道视角方向，因此Epic Games假设视角方向——以及镜面反射方向——总是等于采样方向 $\omega_o$ 。以作进一步近似。翻译成代码如下：

vec3 N = normallize(w_o);
vec3 R = N;
vec3 V = R;

如何计算预滤波环境贴图

首先，对于卷积的每个粗糙度级别，我们将按顺序把模糊后的结果存储在预滤波贴图的mipmap中。

在漫反射辐照度小节中，我们使用球面坐标生成均匀分布在半球 $\Omega$ 上的采样向量，以对环境贴图进行卷积。虽然这个方法非常适用于辐照度，但对于镜面反射效果较差。镜面反射依赖于表面的粗糙度，反射光线可能比较松散，也可能比较紧密，但是一定会围绕着反射向量 $r$ ，除非表面极度粗糙：

在这里插入图片描述

所有可能出射的反射光构成的形状称为镜面波瓣。随着粗糙度的增加，镜面波瓣的大小增加；随着入射光方向不同，形状会发生变化。因此，镜面波瓣的形状高度依赖于材质。在微表面模型里给定入射光方向，则镜面波瓣指向微平面的半向量的反射方向。考虑到大多数光线最终会反射到一个基于半向量的镜面波瓣内，采样时以类似的方式选取采样向量是有意义的，因为大部分其余的向量都被浪费掉了，这个过程称为重要性采样。

对于蒙特卡洛积分，我在光线追踪笔记里已经有过详细介绍。除此之外，生成随机样本的方法也多种多样。默认情况下，每次采样都是我们熟悉的完全(伪)随机，不过利用半随机序列的某些属性，我们可以生成虽然是随机样本但具有一些有趣性质的样本向量。例如，我们可以对一种名为低差异序列的东西进行蒙特卡洛积分，该序列生成的仍然是随机样本，但样本分布更均匀：

在这里插入图片描述

当使用低差异序列生成蒙特卡洛样本向量时，该过程称为拟蒙特卡洛积分。拟蒙特卡洛方法具有更快的收敛速度，这使得它对于性能繁重的应用很有用。

鉴于我们新获得的有关蒙特卡洛（Monte Carlo）和拟蒙特卡洛（Quasi-Monte Carlo）积分的知识，我们可以使用一个有趣的属性来获得更快的收敛速度，这就是重要性采样。我们在前文已经提到过它，但是在镜面反射的情况下，反射的光向量被限制在镜面波瓣中，波瓣的大小取决于表面的粗糙度。既然镜面波瓣外的任何（拟）随机生成的样本与镜面积分无关，因此将样本集中在镜面波瓣内生成是有意义的，但代价是蒙特卡洛估算会产生偏差。

本质上来说，这就是重要性采样的核心：只在某些区域生成采样向量，该区域围绕微表面半向量，受粗糙度限制。通过将拟蒙特卡洛采样与低差异序列相结合，并使用重要性采样偏置样本向量的方法，我们可以获得很高的收敛速度。因为我们求解的速度更快，所以要达到足够的近似度，我们所需要的样本更少。因此，这套组合方法甚至可以允许图形应用程序实时求解镜面积分，虽然比预计算结果还是要慢得多。

低差异序列

具体的理论解释可以参考这篇文章。

GGX重要性采样

参考这篇文章，我们能够了解到NDF有一个归一化条件： $\int_\Omega D(m)cos\theta_md\omega = 1$ ，因此，我们可以根据NDF进行重要性采样，令 $pdf(\omega) = D(m)cos\theta_m$ 。

结合低差异序列和GGX重要性采样，我们可以实现预滤波器的卷积。

预计算BRDF

接下来集中精力求解求和近似的第二个部分：BRDF。

第二部分的方程式为：
$\int_\Omega f_r(p,\omega_i, \omega_o)n\cdot \omega_id\omega_i$
这个积分含有变量 $\omega_o, n, F_0$ 和粗糙度。由于我们使用的GGX是各向同性的，所以我们只需要知道 $\omega_o$ 和 $n$ 的夹角 $\theta$ 就好。

在着色器求解的过程中，我们可以设置 $n$ 为 $(0, 0, 1)$ ，即指向 $z$ 轴正前方，根据设置的 $n$ 以及夹角， $\theta$ ，我们可以求得一个 $\omega_o$ ，因为我们使用的GGX是各向异性的，因此这个 $\omega_o$ 便可代表所有可能的 $\omega_o$ 。

$F_0$ 由金属度和albedo决定。

现在我们还剩三个变量 $\theta$ 、 $F_0$ 和粗糙度。因此我们可以考虑将 $F_0$ 移出积分式。
$\begin{align} \int_\Omega f_r(p,\omega_i, \omega_o)n\cdot \omega_id\omega_i &= \int_\Omega f_r(p,\omega_i,\omega_o)\frac{F(\omega_o,h)}{F(\omega_o,h)}n\cdot \omega_id\omega_i \nonumber \\ &= \int_\Omega \frac{f_r(p,\omega_i,\omega_o)}{F(\omega_o,h)} \left(F_0 + (1 - F_0)(1 - \omega_o \cdot h)^5 \right) n \cdot \omega_id\omega_i \nonumber \\ &= F_0 \int_\Omega \frac{f_r(p,\omega_i, \omega_o)}{F(\omega_o,h)}\left( 1 - (1 - \omega_o \cdot h)^5 \right)n\cdot \omega_id\omega_i \nonumber \\ &+ \int_\Omega \frac{f_r(p,\omega_i,\omega_o)}{F(\omega_o,h)}(1 - \omega_o \cdot h)^5n\cdot \omega_id\omega_i \nonumber \\ &= F_0 * scale + bias \nonumber \end{align}$
$sc a l e$ 和 $bia s$ 可以根据法线分布函数进行重要性采样的蒙特卡洛积分来解决，其中 $sc a l e$ 的推导如下：

$\begin{align} scale &\approx \frac{1}{N}\sum_k^N \frac{ \frac{D(\omega_h^{(k)})G(\omega_o,\omega_i^{(k)})}{4(\omega_o\cdot n)(\omega_i^{(k)}\cdot n)} \left( 1 - (1 - \omega_o \cdot \omega_h^{(k)})^5 \right)(\omega_i^{(k)}\cdot n) }{p(\omega_i^{(k)})} \nonumber \\ &= \frac{1}{N}\sum_k^N \frac{ \frac{D(\omega_h^{(k)})G(\omega_o,\omega_i^{(k)})}{4(\omega_o\cdot n)(\omega_i^{(k)}\cdot n)} \left( 1 - (1 - \omega_o \cdot \omega_h^{(k)})^5 \right)(\omega_i^{(k)}\cdot n) }{\frac{D(\omega_h^{(k)})(\omega_h^{(k)}\cdot n)}{4(\omega_o\cdot \omega_h^{(k)})}} \nonumber \\ &= \frac{1}{N}\sum_k^N \frac{ G(\omega_o,\omega_i^{(k)})(\omega_o \cdot \omega_h^{(k)})\left( 1 - (1 - \omega_o \cdot \omega_h^{(k)})^5 \right)}{(\omega_o\cdot n)(\omega_h^{(k)}\cdot n)} \nonumber \end{align}$

$bia s$ 的推导与上面类似，其结果为：

$\approx \frac{1}{N}\sum_k^N \frac{ G(\omega_o,\omega_i^{(k)})(\omega_o \cdot \omega_h^{(k)})(1 - \omega_o \cdot \omega_h^{(k)})^5}{(\omega_o\cdot n)(\omega_h^{(k)}\cdot n)}$

其中， $p(\omega_i^{(k)}) = \frac{D(\omega_h^{(k)})(\omega_h^{(k)}\cdot n)}{4(\omega_o\cdot \omega_h^{(k)})}$ 为何成立，可以参考这篇博客中的推导。

用着色器来求解预计算BRDF的代码如下：

vec2 IntegrateBRDF(float NdotV, float roughness)
{
    vec3 V;
    V.x = sqrt(1.0 - NdotV*NdotV);
    V.y = 0.0;
    V.z = NdotV;

    float A = 0.0;
    float B = 0.0;

    vec3 N = vec3(0.0, 0.0, 1.0);

    const uint SAMPLE_COUNT = 1024u;
    for(uint i = 0u; i < SAMPLE_COUNT; ++i)
    {
        vec2 Xi = Hammersley(i, SAMPLE_COUNT);
        vec3 H  = ImportanceSampleGGX(Xi, N, roughness);
        vec3 L  = normalize(2.0 * dot(V, H) * H - V);

        float NdotL = max(L.z, 0.0);
        float NdotH = max(H.z, 0.0);
        float VdotH = max(dot(V, H), 0.0);

        if(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);

            A += (1.0 - Fc) * G_Vis;
            B += Fc * G_Vis;
        }
    }
    A /= float(SAMPLE_COUNT);
    B /= float(SAMPLE_COUNT);
    return vec2(A, B);
}
// ---------------------------------------------------------------------------
void main() 
{
    vec2 integratedBRDF = IntegrateBRDF(TexCoords.x, TexCoords.y);
    FragColor = integratedBRDF;
}

可以看到，上述的代码和公式推导是一一对应的。

预计算的BRDF的可视化结果如下，这张图也称为BRDF LUT贴图：

在这里插入图片描述