基于物理的渲染：微平面理论(Cook-Torrance BRDF推导)

(本篇文章同步发表于知乎专栏：https://zhuanlan.zhihu.com/p/152226698 欢迎三连关注)

在前面两篇关于PBR的笔记中，介绍了辐射度量学的一些物理概念及通过这些物理项从而定义出的BRDF与渲染方程，并且通过蒙特卡洛路径追踪的方式成功解出了渲染方程近似解，得到了相当不错的渲染效果。但是对于其中对物体材质起着至关重要作用的BRDF并没有作过多讲解，本次笔记会从微平面理论出发，去推导出广泛使用的Cook-Torrance BRDF。

1 基础概念

在第一章中，先简单得回顾一下辐射度量学的基础概念，以及BRDF的定义，关于这些概念的详细介绍可参考：
基于物理渲染的基础知识(辐射度量学，BRDF和渲染方程)

1.1 辐射能量与辐射功率

辐射能量(Radiant energy), 指电磁辐射的能量，单位为焦耳($J$)，用符号$Q$表示。
辐射功率(Radiant power)，或者称之为辐射通量(Radiant flux)，指每单位时间辐射而出的能量。
$$\Phi \equiv \frac{\mathrm{d}Q}{\mathrm{d}t}$$

1.2 辐照度与辐出度

辐照度(Irradiance)，指每单位面积，单位时间所接收到的能量。
$$E=\frac{d\Phi }{dA}$$
辐出度(Radiance)，指每单位垂直面积，单位立体角，单位时间接受到的能量。
$$L=\frac{d\Phi }{d\omega d{A}^{\perp }}$$
辐照度与辐出度定义上单位面积的差异如下：

1.3 双向反射分布函数

确定入射方向与出射方向，BRDF定义为出射的Radiance与入射的Irradiance的比值，其衡量了物体表面与光线进行反射的特性，表示出物体的材质效果。
$$f(l,v)=\frac{d{L}_o(v)}{dE(l)}$$
其中$l$是入射光方向，$v$是观察方向，也就是反射光方向。
$d{L}_o(v)$表示物体表面反射到方向$v$的微分辐出度，与$L_o(v)$的区别在于，$L_o(v)$计算了物体表面上来自所有方向的入射光的贡献，而$d{L}_o(v)$仅计算物体表面来自入射方向$l$上的贡献。
与之相对的$dE(l)$也是指来自入射方向$l$的辐照度，而$E(l)$则指物体表面接受到所有方向的辐照度。

2 微平面理论

微平面理论(Microfacets Theory)顾名思义，就是从微观的角度去分析，其认为任意物体的表面在到达微观的尺度之后，都是由许许多多的更小的镜面所组成的：

当这些微小的镜面的法线方向较为集中一致的时候，从宏观角度来看，物体的镜面反射会比较明显，反之，当这些微小镜面的法线方向较为杂乱，分布均匀的时候，物体表面则会相对粗糙。如下图所示：

（上方图中，微平面的镜面法线方向比较集中，下方图中微平面的镜面法线方向比较分散）

但对于此，我们仅仅解释了微观角度的镜面反射现象，另外一个比较重要的部分则是漫反射。

对于微观尺度的漫反射可以作如下解释，当光线照射到物体表面时，一部分光线会与表面的许许多多的微小镜面发生镜面反射，另一部分光线则会折射进入物体内部，如下图：

当光线折射进入内部的时候会与物体的微小粒子不断发生碰撞并散射到随机方向，同时才碰撞的过程中一部分光线的能量会被吸收转换为热能，有些光线在多次碰撞之后能量消耗殆尽，则表示该光线完全被物体吸收。还有一部分折射到物体内部的光线会因为散射方向的随机性重新离开表面，而这部分光线就形成了漫反射。
(tips:光线折射进入物体内部之后再被散射出来，离开表面的点可能并不是进入表面的点，但这里并不考虑这种情形，次表面散射(Subsurface Scattering)技术则将这种现象考虑了进去。)

以上，我们从微观的尺度去理解了物体表面镜面反射与漫反射的现象，Cook-Torrance BRDF正是考虑了这两种情形的BRDF，定义如下：
$$f_r=k_d{f}_{lambert}+k_s{f}_{cook?torrance}$$
其中$k_d,k_s$分别指入射光线中被折射部分的能量所占的比率与被反射部分的比率(一般由菲涅尔项决定，下一章会细讲)，而 $f_{lambert}$指漫反射的BRDF，$f_{cook?torrance}$指镜面反射的BRDF。两部分结合得到最终的BRDF。

3 Cook-Torrance BRDF推导

在第二章中得到了Cook-Torrance BRDF的定义:
$$f_r=k_d{f}_{lambert}+k_s{f}_{cook?torrance}$$
首先让我们先看看如何计算出 $f_{lambert}$。

3.1 漫反射项推导

在很久之前谈光栅化的时候就提及过，漫反射会均匀的向每个方向反射

因此，漫反射的BRDF一定是一个常数。假设入射光是均匀且遍布整个半球方向，可以得到如下方程：

（根据定义$dw=sin\theta d\theta d?$， 积分过程省略）
由于假设入射光是均匀且遍布整个半球方向，所以$L_i$与方向无关，且等于$L_o$，最终约去两边$L$，得到结果：
$$f_r=\frac{1}{\pi }$$
但有一点要注意的是，这里并没有考虑能量被吸收，将反射率考虑进来之后，得到最终的漫反射BRDF：
$$f_{lambert}=\frac{\rho }{\pi }$$
这里的$\rho$可以类比光栅化当中的漫反射系数去理解，本质就是一个3维向量，含有物体表面的颜色信息。

3.2 镜面反射项推导

漫反射项的形式简单且比较容易理解，但是对于镜面反射项，其形式则要更加高级一点：
$$f_{cook?torrance}=\frac{F(l,h)G(l,v)D(h)}{4(n?l)(n?v)}$$
看起来确实有些唬人，但是不要急，让我们先简单的解释一下式子当中各项参数，之和再一步步的进行详细阐述：
其中$v$为反射方向(观察方向)，$l$为入射方向，$n$为宏观表面法向，$h$为微平面法向

分子上的$D,F,G$为3个不同的函数：
函数D：法线分布函数(Normal Distribution Function)，其代表了所有微观角度下微小镜面法线的分布情况，粗糙表面法线分布相对均匀，光滑表面法线分布相对集中 (这种解释可能会有些抽象，后面会给出更加直观的物理上的解释)

函数G：几何函数(Geometry Function)，描述了微平面自遮挡的属性。当一个平面相对比较粗糙的时候，平面表面上的微平面有可能挡住其他的微平面从而减少表面所反射的光线。

函数F：菲涅尔方程(Fresnel Rquation)，描述了物体表面在不同观察角度下反射光线所占的比率

可以看出无论是几何函数G，还是菲涅尔方程F，都是由于观察方向，入射方向的不同所导致的只有部分光线能够被反射
$$\begin{gathered} 0<=G(l,v)<=1 \\ 0<=F(l,h)<=1 \end{gathered}$$
因此先假设没有光线会为上述两个函数产生损失（即二者都=1），方便我们进行BRDF推导。

在上文中，函数$D$即法线分布函数，被解释为微观角度下微小镜面法线的分布函数，其函数接受一个参数为$h$，即微平面的法线向量，更直观的来说，D(h)的物理含义为所有法向为$h$的微平面的面积，没错就是面积，你没有看错，只不过我们要说明的是哪个范围里面的法线方向为$h$的微平面面积，需要加上限制：
(1) 第一，我们所关注的区域仅仅是一个着色点，或者说一个微分面积元，因此$D(h)$所指的应该是单位面积下，所有法向为$h$的微平面的面积。
(2) 第二， 微平面的法线$h$是一个离散数值，而微平面法线分布是连续的，因此是几乎不可能精确的找到法线为某个固定值的微平面的，因此为$D(h)$完善第二个限制，即加上每单位立体角，使其变成一个范围，此时$D(h)$则变成了每单位立体角所有法向为$h$的微平面的面积。
(关于第二点限制可能有些不太好理解，其实这可以类比连续型随机变量取一点概率为0，道理是一样的)

好了，加上这两点限制之后，便得到了$D(h)$的真正的物理含义，即每单位面积,每单位立体角所有法向为$h$的微平面的面积。接下来当确定了观察方向$v$和入射方向$l$之后，我们所真正在乎的能对镜面反射做出贡献的只有法线为$h=(l+v)/||l+v||$的微平面，如下图所示：

只有图中红色的微平面才能对镜面反射做出贡献，利用刚刚得到的关于$D(h)$的物理含义，不难计算得到，所有红色微平面面积$dA\left({\omega }_h\right)=D\left({\omega }_h\right)d{\omega }_{h}dA$
(本文以下用${\omega }_i$代替入射方向$l$，${\omega }_o$代替观察方向$v$，${\omega }_h$代替微平面法向$h$)

那么根据该面积可以计算出入射的辐射通量

$$\begin{array}{rl}d{\Phi }_i& =L_i\left({\omega }_i\right)d{\omega }_{i}d{A}^{\perp }\left({\omega }_h\right)\\ & =L_i\left({\omega }_i\right)d{\omega }_i\cos {\theta }_{h}dA\left({\omega }_h\right)\\ & =L_i\left({\omega }_i\right)d{\omega }_i\cos {\theta }_{h}D\left({\omega }_h\right)d{\omega }_{h}dA\end{array}$$
($d{A}^{\perp }({\omega }_h)$是$dA({\omega }_h)$在入射光线方向的投影垂直面积， ${\theta }_h$为入射光线${\omega }_i$与微平面法线${\omega }_h$的夹角)

因为此时不考虑菲涅尔项以及几何函数所带来的光线损失，所以

$$d{\Phi }_o=d{\Phi }_i$$

进一步，根据辐出度radiance定义，计算得

$$d{L}_o\left({\omega }_o\right)=\frac{d{\Phi }_o}{d{\omega }_o\cos {\theta }_{o}dA}=\frac{L_i\left({\omega }_i\right)d{\omega }_i\cos {\theta }_{h}D\left({\omega }_h\right)d{\omega }_{h}dA}{d{\omega }_o\cos {\theta }_{o}dA}$$

最后由BRDF定易得
$$f_{cook?torrance}\left({\omega }_i,{\omega }_o\right)=\frac{d{L}_o\left({\omega }_o\right)}{d{E}_i\left({\omega }_i\right)}=\frac{d{L}_o\left({\omega }_o\right)}{L_i\left({\omega }_i\right)\cos {\theta }_{i}d{\omega }_i}=\frac{\cos {\theta }_{h}D\left({\omega }_h\right)d{\omega }_h}{\cos {\theta }_o\cos {\theta }_{i}d{\omega }_o}$$
其中，${\theta }_i$是入射光线${\omega }_i$与宏观法线方向$n$的夹角，${\theta }_o$是出射光线${\omega }_o$与宏观法线方向$n$的夹角，${\theta }_h$是入射光线${\omega }_i$与微平面法线方向${\omega }_h$的夹角。

最后一步，只需要找到$d{\omega }_h/d{\omega }_o$的比值，就可以成功推导出镜面反射的BRDF了。那么$d{\omega }_h/d{\omega }_o$的从直观角度去看到底是指什么呢？很简单，如下图所示：

首先固定入射方向${\omega }_i$，微分立体角$d{\omega }_h$是一个很小的范围，那么这样一个很小的范围所对应的反射角自然也是一个很小的范围$d{\omega }_o$，形象的来看，就是转动平面得到所有微观法线在$d{\omega }_h$范围之内的平面，通过这些微平面，所得到的反射范围。那么如何求得这样一个比值呢，以入射方向${\omega }_i$为$z$轴负方向，建立球面坐标系如下：

(嗨呀，画的真丑，凑合看看吧)

根据图中的单位球面坐标系以及标注的角度，不难得出${\omega }_h$的球面坐标为(1，${\theta }_h$，$?$)，则不难计算出微分立体角

$d{\omega }_h=sin{\theta }_{\mathrm{h}}\mathrm{d}{\theta }_{\mathrm{h}}\mathrm{d}?$

同时${\omega }_o$的球面坐标为(1，2${\theta }_h$，$?$)，因为$\theta$角变为两倍，那么对于${\omega }_o$来说$\theta$角度的微分变化应该也为两倍，因此计算得到

$d{\omega }_o=sin2{\theta }_h?2d{\theta }_h?d?=4sin{\theta }_{h}cos{\theta }_{h}d{\theta }_{h}d?$

此时已经知道了$d{\omega }_h$与$d{\omega }_o$，那么

$$\frac{d{\omega }_h}{d{\omega }_o}=\frac{1}{4cos{\theta }_h}$$

将此关系代入上面的镜面BRDF得到：
$$f_{cook?torrance}\left({\omega }_i,{\omega }_o\right)=\frac{D\left({\omega }_h\right)}{4\cos {\theta }_o\cos {\theta }_i}$$

最后将我们一开始忽略掉的菲涅尔项与几何函数重新加入
$$f_{cook?torrance}\left({\omega }_i,{\omega }_o\right)=\frac{F_r\left({\omega }_o,{\omega }_h\right)D\left({\omega }_h\right)G\left({\omega }_i,{\omega }_o\right)}{4\cos {\theta }_o\cos {\theta }_i}$$

至此就已经完成了Cook-Torrance BRDF的所有推导了。

tips:
还有一点有关$D(h)$可以注意的是，对于法线分布函数的设计应该满足
$${\int }_{{\mathcal{H}}^2}\cos \theta D\left({\mathbf{\omega }}_h\right)\mathrm{d}{\mathbf{\omega }}_{h}dA=dA$$
（$\theta$为微平面法向${\omega }_h$与宏观法向$n$的夹角）
从几何角度很容易解释，因为$D({\omega }_h)d{\omega }_{h}dA$着色点(微分面积元)内所有法向为$h$的微平面的面积，乘上$cos\theta$则代表将其投影至宏观平面上，将所有的法线方向${\omega }_h$的投影面积都加起来，自然还是原来的微分面积元$dA$保持不变，同时将上式子进一步化简可得：
$${\int }_{{\mathcal{H}}^2}\cos {\theta }_{h}D\left({\mathbf{\omega }}_h\right)\mathrm{d}{\mathbf{\omega }}_h=1$$

可以验证一下广泛使用的GGX (Trowbridge-Reitz)法线分布函数：
$$\begin{gathered} \alpha ={\text{ Roughness }}^2 \\ D\left({\omega }_h\right)=\frac{{\alpha }^2}{\pi {\left({\left(\mathbf{n}?{\mathbf{\omega }}_h\right)}^2\left({\alpha }^2?1\right)+1\right)}^2} \end{gathered}$$
验证如下：

$$\begin{array}{rl}{\int }_{{\mathcal{H}}^2}\cos \theta D\left({\mathbf{\omega }}_h\right)\mathrm{d}{\mathbf{\omega }}_h& ={\int }_0^{2\pi }\mathrm{d}?{\int }_0^{\frac{\pi }{2}}\frac{{\alpha }^2\cos \theta \sin \theta }{\pi {\left(\left({\alpha }^2?1\right){\cos }^2\theta +1\right)}^2}\mathrm{d}\theta \\ & ={\int }_0^{\frac{\pi }{2}}\frac{2{\alpha }^2\cos \theta \sin \theta }{{\left(\left({\alpha }^2?1\right){\cos }^2\theta +1\right)}^2}\mathrm{d}\theta \\ & =?{\alpha }^2{\int }_0^{\frac{\pi }{2}}\frac{1}{{\left(\left({\alpha }^2?1\right){\cos }^2\theta +1\right)}^2}\mathrm{d}{\cos }^2\theta \\ & ={\alpha }^2{\int }_0^1\frac{1}{{\left(\left({\alpha }^2?1\right)t+1\right)}^2}\mathrm{d}t\\ & ={\frac{{\alpha }^2}{1?{\alpha }^2}\frac{1}{\left({\alpha }^2?1\right)t+1}|}_0^1\\ & =1\end{array}$$

3.3 菲涅尔方程与几何函数的补充

最后再简略的讲一下菲尼尔方程和几何函数对光线损失的影响吧。

3.3.1 菲涅尔方程

菲涅尔方程是为了描述物理世界当中，观察角度与法线夹角越大反射程度一般越大的一种情形，但是精确计算计算消耗较大，一般用Fresnel-Schlick近似法求得近似解
$$F_{\text{Schlick}}\left(h,v,F_0\right)=F_0+\left(1?F_0\right)(1?(h?v){)}^5$$
$F_0$表示平面的基础反射率，它是利用所谓折射系数(Indices of Refraction)计算得出的。越是朝球面掠角的方向上看（此时视线和表面法线的夹角接近90度）菲涅尔现象就越明显，反光就越强。但是根据折射系数只能算出电介质(非导体)的$F_0$为了兼容导体，一般会将$F_0$提前计算出来，再代入方程：

关于菲涅尔方程的详细介绍，我在光线追踪笔记中的4.2，4.3节中由详细提及，这里就不再重复了。
计算机图形学十二：Whitted-Style光线追踪原理详解及实现细节

3.3.2几何函数

几何函数G是为了表示微平面的自遮挡从而引起的光线损失，一般会出现如下两种的遮挡情况

左边一幅图中是入射光线无法照射到一些微平面，这种情况称为Shadowing，右边图中是反射光线无法正常到达人眼，称为Masking，而几何函数$G$正是为了模拟出这两种情况所导致的光线所示，在UE4中采用了Schlick-GGX来进行建模：
$$\begin{gathered} k=\frac{(\text{ Roughness }+1{)}^2}{8} \\ {G}_1(\mathbf{v})=\frac{\mathbf{n}?\mathbf{v}}{(\mathbf{n}?\mathbf{v})(1?k)+k} \\ G(1,\mathbf{v},\mathbf{h})=G_1(1)G_1(\mathbf{v}) \end{gathered}$$
(一般来说表面粗糙程度Roughness越大，遮挡程度越大)

4 总结

一开始只看了闫令琪老师的图形课，对于cook-torrance微表面模型确实了解很少，对PBR理解也不是很深刻，这几天趁着有空，找了不少资料看了看，对Cook-Torrance BRDF大概了解了一些，如果有错误希望各位大佬指正。最后列出一些我参考的资料，推荐第4个视频。

Reference

[1] LearnOpenGL PBR Theory
[2] Background: Physics and Math of Shading by Naty Hoffmann
[3] Real shading in Unreal Engine 4
[4] Cook-Torrance BRDF, UC Davis,这个老教授讲的真不错，推荐
[5] Physically Based Rendering, Second Edition 14.5.1 与 8.4.2节，有部分证明
[6] 知乎文章：离线渲染 | 金属工作流
[7] 知乎文章：基于物理着色：BRDF
[8] GAMES101-现代计算机图形学入门-闫令琪