基于物理的实时渲染着色方法

事实上，从IlumEngine能够显示东西之初，我就已经搭建了一套基于金属工作流的延迟PBR渲染管线，但之前对于PBR材质模型的理解还不够深入，很多情况下都是以抄公式为主，这次借实现Kulla-Conty对微表面模型能量不守恒的修正，从理论上来深入理解实时PBR着色方法。

1. 渲染着色的物理原理

从物理角度上看，渲染着色即光与物体的材质表面、材质介质的相互作用后，射入摄像机的结果。

1.1. 光的物理本质

光是具有一定波长的电磁波，限定自由电荷密度$\rho$和自由电流密度$\pmb J$，根据欧姆定律有电流密度正比于电场：
$$
\pmb J=\sigma \pmb E
$$
对于光在线性介质中的传播，满足麦克斯韦方程组：
$$
\begin{aligned}
\nabla \cdot\pmb E &= \frac{1}{\varepsilon}\pmb \rho\\\\
\nabla\cdot\pmb B&=0\\\\
\nabla\times \pmb E&=-\frac{\partial \pmb B}{\partial t}\\\\
\nabla\times \pmb B&=\mu\sigma\pmb E+ \mu_0\varepsilon_0\frac{\partial \pmb E}{\partial t}
\end{aligned}
$$
由自由电荷的连续性方程：
$$
\nabla\cdot\pmb J=-\frac{\partial \rho}{\partial t}
$$
有：
$$
\frac{\partial \rho}{\partial t}=-\sigma(\nabla\cdot\pmb E)=-\frac{\sigma}{\varepsilon}\rho
$$
可以解出：
$$
\rho(t)=\rho(0)e^{-(\sigma/\varepsilon)t}
$$
看到自由电荷密度$\rho$以特征时间$\tau\equiv \varepsilon/\sigma$耗散，

对于理想导体，$\sigma =\infty$，$\tau=0$
对于良导体，$\tau$比电磁波的角周期小得多：$\tau\ll1/\omega$
对于不良导体，则$\tau$比电磁波的角周期大得多：$\tau\gg 1/\omega$

对电场与磁场的旋度再取其旋度，有：
$$
\begin{align}
\nabla\times(\nabla\times \pmb E)&=\nabla(\nabla\cdot\pmb E)-\nabla^2\pmb E=\nabla\times(-\frac{\partial \pmb B}{\partial t})\\\\
&=-\frac{\partial}{\partial t}(\nabla\times\pmb B)=-\mu\sigma\frac{\partial \pmb E}{\partial t}-\mu_0\varepsilon_0\frac{\partial^2\pmb E}{\partial t^2}\\\\
\nabla\times(\nabla\times \pmb B)&=\nabla(\nabla\cdot\pmb B)-\nabla^2\pmb B=\nabla\times(\mu\sigma\pmb E+\mu_0\varepsilon_0\frac{\partial \pmb E}{\partial t})\\\\
&=\mu\sigma(\nabla\times\pmb E) +\mu_0\varepsilon_0\frac{\partial}{\partial t}(\nabla\times\pmb E)=-\mu\sigma\frac{\partial B}{\partial t}-\mu_0\varepsilon_0\frac{\partial^2\pmb B}{\partial t^2}
\end{align}
$$
忽略电磁波在线性介质中的瞬态效应，即自由电荷消失，$\rho=0$，则有：
$$
\begin{align}
&\begin{matrix}
\nabla \cdot \pmb E=0&\nabla\cdot \pmb B=0
\end{matrix}\\\\
\Rightarrow\ \ \ &
\begin{matrix}
\nabla^2\pmb E=\mu\sigma\dfrac{\partial \pmb E}{\partial t}+\mu_0\varepsilon_0\dfrac{\partial^2\pmb E}{\partial t^2}
\\\
\nabla^2\pmb B=\mu\sigma\dfrac{\partial B}{\partial t}+\mu_0\varepsilon_0\dfrac{\partial^2\pmb B}{\partial t^2}
\end{matrix}
\end{align}
$$
可以得到平面波解：
$$
\begin{aligned}
\tilde {\pmb E}(z,t)=\tilde{\pmb E}_0e^{i(\tilde k z-\omega t)}\\\\
\tilde {\pmb B}(z,t)=\tilde{\pmb B}_0e^{i(\tilde k z-\omega t)}
\end{aligned}
$$
其中，波数$\tilde k$为复数：
$$
\tilde k^2=\mu\varepsilon\omega^2+i\mu\sigma\omega
$$
取平方根得：
$$
\tilde k=k+i\kappa
$$
其中，
$$
\begin{align}
k&\equiv \omega\sqrt{\frac{\mu\varepsilon}{2}}\left[\sqrt{1+\left(\frac{\sigma}{\varepsilon\omega}\right)^2}+1 \right]^{1/2}\\\\
\kappa&\equiv \omega\sqrt{\frac{\mu\varepsilon}{2}}\left[\sqrt{1+\left(\frac{\sigma}{\varepsilon\omega}\right)^2}-1 \right]^{1/2}
\end{align}
$$
提取虚部：
$$
\begin{aligned}
\tilde {\pmb E}(z,t)=\tilde{\pmb E}_0e^{-\kappa z} e^{i(k z-\omega t)}\\\\
\tilde {\pmb B}(z,t)=\tilde{\pmb E}_0e^{-\kappa z} e^{i(k z-\omega t)}
\end{aligned}
$$
可以看到：

$\tilde k$的虚部导致波的衰减，振幅随着$z$的增大而减小，度量了波进入导体的深度，也度量了导体阻尼对波能量的吸收（吸收系数$\alpha\equiv 2\kappa$）
$\tilde k$的实部定义了波的属性，如波长$\lambda = \frac{2\pi}{k}$，波速$v=\frac{\omega}{k}$，折射率$n=\frac{ck}{\omega}$

综合起来，材质通过影响光的$\tilde k$，对光的传播效应造成影响，从而带来不同的光照效果。

在成像理论中，光的吸收对视觉效果有直接影响，将降低光的强度，对特定可见波长的吸收将改变光的颜色。

在渲染中，通常将材质与光的相互作用归结为散射和吸收两类：

散射：决定介质的浑浊程度（下图右）
吸收：决定材质的外观颜色（下图左）

每种介质外观均为散射和吸收两种现象的综合结果：

1.2. 渲染中的光学模型

从电磁波角度出发可以很方便我们理解光、物质和成像之间的相互关系，但在渲染领域尤其是实时渲染领域，我们常常对问题进行简化，现在，我们只考虑物体的表面着色以及几何光学模型。

1.2.1. 光与物体表面的交互

当一种介质进入到另一种介质，光线将发生反射、折射及其各种衍生现象：

镜面反射：光线在两种介质交界处发生直接反射，满足反射定律
折射：从表面进入介质的光由于介质折射率的变化，传播方向发生改变，出现折射现象，介质中的光将发生吸收和散射
吸收：由前推导，当介质中的吸收系数$\alpha>0$时，在光传播过程中部分能量将会被介质所吸收
散射：光的传播方向由于折射率的剧烈变化发生改变，分裂为多个方向，但光的总能量保持不变
- 次表面散射：观察像素小于散射距离
  - 透射：入射光进入介质后又穿过该介质出射
- 漫反射：观察像素大于散射距离

1.2.2. 不同物质与光的交互

不同物质对光的反射、折射效应有所不同，根据光学特性大致分为金属与非金属两大类：

金属：金属的外观取决于材质表面的镜面反射，且入射金属的光不存在散射，会被自由粒子完全吸收
非金属：非金属又称电介质，非金属外观取决于吸收和散射特性：
- 均匀介质：主要为透明介质，无折射率变化，光总以直线传播，不存在散射，但存在吸收
- 非均匀介质：通常可建模为具有嵌入散射粒子的均匀介质，具有折射率的变化
  - 浑浊介质：具有弱散射
  - 半透明介质：具有强散射
  - 不透明介质：与半透明介质类似，具有强散射

2. 光线传输基础

2.1. 反射方程

所有的渲染方法其实就是在干一件事，求解渲染方程：
$$
L_o(\pmb p,\pmb \omega_o)=L_e(\pmb p,\pmb \omega_o)+\int_\Omega f_r(\pmb p,\pmb \omega_i\rightarrow\pmb \omega_o)L_i(\pmb p,\pmb \omega_i)\cdot(\pmb \omega_i\cdot\pmb n)\cdot\mathrm d\pmb \omega_i
$$
这里我们不考虑物体的自发光，因此$L_e$项恒为0，渲染方程退化为反射方程：
$$
L_o(\pmb p,\pmb \omega_o)=\int_\Omega f_r(\pmb p,\pmb \omega_i\rightarrow\pmb \omega_o)L_i(\pmb p,\pmb \omega_i)\cdot(\pmb \omega_i\cdot\pmb n)\cdot\mathrm d\pmb \omega_i
$$
其中$f_r(\pmb p,\pmb \omega_i\rightarrow\pmb \omega_o)$定义了入射光与出射光的反射比例，称为BxDF，通常为双向反射分布函数BRDF，在离线渲染中我们还常用到BSDF(=BTDF+BRDF)等：

2.2. BRDF

2.2.1. 辐射度量学的基本概念

BRDF是定义在辐射度量学上的，简单介绍相关概念：

辐射能

光的能量可以由光子动能定义：
$$
Q=hv
$$
其中$h$为普朗克常量，$v$为光子振动频率，由光的波动性，有：
$$
c=\lambda v
$$
因此光的能量也可表示为：
$$
Q=\frac{hc}{\lambda}
$$
其中$\lambda$为光的波长，$c$为光速，能量的量纲为[$J$]，即Joule

辐射通量

辐射通量是单位时间内发射、接收或传输的能量，定义为：
$$
\Phi=\frac{\mathrm dQ}{\mathrm dt}
$$
功率的量纲为[$lm$]即Lumen，

辐射强度Radiant Intensity

定义为单位立体角的光通量：
$$
I=\frac{\mathrm d\Phi}{\mathrm d\Omega}
$$
辐射强度量纲为$[cd]$即candela

辐照度Irradiance

辐照度定义为单位面积的辐射通量：
$$
E=\frac{\mathrm d\Phi}{\mathrm d A}
$$
辐照度量纲为$[lx]$即lux

辐射通量为$\Phi$的点光源在距离$r$处的球面上的辐照度为：
$$
E_e=\frac{\Phi}{4\pi r^2}
$$
辐亮度Radiance

辐亮度定义为单位面积和单位立体角的辐射通量
$$
L=\frac{\mathrm d^2\Phi}{\mathrm d\Omega\mathrm dA^\perp}=\frac{\mathrm d^2\Phi}{\mathrm d\Omega\mathrm dA\cos\theta}
$$
辐亮度的量纲为$[cd/m^2]$或nits

辐亮度是成像理论中最经常使用的物理量，因为它最符合人和摄像机成像的辐射度量。

2.2.2. 双向反射分布函数

对于入射光$L_i(\pmb p,\pmb \omega_i)$，计算它在半球上的辐照度$E$：
$$
E(\pmb p)=\int_\Omega L_i(\pmb p,\pmb \omega_i)\cos\theta_i\mathrm d\pmb \omega_i
$$
即：
$$
\mathrm dE(\pmb p,\pmb \omega_i)=L_i(\pmb p,\pmb \omega_i)\cos\theta_i\mathrm d\pmb \omega_i
$$
而任意方向上，出射光的辐亮度微分应与入射光辐照度成正比，即：
$$
\mathrm dL_o(\pmb p,\pmb \omega_o)\propto\mathrm d E(\pmb p,\pmb \omega_i)
$$
该比例定义为关于$\pmb\omega_i$和$\pmb \omega_o$的双向反射分布函数即BRDF：
$$
f_r(\pmb p, \pmb \omega_o,\pmb \omega_i)=\frac{\mathrm dL_o(\pmb p,\pmb \omega_o)}{\mathrm d E(\pmb p,\pmb \omega_i)}=\frac{\mathrm dL_o(\pmb p,\pmb \omega_o)}{L_i(\pmb p,\pmb \omega_i)\cos\theta_i\mathrm d\pmb \omega_i}
$$
因此在半球上积分有：
$$
L_o(\pmb p,\pmb \omega_o,\pmb \omega_i)=\int_\Omega f_r(\pmb p, \pmb \omega_o,\pmb \omega_i)L_i(\pmb p,\pmb \omega_i)\cos\theta_i\mathrm d\pmb \omega_i
$$
BRDF的性质：

交换律：$f_r(\pmb p,\pmb \omega_o,\pmb \omega_i)=f_r(\pmb p,\pmb \omega_i,\pmb \omega_o)$
能量守恒：$\int_\Omega f_r(\pmb p,\pmb \omega_o,\pmb \omega_i)\cos\theta\mathrm d\omega\leq 1$

2.2.3. Cook-Torrance BRDF

本文将主要介绍目前主流的Cook-Torrance BRDF模型，Cook-Torrance BRDF由漫反射项和镜面反射项组成：
$$
f_r=k_d\cdot f_{\mathrm{diffuse}} +k_s\cdot f_{\mathrm{specular}}
$$
其中漫反射项使用Lambertian光照模型，镜面反射项使用微表面模型进行建模，随后将具体介绍。

3. Lambertian漫反射材质模型

Lambertian漫反射模型假定出射光线在各个方向上都是均匀的（强度相等），由反射方程：
$$
\begin{aligned}
L_o(\pmb\omega_o)&=\int_\Omega f_rL_i(\pmb\omega_i)\cos\theta_i\mathrm d\pmb \omega_i\\\\
&=f_rL_i\int_\Omega \cos\theta_i\mathrm d\pmb \omega_i\\\\
&=\pi f_rL_i
\end{aligned}
$$
为了保证能量守恒，BRDF选择：
$$
f_r=\frac{\rho}{\pi}
$$
其中$\rho$为材质颜色Albedo参数，通常有RGB三个分量，当$\rho=(1,1,1)$时，说明材质为白色，没有吸收任何色光，此时$L_o=L_i$。

4. 镜面反射模型

4.1. 微表面理论

真实世界中不存在绝对光滑的物体表面，即使是看起来很光滑的表面也具有比光的波长大得多的微尺度不规则性。这种微表面现象将导致每个表面点反射和折射不同方向的光，带来不同的材质外观。

同时，在微表面作用的光还会因为自遮挡带来阴影和遮罩：

微表面BRDF很大程度由材质的粗糙度所影响，粗糙度度量了微表面材质表面的粗糙程度，如下图从左到右粗糙度依次降低：

4.2. 法线分布函数（Normal Distribution Function）

法线分布是描述微表面的一个重要几何性质，通常物体表面会有一个法向量$\pmb n$，可以认为是物体表面绝对光滑时具有的法向量，但由于微表面具有的各种微小几何结构，局部的表面法向量可能与法向量$\pmb n$不同，这种局部法向量称为微表面法向量$\pmb m$，而法线分布函数$D(\pmb m)$则是描述了微表面法线的分布概率密度。

在渲染中，常用宏观表面的半矢量$\pmb h$来表示微表面法向，因为只有$\pmb m=\pmb h$的表面点能够将入射光反射到视线方向，而其他朝向的表面点对最终的渲染结果没有贡献

在整个微表面法线上积分$D(\pmb m)$将得到微表面的面积，不过更常用的是对$D(\pmb m)(\pmb n\cdot\pmb m)$进行积分，将$D(\pmb m)$投影到宏观表面平面上，将得到宏观表面的元面积，规定为1，如下图所示：

又称为投影$D(\pmb m)(\pmb n\cdot\pmb m)$的归一化：
$$
\int_{\pmb m\in\Theta}D(\pmb m)(\pmb n\cdot\pmb m)\mathrm d{\pmb m}=1
$$
上式中的$\Theta$表示在整个球体上进行积分，虽然在图形学中的微结构模型多为高度场，即在以$\pmb n$为中心的半球$\Omega$以外的所有方向有$D(\pmb m)\equiv 0$。

更一般地，微表面和宏观表面在垂直于任何视图方向$\pmb v$的平面上的投影总是相等的：

即：
$$
\int_{\pmb m\in\Theta}D(\pmb m)(\pmb v\cdot\pmb m)\mathrm d{\pmb m}=\pmb v\cdot\pmb n
$$
上式定义了一个合理的法线分布函数应当具有的约束性质，一种直观的理解是一片宏观表面的法向等于其微表面所有法向的均值，至于实际法线分布函数应长什么样，那就有特别多的建模手段得到了，例如，对于绝对光滑的表面，显然$\pmb m$与$\pmb n$恒相等，法线分布函数则可以由一个冲激函数进行表示：
$$
D(\pmb m)=\delta(\pmb m-\pmb n)
$$
注意到式子$\int_{\pmb m\in\Theta}D(\pmb m)(\pmb v\cdot\pmb m)\mathrm d{\pmb m}=\pmb v\cdot\pmb n$中，$\pmb v\cdot\pmb m$的正负不受限制，这是由于投影会产生正负造成抵消，如前图所示，但在实际渲染中，我们常常只关注可见的微表面，即离相机最近的那个微表面，可以通过定义几何函数$G_1(\pmb m,\pmb v)$来给出沿着视图向量$\pmb v$可见的具有法向$\pmb m$的微平面比例，则原式可改写为：
$$
\int_{m\in\Theta}G_1(\pmb m, \pmb v)D(\pmb m)(\pmb v\cdot\pmb m)^+\mathrm d\pmb m=\pmb v\cdot \pmb m
$$
如下图所示，而几何遮蔽函数的具体内容待会将作介绍。

4.2.1. 法线分布函数的基本性质

法线分布函数是一种概率密度函数，很自然地有如下性质：

非负性
$$
0\leq D(\pmb m)\leq \infty
$$
法线分布函数的积分表征了微表面的面积，显然微表面面积应大于等于宏观表面面积
$$
\int_{\pmb m\in \Theta}D(\pmb m)\mathrm d\pmb m\geq 1
$$
微表面和宏观表面在垂直于任何视图方向$\pmb v$的平面上的投影总是相等的
$$
\int_{\pmb m\in\Theta}D(\pmb m)(\pmb v\cdot\pmb m)\mathrm d{\pmb m}=\pmb v\cdot\pmb n
$$
特别地，当视图法向沿着宏观法向时，积分归一化
$$
\int_{\pmb m\in\Theta}D(\pmb m)(\pmb n\cdot\pmb m)\mathrm d{\pmb m}=1
$$

4.2.2. 法线分布函数的形状不变性

形状不变性（shape-invariant）是设计合理法线分布函数的一个重要依据，具有形状不变性的法线分布函数，能够推导该函数归一化的各向异性版本，且可以很方便地推导相应的几何函数项$G$

Heitz定义各向同性的法线分布函数具有形状不变性即材质粗糙度的影响等价于微表面的拉伸，具有形状不变性的微表面法线分布函数可以表示为以下形式：
$$
D(\pmb m)=\frac{\chi^+(\pmb n\cdot\pmb m)}{\alpha^2(\pmb n\cdot\pmb m)^4}g\left(\frac{\sqrt{1-(\pmb n\cdot\pmb m)^2}}{\alpha(\pmb n\cdot\pmb m)}\right)
$$
其中$g(\cdot)$代表一个表示法线分布函数形状的一维函数

4.2.3. 常用的法线分布函数

Blinn-Phong分布

Blinn-Phong分布函数具有如下形式：
$$
D(\pmb m)=\frac{\alpha_p+2}{2\pi}(\pmb n\cdot\pmb m)^{\alpha_p}
$$
其中，幂$\alpha_p$为Blinn-Phong法线分布函数的粗糙度参数，$\alpha_p$越高，表示表面越光滑，$\alpha_p$的取值可以是$[0,\infty)$，$\alpha_p$参数不便于艺术家进行调节，UE4中采用映射$\alpha_p=2\alpha^{-2}-2$，使用材质粗糙度$\alpha$进行控制，得到的Blinn-Phong法线分布函数如下：
$$
D(\pmb m)=\frac{1}{\pi\alpha^2}(\pmb n\cdot\pmb m)^{(\frac{2}{\alpha^2}-2)}
$$
Blinn-Phong分布函数不具有形状不变性。

Beckmann分布

Beckmann分布函数具有如下形式：
$$
D(\pmb m)=\frac{1}{\pi\alpha^2(\pmb n\cdot\pmb m)^4}\exp(\frac{(\pmb n\cdot\pmb m)^2-1}{\alpha^2(\pmb n\cdot\pmb m)^2})
$$

Beckmann分布函数呈高斯函数形态，在对称轴$\pmb n\cdot\pmb m$处取得最值，如下图所示呈现一种中间亮两边暗的渲染效果。

Beckmann分布与Blinn-Phong分布可以通过关系式$\alpha_p=2\alpha_b^{-2}-2$进行等效，如下图所示，蓝色虚线为Blinn-Phong分布、绿色实线为Beckmann分布：

Beckmann分布函数具有形状不变性。

GGX（Trowbridge-Reitz）分布

GGX分布函数具有如下形式：
$$
D(\pmb m)=\frac{\alpha^2}{\pi((\pmb n\cdot \pmb m)^2(\alpha^2-1)+1)^2}
$$
GGX是目前业界最流行的微表面法线分布函数模型，因为它在主流模型中拥有最长的尾部，如下图所示，Beckmann分布函数会很快地逼近于0，而GGX能维持较长的一段尾部：

GGX分布函数具有形状不变性。

GTR（Generalized-Trowbridge-Reitz）分布

广义的Trowbridge-Reitz分布，具有如下形式：
$$
D(\pmb m)=\frac{c}{(1+(\pmb n\cdot\pmb m)^2(\alpha^2-1))^\gamma}
$$
其中，参数$\gamma$用于控制函数的尾部形状，当$\gamma=2$时，则退化为GGX分布。各个$\gamma$取值与GTR分布曲线的关系如下：

GTR分布函数不具有形状不变性。

4.2.4. 各项异性的法线分布函数

我们已知对于一个各向同性的具有形状不变性的法线分布函数，可以用下述形式所表示（考虑定义在正半球上）：

$$
D(\pmb m)=\frac{1}{\alpha^2(\pmb n\cdot\pmb m)^4}g\left(\frac{\sqrt{1-(\pmb n\cdot\pmb m)^2}}{\alpha(\pmb n\cdot\pmb m)}\right)
$$

则对应的各项异性法线分布函数则为：
$$
D(\pmb m)=\frac{1}{\alpha_x\alpha_y(\pmb n\cdot\pmb m)^4}g\left(\frac{
\sqrt{\frac{(\pmb t\cdot\pmb m)^2}{\alpha_x^2}+\frac{(\pmb b\cdot\pmb m)^2}{\alpha_y^2}}}{\pmb n\cdot\pmb m}\right)
$$
其中，$\alpha_x$和$\alpha_y$分别表示沿切线$\pmb t$方向和副法线$\pmb b$方向的粗糙度，当$\alpha_x=\alpha_y$时，则退化为各向同性的形式。

各向异性的Beckmann分布
$$
D(\pmb m)=\frac{1}{\pi\alpha_x\alpha_y(\pmb n\cdot\pmb m)^4}\exp\left(
-\frac{\frac{(\pmb t\cdot\pmb m)^2}{\alpha_x^2}+\frac{(\pmb b\cdot\pmb m)^2}{\alpha_y^2}}{(\pmb n\cdot\pmb m)^2}\right)
$$
各项异性的GGX分布
$$
D(\pmb m)=\frac{1}{\pi\alpha_x\alpha_y}\frac{1}{\left(\frac{(\pmb t\cdot\pmb m)^2}{\alpha_x^2}+\frac{(\pmb b\cdot\pmb m)^2}{\alpha_y^2}+(\pmb n\cdot\pmb m)^2\right)^2}
$$

4.2.5. Beckmann与GGX法线分布函数效果比较

各向同性比较：

各向异性比较（上一行为Beckmann分布，下一行为GGX分布，$\alpha_x$从左到右依次递增，$\alpha_y$保持不变）：

4.3. 几何函数（Geometry Function）

从法线分布函数一节中我们从微表面模型导出了几何函数的基本概念：
$$
\int_{m\in\Theta}G_1(\pmb m, \pmb v)D(\pmb m)(\pmb v\cdot\pmb m)^+\mathrm d\pmb m=\pmb v\cdot \pmb m
$$
几何函数$G$作为一个取值范围为$[0,1]$的标量，描述了微表面的自阴影的属性，表示了具有半矢量的微表面中，同时被入射方向和反射方向可见的比例

几何函数具有两种主要的形式：

$G_1$：微表面在单个方向（光照方向或观察方向）上的可见比例，一般表示阴影函数或遮蔽函数
$G_2$：微表面在光照方向和观察方向上均可见的比例，一般表示联合遮蔽阴影函数

一般默认情况下，几何函数指代$G_2$。

几何函数的设计也有很多，其中基于物理的相对经典的模型为Smith遮蔽函数和V-cavity遮蔽函数

Smith模型具有不相关的表面，即每个微表面与其邻域不相关，与真实世界的连续微表面对比如下（右图为Smith模型）：

V-cavity遮蔽函数计算单独微表面上的散射并混合计算结果：

Heitz证明了Smith遮蔽函数是唯一既遵循公式：
$$
\int_{m\in\Theta}G_1(\pmb m, \pmb v)D(\pmb m)(\pmb v\cdot\pmb m)^+\mathrm d\pmb m=\pmb v\cdot \pmb m
$$
又具有法线遮蔽独立性的便利特性的函数，因此业界更青睐于使用Smith模型，Smith模型也比V-cavity模型能够更好地拟合真实世界的反射现象：

4.3.1. Smith遮蔽函数

下面我们将从各项同性的法线分布函数出发，详细推导Smith几何函数模型。

在前面的介绍中我们已经得到微表面模型的基本方程：
$$
\int_{m\in\Theta}G_1(\pmb m, \pmb v)D(\pmb m)(\pmb v\cdot\pmb m)^+\mathrm d\pmb m=\pmb v\cdot \pmb m
$$
其中的$G_1(\pmb m,\pmb v)$项则是我们要求的几何函数项，这里我们仅考虑以高度场建模的微表面，由于几何函数$G_1$将筛选出在单个方向上（即$\pmb m\cdot\pmb v>0$）可见的光，我们可以将$G_1$项分离为：
$$
G_1(\pmb m,\pmb v)=\chi^+(\pmb m,\pmb v)G_1’(\pmb v)
$$
原方程可化为：
$$
\int_{m\in\Theta}\chi^+(\pmb m,\pmb v)G_1’(\pmb v)D(\pmb m)(\pmb v\cdot\pmb m)\mathrm d\pmb m=\pmb v\cdot \pmb m
$$
提取出与$\pmb m$无关的$G’(\pmb v)$，有：
$$
\pmb v\cdot \pmb m=G_1’(\pmb v)\int_{m\in\Theta}\chi^+(\pmb m,\pmb v)D(\pmb m)(\pmb v\cdot\pmb m)\mathrm d\pmb m
$$
定义与微表面法向$\pmb m(x_m,y_m,z_m)$相关的微表面坡度：
$$
\tilde m(\pmb m)=(x_{\tilde m}, y_{\tilde m})=(-x_m/z_m,-y_m/z_m)
$$
反过来，微表面法向也可表示为：
$$
\pmb m(\tilde m)=(x_m,y_m,z_m)=\frac{1}{\sqrt{x_{\tilde m}^2+ y_{\tilde m}^2+1}}(-x_{\tilde m}, -y_{\tilde m},1)
$$
定义微表面坡度的分布函数$P^{22}$，由$P^{22}$所应满足的归一性：
$$
\int_{-\infty}^\infty\int_{-\infty}^\infty P^{22}(\tilde m)\mathrm d\tilde m=1
$$
以及法线分布函数的归一性（高度场，仅考虑半球）：
$$
\int_{\Omega}D(\pmb m)(\pmb n\cdot\pmb m)\mathrm d{\pmb m}=1
$$
建立关联：
$$
P^{22}(\tilde m)\mathrm d\tilde m=D(\pmb m)(\pmb n\cdot\pmb m)\mathrm d{\pmb m}
$$
整理后积分有：
$$
\int_{\Omega}\chi^+(\pmb m,\pmb v)D(\pmb m)(\pmb v\cdot\pmb m)\mathrm d\pmb m=\int_{-\infty}^\infty\int_{-\infty}^\infty \chi^+(\pmb m,\pmb v)\frac{\pmb v\cdot\pmb m(\tilde m)}{\pmb n\cdot\pmb m(\tilde m)}P^{22}(\tilde m) \mathrm d\tilde m
$$
由于$\pmb n=(0,0,1)$，因此
$$
\pmb n\cdot\pmb m(\tilde m)=\frac{1}{\sqrt{x_{\tilde m}^2+ y_{\tilde m}^2+1}}
$$
与$\chi^+(\pmb m,\pmb v)$项结合起来有：
$$
\chi^+(\pmb m,\pmb v)\pmb v\cdot\pmb m(\tilde m)=\frac{\chi^+(-x_0x_{\tilde m}-y_0y_{\tilde m}+z_0)(-x_0x_{\tilde m}-y_0y_{\tilde m}+z_0)}{\sqrt{x_{\tilde m}^2+ y_{\tilde m}^2+1}}
$$
因此积分可以写为：
$$
\int_{\Omega}\chi^+(\pmb m,\pmb v)D(\pmb m)(\pmb v\cdot\pmb m)\mathrm d\pmb m=
\int_{-\infty}^\infty\int_{-\infty}^\infty \chi^+(-x_0x_{\tilde m}-y_0y_{\tilde m}+z_0)(-x_0x_{\tilde m}-y_0y_{\tilde m}+z_0) P^{22}(x_{\tilde m},y_{\tilde m})\mathrm dx_{\tilde m}\mathrm dy_{\tilde m}
$$
不失一般性地，我们可以假设视角方向与$x$轴对齐，即$\pmb v=(\sin\theta_o,0,\cos\theta_o)$，则积分表示为：
$$
\begin{align}
&\int_{\Omega}\chi^+(\pmb m,\pmb v)D(\pmb m)(\pmb v\cdot\pmb m)\mathrm d\pmb m\\\\
=&\int_{-\infty}^\infty\int_{-\infty}^\infty\chi^+(-\sin\theta_ox_{\tilde m}+\cos\theta_o)(-\sin\theta_ox_{\tilde m}+\cos\theta_o)P^{22}(x_{\tilde m},y_{\tilde m})\mathrm dx_{\tilde m}\mathrm dy_{\tilde m}\\\\
=&\int_{-\infty}^\infty\chi^+(-\sin\theta_ox_{\tilde m}+\cos\theta_o)(-\sin\theta_ox_{\tilde m}+\cos\theta_o)\left(\int_{-\infty}^\infty P^{22}(x_{\tilde m},y_{\tilde m})\mathrm dy_{\tilde m}\right)\mathrm dx_{\tilde m}\\\\
=&\int_{-\infty}^\infty \chi^+(-\sin\theta_ox_{\tilde m}+\cos\theta_o)(-\sin\theta_ox_{\tilde m}+\cos\theta_o)P^2(x_{\tilde m})\mathrm dx_{\tilde m}
\end{align}
$$
其中，$P^2(x_{\tilde m})$为沿着视角法向的一维坡度分布函数，定义为：
$$
P^2(x_{\tilde m})=\int_{-\infty}^{+\infty}P^{22}(x_{\tilde m},y_{\tilde m})\mathrm dy_{\tilde m}
$$
又：
$$
-\sin\theta_ox_{\tilde m}+\cos\theta_o>0\Rightarrow x_{\tilde m}<\cot\theta_o
$$

因此原积分又可化为：
$$
\int_{-\infty}^\infty \chi^+(-\sin\theta_ox_{\tilde m}+\cos\theta_o)(-\sin\theta_ox_{\tilde m}+\cos\theta_o)P^2(x_{\tilde m})\mathrm dx_{\tilde m}=
\int_{-\infty}^{\cot\theta_o} (-\sin\theta_ox_{\tilde m}+\cos\theta_o)P^2(x_{\tilde m})\mathrm dx_{\tilde m}
$$
代入到方程：
$$
\pmb v\cdot \pmb m=G_1’(\pmb v)\int_{m\in\Theta}\chi^+(\pmb m,\pmb v)D(\pmb m)(\pmb v\cdot\pmb m)\mathrm d\pmb m
$$
可以得到：
$$
\cos\theta_o=G_1’(\pmb v)\int_{-\infty}^{\cot\theta_o} (-\sin\theta_ox_{\tilde m}+\cos\theta_o)P^2(x_{\tilde m})\mathrm dx_{\tilde m}
$$
两边同除以$\sin\theta_o$可得：
$$
\cot\theta_o=G_1’(\pmb v)\int_{-\infty}^{\cot\theta_o} (-x_{\tilde m}+\cot\theta_o)P^2(x_{\tilde m})\mathrm dx_{\tilde m}
$$
由于微表面分布是中心对称的，且任意方向的平均坡度为零（正负抵消，分布的均值为0），即
$$
\int_{-\infty}^{+\infty} x_{\tilde m}P^2(x_{\tilde m})\mathrm dx_{\tilde m}=0
$$
引入该项有：
$$
\cot\theta_o=G_1’(\pmb v)\int_{-\infty}^{+\infty}x_{\tilde m}P^2(x_{\tilde m})\mathrm dx_{\tilde m}+ G_1’(\pmb v)\int_{-\infty}^{\cot\theta_o} (-x_{\tilde m}+\cot\theta_o)P^2(x_{\tilde m})\mathrm dx_{\tilde m}
$$
利用$\cot\theta_o=(1-G_1’(\pmb v))\cot\theta_o+G_1’(\pmb v)\cot\theta_o$，代入有：
$$
\begin{align}
(1-G_1’(\pmb v))\cot\theta_o+G_1’(\pmb v)\cot\theta_o
=&G_1’(\pmb v)\int_{-\infty}^{+\infty}x_{\tilde m}P^2(x_{\tilde m})\mathrm dx_{\tilde m}+ G_1’(\pmb v)\int_{-\infty}^{\cot\theta_o} (-x_{\tilde m}+\cot\theta_o)P^2(x_{\tilde m})\mathrm dx_{\tilde m}\\\\
(1-G_1’(\pmb v))\cot\theta_o
=&G_1’(\pmb v)\int_{-\infty}^{+\infty}x_{\tilde m}P^2(x_{\tilde m})\mathrm dx_{\tilde m}+ G_1’(\pmb v)\int_{-\infty}^{\cot\theta_o} (-x_{\tilde m}+\cot\theta_o)P^2(x_{\tilde m})\mathrm dx_{\tilde m}\\\\
&-G_1’(\pmb v)\cot\theta_o
\end{align}
$$
由于概率密度函数的归一性：
$$
\int_{-\infty}^{+\infty} P^2(x_{\tilde m})\mathrm dx_{\tilde m}=1
$$
因此有
$$
G_1’(\pmb v)\cot\theta_o=G_1’(\pmb v) \int_{-\infty}^{+\infty} \cot\theta_oP^2(x_{\tilde m})\mathrm dx_{\tilde m}
$$
代入得：
$$
\begin{align}
(1-G_1’(\pmb v))\cot\theta_o=&G_1’(\pmb v)\int_{-\infty}^{+\infty}x_{\tilde m}P^2(x_{\tilde m})\mathrm dx_{\tilde m}+ G_1’(\pmb v)\int_{-\infty}^{\cot\theta_o} (-x_{\tilde m}+\cot\theta_o)P^2(x_{\tilde m})\mathrm dx_{\tilde m}\\\\
&-G_1’(\pmb v) \int_{-\infty}^{+\infty} \cot\theta_oP^2(x_{\tilde m})\mathrm dx_{\tilde m}\\\\
=&G_1’(\pmb v)\left(\int_{-\infty}^{+\infty}x_{\tilde m}P^2(x_{\tilde m})\mathrm dx_{\tilde m}-\int_{-\infty}^{\cot\theta_o} x_{\tilde m}P^2(x_{\tilde m})\mathrm dx_{\tilde m}\right)\\\\
&+G_1’(\pmb v)\left(\int_{-\infty}^{\cot\theta_o}\cot\theta_oP^2(x_{\tilde m})\mathrm dx_{\tilde m}-\int_{-\infty}^{\infty}\cot\theta_oP^2(x_{\tilde m})\mathrm dx_{\tilde m} \right)\\\\
=&G_1’(\pmb v)\int_{\cot\theta_o}^{+\infty}x_{\tilde m}P^2(x_{\tilde m})\mathrm dx_{\tilde m}-G_1’(\pmb v)\int_{\cot\theta_o}^{+\infty}\cot\theta_oP^2(x_{\tilde m})\mathrm dx_{\tilde m}\\\\
=&G_1’(\pmb v)\int_{\cot\theta_o}^{+\infty}(x_{\tilde m}-\cot\theta_o)P^2(x_{\tilde m})\mathrm dx_{\tilde m}
\end{align}
$$
等式变换一下有：
$$
\frac{1-G’1(\pmb v)}{G’1(\pmb v)}=\frac{1}{\cot\theta_o}\int{\cot\theta_o}^{+\infty}(x{\tilde m}-\cot\theta_o)P^2(x_{\tilde m})\mathrm dx_{\tilde m}
$$
于是得到最终形式：
$$
G_1’(\pmb v)=\frac{1}{1+\Lambda(\pmb v)}
$$
其中，$\Lambda(\cdot)$定义为：
$$
\Lambda(\pmb v)=\frac{1}{\cot\theta_o}\int_{\cot\theta_o}^{+\infty}(x_{\tilde m}-\cot\theta_o)P^2(x_{\tilde m})\mathrm dx_{\tilde m}
$$
这是一个与法线分布函数$D(\pmb m)$相关的函数，每个$D(\pmb m)$都有唯一的$\Lambda(\pmb v)$与之对应。

综合起来可得Smith几何函数$G_1$的具体形式：
$$
G_1(\pmb m, \pmb v)=\frac{\chi^+(\pmb m\cdot\pmb v)}{1+\Lambda(\pmb v)}
$$
具有形状不变性的法线分布函数能够得到$\Lambda(\pmb v)$的解析解，如：

Beckmann分布函数
$$
\begin{align}
P^{22}(x_{\tilde m},y_{\tilde m})&=\frac{1}{\pi\alpha^2}\exp\left(-\frac{x^2_{\tilde m}+y^2_{\tilde m}}{\alpha^2}\right)\\\\
D(\pmb m)&=\frac{1}{\pi\alpha^2(\pmb n\cdot\pmb m)^4}\exp(\frac{(\pmb n\cdot\pmb m)^2-1}{\alpha^2(\pmb n\cdot\pmb m)^2})\\\\
\Lambda(a)&=\frac{\mathrm{erf}(a)-1}{2}+\frac{1}{2a\sqrt{\pi}}\exp(-a^2)
\end{align}
$$
其中，
$$
\begin{align}
a&=\frac{\pmb n\cdot \pmb v}{a\sqrt{1-(\pmb n\cdot\pmb v)^2}}\\\\
\mathrm{erf}(x)&=\frac{2}{\sqrt{\pi}}\int_0^xe^{-x’^2}\mathrm dx’^2
\end{align}
$$
当然上述公式的计算开销很大，Walter等人给出一种近似逼近：
$$
\Lambda(a)=\begin{cases}
\frac{1-1.259a+0.396a^2}{3.535a+2.181a^2},&a<1.6\\\\
0,&a\geq 1.6
\end{cases}
$$
GGX分布函数
$$
\begin{aligned}
P^{22}(x_{\tilde m},y_{\tilde m})&=\frac{1}{\pi\alpha^2\left(1+\frac{x^2_{\tilde m}+y_{\tilde m}^2}{\alpha^2}\right)^2}\\\\
D(\pmb m)&=\frac{\alpha^2}{\pi((\pmb n\cdot \pmb m)^2(\alpha^2-1)+1)^2}\\\\
\Lambda(a)&=\frac{-1+\sqrt{1+\frac{1}{a^2}}}{2}
\end{aligned}
$$
其中，
$$
a=\frac{\pmb n\cdot \pmb v}{a\sqrt{1-(\pmb n\cdot\pmb v)^2}}
$$
同时，Karis也给出了GGX的Smith $G_1$计算近似公式：
$$
G_1(\pmb v)\approx \frac{2(\pmb n\cdot\pmb v)}{(\pmb n\cdot\pmb v)(2-\alpha)+\alpha}
$$
可以很方便地用在实时渲染中。

对于各向异性的法线分布函数，可以证明$\Lambda(\cdot)$与各向同性版本的一致。

4.3.2. Smith联合遮蔽函数-阴影函数

分离的遮蔽阴影函数

最简单和广泛使用的遮蔽阴影函数，认为遮蔽和阴影是独立的，可以分别计算并进行相乘：
$$
\begin{align}
G_2(\pmb v,\pmb l,\pmb m)&=G_1(\pmb v,\pmb m)G_1(\pmb l,\pmb m)\\\\
&=\frac{\chi^+(\pmb v\cdot\pmb m)}{1+\Lambda(\pmb v)}\frac{\chi^+(\pmb l\cdot\pmb m)}{1+\Lambda(\pmb l)}
\end{align}
$$
这种形式不模拟遮蔽和阴影之间的相关性，因此总会多估算阴影，因为一些相关性总是存在的

高度相关的遮蔽阴影函数

这种形式的遮蔽阴影函数更好地模拟了由于微表面高度引起的遮蔽和阴影之间的相关性。直观地说，微平面在微表面升高得越多，对于出射方向未被遮蔽和入射方向未被掩蔽的可见概率会同时增加，形式如下：
$$
G_2(\pmb v,\pmb l,\pmb m)=\frac{\chi^+(\pmb v\cdot\pmb m)\chi^+(\pmb l\cdot\pmb m)}{1+\Lambda(\pmb v)+\Lambda(\pmb l)}
$$
方向相关的遮蔽阴影函数

通过混合可分离的遮蔽阴影函数与两个方向完全相关的情形来表达方向相关：
$$
G_2(\pmb v,\pmb l,\pmb m)=\lambda(\phi)G_1(\pmb v,\pmb m)G_1(\pmb l,\pmb m)+(1-\lambda(\phi))\min(G_1(\pmb v,\pmb m),G_1(\pmb l,\pmb m))
$$
其中$\lambda(\phi)$是一个经验公式，随着角度$\phi$的增大，$\lambda$从0增大到1，Ashikhmin等人建议使用一个高斯函数进行处理：
$$
\lambda(\phi)=1-e^{-7.3\phi^2}
$$
van Ginneken则提出另外一种计算方法：
$$
\lambda(\phi)=\frac{4.41\phi}{4.41\phi+1}
$$
高度方向相关的遮蔽阴影函数

结合高度相关形式和方向相关形式，可以得到：
$$
G_2(\pmb v, \pmb l,\pmb m)=\frac{\chi^+(\pmb v\cdot\pmb m)\chi^+(\pmb l\cdot\pmb m)}{1+\max(\Lambda(\pmb v),\Lambda(\pmb l))+\lambda(\pmb v\cdot\pmb l)\min(\Lambda(\pmb v),\Lambda(\pmb l))}
$$

当出射方向与入射方向平行且$\lambda=0$时，遮蔽和阴影完全相关
相关性随着方向之间的夹角增加而减小，当$\lambda=1$时，遮蔽和阴影不再方向相关，而是退化为高度相关的形式
$\lambda$的计算同样可以通过：
$$
\lambda(\phi)=\frac{4.41\phi}{4.41\phi+1}
$$
$\phi$为$\pmb v$和$\pmb l$之间的方位角

4.4. 菲涅尔反射

我们已经知道光在两介质的交界表面会发生反射和折射现象，在不考虑次表面散射和透射的情况下我们常常只关注反射光（只有反射光对BRDF有贡献），由费马原理（光路最短原理）我们能够从几何上求出入射角和折射角的关系，但无法确定入射光和折射光、反射光的比例关系。而菲涅尔方程则能够描述了光在介质表面传播的能量关系。

4.4.1. 菲涅尔方程的推导

要推导出菲涅尔方程，我们需要重新回到麦克斯韦方程组，假设在均匀线性物质中，不存在自由电荷和电流，则麦克斯韦方程组变为：
$$
\begin{aligned}
\nabla \cdot\pmb E &= 0\\\\
\nabla\cdot\pmb B&=0\\\\
\nabla\times \pmb E&=-\frac{\partial \pmb B}{\partial t}\\\\
\nabla\times \pmb B&=\varepsilon\mu\frac{\partial \pmb E}{\partial t}
\end{aligned}
$$
由前述推导中，我们可以很容易求得：
$$
\begin{matrix}
\nabla^2\pmb E=\mu\varepsilon\dfrac{\partial^2\pmb E}{\partial t^2}&\nabla^2\pmb B=\mu\varepsilon\dfrac{\partial^2\pmb B}{\partial t^2}
\end{matrix}
$$
可以看到$\pmb E$和$\pmb B$的每个直角坐标分量都满足三维波方程的形式：
$$
\nabla^2f=\frac{1}{v^2}\frac{\partial^2 f}{\partial t^2}
$$
波速：
$$
v=\frac{1}{\sqrt{\mu\varepsilon}}
$$
这个波速也为光速，在真空中$v=1/\sqrt{\mu_0\varepsilon_0}=3.00\times 10^8m/s$

定义折射率
$$
n=\sqrt{\frac{\varepsilon\mu}{\varepsilon_0\mu_0}}
$$
则均匀线性物质中的光速为$v=c/n$

能量密度为：
$$
u=\frac{1}{2}\left(\varepsilon E^2+\frac{1}{\mu}B^2 \right)
$$
对于单色平面波，波的强度为：
$$
I=\frac{1}{2}\varepsilon vE_0^2
$$
边界条件

根据麦克斯韦的积分形式，我们可以推导出电磁场在无自由电荷或自由电流的线性介质的边界条件：
$$
\begin{aligned}
\varepsilon_1\pmb E_1^\perp-\varepsilon_2\pmb E_2^\perp&=0\\\\
\pmb E_1^{\parallel}-\pmb E_2^{\parallel}&=0\\\\
\pmb B_1^\perp-\pmb B_2^\perp&=0\\\\
\frac{1}{\mu_1}\pmb B_1^\parallel-\frac{1}{\mu_2}\pmb B_2^\parallel&=0
\end{aligned}
$$
垂直入射的反射与折射

如上图所示，假设$xy$平面为两线性介质的界面，角频率为$\omega$，沿$z$方向传播和沿$x$方向偏振的平面波从左边入射到界面：
$$
\begin{align}
\tilde{\pmb E}_ I(z,t)&=\tilde E_{0I}e^{i(k_1z-\omega t)}\hat{\pmb x}\\\\
\tilde{\pmb B}_ I(z,t)&=\frac{1}{v_I}\tilde E_{0I}e^{i(k_1z-\omega t)}\hat{\pmb y}
\end{align}
$$
它将产生一个反射波：
$$
\begin{align}
\tilde{\pmb E}_ R(z,t)&=\tilde E_{0R}e^{i(-k_1z-\omega t)}\hat{\pmb x}\\\\
\tilde{\pmb B}_ R(z,t)&=-\frac{1}{v_1}\tilde E_{0R}e^{i(-k_1z-\omega t)}\hat{\pmb y}
\end{align}
$$
和一个透射波：
$$
\begin{align}
\tilde{\pmb E}_ T(z,t)&=\tilde E_{0T}e^{i(k_2z-\omega t)}\hat{\pmb x}\\\\
\tilde{\pmb B}_ T(z,t)&=\frac{1}{v_2}\tilde E_{0T}e^{i(k_2z-\omega t)}\hat{\pmb y}
\end{align}
$$
由边界条件可以得到：
$$
\begin{align}
\tilde E_{0I}+\tilde E_{0R}&=\tilde E_{0T}\\\\
\frac{1}{\mu_1}\left(\frac{1}{v_1}\tilde E_{0I}-\frac{1}{v_1}\tilde E_{0R} \right)&=\frac{1}{\mu_2}\left(\frac{1}{v_2}\tilde E_{0T}\right)
\end{align}
$$
解得：
$$
\begin{aligned}
\tilde E_{0R}=\left(\frac{1-\beta}{1+\beta}\right)\tilde E_{0I}\\\\
\tilde E_{0T}=\left(\frac{2}{1+\beta}\right)\tilde E_{0I}
\end{aligned}
$$
其中，
$$
\beta\equiv \frac{\mu_1v_1}{\mu_2v_2}=\frac{\mu_1n_2}{\mu_2n_1}
$$
若介质的磁导率$\mu$与真空中的接近（大多数介质是这样），则$\beta\approx v_1/v_2$，有：
$$
\begin{aligned}
\tilde E_{0R}=\left(\frac{v_2-v_1}{v_2+v_1}\right)\tilde E_{0I}\\\\
\tilde E_{0T}=\left(\frac{2v_2}{v_2+v_1}\right)\tilde E_{0I}
\end{aligned}
$$
实振幅间有以下关系：
$$
\begin{aligned}
E_{0R}=\left|\frac{v_2-v_1}{v_2+v_1}\right|E_{0I}=\left|\frac{n_1-n_2}{n_1+n_2}\right|E_{0I}\\\\
E_{0T}=\left|\frac{2v_2}{v_2+v_1}\right|E_{0I}=\left|\frac{2n_1}{n_1+n_2}\right|E_{0I}
\end{aligned}
$$
从能量角度看，反射波与入射波强度之比为：
$$
R\equiv\frac{I_R}{I_I}=\left(\frac{E_{0R}}{E_{0I}}\right)^2=\left(\frac{n_1-n_2}{n_1+n_2}\right)^2
$$
透射波与入射波强度之比为：
$$
T\equiv\frac{I_T}{I_I}=\frac{\varepsilon_2v_2}{\varepsilon_1v_1}\left(\frac{E_{0T}}{E_{0I}}\right)^2=\frac{4n_1n_2}{(n_1+n_2)^2}
$$
其中，$R$称为反射系数，$T$称为透射系数，注意到，$R+T=1$，也是一种能量守恒。

倾斜入射的反射与折射

如上图所示，考虑更一般的入射情况，入射角度为$\theta_I$，假设一个单色平面波：
$$
\begin{align}
\tilde{\pmb E}_ I(\pmb r,t)&=\tilde E_{0I}e^{i(\pmb k_I\cdot\pmb r-\omega t)}\\\\
\tilde{\pmb B}_ I(\pmb r,t)&=\frac{1}{v_1}(\tilde{\pmb k_I}\times\tilde{\pmb E_I})
\end{align}
$$
从左边入射，得到一个反射波：
$$
\begin{align}
\tilde{\pmb E}_ R(\pmb r,t)&=\tilde {\pmb E}{0R}e^{i(\pmb k_R\cdot\pmb r-\omega t)}\\\\
\tilde{\pmb B} R(\pmb r,t)&=\frac{1}{v_1}(\tilde{\pmb k_R}\times\tilde{\pmb E_R})
\end{align}
$$
和一个透射波：
$$
\begin{align}
\tilde{\pmb E}_ T(\pmb r,t)&=\tilde {\pmb E}{0T}e^{i(\pmb k_T\cdot\pmb r-\omega t)}\\\\
\tilde{\pmb B} T(\pmb r,t)&=\frac{1}{v_2}(\tilde{\pmb k_T}\times\tilde{\pmb E_T})
\end{align}
$$
所有三个波具有相同的频率$\omega$，该物理量由光源所决定，因此三个波的波数之间存在关系：
$$
k_Iv_1=k_Rv_1=k_Tv_2=\omega
$$
由边界条件，介质(1)中的合场强$\tilde{\pmb E}_I+\tilde{\pmb E}_R$和$\tilde{\pmb B}_I+\tilde{\pmb B}_R$必须拟合介质(2)中的合场强$\tilde{\pmb E_T}+\tilde{\pmb B}_T$，因此当$z=0$时，应有：
$$
\pmb k_I\cdot \pmb r=\pmb k_R\cdot\pmb r=\pmb k_T\cdot\pmb r
$$
保证指数项相等，因此，对所有的$x$和$y$有：
$$
x(k_I)_x+y(k_I)_y=
x(k_R)_x+y(k_R)_y=
x(k_T)_x+y(k_T)_y
$$
对$x=0$，有
$$
(k_I)_y=(k_R)_y=(k_T)_y
$$
对$y=0$，有
$$
(k_I)_x=(k_R)_x=(k_T)_x
$$
因此，通过改变坐标轴的方向，我们可以使得$\pmb k_I$、$\pmb k_R$和$\pmb k_T$均在xz平面内，至此，我们可以得出几何光学的三个基本定律：

入射光、反射光和折射光矢量在同一个平面内（称为入射面），入射面法线也在入射面内，同时满足：
$$
k_I\sin\theta_I=k_R\sin\theta_R=k_T\sin\theta_T
$$
式中，$\theta_I$为入射角，$\theta_R$为反射角，$\theta_T$为折射角，都是相对于法线方向考虑
由于$k_1v_1=k_Rv_1=\omega$，因此$k_I=k_R$，于是有入射角等于反射角：
$$
\theta_I=\theta_R
$$
即反射定律
由$k_Iv_1=k_Tv_2=\omega$，可得
$$
k_I=\frac{n_1}{n_2}k_T
$$
因此有
$$
\frac{\sin\theta_T}{\sin\theta_I}=\frac{n_1}{n_2}
$$
即折射定律，又称斯涅尔定律

回到边界条件问题，我们已经可以将指数项进行抵消，则可以得到下述新的边界条件方程：
$$
\begin{aligned}
\varepsilon_1(\tilde {\pmb E}_ {0I}+\tilde{\pmb E}_ {0R})z&=\varepsilon_2(\tilde {\pmb E} {0T})z\\\\
(\tilde {\pmb B}{0I}+\tilde{\pmb B}_ {0R})z&=(\tilde {\pmb B} {0T})z\\\\
(\tilde{\pmb E}{0T}+\tilde {\pmb E}_ {0R}){x,y}&=(\tilde {\pmb E} {0T}){x,y}\\\\
\frac{1}{\mu_1}(\tilde{\pmb B} {0T}+\tilde {\pmb B}_ {0R}){x,y}&=\frac{1}{\mu_2}(\tilde {\pmb B} {0T})_{x,y}
\end{aligned}
$$
其中，$\tilde{\pmb B}_0=(1/v)\hat{\pmb k}\times\tilde{\pmb E}_0$

为了简化问题，假设入射波的电场的偏振方向平行于入射面，即图中的$xz$平面，这样一来反射和透射波的电场的偏振方向也在这个面内，则边界条件方程可表示为：
$$
\begin{align}
\varepsilon_1(-\tilde E_{0I}\sin\theta_I+\tilde E_{0R}\sin\theta_R)&=\varepsilon_2(-\tilde E_{0T}\sin\theta_T)\\\\
\tilde E_{0I}\cos\theta_I+\tilde E_{0R}\cos\theta_R&=\tilde{E}{0T}\cos\theta_T\\\\
\frac{1}{\mu_1v_1}(\tilde E{0I}-\tilde E_{0R})&=\frac{1}{\mu_2v_2}\tilde E_{0T}
\end{align}
$$
由反射定律和折射定律，上式可简化为：
$$
\begin{aligned}
\tilde E_{0I}-\tilde E_{0R}&=\beta\tilde E_{0T}\\\\
\tilde E_{0I}+\tilde E_{0R}&=\alpha\tilde E_{0T}
\end{aligned}
$$
其中，
$$
\begin{align}
\alpha&\equiv\frac{\cos\theta_T}{\cos\theta_I}\\\\
\beta&\equiv\frac{\mu_1v_1}{\mu_2v_2}=\frac{\mu_1n_2}{\mu_2n_1}
\end{align}
$$
解得：
$$
\begin{aligned}
\tilde E_{0R}&=\left(\frac{\alpha-\beta}{\alpha+\beta}\right)\tilde E_{0I}\\\\
\tilde E_{0T}&=\left(\frac{2}{\alpha+\beta}\right)\tilde E_{0I}
\end{aligned}
$$

上式则是大名鼎鼎的菲涅尔方程了，可以看到透射波和反射波的振幅受入射角度影响，因为$\alpha$是关于$\theta_I$的函数：
$$
\alpha=\frac{\sqrt{1-\sin^2\theta_T}}{\cos\theta_I}=\frac{\sqrt{1-[(n_1/n_2)\sin\theta_I]^2}}{\cos\theta_I}
$$
反射率和透射率随入射角度的变化如下图所示：

而且可以看到，当$\alpha=\beta$时，即
$$
\sin^2\theta_B=\frac{1-\beta^2}{(n_1/n_2)^2-\beta^2}
$$
反射波完全消失，$\theta_B$成为布鲁斯特角。而当入射角$\theta_I=90°$时，$\tilde E_{0R}=\tilde E_{0I}$，入射光全部反射，因此当光打在一个球体上，球面掠角方向应会有反光现象：

再从能量角度看，单位入射面上的功率为$\pmb S\cdot\hat{\pmb z}$，入射强度为：
$$
I_I=\frac{1}{2}\varepsilon_1v_1E_{0I}^2\cos\theta_I
$$
反射和折射强度分别为：
$$
\begin{align}
I_R&=\frac{1}{2}\varepsilon_1v_1E_ {0R}^2\cos\theta_R\\\\
I_T&=\frac{1}{2}\varepsilon_2v_2E_ {0T}^2\cos\theta_T
\end{align}
$$
反射系数和透射系数为：
$$
\begin{aligned}
R&\equiv\frac{I_R}{I_I}=\left(\frac{E_{0R}}{E_{0I}}\right)^2=\left(\frac{\alpha-\beta}{\alpha+\beta}\right)^2\\\\
I_T&\equiv\frac{I_T}{I_I}=\frac{\varepsilon_2v_2}{\varepsilon_1v_1}\left(\frac{E_{0T}}{E_{0I}}\right)^2\frac{\cos\theta_T}{\cos\theta_I}=\alpha\beta\left(\frac{2}{\alpha+\beta}\right)^2
\end{aligned}
$$

4.4.2. Schlick近似

在实际渲染中，对菲涅尔项的处理经常采用Schlick近似方法进行逼近：
$$
R(\theta)=R_0+(1-R_0)(1-\cos\theta)^5
$$
其中，
$$
R_0=\left(\frac{n_1-n_2}{n_1+n_2}\right)^2
$$

4.5. 镜面反射项

至此，我们已经详细推导了基于物理的BRDF中镜面反射项的三个重要组成部分：

法线分布函数$D(\pmb m)$
几何函数$G(\pmb m,\pmb v)$
菲涅尔项$F(\pmb v,\pmb m)$

现在来研究如何将这三项整合到我们的PBR着色模型中。

4.5.1. 半矢量

前面已经简单介绍过半矢量$\pmb h$的概念，这里再提一下。如上图所示，只有当入射光通过微表面反射后能够沿着视角方向的反射光才能对BRDF有贡献，因此我们定义使得反射光能够沿视角方向的法向量为半矢量：
$$
\pmb h=\widehat{\pmb v+\pmb l}
$$
在渲染中，我们常用半矢量$\pmb h$代替微表面法向$\pmb m$。

下面推导中，我们采用的符号意义如上图所示

4.5.2. Torrance-Sparrow模型

Torrance-Sparrow模型是微表面BRDF的一个经典模型，其镜面反射项与IlumEngine中所用的Cook-Torrance BRDF几乎一致，这里以Torrance-Sparrow模型为例推导从微表面模型到基于物理的BRDF。

微表面上方向为$\pmb \omega_h$的微元面积为：
$$
\mathrm dA(\pmb \omega_h)=D(\pmb \omega_h)\mathrm d\pmb \omega_h\mathrm dA
$$
因此微元面积上的辐射通量为：
$$
\begin{align}
\mathrm d\Phi_h&=L_i(\pmb \omega_i)\mathrm d\pmb \omega_i\mathrm dA^\perp(\pmb \omega_h)\\\\
&=L_i(\pmb \omega_i)\mathrm d\pmb \omega_i\cos\theta_h\mathrm dA(\pmb \omega_h)\\\\
&=L_i(\pmb \omega_i)\mathrm d\pmb \omega_i\cos\theta_hD(\pmb \omega_h)\mathrm d\pmb \omega_h\mathrm dA
\end{align}
$$
又由菲涅尔定律，出射辐射通量有：
$$
\mathrm d\Phi_o=F(\pmb \omega_o)\mathrm d\Phi_h
$$
又根据辐亮度定义：
$$
L(\pmb \omega_o)=\frac{\mathrm d\Phi_o}{\mathrm d\pmb \omega_o\cos\theta_o\mathrm dA}
$$
联立有：
$$
L(\pmb \omega_o)=\frac{F(\pmb \omega_o)L_i(\pmb \omega_i)\mathrm d\pmb \omega_iD(\pmb \omega_h)\mathrm d\pmb \omega_h\mathrm dA\cos\theta_h}{\mathrm d\pmb \omega_o\mathrm dA\cos\theta_o}
$$
但是对半矢量$\pmb \omega_h$是困难的，我们将其转化为简单地对$\pmb \omega_o$采样，由立体角公式，我们有：
$$
\frac{\mathrm d\pmb \omega_h}{\mathrm d\pmb \omega_i}=\frac{\sin\theta_h\mathrm d\theta_h\mathrm d\phi_h}{\sin\theta_i\mathrm d\theta_i\mathrm d\phi_i}
$$
如下图所示，

由反射定律，显然有$\theta_i=2\theta_h$，又由反射面性质，有$\phi_i=\phi_h$，因此，
$$
\begin{align}
\frac{\mathrm d\pmb \omega_h}{\mathrm d\pmb \omega_i}=\frac{\mathrm d\pmb \omega_h}{\mathrm d\pmb \omega_o}
&=\frac{\sin\theta_h\mathrm d\theta_h\mathrm d\phi_h}{\sin\theta_i\mathrm d\theta_i\mathrm d\phi_i}\\\\
&=\frac{\sin\theta_h\mathrm d\theta_h\mathrm d\phi_h}{2\sin2\theta_h\mathrm d\theta_h\mathrm d\phi_h}\\\\
&=\frac{1}{4\cos\theta_h}\\\\
&=\frac{1}{4(\pmb \omega_i\cdot\pmb \omega_h)}=\frac{1}{4(\pmb \omega_o\cdot\pmb \omega_h)}
\end{align}
$$
代入到辐亮度求算公式，有：
$$
L(\pmb \omega_o)=\frac{F(\pmb \omega_o)L_i(\pmb \omega_i)D(\pmb \omega_h)\mathrm d\pmb \omega_i}{4\cos\theta_o}
$$
进一步考虑几何函数，则BRDF可以写为：
$$
f_r(\pmb \omega_o,\pmb \omega_i)=\frac{F(\pmb \omega_o)G(\pmb \omega_o,\pmb \omega_i)D(\pmb \omega_h)}{4\cos\theta_o\cos\theta_i}=\frac{F(\pmb \omega_o)G(\pmb \omega_o,\pmb \omega_i)D(\pmb \omega_h)}{4(\pmb \omega_o\cdot \pmb n)(\pmb \omega_i\cdot \pmb n)}
$$

至此，我们的基于物理BRDF的镜面反射项推导完成。

5. PBR着色模型

5.1. Cook-Torrance BRDF

综合前述推导，我们得到Cook-Torrance BRDF的最终形式：
$$
f_r(\pmb \omega_o,\pmb \omega_i)=(1-F(\pmb\omega_o))\frac{\rho}{\pi}+\frac{F(\pmb \omega_o)G(\pmb \omega_o,\pmb \omega_i)D(\pmb \omega_h)}{4(\pmb \omega_o\cdot \pmb n)(\pmb \omega_i\cdot \pmb n)}
$$
其中，

$\rho$为Albedo，材质的颜色
$F(\pmb\omega_o)$为菲涅尔系数，选用Schlick近似
$$
\begin{align}
F(\theta)&=F_0+(1-F_0)(1-\cos\theta)^5\\\\
F_0&=\left(\frac{n_1-n_2}{n_1+n_2}\right)^2
\end{align}
$$
$D(\pmb \omega_h)$为法线分布函数，选用GGX分布
$$
D(\pmb \omega_h)=\frac{\alpha^2}{\pi((\pmb n\cdot \pmb \omega_h)^2(\alpha^2-1)+1)^2}
$$
$G(\pmb\omega_o,\pmb\omega_i)$为几何函数，选用Smith分布的Schlick近似
$$
G(\pmb \omega_o,\pmb\omega_i)=\frac{\pmb \omega_o\cdot\pmb\omega_h}{(\pmb \omega_o\cdot\pmb\omega_h)(1-k)+k}
$$
其中，$k=\frac{(\alpha+1)^2}{8}$

5.2. 金属工作流

Cook-Torrance BRDF中的Schlick近似的菲涅尔公式只能用于描述电介质材质（水、玻璃），对于导体（如金属）则需要使用另外一个菲涅尔公式，对材质设计造成不变，金属工作流则使用金属度参数进行插值金属属性，使得菲涅尔公式同时可以用于描述电介质和导体，非金属的$F_0$默认为0.04，而金属的$F_0$则为表面颜色，通过金属度插值$F_0$的GLSL代码如下：

1 2	vec3 F0 = vec3(0.04); F0 = mix(F0, albedo, metallic);

不同金属度和粗糙度的着色结果：

	Roughness = 0.1	Roughness = 0.4	Roughness = 0.7	Roughness = 1.0
Metallic = 0.1
Metallic = 0.4
Metallic = 0.7
Metallic = 1.0

5.3. 多次弹射Kulla-Conty能量补偿

前述的微表面模型仅考虑光线与微表面的单次反射，而不考虑光线多次弹射的结果，这将带来能量损失的问题，具体表现为当材质粗糙度越大时，物体的亮度越低，如下图所示：

Heitz于2016年给出了能量补偿问题的精确解，但由于计算量过大并不适合于实时渲染，Kulla和Conty则给出了一种近似的能量补偿算法。

首先，定义给定出射方向$\mu_o$二维BRDF lobe的总能量，这个能量表示了单次弹射的能量占比：
$$
E(\mu_o)=\int_0^{2\pi}\int_0^1f(\mu_o,\mu_i,\phi)\mu_i\mathrm d\mu_i\mathrm d\phi
$$
在半球上进行积分，其中$\mu=\sin\theta$

Kulla-Conty方法的核心思想是设计另外一个BRDF，满足积分为$1-E(\mu_0)$，则这个能量将作为能量补偿项。

通过经验性的拼凑，凑出了一个满足要求的BRDF函数：
$$
f_{ms}(\mu_o,\mu_i)=\frac{(1-E(\mu_o))(1-E(\mu_i))}{\pi(1-E_{avg})}
$$
其中，
$$
E_{\mathrm{avg}}=2\int_0^1E(\mu)\mu\mathrm d\mu
$$

验证一下$f_{ms}(\mu_o,\mu_i)$是否满足要求：
$$
\begin{align}
E_{ms}(\mu_o)
&=\int_0^{2\pi}\int_0^1f_{ms}(\mu_o,\mu_i,\phi)\mu_i\mathrm d\mu_i\mathrm d\phi\\\\
&=2\pi\int_0^1\frac{(1-E(\mu_o))(1-E(\mu_i))}{\pi(1-E_{\mathrm{avg}})}\mu_i\mathrm d\mu_i\\\\
&=2\frac{1-E(\mu_o)}{1-E_{\mathrm {avg}}}\int_0^1(1-E(\mu_i))\mu_i\mathrm d\mu_i\\\\
&=\frac{1-E(\mu_o)}{1-E_{\mathrm {avg}}}(1-E_{\mathrm {avg}})\\\\
&=1-E(\mu_o)
\end{align}
$$
确实是可以的。

然而$E(\mu)$和$E_{\mathrm {avg}}$都没有解析表达，因此在实际开发中应使用预计算的方法，用贴图预存积分结果。

对于$E(\mu)$，采用一张单通道二维贴图进行存储，从上到下粗糙度依次增加，从左到右$\mu$依次增加

对于$E_{\mathrm{avg}}$，采用一张单通道一维贴图进行存储，为了可视化方便，这里仍采用二维贴图进行展示，从上到下粗糙度依次增加

从上图中可以看到的能量，在粗糙度低时微表面能量损失少，能量多存储结果较大；在粗糙度高时微表面能量损失多，能量少结果偏黑，与我们的认知是一致的。

截至目前，我们仅讨论没有颜色的BRDF的能量补偿情况，对于具有颜色的BRDF，意味着会发生能量的吸收，能量将发生损失，因此我们仍需要计算损失的能量。

定义平均菲涅尔系数（表示了被反射的能量的多少）：
$$
E_{\mathrm{avg}}=\frac{\int_0^1F(\mu)\mu\mathrm d\mu}{\int_0^1\mu\mathrm d\mu}=2\int_0^1F(\mu)\mu\mathrm d\mu
$$
因此带有颜色的反射能量为：$F_{\mathrm{avg}}E_{\mathrm{avg}}$

一次弹射之后：$F_{\mathrm{avg}}(1-E_{\mathrm{avg}})\cdot F_{\mathrm{avg}}E_{\mathrm{avg}}$
……
$k$次弹射之后：$F_{\mathrm{avg}}^k(1-E_{\mathrm{avg}})^k\cdot F_{\mathrm{avg}}E_{\mathrm{avg}}$

等比数列求和得到颜色项：
$$
\sum_{k=0}^\infty F_{\mathrm{avg}}^k(1-E_{\mathrm{avg}})^k\cdot F_{\mathrm{avg}}E_{\mathrm{avg}}=\frac{F_{\mathrm{avg}}E_{\mathrm{avg}}}{1-F_{\mathrm{avg}}(1-E_{\mathrm{avg}})}
$$
其中$F_{\mathrm{avg}}$也没有很好求，不过有一种近似的计算方法：
$$
\begin{aligned}
F_{\mathrm{avg}}(r,g)&\approx 0.087237 + 0.0230685g - 0.0864902g^2 + 0.0774594g^3\\\\
&+ 0.782654r - 0.136432r^2 + 0.278708r^3\\\\
&+ 0.19744gr + 0.0360605g^2r - 0.2586gr^2;
\end{aligned}
$$
其中$r$为反射率，$g$为边缘色调，通常取$r$为Albedo颜色，$g=(0.827, 0.792, 0.678)$即可

下图展示了Kulla-Conty能量补偿的结果，其中上面一行为进行能量补偿的结果，下一行为未进行能量补偿的结果