图形学基础|PBR回顾

图形学基础|PBR回顾

一、前言

在辐射度量学和渲染方程中简单介绍了辐射度量学以及提及了渲染方程。

实时渲染，其本质就是模拟光线与各类材质进行交互的过程，在有限的时间内尽可能准确地近似求解渲染方程。

这就不得不提到基于物理的渲染（Physically Based Rendering，PBR）。

在之前学习OpenGL的过程中，也曾学习过一些PBR相关的基础，大多来自Learn OpenGL PBR Theory。

图形学基础 | 基于物理的渲染理论（PBR)

图形学基础 | 基于物理的渲染PBR之直接光照

图形学基础 | 基于物理的渲染PBR之基于图像的光照IBL(Diffuse篇)

图形学基础 | 基于物理的渲染PBR之基于图像的光照IBL(Specular篇上)

当时对上述的知识囫囵吞枣，没有进行比较深刻的比较理解，在实现方面也比较简单，有一些预计算的贴图都是从网络上截取的资源。

这些过程的简化，对自己的图形渲染水平的提高都是不利的。因而，找个机会对这些知识再一次温习总结并实践理解。

本文的主要内容为：

• 概率论相关的基础知识，这部分是预计算时进行的采样和积分计算的基础知识。
• 蒙特卡洛积分，结合不同的采样给出如何对定积分进行近似求解。
• 反射方程，结合UE4的ShadingModel介绍反射方程中的各个项。
• IBL，以公式推导结合代码的方式对环境光的实现进行介绍。

二、概率论基础

2.1 随机变量

随机试验（E）

1. 可以在相同的条件下重复地进行
2. 每次试验的可能结果不止一个，并且能事先明确试验的所有可能结果
3. 进行一次试验之前不确定哪一个结果会出现

例子：抛一枚硬币，观察正面，反面出现的情况。

样本空间（S）

随机试验所有可以能的结果组成的集合。

样本点

样本空间的元素，即每个可能的结果。

随机事件

随机试验E的样本空间S的子集称为随机事件。

在随机现象中，有很多问题与数值有密切的关系，这时引入了随机变量的概念。

将定义在样本空间$S$上的函数称之为随机变量。这种函数的取值是不可预先确定的，是随机的。

随机变量可以分为：离散型和非离散型。

其中非离散型主要以连续型随机变量为主。离散型研究其分布律，而连续型则研究其概率密度（某一个点的改变不影响总体，因此求积分时开区间和闭区间没有太大区别）。

离散型随机变量：随机变量可能取到的值时有限个数或可列无限多个。

连续型随机变量：随机变量可能取到的值时无限个数。

这里举一个离散型随机变量的例子：

抛一枚骰子，观察其正面向上的数值，则有1，2，3，4，5，6这些结果。

2.2 分布函数（CDF）和概率密度函数（PDF）

随机变量的分布函数定义如下：（设$X$为随机变量，x为任意实数）

$$F(x)=P(X\le x),x\in (-\infty ,\infty )$$

2.2.1 离散型变量

对于离散型随机变量$X$，可以取的值有$x1,...,xi,...,xn$ ，对应的概率为$P(x1),...,P(xi),...,P(xn)$。

那么，其分布函数则为：

$$F(x)=P(X\le x)=\sum _{xi<x}P(x_i)$$

对满足条件的随机变量对应的概率进行累加。

其中有很多常用的离散型随机分布，例如：

• 伯努利分布、二值分布、泊松分布。

2.2.2 连续型变量

对于连续型随机变量$X$，$f(x)$为x的概率密度函数（PDF），简称为概率密度。

概率密度表示：随机变量取到$x$值的概率为$f(x)$。

分布函数定义如下：

$$F(x)={\int }_{-\infty }^{x}f(t)dt$$

概率密度函数的积分，即围成的面积，为随机变量落入某一区间的概率：

$$F(x1\le x\le x_2)={\int }_{x_1}^{x_2}f(t)dt$$

其中有很多常用的随机型随机分布，例如：

• 均匀分布、正态分布、指数分布

更多详细细节可以参考 随机变量及其分布

2.3 期望和方差

数学期望，简称期望，又称为均值。数学期望完全由随机变量的分布所确定，若$X$服从某一分布，也称$E(X)$是这一分布的数学期望。

离散型随机变量：

$$E(X)=\sum _{k=1}^{k=\infty }{x}_k{p}_k$$

连续型随机变量：

$$E(X)={\int }_{-\infty }^{\infty }xf(x)dx$$

方差，表达了随机变量X的取值与其数学期望的偏离程度。

$$\begin{array}{rl}D(X)=& Var(X)=E[(X-E[X]{)}^2]\\ \sigma （X)=& \sqrt{D(X)}\end{array}$$

其中，

由定义可知，方差实际上就是随机变量X的函数$g(x)$的数学期望。

$$g(X)=(X-E(X){)}^2$$

离散型随机变量的方差：

$$D(X)=\sum _{k=1}^{\infty }[x_k-E(X){]}^2{p}_k$$

连续型随机变量的方差：

$$D(X)={\int }_{-\infty }^{\infty }(x-E(X){)}^{2}f(x)dx$$

2.4 二维随机变量

对于二维随机变量，而需要了解以下的概念。

设现在有二维连续型随机变量$(X,Y)$，二维随机变量的联合概率密度函数为$f(x,y)$，那么其联合分布函数为：

F ( x , y ) = P { X < x , Y < y } = ∫ − ∞ x ∫ − ∞ y f ( u , v ) d u d v F(x,y)=P\{X<x,Y<y\}=\int_{-\infty}^{x}\int_{-\infty}^{y}f(u,v)dudv F(x,y)=P{X<x,Y<y}=∫−∞x​∫−∞y​f(u,v)dudv

联合概率密度满足：

$$f(x,y)\ge 0,{\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }f(x,y)dxdy=1$$

$X$的边缘概率密度为：

$$f_X(x)={\int }_{-\infty }^{\infty }(x,y)dy$$

在$X=x$时，$Y$的条件概率密度为：

$$f_{Y|X}(y)=\frac{f(x,y)}{f_X(x)}$$

2.5 概率分布变换

已知一个随机变量$X$的概率密度为 $f_X(x)$，现在我们令随机变量$Y=T(X)$。

尝试求出$Y$的概率密度函数$f_Y(y)$。为了计算方便，仅考虑T为单调递增的情况。

由分布函数定义已知：

P { Y < T ( x ) } = P { X < x } P\{Y<T(x)\} = P\{X<x\} P{Y<T(x)}=P{X<x}

因此有：

$$f_Y(y)=f_Y(T(x))=f_X(x)$$

对两边一起求导：

$$\begin{array}{rl}{f}_y(y)\frac{dy}{dx}=& f_X(x)\\ f_y(y)=& (\frac{dy}{dx}{)}^{-1}{f}_X(x)\end{array}$$

这样，即可求得Y的概率密度函数。

而对于多维随机变量，设$X$和$Y$都是n维随机变量，$X$和$Y$有以下转移关系：$Y=M(X)$。

有以下结论：

$$f_Y(y)=f_Y(M(x))=\frac{fX(x)}{|J_T(x)|}$$

其中，$|J_T(x)|$表示M的雅克比矩阵的行列式的绝对值。M的雅克比矩阵为：

$$\left(\begin{array}{ccc}\partial M_1/\partial x_1& \cdots & \partial M_1/\partial x_n\\ ⋮& \ddots & ⋮\\ \partial M_n/\partial x_1& \cdots & \partial M_n/\partial x_n\end{array}\right)$$

在图形渲染中，使用得最多的有以下的情况：

球坐标系

单位球上的坐标$(\theta ,\varphi )$和笛卡尔坐标系$(x,y,z)$的转换如下：

$$\begin{array}{rl}x=& cos\theta cos\varphi \\ y=& sin\theta sin\varphi \\ z=& cos\theta \end{array}$$

根据推导有，相应的概率密度为：

$$f(\theta ,\varphi )=sin\theta f(x,y,z)$$

在球坐标系中，立体角的定义如下：

$$dw=sin\theta d\theta d\varphi$$

立体角在某个$\Omega$范围内的概率为：

$$P(w\in \Omega )={\int }_{\Omega }f(w)dw$$

又$f(\theta ,\varphi )d\theta d\varphi =f(w)dw$可以推导得到：

$$f(\theta ,\varphi )=sin(\theta )f(w)$$

注！这个概率密度公式在后续大有作为！

三、蒙特卡洛积分和重要性采样

蒙特卡洛方法是一种采用统计思路，用于求解难以直接求出具体值问题的方法。

在图形学中，蒙特卡洛方法最主要的应用就是帮助我们快速的求解积分。

什么是蒙特卡洛积分呢？

简而言之就是，在求积分时，如果找不到被积函数的原函数，那么利用经典积分方法是得不到积分结果的。

但是蒙特卡洛积分方法告诉我们，利用一个随机变量对被积函数进行采样，并将采样值进行一定的处理，那么当采样数量达到一定时，得到的结果可以很好的近似原积分的结果。这样一来，我们就不用去求原函数的形式，就能求得积分的近似结果。

而重要性采样和蒙特卡洛积分二者是两个不同的概念，仅仅是因为经常一起使用。

3.1 黎曼和

在介绍蒙特卡洛积分之前，让我们先来了解一下黎曼和（Riemann Sum）。

黎曼和是一种很简单的估计定积分的数值方法。

给定我们要求解的定积分：${\int }_a^{b}f(x)dx$

如图所示，对于一个这样不规则的图形，很难直接求解出它的面积。

不过根据积分的第一中值定理，在积分区间$[a,b]$一定可以找到存在一点$\epsilon$使得：

$${\int }_a^{b}f(x)dx=f(\epsilon )(b-a)$$

即可以转换为一个长度为$b-a$，高度为$f(\epsilon )$的矩形的面积。

那么，如何求解这个高度$f(\epsilon )$呢？

Riemann Sum是这样来求解的：

• 设定一个步进值step，从函数的定义域$[a,b]$起点出发，一步一步的往终点去移动，计算每一个函数的值，然后将这些值求平均数（除以步进的步数）。

$\frac{1}{N}{\sum }_{i=1}^{N}f(x_i)$

则该定积分的估计为：

$${\int }_a^{b}f(x)dx\approx \frac{b-a}{N}\sum _{i=1}^{N}f(x_i)$$

3.2 蒙特卡洛积分

同样的，对于定积分${\int }_a^{b}f(x)dx$，蒙特卡洛积分的公式如下：

$${\int }_a^{b}f(x)dx\approx \frac{1}{N}\sum _{i=1}^N\frac{f(x_i)}{pdf(x_i)}$$

其中，$pdf(x)$表示概率密度分布函数。

我们对这个公式求期望E：

$$\begin{array}{rl}E[\frac{1}{N}\sum _{i=1}^N\frac{f(x_i)}{p(x_i)}]=& \frac{1}{N}\sum _{i=1}^{N}E[\frac{f(x_i)}{p(x_i)}]\\ =& \frac{1}{N}\sum _{i=1}^N{\int }_a^b\frac{f(x)}{p(x)}p(x)dx\\ =& \frac{1}{N}\sum _{i=1}^N{\int }_a^{b}f(x)dx\\ =& {\int }_a^{b}f(x)dx\end{array}$$

蒙特卡洛积分的数学期望与所要求的定积分是完全相同的。

pdf的准确度会大大的影响我们的收敛速度（也就是N的大小）。

当pdf较为准确时，那么我们只需少数的样本就对我们的积分值有一个较为准确的估计。

当pdf不够准确时，那么我们只需很多的样本才能对我们的积分值有一个较为准确的估计。

3.3 采样

在有了蒙特卡洛积分之后，我们就需要确定随机变量$X$的概率分布了。

由于我们并不知道随机变量$X$的概率密度函数$p(X)$，所以我们只能进行猜测。

其中的一种方案就是假设$X$在样本空间是均分分布的，那么$p(X)$即为一个常数。

3.3.1 均匀采样

均匀采样是图形学中常见的一个问题和方法。

这类问题的描述如下：

如何用[0,1]之间的均匀采样获得在XXX上的均匀采样？

其中，XXX可以是圆盘，球面，球体等。我们可以利用的只有[0,1]之间的均匀采样（可以通过程序的伪随机实现，后续再单独进行介绍）。

下面我们将对一些均匀采样的获得进行推导介绍。

单位半球面均匀采样

对于半球面的均匀采样而言，采样的是面积，半球面的面积可以用立体角$\omega$来表示，即$p(\omega )$为常数，根据半球面积可求得：

$$p(w)=\frac{1}{2\pi }$$

再结合2.4（概率分布变换）中所得的结论，则可以得到球坐标下的概率密度：

$$p(\theta ,\varphi )=sin\theta f(w)=\frac{sin\theta }{2\pi }$$

有了pdf，则可以求得边缘pdf以及条件pdf：

$$p(\theta )={\int }_0^{2\pi }p(\theta ,\varphi )d\varphi =sin\theta$$
$$p(\varphi |\theta )=\frac{p(\theta ,\varphi )}{p(\theta )}=\frac{1}{2\pi }$$

接下来可以求解CDF（概率分布函数）：
$$P(\theta )={\int }_0^{\theta }sin{\theta }^{\prime }d{\theta }^{\prime }=1-cos\theta$$

$$P(\varphi )={\int }_0^{\varphi }\frac{1}{2\pi }d{\varphi }^{\prime }=\frac{\varphi }{2\pi }$$

根据生成的[0,1]之间的均匀采样，则可以反求出对应的$(\theta ,\varphi )$：

$$\theta =arccos(1-{\xi }_1)$$
$$\varphi =2\pi {\xi }_2$$

其中，${\xi }_1$和${\xi }_2$则是通过[0,1]均匀分布生成的。

单位半球面均匀采样，是在图形渲染中经常使用到的一个采样方式。

单位半球面余弦权重采样

余弦采样指的是采样点在半球面分布的概率$p(w)$与角度$\theta$的余弦$cos\theta$正相关，概率密度的求解过程如下：

$$\begin{array}{rl}{\int }_{\Omega }p=& 1\\ {\int }_0^{2\pi }{\int }_0^{\pi /2}c\ast cos\theta sin\theta d\theta d\varphi =& 1\\ c=& \frac{1}{\pi }\\ p(\theta ,\varphi )=& \frac{1}{\pi }cos\theta sin\theta \end{array}$$

有了pdf，则可以求得边缘pdf以及条件pdf：

$$\begin{array}{rl}p(\theta )=& sin2\theta \\ p(\varphi |\theta )=& \frac{1}{2\pi }\end{array}$$

再根据pdf求解CDF：

$$\begin{array}{rl}P(\theta )=& \frac{1-cos2\theta }{2}=si{n}^2\theta \\ P(\varphi )=& \frac{\varphi }{2\pi }\end{array}$$

根据生成的[0,1]之间的均匀采样，则可以反求出对应的$(\theta ,\varphi )$：

$$\begin{array}{rl}\theta =& \sqrt{1-{\xi }_1^2}\\ \varphi =& 2\pi {\xi }_2\end{array}$$

通过上述两个例子，可以小结一下求解思路：

1. 用PDF来描述题目要求，求出所有变量的联合PDF；
2. 求出第一个变量的边缘PDF
3. 根据已知的边缘PDF或条件PDF，依次求其余变量的条件PDF
4. 求所有变量的CDf，并令其分别等于${\xi }_i$（来自于[0,1]均匀分布的采样值）
5. 反求各个变量和${\xi }_i$的关系

备注：这里我们来看一种均匀分布的2D随机采样。

他将十进制转换成二进制，再将二进制转换到[0,1]之间的小数，这一过程被称作 Radical Inverse。具体的方法如下图。

Hammersley采样点集合：

$(x_i,y_i)=(i/N,\varphi (i))$

其中，N为采样点总数，$\varphi (i)$为RadicalInverse之后得到的小数。

代码实现如下：

float RadicalInverse( uint bits ){
          //reverse bit
          //高低16位换位置
          bits = (bits << 16u) | (bits >> 16u); 
          //A是5的按位取反
          bits = ((bits & 0x55555555) << 1u) | ((bits & 0xAAAAAAAA) >> 1u);
          //C是3的按位取反
          bits = ((bits & 0x33333333) << 2u) | ((bits & 0xCCCCCCCC) >> 2u);
          bits = ((bits & 0x0F0F0F0F) << 4u) | ((bits & 0xF0F0F0F0) >> 4u);
          bits = ((bits & 0x00FF00FF) << 8u) | ((bits & 0xFF00FF00) >> 8u);
          return  float(bits) * 2.3283064365386963e-10;
}

vec2 Hammersley(uint i,uint N){
          return vec2(float(i) / float(N), RadicalInverse(i))
}

3.3.2 重要性采样

为什么要使用重要性采样？

比如：

• 下图表示一束光线沿着$w_i$方向，遇到一个微平面的反射情况，反射的光线呈现一个椭圆。
• 如果我们想要计算所有反射光线总和，可以对每一个方向都进行计算然后累加即可。
• 但我们可以知道如果我们的采样的光线超过这个椭圆的范围，那么贡献为0，对最后的结果不会有任何的影响。
• 因此，我们需要尽量采样落在椭圆范围内的光线，使用少量的采样获得较好的增益。

因此，重要性采样的目的就是有针对性的进行采样计算，而不是完全的进行随机采样，从而减少一些无效值的计算。

重要性采样小结这篇文章给出了如何根据NDF得到相应的重要性采样方法。

例如通过GGX的NDF推导出了类似于我们上述均匀采样所得到的$(\theta ,\varphi )$的计算公式。

四、反射方程

基于物理的渲染遵循的是反射方程（The Reflectance Equation），这是渲染方程特化版本。

其公式如下：

$$L_o(p,w_o)={\int }_{\Omega }{f}_r(p,w_i,w_o)(w_i\cdot n)d{w}_i$$

其中，$f_r$ 为BRDF（双向反射分布函数）。

下面将介绍一下《Real Shading in Unreal Engine 4》中提及的虚幻的着色模型（Shading Model）。

$$f_r=k_d{f}_{lambert}+k_s{f}_{cook-torrance}$$

其中，$k_s$表示入射光被反射部分的比例，$k_d$表示被折射的比例。

4.1 Diffuse BRDF

Unreal在评估了Burley的漫反射模型，只观察到和与Lambertain模型相比轻微的差别。

除此之外，任何更复杂的漫反射模型很难高效的使用基于图像光照或球谐光照（IBL和SH）。

UE在Diffuse BRDF选择了Lambertain模型，公式如下：

$$f(l,v)=\frac{c_{diff}}{\pi }$$

其中，$c_{diff}$为材质的漫反射率（diffuse albedo）。

float3 Diffuse_Lambert( float3 DiffuseColor ) { 
	return DiffuseColor * (1 / PI); 
}

4.2 Microfacet Specular BRDF

常规的Cook-Torrance微表面镜面反射着色模型如下：

$$f(l,v)=\frac{D(h)F(v,h)G(l,v,h)}{4(n\cdot l)(n\cdot v)}$$

其中，各个项都具有其特定的含义。

• D，法线分布函数 Normal Distribution Function（NDF）。
• G，几何函数。
• F，菲涅尔方程。

4.2.1 法线分布函数 D

D，法线分布函数 Normal Distribution Function（NDF），描述的是微表面法线的分布。

当越多的微表面法线与宏观表面法线相近的时候，高光区域会越亮，区域越小，能量越集中。

反之，高光区域会越暗，区域越大。

举例来说，假设给定一个向量v，如果我们的微平面中有35%与向量v取向一致，则法线分布函数（NDF）将会返回0.35。

UE这里选择了GGX法线分布函数，原因如下：

• 相比于Blinn-Phong的公式，GGX的额外消耗更小，并提供了更好的”长尾效果”，这个长尾效果的表现吸引了艺术家。

同时，其采用了迪士尼重新定义的参数$\alpha =Rougnes{s}^2$。

D项的公式如下：

$$D(h)=\frac{a^2}{\pi ((n\cdot h{)}^2(a^2-1)+1{)}^2}$$

该函数返回的是一个标量。

主流NDF函数高光长尾效果的对比，GGX拥有最好的长尾表现，效果如下：

UE4中着色器代码的D项实现：

float D_GGX( float Roughness, float NoH )
{
	float a = Roughness * Roughness;
	float a2 = a * a;
	float d = ( NoH * a2 - NoH ) * NoH + 1;	// 2 mad
	return a2 / ( PI*d*d );			// 4 mul, 1 rcp
}

4.2.2 几何函数 G

G，几何函数，从统计学上近似地计算微平面间相互遮蔽的比率。这种相互遮蔽会损耗光线的能量。

几何函数采用一个材料的粗糙度参数作为输入参数，粗糙度较高的表面其微平面间相互遮蔽的概率就越高。

UE评估了几何衰减项的多种选择，最后选择使用Schlick模型。但$k=\frac{\alpha }{2}$，以更好地拟合Smith模型GGX。

通过这种修改，Schlick模型与$\alpha =1$的Smith完全匹配，并且相当$[0,1]$范围内的近似逼近。如下图所示。

公式如下：

$$\begin{array}{rl}k=& \frac{(roughness+1{)}^2}{8}\\ G_l(v)=& \frac{n\cdot v}{(n\cdot v)(1-k)+k}\\ G(l,v,h)=& G_l(l)\ast G_l(v)\end{array}$$

几何函数（Geometry Function）是一个0到1之间的标量，描述了微平面自阴影的属性，表示了具有半矢量法线的微平面（microfacet）中，同时被入射方向和反射方向可见（没有被遮挡的）的比例，即未被遮挡的m= h微表面的百分比。

几何函数（Geometry Function）即是对能顺利完成对光线的入射和出射交互的微平面概率进行建模的函数。

• 观察方向，几何遮蔽（Geometry Obstruction）
• 光线入射方向，几何阴影（Geometry Shadowing）

之前的Schlick模型的实现代码如下，但目前UE中并没有使用这个。

// Tuned to match behavior of Vis_Smith
// [Schlick 1994, "An Inexpensive BRDF Model for Physically-Based Rendering"]
float Vis_Schlick( float a2, float NoV, float NoL )
{
	float k = sqrt(a2) * 0.5;
	float Vis_SchlickV = NoV * (1 - k) + k;
	float Vis_SchlickL = NoL * (1 - k) + k;
	return 0.25 / ( Vis_SchlickV * Vis_SchlickL );
}

当前UE使用的几何函数G为：

// Appoximation of joint Smith term for GGX
// [Heitz 2014, "Understanding the Masking-Shadowing Function in Microfacet-Based BRDFs"]
float Vis_SmithJointApprox(float a2, float NoV, float NoL)
{
    float a = sqrt(a2);
    float Vis_SmithV = NoL * (NoV * (1 - a) + a);
    float Vis_SmithL = NoV * (NoL * (1 - a) + a);
    return 0.5f * 1.0f /(Vis_SmithV + Vis_SmithL);
}

4.2.3 菲涅尔方程 F

F，菲涅尔方程，表示的是入射光线被表面反射所占的比例。

这个比例会随着我们观察的角度的不同而不同。

经典的Schlick近似公式如下：

$$F(h,v,F_0)=F_0+(1-F_0)(1-h\cdot v{)}^5$$

其中，$F0$为表面的基本反射属性。

F项描述的是镜面反射部分的比例，那么漫反射部分的比例自然就是1-F。

对于金属材质而言，只有镜面反射，没有漫反射。则对于F项而言，需要区别对待金属材质和非金属材质。

因而，材质的属性输入存在一个Metallic金属项。

前面所提到的$k_s$其本质就是F项，因而可以求解出$k_d$，代码如下：

vec3 ks = F;
vec3 kd = vec3(1.0f)-ks;
kd *= 1.0f - Metallic; // 对于金属而言，只有镜面反射（Metallic=1，kd=0）

虚幻采用对Schlick的近似作为F项，即使用了球面高斯Spherical Gaussian近似代替power计算，以此提高计算效率，并且差异微不可察，公式为：

$$F(h,v,F_0)=F_0+(1-F_0)2^{(-5.55473(v\cdot h)-6.98316)(v\cdot h)}$$

采用球面高斯的近似来代替Power的选择，主要是考虑在GPU上运算的效率来考虑（GPU指令减少了50%）。

然而，实际F项并没有使用球面高斯近似，而是采用一种更高效的方式：

• Fc + (1-Fc) * SpecularColor，SpecularColor本质就是$F_0$。

// [Schlick 1994, "An Inexpensive BRDF Model for Physically-Based Rendering"]
float3 F_Schlick( float3 SpecularColor, float VoH )
{
	float Fc = Pow5( 1 - VoH );					// 1 sub, 3 mul
	//return Fc + (1 - Fc) * SpecularColor;		// 1 add, 3 mad
	
	// Anything less than 2% is physically impossible and is instead considered to be shadowing
    // 50的意义：F0小于0.02在物理上是没有意义的，小于0.02的值被当做阴影处理。
	return saturate( 50.0 * SpecularColor.g ) * Fc + (1 - Fc) * SpecularColor;
}

至此，对于UE的PBR模型的各项都有了基本的了解，在直接光照中，可以依据材质的属性对各个项进行计算，最后获得光照结果。

五、IBL

PBR 其本质是一个光照模型，相比于经典的光照模型：

• 漫反射采用 Lambert

• 镜面反射采用 blinn-Phong

• 环境光采用 C*(1-AO)

在PBR中：

• 漫反射为：Lambert/PI
• 镜面反射采用 Cook-Torrance 模型；
• 环境光则通过 IBL（Image-Based Lighting）实现。

对于直接光照，漫反射和镜面反射的计算在前面已经介绍过光照模型了。

对于一个着色点而言，其所受的环境光来自于四面八方（各个方向），因而无法实时地对每个着色点的环境光进行计算，需要采用一系列算法进行预计算，从而在实时中可以直接获取或经过简单的计算进行实现环境光效果。

一般而言，可以通过一张环境贴图（Environment Map）来保存一个物体周围的环境信息。然后通过某种处理，来实现丰富的环境光照效果。

漫反射环境光照部分，一般采用传统IBL中Irradiance Environment Mapping技术，也可以采用SH进行存储和重建（可参考另一篇博文球谐光照（Spherical Harmonics Lighting））。

镜面反射环境光照部分，一般则采用了UE提供的Split Sum Approximation分解求和近似。

下面将分别进行介绍和提供实现代码。

5.1 漫反射环境光照

在这里，我们介绍IrradianceMap。

首先，看一下反射方程的漫反射部分。

$$\begin{array}{rl}{L}_{diff}=& k_d{\int }_{\Omega }{f}_{d}Li(p,w_i)(n\cdot w_i)d{w}_i\\ f_d=& \frac{c_{diff}}{\pi }\end{array}$$

其中，$k_d$表示环境光中漫反射的比例。又由于$f_d=\frac{c_{diff}}{\pi }$为常数。可以提取到积分符号之外，则有：

$$\begin{array}{rl}{L}_{diff}=& k_d\frac{c_{diff}}{\pi }{\int }_{\Omega }Li(p,w_i)(n\cdot w_i)d{w}_i\\ =& k_d\frac{c_{diff}}{\pi }{\int }_0^{2\pi }{\int }_0^{\frac{\pi }{2}}{L}_i(p,w_i)cos\theta sin\theta d\theta d\varphi \end{array}$$

现在，积分项只取决于法线$n$，与观察方向（相机方向）$w_0$无关，$\Omega$是指法线方向所在的上半球。

因而，积分项是可以提前进行计算的！在实时渲染时仅需用法线$n$进行采样即可。

这个思路就是IrradianceMap。

要求解定积分的数值可以采用蒙特卡洛积分。但不同的采样方式获得求解公式亦不同。

5.1.1 均匀采样

这里的均匀采样指的是对角度$(\theta ,\varphi )$进行均匀采样。

相当于将角度看作一个二维平面，进行均匀采样。

则概率密度函数为：

$$p({\theta }_i,{\varphi }_j)=\frac{1}{2\pi }\ast \frac{1}{0.5\pi }$$

则积分项可以进行如下近似：

$$\begin{array}{rl}{\int }_{\Omega }Li(p,w_i)(n\cdot w_i)d{w}_i=& {\int }_0^{2\pi }{\int }_0^{\frac{\pi }{2}}{L}_i(p,w_i)cos\theta sin\theta d\theta d\varphi \\ =& \frac{1}{N_1{N}_2}\sum _i^{N_1}\sum _j^{N_2}\frac{L_i(p,{\theta }_i,{\varphi }_j)cos{\theta }_{i}sin{\theta }_i}{p({\theta }_i,{\varphi }_j)}\\ =& \frac{1}{N_1{N}_2}\sum _i^{N_1}\sum _j^{N_2}\frac{L_i(p,{\theta }_i,{\varphi }_j)cos{\theta }_{i}sin{\theta }_i}{\frac{1}{2\pi }\frac{1}{0.5\pi }}\\ =& \frac{{\pi }^2}{N_1{N}_2}\sum _i^{N_1}\sum _j^{N_2}{L}_i(p,{\theta }_i,{\varphi }_j)cos{\theta }_{i}sin{\theta }_i\end{array}$$

代码如下：

// NumSamplesPerDir 每个方向采样的数量
// sampleDelta 采样的角度间隔

float sampleDelta = 1.0 / NumSamplesPerDir;
float NumSamples = 0.0;
for (float phi = 0.0; phi < 2.0 * PI; phi += sampleDelta)
{
    for (float theta = 0.0; theta < 0.5 * PI; theta += sampleDelta)
    {
        // spherical to cartesian (in tangent space)
        float sintheta = sin(theta);
        float costheta = cos(theta);
        float3 tangentSample = float3(sintheta * cos(phi),  sintheta * sin(phi), costheta);
        // tangent space to world
        float3 sampleVec = TangentToWorld(tangentSample, Normal);
        // 采样颜色
        float3 sampleColor = CubeEnvironment.SampleLevel(LinearSampler, sampleVec, 0).rgb;
		// 累加Radiance
        Irradiance += sampleColor * costheta * sintheta;
        // 采样数量++
        NumSamples++;
    }
}
// 除以采样数量
Irradiance = PI*Irradiance / NumSamples;

5.1.2 余弦权重的半球采样

将该概率密度函数应用到Irradiancemap求解的积分项中则有：

$$\begin{array}{rl}{\int }_{\Omega }Li(p,w_i)(n\cdot w_i)d{w}_i=& {\int }_0^{2\pi }{\int }_0^{\frac{\pi }{2}}{L}_i(p,w_i)cos\theta sin\theta d\theta d\varphi \\ =& \frac{1}{N}\sum _i^N\frac{L_i(p,\theta ,\varphi )cos\theta sin\theta }{p(\theta ,\varphi )}\\ =& \frac{1}{N}\sum _i^N\frac{L_i(p,\theta ,\varphi )cos\theta sin\theta }{\frac{1}{\pi }cos\theta sin\theta }\\ =& \frac{\pi }{N}\sum _i^N{L}_i(p,{\omega }_i)\end{array}$$

代码如下：

float3 Irradiance = 0;
// 一个采样数量值
const uint NumSamples = 400;
for (uint i = 0; i < NumSamples; i++)
{
    float2 E = Hammersley(i, NumSamples, Random);
    float3 L = TangentToWorld(CosineSampleHemisphere(E).xyz, N);
    float NoL = saturate(dot(N, L));
    float Factor = NoL > 0 ? 1.0 : 0.0;
    Irradiance += CubeEnvironment.SampleLevel(LinearSampler, L, 0).rgb * Factor;
}
Irradiance /= NumSamples;
return float4(Irradiance, 1.0);

至此，我们完成了Irradiancemap的预计算，从上述提供的代码当中可以看出，我们将系数${\pi }^2$或$\pi$都不包含在内，后续在使用的过程要注意这个问题。

5.1.3 One More Thing

在实践过程中发现，如果采用过程中不配合LOD进行的话，会出现噪点。例如对于下图，仅采样LOD0生成IrradianceMap。

而采用了LOD则有更好的效果，效果如下：

结合LOD采样的实现代码如下所示：

float sampleDelta = 1.0 / NumSamplesPerDir;

// 计算LOD级别
uint2 Dimension;
CubeEnvironment.GetDimensions(Dimension.x, Dimension.y);
float lod = max(log2(Dimension.x / float(NumSamplesPerDir)) + 1.0, 0.0);

float NumSamples = 0.0;
for (float phi = 0.0; phi < 2.0 * PI; phi += sampleDelta)
{
    for (float theta = 0.0; theta < 0.5 * PI; theta += sampleDelta)
    {
        // spherical to cartesian (in tangent space)
        float sintheta = sin(theta);
        float costheta = cos(theta);
        float3 tangentSample = float3(sintheta * cos(phi),  sintheta * sin(phi), costheta);
        // tangent space to world
        float3 sampleVec = TangentToWorld(tangentSample, Normal);
        // 采样LOD级别的环境贴图
        float3 sampleColor = CubeEnvironment.SampleLevel(LinearSampler, sampleVec, lod).rgb;

        Irradiance += sampleColor * costheta * sintheta;
        NumSamples++;
    }
}

Irradiance = PI*Irradiance / NumSamples;

5.2 镜面反射环境光

让我们来看一下反射方程的镜面反射部分。

$$\begin{array}{rl}{L}_{Specular}=& {\int }_{\Omega }{f}_{s}Li(w_i)(n\cdot w_i)d{w}_i\\ f_s=& \frac{D(h)F(w_o,w_h)G(w_i,w_o,w_h)}{4(n\cdot w_i)(n\cdot w_o)}\end{array}$$

镜面反射部分，取决于$n$、$w_o$、还有 BRDF 中的粗糙度Roughness、$F_0$（金属度和albedo决定），总共是9个维度。

方向向量用球面坐标表示的话就是 2 个维度。

因此，不能仅存储一张CubeMap就解决问题。

考虑暴力地实时求解方法，采用重要性采样，将上述渲染方程转换为重要性采样之后的样式：

$$pdf=p(w_i)=\frac{D\cdot (n\cdot h)}{4(v\cdot h)}$$

$$L_{Specular}={\int }_{\Omega }{f}_{s}Li(w_i)(n\cdot w_i)d{w}_i\approx \frac{1}{N}\sum _{i=1}^N\frac{f_s{L}_i(n\cdot w_i)}{p(w_i)}$$

其中，N为采样的总样本数量。$p(w_i)$为重要采样的pdf（概率密度函数）。

这种方法虽然可以求解出来，但为了快速计算，不能进行过多的采样。

采样的数量少带来的误差就比较大，因而导致效果较差。

为了加速对该公式的计算，UE4将上述的公式划分为两个不同的求和部分！而这两个求和部分能够分别通过预计算得到结果。

推导的公式如下：

$$\begin{array}{rl}{f}_s=& \frac{D(h)F(w_o,w_h)G(w_i,w_o,w_h)}{4(n\cdot w_i)(n\cdot w_o)}\\ F=& F_0+(1-F_0)(1-(w_o\cdot w_h{)}^5\\ L_{Specular}=& {\int }_{\Omega }{f}_{s}Li(w_i)(n\cdot w_i)d{w}_i\\ =& {\int }_{\Omega }Li(w_i)d{w}_i\ast {\int }_{\Omega }{f}_s(n\cdot w_i)d{w}_i\end{array}$$

可以看出：

• 第一个积分项，仅与环境贴图的纹理颜色有关；被称为pre-filter environment map。

5.2.1 BRDF-LUT

首先，来看第二个积分项！

$${\int }_{\Omega }f(p,w_i,w_o)(w_i\cdot n)d{w}_i$$

公式中 $f$ 包含了 $F_0$，将其进行分离：

$$f(p,w_i,w_0)=\frac{DGF}{4(w_i\cdot n)(w_i\cdot n)}$$

$$F=F_0+(1-F_0)(1-v\cdot h{)}^5$$

公式就变成以下：

$$\begin{array}{r}\int (F_0+(1-F_0)(1-w_o\cdot h{)}^5)\frac{DG}{4(w_i\cdot n)(w_i\cdot n)}(w_i\cdot n)d{w}_i\\ =F_0\ast \int (1-(1-w_o\cdot h{)}^5)\frac{DG}{4(w_o\cdot n)}d{w}_i+\int (1-v\cdot h{)}^5\frac{DG}{4(w_o\cdot n)}d{w}_i\end{array}$$

可以将公式理解成为：

$$F_0\ast Scale+Bias$$

从而将Scale和Bias进行预计算存储。

可以看出这Scale、Bias两项，和$n$、$w_o$、粗糙度Roughness三个变量有关！

不过可以注意到在粗糙度相同的情况下，对于任意的法线$n$，只要观察方向和法线之间的夹角相同，则Scale和Bias就保持不变。

如此一来，我们仅需要关心法线$n$和观察法向$w_o$的夹角即可，这样的好处是根据$w_o$和法线的点积（余弦值）作为一个自变量，同时完全可以根据这个余弦值还原出观察方向$w_o$，此时法线为固定值$(0,0,1)$。则$w_o$在XZ平面上变化：

通过上图的关系，可以求得$w_0$，$x=\sqrt{1-co{s}^2\theta }$，$z=cos\theta$，则$w_0$为$(1-co{s}^2\theta ,0,cos\theta )$。

通过上述这种方式就把三个自变量降维简化为两个了！即可以预计算一张二维纹理了！

• 例如水平方向$NdotV$，竖直方向为$Roughness$。

/**
* ┌---→ NdotV
* |
* ↓
* Roughness
*/

计算则采用蒙特卡洛重要性采样对积分进行近似求解，选取 $pdf=\frac{D(w_i\cdot n)}{4(w_o\cdot w_h)}$。有：

$$\begin{array}{rl}& \int (F_0+(1-F_0)(1-w_o\cdot h{)}^5)\frac{DG}{4(w_i\cdot n)(w_i\cdot n)}(w_i\cdot n)d{w}_i\\ =& F_0\ast \int (1-(1-w_o\cdot h{)}^5)\frac{DG}{4(w_o\cdot n)}d{w}_i+\int (1-v\cdot h{)}^5\frac{DG}{4(w_o\cdot n)}d{w}_i\\ \approx & F_0\ast \frac{1}{N}\sum _{i=1}^N(1-(1-w_o\cdot h{)}^5)\frac{G(w_o\cdot h)}{(w_i\cdot n)(w_o\cdot n)}+\frac{1}{N}\sum _{i=1}^N(1-w_o\cdot h{)}^5\frac{G(w_o\cdot h)}{(w_i\cdot n)(w_o\cdot n)}\end{array}$$

可以将除了$F_0$外的两项分别进行预计算存储。

代码如下：

const float PI = 3.14159265358979323846264338327950288;

float RadicalInverse_VdC(unsigned int bits)
{
    bits = (bits << 16u) | (bits >> 16u);
    bits = ((bits & 0x55555555u) << 1u) | ((bits & 0xAAAAAAAAu) >> 1u);
    bits = ((bits & 0x33333333u) << 2u) | ((bits & 0xCCCCCCCCu) >> 2u);
    bits = ((bits & 0x0F0F0F0Fu) << 4u) | ((bits & 0xF0F0F0F0u) >> 4u);
    bits = ((bits & 0x00FF00FFu) << 8u) | ((bits & 0xFF00FF00u) >> 8u);
    return float(bits) * 2.3283064365386963e-10;
}

Vector2f Hammersley(unsigned int i, unsigned int N)
{
    return Vector2f(float(i) / float(N), RadicalInverse_VdC(i));
}

Vector3f ImportanceSampleGGX(Vector2f Xi, float roughness, Vector3f N)
{
    float a = roughness * roughness;

    float phi = 2.0 * PI * Xi.x;
    float cosTheta = sqrt((1.0 - Xi.y) / (1.0 + (a * a - 1.0) * Xi.y));
    float sinTheta = sqrt(1.0 - cosTheta * cosTheta);

    // 这里采样的球是一个切线空间的球
    // from spherical coordinates to cartesian coordinates
    Vector3f H;
    H.x = cos(phi) * sinTheta;
    H.y = sin(phi) * sinTheta;
    H.z = cosTheta;

    // from tangent-space vector to world-space sample vector
    Vector3f up = abs(N.z) < 0.999 ? Vector3f(0.0, 0.0, 1.0) : Vector3f(1.0, 0.0, 0.0);
    Vector3f tangent = Cross(up, N).Normalize();
    Vector3f bitangent = Cross(N, tangent);

    Vector3f sampleVec = tangent * H.x + bitangent * H.y + N * H.z;
    return sampleVec.Normalize();
}

float GeometrySchlickGGX(float NdotV, float roughness)
{
    float a = roughness;
    float k = (a * a) / 2.0;

    float nom = NdotV;
    float denom = NdotV * (1.0 - k) + k;

    return nom / denom;
}

float GeometrySmith(float roughness, float NoV, float NoL)
{
    float ggx2 = GeometrySchlickGGX(NoV, roughness);
    float ggx1 = GeometrySchlickGGX(NoL, roughness);

    return ggx1 * ggx2;
}

Vector2f IntegrateBRDF(float NdotV, float roughness, unsigned int samples = 128)
{
    Vector2f Eval;

    // View
    Vector3f V;
    V.x = sqrt(1 - NdotV * NdotV);	// sin
    V.y = 0;
    V.z = NdotV;					// cos

    float A = 0;
    float B = 0;

    // Normal
    Vector3f N = Vector3f(0, 0, 1);

    for (unsigned int i = 0; i < samples; ++i)
    {
        // Hammersley产生均匀分布的2D随机采样
        Vector2f Xi = Hammersley(i, samples);
        // GGX重要性采样，采样得到的是半程向量H
        Vector3f H = ImportanceSampleGGX(Xi, roughness, N);
        // H确定了，可以计算出L
        Vector3f L = Vector3f(2.0f * V.Dot(H) * H - V).Normalize();

        float NoL = std::max(L.z, 0.0f);
        float NoH = std::max(H.z, 0.0f);
        float VoH = std::max(V.Dot(H), 0.0f);
        float NoV = std::max(N.Dot(V), 0.0f);

        // 上半球
        if (NoL > 0.0)
        {
            float G = GeometrySmith(roughness, NoV, NoL);
            float G_Vis = (G * VoH) / (NoH * NoV);
            float Fc = pow(1.0 - VoH, 5.0);

            A += (1.0 - Fc) * G_Vis;
            B += Fc * G_Vis;
        }
    }

    Eval.x = A / (1.0f * samples);
    Eval.y = B / (1.0f * samples);
    return Eval;
}

void GenerateBRDFLut()
{
    // 分辨率
    int N = 128;
    DirectX::ScratchImage BRDFLut;
    BRDFLut.Initialize2D(DXGI_FORMAT_R32G32_FLOAT, N, N, 1, 1);
    const DirectX::Image img = BRDFLut.GetImages()[0];
    float* dst = (float*)(img.pixels);
    size_t rowPitch = 0;
    size_t slicePitch = 0;
    DirectX::ComputePitch(img.format, img.width, img.height, rowPitch, slicePitch);

    /**
		* ┌---→ NdotV
		* |
		* ↓
		* Roughness
		*/
    for (int rowIndex = 0; rowIndex < N; ++rowIndex)
    {
        for (int colInddex = 0; colInddex < N; ++colInddex)
        {
            float roughness = (rowIndex + 0.5f) * (1.0f / N);
            float NoV = (colInddex + 0.5f) * (1.0f / N);
            Vector2f eval = IntegrateBRDF(NoV, roughness);

            // write pixel value
            float* dst = (float*)(img.pixels + rowPitch * (rowIndex));
            dst[colInddex * 2 + 0] = eval.x;
            dst[colInddex * 2 + 1] = eval.y;
        }
    }

    // 可压缩成为R16G16_FLOAT
    // BRDFLut.dds
    HRESULT hr = DirectX::SaveToDDSFile(BRDFLut.GetImages(), BRDFLut.GetImageCount(), BRDFLut.GetMetadata(), DirectX::DDS_FLAGS_NONE, L"dx12_brdflut.dds");
    BRDFLut.Release();
}

而虚幻中的LUT图，宽高比并不为1，其分辨率为128x32，其对粗糙度的精度要求不如NdotV，如下所示：

5.2.2 Prefilter EnvironmentMap

再来看第一个积分项：

$${\int }_{\Omega }Li(w_i)d{w}_i$$

这个积分项看似和漫反射环境光的积分类似，可以通过一张CubeMap纹理存储即可。

但其真实的物理含义却告诉你，不可以如此处理。因为在镜面反射中，光照是与视角有关的（$V\cdot R$），并非像漫反射那样结果都一致（$n\cdot l$)。

如果，对于每个视角方向$w_o$和每一个法线方向$n$都进行预计算，那么需要N张纹理，完全不能接受。

这里的优化方式是假设：

• 视角方向$w_o$ = 法线方向$n$ = 反射方向$R$

采样$L_i$仅需要一个$l$方向。

采用中程向量的计算办法也可以很轻松的计算出入射光方向。计算方法如下：

$$l=2\ast (H\cdot w_o)H-w_o$$

其中，这里的$H$和$w_o$都是单位向量。中程向量就是我们需要进行随机生成了。

伪代码如下：

vec3 R = n;
vec3 V = w_o;
for(int i = 0; i < N; i++)
{
    vec3 randNum = randFunction();
    vec3 H = importanceSamlpe(randNum,n);
    vec3 L = normalize(2.0 * dot(V, H) * H - V);
    color += texture(cubeMap,L).xyz;
}
color /= N;

再考虑加权平均，对于不同角度的光线，贡献应该不同。

当NdotL较小时，期待其对最终的结果影响较小，因而赋予较小的权重！

• 入射光线方向越接近法线方向，那么他的光照贡献应该更大！

伪代码如下：

vec3 R = n;
vec3 V = w_o;
for(int i = 0; i < N; i++)
{
    vec3 randNum = randFunction();
    vec3 H = importanceSamlpe(randNum,n);
    vec3 L = normalize(2.0 * dot(V, H) * H - V);
    // dot(n,L)权重
    color += texture(cubeMap,L).xyz * dot(n,L);
    weight += dot(n,L);
}
color /= weight;

同时，当粗糙度不同的时候，CubeMap上影响的着色点光照范围也不同，重要性采样还要加入粗糙度的影响。

最终的伪代码如下：

vec3 R = n;
vec3 V = w_o;
for(int i = 0; i < N; i++)
{
    vec3 randNum = randFunction();
    // 重要性采样加入粗糙度影响 Roughness!
    vec3 H = importanceSamlpe(randNum,n,roughness);
    vec3 L = normalize(2.0 * dot(V, H) * H - V);
    // dot(n,L)权重
    color += texture(cubeMap,L).xyz * dot(n,L);
    weight += dot(n,L);
}
color /= weight;

上述的描述可知，需要给每层Mipmap一个粗糙度值。

在虚幻中是这么计算的：

// Mip，当前CubeMap的Mip等级
// CubemapMaxMip，最大的MipMap等级
// 通过以上公式计算得到一个粗糙度值
float ComputeReflectionCaptureRoughnessFromMip(float Mip, half CubemapMaxMip)
{
	float LevelFrom1x1 = CubemapMaxMip - 1 - Mip;
	return exp2((1.0 - LevelFrom1x1) / 1.2);
}

那么就有对应的公式根据粗糙度计算当前的Mipmap等级：

half ComputeReflectionCaptureMipFromRoughness(float Roughness, half CubemapMaxMip)
{
	half LevelFrom1x1 = 1.0 - 1.2 * log2(max(Roughness, 0.001));
	return CubemapMaxMip - 1 - LevelFrom1x1;
}

虚幻中预过滤环境贴图的着色器代码可参考：ReflectionEnvironmentShaders.usf

void FilterPS(
	FScreenVertexOutput Input,
	out float4 OutColor : SV_Target0
	)
{ 
	float2 ScaledUVs = Input.UV * 2 - 1;
	const int SelectedCubeFace = CubeFace;

    float3 CubeCoordinates = GetCubemapVector(ScaledUVs, SelectedCubeFace);

	float3 N = normalize(CubeCoordinates);
	float3x3 TangentToWorld = GetTangentBasis( N );

    // 计算粗糙度
	float Roughness = ComputeReflectionCaptureRoughnessFromMip( MipIndex, NumMips - 1 );

	if( Roughness < 0.01 )
	{
		OutColor = SourceCubemapTexture.SampleLevel(SourceCubemapSampler, CubeCoordinates, 0 );
		return;
	}

	uint CubeSize = 1 << ( NumMips - 1 );
	const float SolidAngleTexel = 4*PI / ( 6 * CubeSize * CubeSize ) * 2;

	//const uint NumSamples = 1024;
    // 采样数量
	const uint NumSamples = Roughness < 0.1 ? 32 : 64;
	
	float4 FilteredColor = 0;
	BRANCH
	if( Roughness > 0.99 )
	{
		// Roughness=1, GGX is constant. Use cosine distribution instead
		// 当粗糙度为1时，GGX为常数，使用余弦分布代替
		LOOP
		for( uint i = 0; i < NumSamples; i++ )
		{
			float2 E = Hammersley( i, NumSamples, 0 );

			float3 L = CosineSampleHemisphere( E ).xyz;

			float NoL = L.z;

			float PDF = NoL / PI;
			float SolidAngleSample = 1.0 / ( NumSamples * PDF );
			float Mip = 0.5 * log2( SolidAngleSample / SolidAngleTexel );

			L = mul( L, TangentToWorld );
			FilteredColor += SourceCubemapTexture.SampleLevel(SourceCubemapSampler, L, Mip );
		}

		OutColor = FilteredColor / NumSamples;
	}
	else
	{
        
		float Weight = 0;
		LOOP
		for( uint i = 0; i < NumSamples; i++ )
		{
			float2 E = Hammersley( i, NumSamples, 0 );
			// 6x6 Offset rows. Forms uniform star pattern
			//uint2 Index = uint2( i % 6, i / 6 );
			//float2 E = ( Index + 0.5 ) / 5.8;
			//E.x = frac( E.x + (Index.y & 1) * (0.5 / 6.0) );
            
			E.y *= 0.995;
			// GGX重要性采样
			float3 H = ImportanceSampleGGX( E, Pow4(Roughness) ).xyz;
			float3 L = 2 * H.z * H - float3(0,0,1);

			float NoL = L.z;
			float NoH = H.z;

			if( NoL > 0 )
			{
				//float TexelWeight = CubeTexelWeight( L );
				//float SolidAngleTexel = SolidAngleAvgTexel * TexelWeight;

				//float PDF = D_GGX( Pow4(Roughness), NoH ) * NoH / (4 * VoH);
                // 概率密度函数
				float PDF = D_GGX( Pow4(Roughness), NoH ) * 0.25;
				// mipmap等级
                float SolidAngleSample = 1.0 / ( NumSamples * PDF );
				float Mip = 0.5 * log2( SolidAngleSample / SolidAngleTexel );
		
                // 切线空间到世界空间
				L = mul( L, TangentToWorld );
                // 采样并乘以NoL权重
				FilteredColor += SourceCubemapTexture.SampleLevel(SourceCubemapSampler, L, Mip ) * NoL;
				Weight += NoL;
			}
		}

		OutColor = FilteredColor / Weight;
	}
}

参考文献

• 概率论基础

• 看懂蒙特卡洛积分(一) 概率分布变换与随机采样

• 看懂蒙特卡洛积分(二) 蒙特卡洛估计与重要性采样

• 看懂蒙特卡洛积分(三) 低差异采样序列

• 均匀采样问题总结

• 关于Hammersley采样，球面上的均匀采样

• 蒙特卡洛积分与球面上的均匀采样

• 图形学数学基础之Hammersley采样

• GGX重要性采樣

• 彻底看懂PBR/BRDF方程

• 如何看懂这些"该死的"图形学公式

• PBR（基于物理的渲染）学习笔记2

• 【知识补充】基于图像的光照 IBL

• 无偏估计与有偏估计

• 蒙特卡洛积分