计算机图形学（二）——微表面模型

计算机图形学中基于物理建模的渲染技术之所以能给人极佳的视觉体验，是因为利用这些渲染技术能够很真实的反映出每种物体独有的“质感”。我们能通过人眼观察来感受物体表面“质感”的原因，也是因为物体表面反射周围环境的特性不同而造成的，因此对物体表面的物理建模对于其表面本身的质感表现至关重要。对物体表面的建模，最简单的是镜面模型。利用镜面模型渲染出的物体具有十分光滑的感觉。然而现实生活中很多物体表面一般是粗糙的，因此为了对上述的一般表面进行物理建模并将该物理模型应用到实际的渲染技术当中去，我们首先要引入微表面模型。微表面模型假定物体表面是由一系列朝向各异的微小镜面组成的，各个微小镜面的朝向分布越分散，则宏观的物体表面就会越粗糙；若越集中，则物体表面就会越光滑。在这一部分中，我们首先详细介绍法线分布函数(Microfacet Distribution Function)和自遮挡函数(Masking and Shadowing Function)，随后基于上述两部分内容，详细地推导Torrance-Sparrow微表面模型的BRDF解析形式。

 

2.1 法线分布函数(Microfacet Distribution Function)

 

 法线分布函数类似概率论中的概率密度分布函数，它描述了每一个微表面法线朝向一个特定方向的概率。该函数写为$D\left({{\bf{\omega}}_h}\right)$，这里的${\bf{\omega}}_h$是每一个微小面元(microfacet)的朝向，即法线方向。如图4所示，我们假设微表面(蓝色部分)是原物体光滑表面(黑色部分)沿法线${\bf{n}}$做随机扰动后所形成的不平整表面，其中的红色部分是朝向参考方向${\bf{\omega}}_h$附近的微表面集合，我们设扰动前光滑表面(黑色部分)的表面积总和为$A$，经扰动之后的微表面(蓝色部分)的表面积总和为$A'$，并假定$A$的面积足够小。于是$D\left( {\bf{\omega}}_h \right){\rm{d}}{\bf{\omega}}_h $被定义为：图4中所有朝向位于${\bf{\omega}}_h$附近的红色微小面元面积之和占$A$的比例(注意是占$A$的比例，而不是占$A'$的比例。这个定义似乎有些反常，从直觉上来看，我们求红色面积对$A'$的比例更有意义，但由于微表面的总表面积$A'$难于求得，且$A'$和物体表面的粗糙度有关，粗糙度越大则相应的$A'$越大，但我们发现$A$却与物体表面的粗糙度无关，这也是不得已而做出的权衡)。因此，表达式${A}{\cdot}D({\bf{\omega}}_h){\rm{d}}{{\bf{\omega}}_h}$即为朝向${\bf{\omega}}_h$附近的所有微表面的面积总和。由于在物体表面任意位置$x$处形成的半球面$\Omega_x$包含了所有$\bf{\omega}_{h}$的可能取值，根据上述的叙述过程，所以该处微表面的表面积总和$A' = \int\limits_{\Omega _x } {A} \cdot {D\left( {{\bf{\omega}} }_h \right){\rm{d}}{{\bf{\omega}} }_h}$。

 

 

根据以上的陈述，我们可以得到有关函数$D({\bf{\omega}}_h)$的一些性质如下：

1） 非负性。

\begin{equation}
D({\bf{\omega}}_h) \ge 0
\end{equation}

2） 在物体表面任意位置$x$处形成的半球面$\Omega_x$内，微表面的面积$A'$应不小于$A$。即$A' = \int\limits_{\Omega _x } {A} \cdot {D\left( {{\bf{\omega}} }_h \right){\rm{d}}{{\bf{\omega}} }_h} \ge A $。约去两边的$A$，有

\begin{equation}
\int\limits_{\Omega _x } {D\left( {{\bf{\omega}} }_h \right){\rm{d}}{{\bf{\omega}} }_h } \ge 1
\end{equation}
这样看来，$D(\bf{\omega}_{h})$并不能算是严格意义上的概率密度函数，因为原始的$D(\bf{\omega}_{h})$在半球面域上的积分不一定为$1$。这是因为前文所述$D({\bf{\omega}}_{h}) {\rm{d}} {\bf{\omega}_{h}}$的含义较为奇怪所造成的。

3） 原光滑表面在任意方向$\bf{v}$上的投影面积，也等于微表面在该方向上的投影面积。
这一点可以通过图5说明。注意朝向与投影方向相反的微表面，其投影面积按负值计算，体现为图中的红色部分。设每个微表面的法线方向为$\bf{h}$，原表面法线为$\bf{n}$，则有$\int\limits_{\Omega _x } {A \cdot D\left( {{\bf{\omega }}_h } \right)\left( {{\bf{v}} \cdot {\bf{h}}} \right){\rm{d}}{\bf{\omega }}_h } = A\left( {{\bf{v}} \cdot {\bf{n}}} \right)$。即
\begin{equation}
\label{eq:microfacet-projection}
\int\limits_{\Omega _x } { D\left( {{\bf{\omega }}_h } \right)\left( {{\bf{v}} \cdot {\bf{h}}} \right){\rm{d}}{\bf{\omega }}_h } = \left( {{\bf{v}} \cdot {\bf{n}}} \right)
\end{equation}
特殊地，若取$\bf{v}=\bf{n}$，则式(\ref{eq:microfacet-projection})变为
\begin{equation}
\label{eq:microfacet-projection-2}
\int\limits_{\Omega _x } { D\left( {{\bf{\omega }}_h } \right) {\rm{cos} {\theta}} {\rm{d}}{\bf{\omega }}_h } = 1
\end{equation}
式(\ref{eq:microfacet-projection-2})中$\theta$是每个微表面的法线$\bf{h}$与原表面法线$\bf{n}$的夹角。该式实质上定义了法线分布函数$D(\bf{\omega}_h)$的归一化条件。

 

以上内容详细介绍了有关法线分布函数的特性，而针对法线分布函数的具体表现形式，有如下几种典型的实现。接下来主要介绍Beckmann–Spizzichino分布函数(后文简称Beckmann分布)和Walter GGX分布函数(后文简称GGX分布)。

 

2.1.1 Beckmann分布

Beckmann和Spizzichino两人在1963年提出了$D({\bf{\omega}}_{h})$的一种实现形式如下
\begin{equation}
\label{eq:beckmann}
D\left( {{\bf{\omega }}_h } \right) = \frac{{e^{ - \tan ^2 \theta _h /\alpha ^2 } }}{{\pi \alpha ^2 \cos ^4 \theta _h }}
\end{equation}
其中$\theta_h$是${\bf{\omega}}_{h}$在球坐标系表示下的极角。

式(\ref{eq:beckmann})可以用于各向同性表面，针对各向异性的表面，可以将$\alpha$分解为$\alpha_x$和$\alpha_y$分量，并模仿椭圆解析式的形式对该式进行改进，得到如下所示的各向异性表达式

\begin{equation}
\label{eq:beckmann-aniso}
D\left( {{\bf{\omega }}_h } \right) = \frac{{e^{ - \tan ^2 \theta _h \left( {\cos ^2 \varphi _h /\alpha _x^2 + \sin ^2 \varphi _h /\alpha _y^2 } \right)} }}{{\pi \alpha _x \alpha _y \cos ^4 \theta _h }}
\end{equation}
其中$\theta_h$、$\varphi_h$分别是${\bf{\omega}}_{h}$在球坐标系表示下的极角和方位角。各向同性表达式仅是各向异性在$\alpha_x=\alpha_y=\alpha$的情况，其分布特点如图6所示，每幅子图沿径向外移对应球坐标系的极角$\theta$增大，球坐标系的方位角$\varphi$变化一周，对应每幅子图等半径的同心圆。

 

 

2.1.2 GGX分布
Walter在2007年基于Trowbridge, T. S.和Reitz, K. P.两人提出的分布函数，进一步推导出了GGX分布。利用这种分布构造的$D({\bf{\omega}} _ {h})$，具有较长的拖尾，可以很好的模拟普通物体的高光特性。这里直接给出其各向异性的$D({\bf{\omega}} _ {h})$形式如下。各向同性形式仅需令$\alpha_x=\alpha_y$即可。其分布特点如图7所示。

\begin{equation}
\label{eq:trowbridge-aniso}
D\left( { {\bf{\omega}} _h } \right) = \frac{1}{{\pi \alpha _x \alpha _y \cos ^4 \theta _h \left( {1 + \tan ^2 \theta _h \left( {\frac{{\cos ^2 \varphi _h }}{{\alpha _x^2 }} + \frac{{\sin ^2 \varphi _h }}{{\alpha _y^2 }}} \right)} \right)^2 }}
\end{equation}

 

2.2 自遮挡函数(Masking and Shadowing Function)

仅定义法线分布函数，是不足以描述微表面本身特征的，我们还要再引入自遮挡函数(Masking and Shadowing Function)，这种函数意在模拟随光照方向发生变化时，微表面自身的遮挡特性。或者说，当观察视角发生变化时，某些微表面会被其它微表面遮挡，我们需要一个函数来阐述这种遮挡关系。这种函数的原型为$G_1( { \bf{\omega},\bf{\omega}_{h} } )$，其含义为：沿视线方向$\bf{\omega}$看去，朝向$\bf{\omega}_{h}$附近的微表面中能被人眼看到的比例(图8中蓝色部分与绿色部分之和的比例)。图形学中为了简化模型，常假设$G_1$与微表面自身朝向$\bf{\omega}_{h}$无关，此时$G_1$常写为$G_1( \bf{\omega} )$。

与式(\ref{eq:microfacet-projection-2})中描述的法线分布函数的归一化约束条件类似，这里我们开始推导自遮挡函数$G_1( { \bf{\omega},\bf{\omega}_{h} } )$的约束条件。在介绍法线分布函数部分中我们曾提及，朝向$\bf{\omega}_{h}$附近的微表面的总表面积为${A \cdot D\left( {{\bf{ \bf{\omega} }}_h } \right){\rm{d}}{\bf{ \bf{\omega} }}_h }$。如果我们沿着视线方向$\bf{\omega}$看去，则朝向$\bf{\omega}_{h}$附近的微表面沿着视线方向的投影面积为
$$
{\max \left( {0,{\bf{\omega }} \cdot {\bf{\omega }}_h } \right) \cdot A \cdot D\left( {{\bf{\omega }}_h } \right){\rm{d}}{\bf{\omega }}_h }
$$
上式用$\max$来过滤掉背向视线的微表面投影面积，使得当微表面朝向${\bf{\omega}_{h}}$与视线${\bf{\omega}}$成角大于$\frac{\pi}{2}$时，它们的投影(负值)不会对最终投影结果产生抵消作用(因为我们从任意方向观察时本身就看不到背离视线方向的微表面，它们对可见面积没有贡献)。注意这里即便是我们得到了所有朝向${\bf{\omega}_{h}}$的微表面投影，它们自身也可能发生遮挡(图8中沿视线方向$\bf{\omega}$看去，绿色部分之间有相互重叠的部分)，使得它们的实际可见面积小于投影面积，而$G_1$函数可以得到其可见面积与实际面积之比。根据前文对$G_1$函数的描述，我们可以得到实际可见的朝向${\bf{\omega}_{h}}$附近的微表面投影面积为
$$
{ G_1({\bf{\omega}},{\bf{\omega}}_{h} ) \cdot \max \left( {0,{\bf{\omega }} \cdot {\bf{\omega }}_h } \right) \cdot A \cdot D\left( {{\bf{\omega }}_h } \right){\rm{d}}{\bf{\omega }}_h }
$$
上述面积对位于点$x$处附近半球面$\Omega_x$内所有$\bf{\omega}_{h}$进行积分，我们就得到了所有微表面沿视线$\bf{\omega}$投影后的可见面积为
$$
\int\limits_{\Omega _x } {G_1 \left( {\bf{\omega} ,\bf{\omega} _h } \right) \cdot \max \left( {0,\bf{\omega} \cdot \bf{\omega} _h } \right) \cdot A \cdot D\left( {\bf{\omega} _h } \right){\rm{d}}\bf{\omega} _h }
$$
而这个面积又一定与$A$沿着视线方向的投影面积严格相等，即
$$
A \cos \theta = \int\limits_{\Omega _x } {G_1 \left( {\bf{\omega} ,\bf{\omega} _h } \right) \cdot \max \left( {0,\bf{\omega} \cdot \bf{\omega} _h } \right) \cdot A \cdot D\left( {\bf{\omega} _h } \right){\rm{d}}\bf{\omega} _h }
$$

 

由于$A$为常数，两边约去重复项$A$，得到
\begin{equation}
\label{eq:g1-constraint}
\cos \theta = \int\limits_{\Omega _x } {G_1 \left( {\bf{\omega} ,\bf{\omega} _h } \right) \cdot \max \left( {0,\bf{\omega} \cdot \bf{\omega} _h } \right) \cdot D\left( {\bf{\omega} _h } \right){\rm{d}}\bf{\omega} _h }
\end{equation}
式(\ref{eq:g1-constraint})即为$G_1$和$D$的共同约束条件，当$D$的形式给定且已满足式(\ref{eq:microfacet-projection-2})的归一化条件时，则进一步化为$G_1$自身的约束条件。此条件可以用于后续对$G_1$解析形式的求解。

 

基于上述讨论，我们只定性分析了$G_1$的形式与作用，这里我们继续推导$G_1$的解析形式。为了简化$G_1$的形式，我们假定假设$G_1$与微表面自身朝向$\bf{\omega}_{h}$无关。此时简化后的$G_1$如图9所示。注意图9与图8的区别，图8这里由于已经假设$G_1$仅与${\bf{\omega}}$有关，所以在投影的时候，我们把所有与视线成角小于$\frac{\pi}{2}$的微表面一并投影，得到的绿色面积之和设为$A^+({\bf{\omega}})$，把所有与视线成角大于$\frac{\pi}{2}$的微表面一并投影，得到的红色面积之和设为$A^-({\bf{\omega}})$。这样我们可以轻松得到
$$
A \cos \theta = A^+({\bf{\omega}}) - A^-({\bf{\omega}})
$$
则
\begin{equation}
\label{eq:g1-2}
G_1({\bf{\omega}})= \frac{A \cos \theta}{A^+({\bf{\omega}})} = \frac{A^+({\bf{\omega}}) - A^-({\bf{\omega}})}{A^+({\bf{\omega}})}
\end{equation}
此时构造辅助函数
$$
\Lambda \left( {\bf{\omega}} \right) = \frac{{A^ - \left( {\bf{\omega}} \right)}}{{A^ + \left( {\bf{\omega}} \right) - A^ - \left( {\bf{\omega}} \right)}}
$$
则
$$
G_1 \left( {\bf{\omega}} \right) = \frac{1}{{1 + \Lambda \left( {\bf{\omega}} \right)}}
$$
如果我们假设，微表面之间的高度互不相关(很明显这个假设并不正确，因为微表面是连续的，相邻微表面的高度一定相关，这里做出假设是因为要根据$D({\bf{\omega}}_{h})$得出$\Lambda$的确解，否则将会有很多$\Lambda$的解满足式(\ref{eq:g1-constraint})。实际应用表明，这种假设下得到的$\Lambda$仍能对生活中的物体表面做出很好的近似)，则对于Beckmann分布，有
\begin{equation}
\label{eq:beckmann-lambda}
\Lambda \left( {\bf{\omega}} \right) = \frac{1}{2}\left( { {\rm{erf}} \left( a \right) - 1 + \frac{{e^{ - a^2 } }}{{a\sqrt \pi }}} \right)
\end{equation}
其中$a=\frac{1}{\alpha \tan \theta}$，$ {\rm{erf}} \left( x \right) = \frac{2}{{\sqrt \pi }}\int_0^x {e^{ - t^2 } {\rm{d}}t} $。对于各向异性表面，取
$$
\alpha = \sqrt {\alpha _x^2 \cos ^2 \varphi + \alpha _y^2 \sin ^2 \varphi }
$$
对于GGX分布，有
\begin{equation}
\label{eq:GGX-lambda}
\Lambda \left( {\bf{\omega}} \right) = \frac{{\sqrt {1 + \alpha ^2 \tan ^2 \theta } - 1}}{2}
\end{equation}

 

以上详细介绍了$G_1$函数的表现形式，然而在实际渲染中，$G_1$并不能直接套用在渲染方程里，我们还需要定义$G({\bf{\omega}_{i}},{\bf{\omega}_{o}})$，$G$的含义与$G_1$类似，$G_1$的分子是仅从$\bf{\omega}$方向能看到的微表面投影面积，$G$的分子给出的是在方向${\bf{\omega}_{i}}$和${\bf{\omega}_{o}}$上都能看见的微表面的投影面积。如果假设微表面同时被${\bf{\omega}_{i}}$和${\bf{\omega}_{o}}$看见的概率正比于微表面高出物体原表面的高度，则可按照该假设构造出$G$的经验公式

\begin{equation}
\label{eq:G_wi_wo}
G({\bf{\omega}_{i}},{\bf{\omega}_{o}}) = \frac{1}{1+{\Lambda(\bf{{\omega}}_{i})} + {\Lambda(\bf{{\omega}}_{o})}}
\end{equation}

实践中表明，该经验公式与实际测量的参数拟合的很好。注意式(\ref{eq:G_wi_wo})成立的条件实质上是基于三个假设：1)自遮挡函数$G_1$与${\bf{\omega}}_h$无关；2)微表面之间的高度互不相关；3)微表面同时被${\bf{\omega}_{i}}$和${\bf{\omega}_{o}}$看见的概率正比于微表面高出物体原表面的高度。

 

2.3 Torrance–Sparrow 模型

该模型是Torrance和Sparrow两人在1967年为了对粗糙金属表面进行物理建模时提出的，该模型假设金属物体表面是由许多具有理想镜面特性的微表面组成的。下面我们就来详细推导Torrance–Sparrow模型的解析形式。

 

如图10所示，设${\bf{\omega}_{h}}$为某一个微表面的法向方向，且为入射方向${\bf{\omega}_{i}}$和出射方向${\bf{\omega}_{o}}$的中间向量，$A$为物体表面处的微小面元，$\theta_o$为${\bf{\omega}}_{o}$和物体表面法向$\bf{n}$的夹角，$\theta_h$为${\bf{\omega}}_{i}$和${\bf{\omega}}_{h}$的夹角(也即${\bf{\omega}}_{o}$和${\bf{\omega}}_{h}$的夹角)。此时${\bf{\omega}_{i}}$和${\bf{\omega}_{o}}$呈镜面反射关系。根据辐射通量的定义和在“计算机图形学（一）——辐照度学概述”中所推导得到的式子

 

有
$$
{\rm{d}}\Phi _h = L_i \left( {{\bf{\omega}} _i } \right)\cos \theta_h {\rm{d}}{\bf{\omega}}_{i} {\rm{d}}A\left( {{\bf{\omega}} _h } \right)
$$
由2.1部分有关微表面表面积积分的内容，可以得到
$$
{\rm{d}}A\left( {{\bf{\omega}} _h } \right) = A \cdot D\left( {{\bf{\omega}} _h } \right){\rm{d}}{\bf{\omega}} _h
$$
所以
\begin{equation}
\label{eq:dphih}
{\rm{d}}\Phi _h = L_i \left( {{\bf{\omega}} _i } \right)A \cdot D\left( {{\bf{\omega}} _h } \right)\cos \theta_h {\rm{d}}{\bf{\omega}}_{i} {\rm{d}}{\bf{\omega}} _h
\end{equation}
设$F_r({\bf{\omega}}_{o})$为沿着入射方向${\bf{\omega}_{i}}$朝出射方向${\bf{\omega}_{o}}$的菲涅尔反射系数，$F_r({\bf{\omega}}_{o}) \in [0,1] $，则沿出射方向${\bf{\omega}_{o}}$的辐射通量${{\rm{d}}\Phi_o}$为
\begin{equation}
\label{eq:dphio}
{{\rm{d}}\Phi_o} = F_r({\bf{\omega}}_{o}){{\rm{d}}\Phi_h}
\end{equation}

根据辐射率的定义，又很容易得到出射辐射率$L_o$为
\begin{equation}
\label{eq:lowo}
{L_o}\left( {{{\bf{\omega }}_o}} \right) = \frac{{{\rm{d}}{\Phi _o}}}{{ A \cos {\theta_o }{\rm{d}}{ {\bf{\omega }} _o}}}
\end{equation}

将式(\ref{eq:dphih})、式(\ref{eq:dphio})代入式(\ref{eq:lowo})，得到
\begin{equation}
\label{eq:lo2}
L_o \left( {{\bf{\omega}} _o } \right) = \frac{{F_r \left( {{\bf{\omega}} _o } \right)L_i \left( {{\bf{\omega}} _i } \right)D\left( {{\bf{\omega}} _h } \right)\cos \theta _h {\rm{d}}{\bf{\omega}}_{i} {\rm{d}}{\bf{\omega}} _h }}{{\cos \theta _o {\rm{d}}{\bf{\omega}} _o }}
\end{equation}

 

假设图11中，入射光线沿着$IO$方向射入表面处$O$点，并沿着镜面反射方向$OR$出射，微表面局部法线为$\bf{h}$。由于$IP=PR$，因此${\rm{d}}A_R=4 \cdot {\rm{d}}A_P$。假设该半球面为单位半球，即$r=|OQ|=1$，则$OP=\cos \theta_i$。所以又有${\rm{d}}A_P=\cos^2 \theta_i \cdot {\rm{d}}A_Q'$。这造成${\rm{d}}A_R=4 \cos^2 \theta_i \cdot {\rm{d}}A_Q'$。注意到${\rm{d}}A_Q$和${\rm{d}}A_Q'$之间满足的投影关系，${\rm{d}}A_Q= \cos \theta_i \cdot {\rm{d}}A_Q'$，整理上述关系，很容易得到
$$
{\rm{d}}A_R= 4 \cos \theta_i \cdot {\rm{d}}A_Q
$$

由于球面是单位球面，根据立体角的定义，${\rm{d}}A_R={\rm{d}}{\bf{\omega}}_o$，${\rm{d}}A_Q={\rm{d}}{\bf{\omega}}_h$；由于图11中$\theta_i$和图10中$\theta_h$等价，$\theta_i=\theta_h$。因此
\begin{equation}
\label{eq:wo=4coswh}
{\rm{d}}{\bf{\omega}}_o = 4 \cos \theta_h \cdot {\rm{d}}{\bf{\omega}}_h
\end{equation}

将式(\ref{eq:wo=4coswh})代入式(\ref{eq:lo2})，得到
\begin{equation}
\label{eq:lo3}
L_o \left( {{\bf{\omega}} _o } \right) = \frac{{F_r \left( {{\bf{\omega}} _o } \right)L_i \left( {{\bf{\omega}} _i } \right)D\left( {{\bf{\omega}} _h } \right){\rm{d}}{\bf{\omega}}_{i} }}{{4 \cos \theta _o }}
\end{equation}

上式是仅在考虑单个微表面时的情况。对于由大量微表面组成的一般材质表面，一定会产生微表面间的遮挡，这时候就要引入上文提及的自遮挡函数$G$。修正后的式(\ref{eq:lo3})变为
\begin{equation}
\label{eq:lo4}
L_o \left( {{\bf{\omega}} _o } \right) = \frac{{F_r \left( {{\bf{\omega}} _o } \right)L_i \left( {{\bf{\omega}} _i } \right)D\left( {{\bf{\omega}} _h } \right)G\left( {{\bf{\omega}} _i ,{\bf{\omega}} _o } \right){\rm{d}}{\bf{\omega}}_{i} }}{{4\cos \theta _o }}
\end{equation}
根据BRDF的定义式
$$
f_r \left( {{\bf{\omega}} _i ,{\bf{\omega}} _o } \right) = \frac{{L_o \left( {{\bf{\omega}} _o } \right)}}{{L_i \left( {{\bf{\omega}} _i } \right)\cos \theta _i {\rm{d}}{\bf{\omega}} _i }}
$$
得到最终的Torrance-Sparrow BRDF解析形式为
\begin{equation}
\label{eq:torrance-sparrow-brdf}
f_r \left( {{\bf{\omega}} _i ,{\bf{\omega}} _o } \right) = \frac{{F_r \left( {{\bf{\omega}} _o } \right)D\left( {{\bf{\omega}} _h } \right)G\left( {{\bf{\omega}} _i ,{\bf{\omega}} _o } \right)}}{{4\cos \theta _o \cos \theta _i }}
\end{equation}

观察上述推导过程，我们发现Torrance-Sparrow BRDF模型并没有依赖于任何的微表面分布特征，因此该BRDF具有一定的普适性。但是该BRDF的推导仍基于一个基本假设：光线在射入微表面时满足镜面反射规律。