基于微平面理论的BRDF模型（PBR材质、Cook-Torrance、GGX）

前言

因为前边讲了在Unity中如何实现一个BRDF——简单来说就是把公式抄到Shader里，于是这篇文章则主要从原理角度来讲一讲基于物理的BRDF公式为什么长那个样子。本篇文章主要是整理一下去年（2022）十二月底写的关于基于微平面理论的BRDF的一些东西，主要依据的是两篇经典论文：一篇是Cook-Torrance模型的A Reflectance Model for Computer Graphics[1982]，另一篇则是GGX模型Microfacet Models for Refraction through Rough Surfaces[2007]（这篇论文的作者里也有我的导师…真要说的话本篇文章内容的由来还有一些故事）。

原本计划中应该还有两三篇文章，但这可能是这几天（或者很长一段时间）的最后一篇了，连更几篇确实把我写累了（倒），而且因为一直忙着整理以前的东西，导致我这几天的作业也还没写，导师的项目也几天没动了（其实也是因为还有一些进度存货）。

BRDF

双向反射分布函数（Bidirectional Reflectance Distribution Function，BRDF）描述了物体表面入射光与反射光的关系，具体从定义上来说是出射方向Radiance与入射方向Irradiance的比值，这里也简单补上辐射度量学里的相关定义：

辐射能量(Radiant Energy)：电磁辐射的能量， Q Q Q
辐射通量/功率(Radiant Flux/Power)：单位时间内通过某一区域的能量， Φ = d Q d t \Phi = \frac{dQ}{dt} Φ=dtdQ
辐射强度(Radiant Intensity)：单位立体角内的辐射通量， I ( ω ) = d Φ d ω I(\omega) = \frac{d\Phi}{d\omega} I(ω)=dωdΦ
辐射照度(Irradiance)：单位面积内的辐射通量， E ( x ) = d Φ ( x ) d A E(x) = \frac{d\Phi(x)}{dA} E(x)=dAdΦ(x)
辐射亮度(Radiance)：单位立体角、单位投影面积内的辐射通量， L ( p , ω ) = d 2 Φ ( p , ω ) d ω d A c o s θ L\left( {p,\omega} \right) = \frac{d^{2}\Phi(p,\omega)}{d\omega dAcos\theta} L(p,ω)=dωdAcosθd2Φ(p,ω)

这时候我们再来看渲染方程， L i ( p , ω i ) n ⋅ ω i d ω i L_{i}\left( {p,\omega_{i}} \right)n \cdot \omega_{i}\mathbb{d}\omega_{i} Li(p,ωi)n⋅ωidωi其实算的就是入射方向的Irradiance d E ( ω i ) {d}E(\omega_{i}) dE(ωi)，中间的点乘是因为入射的角度不同接收到的能量就不同。

而 f r ( p , ω i , ω o ) L i ( p , ω i ) n ⋅ ω i d ω i f_{r}\left( {p,\omega_{i},\omega_{o}} \right)L_{i}\left( {p,\omega_{i}} \right)n \cdot \omega_{i}\mathbb{d}\omega_{i} fr(p,ωi,ωo)Li(p,ωi)n⋅ωidωi就是从这个方向入射的Radiance经过物体表面后又有多少被分配至出射方向，于是这么一积分再加上自身自发光的Radiance便是我们要的出射方向的Radiance，这也就是BRDF的作用。

L o ( p , ω o ) = L e ( p , ω o ) + ∫ Ω f r ( p , ω i , ω o ) L i ( p , ω i ) n ⋅ ω i d ω i L_{o}\left( {p,\omega_{o}} \right) = L_{e}\left( {p,\omega_{o}} \right) + \int_{\Omega}^{~}{f_{r}\left( {p,\omega_{i},\omega_{o}} \right)L_{i}\left( {p,\omega_{i}} \right)n \cdot \omega_{i}\mathbb{d}\omega_{i}} Lo(p,ωo)=Le(p,ωo)+∫Ω fr(p,ωi,ωo)Li(p,ωi)n⋅ωidωi

微平面理论与Cook-Torrance

回到BRDF本身，基于物理的BRDF基本理论便是微平面理论，其认为物体的宏观表面由一个个微小的表面组成，光线在微平面上发生理想镜面反射或折射，微平面分布越粗糙，物体高光越模糊：

宏观表面法线方向为n，微观表面法线方向为m

一个经典的基于微平面理论提出的BRDF模型就是Cook-Torrance模型，其BRDF项包含了漫反射 R d R_{d} Rd和镜面反射 R s R_{s} Rs两部分，可以看到公式结构与之前实现里提到的GGX公式基本相同，不同的地方在于 R s R_{s} Rs项Cook-Torrance分母中为 π \pi π而GGX为 4 4 4，之后普遍认为这里取 4 4 4更加合理：
R = s R s + d R d R = sR_{s} + dR_{d} R=sRs+dRd

R s = F π D G ( N ⋅ L ) ( N ⋅ V ) R_{s} = \frac{F}{\pi}\frac{DG}{({\mathbf{N}} \cdot {\mathbf{L}})({\mathbf{N}} \cdot {\mathbf{V}})} Rs=πF(N⋅L)(N⋅V)DG

其中镜面反射描述了从物体表面反射出去的光，而漫反射本质上与次表面散射相同（在Cook-Torrance的论文中，认为漫反射也包含了当表面足够粗糙时光在物体微表面的多次弹射，因为从下文就可以看到，镜面反射中没有考虑到多次弹射的反射现象。而现在似乎有单独对其进行的研究，粗糙度改变时能量应该守恒），都是在进入物体表面后散射出去的光,不同的是相对于渲染像素的距离：如果光线在同一像素内进入并离开则被认为是漫反射，不同像素则认为是次表面散射。

对于漫反射，通常都采用归一化（归一化是为了保证能量守恒，当漫反射强度系数为1时，进入的能量等于出去的能量）的Lambert漫反射模型，模型虽然简单同时也没有考虑菲涅尔反射系数对入射光数量的影响，但效果很好，其公式如下：

R d = 1 π R_{d} = \frac{1}{\pi} Rd=π1

菲涅尔项F

菲涅尔项描述了一个现象，视线与物体表面夹角越小（越接近掠射角）反射越强，如近处的海面可以看到水面以下的东西，而远处的海面则反射出天空的颜色（实质上菲涅尔反射系数是描述了对于一份光，有多少进入了物体表面，又有多少在物体表面被直接反射）。对于同一种材质，在不同波长下，其强度可能不同，这也是为什么金属的高光会有颜色。

通常我们采用一个近似公式来进行描述，其中 F 0 F_{0} F0为入射角为0时的菲涅尔反射率，通常金属较高，而非金属较低。注意宏观表面的菲涅尔反射效应是由微观表面菲涅尔反射产生的，在公式的表现上就是公式描述了微观表面（理想表面）下的菲涅尔反射（取 m \mathbf{m} m的原因），并通过推导成为了宏观BRDF公式的一部分：
F S c h l i c k ( i , m ) = F 0 + ( 1 − F 0 ) ( 1 − ( i ⋅ m ) ) 5 F_{Schlick}\left( {{\mathbf{i}},{\mathbf{m}}} \right) = F_{0} + \left( 1 - F_{0} \right)\left( 1 - ({\mathbf{i}} \cdot {\mathbf{m}}) \right)^{5} FSchlick(i,m)=F0+(1−F0)(1−(i⋅m))5

法线分布函数D

法线分布函数描述了微平面法线 m \mathbf{m} m的统计分布。 D ( m ) D({\mathbf{m}}) D(m)是一个密度函数，其定义如下：以 d ω m \mathbb{d}\omega_{m} dωm表示以m为中心的无限小的立体角， d A \mathbb{d}A dA表示一个无限小的宏观面积，则 D ( m ) d ω m d A D({\mathbf{m}})d\omega_{m}dA D(m)dωmdA表示法线位于立体角内的对应微表面的总面积。下边是一些性质：

0 ≤ D ( m ) ≤ ∞ 0 \leq D({\mathbf{m}}) \leq \infty 0≤D(m)≤∞
1 ≤ ∫ D ( m ) d ω m 1 \leq \int{D({\mathbf{m}})d\omega_{m}} 1≤∫D(m)dωm，微表面面积大于等于对应的宏观表面面积
( v ⋅ n ) = ∫ D ( m ) ( v ⋅ m ) d ω m \left( {{\mathbf{v}} \cdot {\mathbf{n}}} \right) = \int{D({\mathbf{m}})\left( {{\mathbf{v}} \cdot {\mathbf{m}}} \right)d\omega_{m}} (v⋅n)=∫D(m)(v⋅m)dωm，两者投影面积相同
∫ D ( m ) ( n ⋅ m ) d ω m = 1 \int{D({\mathbf{m}})\left( {{\mathbf{n}} \cdot {\mathbf{m}}} \right)d\omega_{m}} = 1 ∫D(m)(n⋅m)dωm=1，即上式中 v \mathbf{v} v与 n \mathbf{n} n重合的特殊情况，也就是归一化，保证能量守恒

投影面积相同

常见如Phong分布、Beckmann分布、GGX（Trowbridge-Reitz）分布，其实都是在尝试用函数来近似拟合物体表面的法线分布情况，其公式分别如下（均做了归一化处理，主要就是为了保证能量守恒，其中Cook-Torranc采用的便是未归一化的Beckmann分布）：

Phong分布： D ( m ) = α + 2 2 π c o s α θ D({\mathbf{m}}) = \frac{\alpha + 2}{2\pi}{cos}^{\alpha}\theta D(m)=2πα+2cosαθ， c o s θ = ( n ⋅ m ) cos\theta = ({\mathbf{n}} \cdot {\mathbf{m}}) cosθ=(n⋅m)
Beckmann分布： D ( m ) = 1 π α 2 c o s 4 θ e − ( t a n θ α ) 2 D({\mathbf{m}}) = \frac{1}{\pi\alpha^{2}{cos}^{4}\theta}e^{- {(\frac{tan\theta}{\alpha})}^{2}} D(m)=πα2cos4θ1e−(αtanθ)2， c o s θ = ( n ⋅ m ) cos\theta = ({\mathbf{n}} \cdot {\mathbf{m}}) cosθ=(n⋅m)
GGX分布： D ( m ) = α 2 π c o s 4 θ ( α 2 + t a n 2 θ ) 2 D({\mathbf{m}}) = \frac{\alpha^{2}}{\pi{cos}^{4}\theta\left( \alpha^{2} + {tan}^{2}\theta \right)^{2}} D(m)=πcos4θ(α2+tan2θ)2α2， c o s θ = ( n ⋅ m ) cos\theta = ({\mathbf{n}} \cdot {\mathbf{m}}) cosθ=(n⋅m)

不同分布函数对比，其中绿色为GGX，明显具有更长的尾部

阴影遮蔽函数G

阴影遮蔽项描述了光线被物体微平面遮挡的情况。阴影指入射光方向的遮挡，遮蔽指出射光方向的遮挡，其具有的一些性质如下：

0 ≤ G ( i , o , m ) ≤ 1 0 \leq G({\mathbf{i}},{\mathbf{o}},{\mathbf{m}}) \leq 1 0≤G(i,o,m)≤1，小于等于1同样是保证能量守恒
G ( i , o , m ) = G ( o , i , m ) G\left( {{\mathbf{i}},{\mathbf{o}},{\mathbf{m}}} \right) = G({\mathbf{o}},{\mathbf{i}},{\mathbf{m}}) G(i,o,m)=G(o,i,m)，交换入射与出射方向遮挡系数不变
G ( i , o , m ) = 0 G\left( {{\mathbf{i}},{\mathbf{o}},{\mathbf{m}}} \right) = 0 G(i,o,m)=0，当 ( i ⋅ m ) ( i ⋅ n ) ≤ 0 ({\mathbf{i}} \cdot {\mathbf{m}})({\mathbf{i}} \cdot {\mathbf{n}}) \leq 0 (i⋅m)(i⋅n)≤0或 ( o ⋅ m ) ( o ⋅ n ) ≤ 0 ({\mathbf{o}} \cdot {\mathbf{m}})({\mathbf{o}} \cdot {\mathbf{n}}) \leq 0 (o⋅m)(o⋅n)≤0，即微表面的背面在宏观表面正方向上不可见

GGX中采用了Smith阴影遮蔽函数，将双向阴影遮蔽 G G G近似为两个单向阴影项 G 1 G_1 G1的可分离乘积：
G ≈ G 1 ( i , m ) G 1 ( o , m ) G \approx G_{1}({\mathbf{i}},{\mathbf{m}})G_{1}({\mathbf{o}},{\mathbf{m}}) G≈G1(i,m)G1(o,m)
G 1 G_1 G1是从法线分布函数D推导得出的，GGX分布推导出Smith G 1 G_1 G1如下：
G 1 ( v , m ) = 2 1 + 1 + α 2 t a n 2 θ ， c o s θ = ( v ⋅ n ) G_{1}\left( {{\mathbf{v}},{\mathbf{m}}} \right) = \frac{2}{1 + \sqrt{1 + \alpha^{2}{tan}^{2}\theta}}，cos\theta = ({\mathbf{v}} \cdot {\mathbf{n}}) G1(v,m)=1+1+α2tan2θ 2，cosθ=(v⋅n)

GGX公式推导

有了上边这些项的定义，就可以利用微平面的理想镜面反射推导宏观粗糙表面BRDF，通过积分对应微面元的贡献计算宏观BRDF（当然，GGX其实提出的是一个完整的BSDF模型，这里只讲其BRDF部分）：
f r ( i , o , n ) = ∫ ∣ i ⋅ m i ⋅ n ∣ f r m ( i , o , m ) ∣ o ⋅ m o ⋅ n ∣ G ( i , o , m ) D ( m ) d ω m f_{r}\left( {{\mathbf{i}},{\mathbf{o}},{\mathbf{n}}} \right) = \int{\left| \frac{{\mathbf{i}} \cdot {\mathbf{m}}}{{\mathbf{i}} \cdot {\mathbf{n}}} \right|f_{r}^{m}\left( {{\mathbf{i}},{\mathbf{o}},{\mathbf{m}}} \right)\left| \frac{{\mathbf{o}} \cdot {\mathbf{m}}}{{\mathbf{o}} \cdot {\mathbf{n}}} \right|G\left( {{\mathbf{i}},{\mathbf{o}},{\mathbf{m}}} \right)D({\mathbf{m}})d\omega_{m}} fr(i,o,n)=∫ i⋅ni⋅m frm(i,o,m) o⋅no⋅m G(i,o,m)D(m)dωm

其中 ∣ i ⋅ m i ⋅ n ∣ \left| \frac{{\mathbf{i}} \cdot {\mathbf{m}}}{{\mathbf{i}} \cdot {\mathbf{n}}} \right| i⋅ni⋅m ， ∣ o ⋅ m o ⋅ n ∣ \left| \frac{{\mathbf{o}} \cdot {\mathbf{m}}}{{\mathbf{o}} \cdot {\mathbf{n}}} \right| o⋅no⋅m 为宏观表面与微面元之间转换的系数，这与之前提到的辐射度量学的相关定义有关，因为微平面的法线与宏观表面的法线可能并不相同。其中：
f r m ( i , o , m ) = F ( i , m ) δ ω m ( h , m ) 4 ( i ⋅ h ) 2 f_{r}^{m}\left( {{\mathbf{i}},{\mathbf{o}},{\mathbf{m}}} \right) = F({\mathbf{i}},{\mathbf{m}})\frac{\delta_{\omega_{m}}\left( {\mathbf{h}},{\mathbf{m}} \right)}{4\left( {\mathbf{i}} \cdot {\mathbf{h}} \right)^{2}} frm(i,o,m)=F(i,m)4(i⋅h)2δωm(h,m)

δ ω m ( h r , m ) \delta_{\omega_{m}}\left( {\mathbf{h}}_{r},{\mathbf{m}} \right) δωm(hr,m)为狄拉克函数，当 h = m {\mathbf{h}} = {\mathbf{m}} h=m时为无穷大，否则为0，也就是只有镜面反射方向不为0，或者说只有 m \mathbf{m} m与 h \mathbf{h} h重合的微平面才会对宏观BRDF产生贡献（也能用来理解为什么 m \mathbf{m} m最后被换成了 h \mathbf{h} h），其性质如下（g为任意函数）：

∫ Ω g ( o ) δ ω o ( s , o ) d ω o = { g ( s ) , i f s ∈ Ω 0 , 其他 \int_{\Omega}^{~}{g({\mathbf{o}})\delta_{\omega_{o}}\left( {{\mathbf{s}},{\mathbf{o}}} \right)d\omega_{o} = \left\{ \begin{matrix}{g({\mathbf{s}}),~if{\mathbf{s}} \in \Omega} \\ {0,其他} \\ \end{matrix} \right.} ∫Ω g(o)δωo(s,o)dωo={g(s), ifs∈Ω0,其他

这个公式实际上是经过一定推导得来的，因为考虑镜面反射，所以 i ⋅ h \mathbf{i} \cdot \mathbf{h} i⋅h和 o ⋅ h \mathbf{o} \cdot \mathbf{h} o⋅h相同，以及狄拉克函数为了进行积分然后换元时产生的雅可比行列式，于是才有了分母中的式子，对详细的数学推导感兴趣的话可以看原论文。这里补充上狄拉克函数的

现在我们将刚刚提到的式子代入上方进行积分，就得到了我们常见的这个BRDF公式：
f r ( i , o , n ) = F ( i , h ) G ( i , o , h ) D ( h ) 4 ∣ i ⋅ n ∣ | o ⋅ n ∣ f_{r}\left( {{\mathbf{i}},{\mathbf{o}},{\mathbf{n}}} \right) = \frac{F\left( {{\mathbf{i}},{\mathbf{h}}} \right)G\left( {{\mathbf{i}},{\mathbf{o}},{\mathbf{h}}} \right)D({\mathbf{h}})}{\left. 4\left| {{\mathbf{i}} \cdot {\mathbf{n}}} \right| \middle| {\mathbf{o}} \cdot {\mathbf{n}} \right|} fr(i,o,n)=4∣i⋅n∣∣o⋅n∣F(i,h)G(i,o,h)D(h)