着色模型简介和实现（上）

摘要

着色模型（shading model）描述的是材质表面对光线的作用，它的输入是光线、入射角、反射角、材质的属性（例如粗糙度、金属度等），输出是材质的颜色。前面提到的BRDF模型就是shading model的一种，但是BRDF模型过于简单，很多复杂的物理现象无法描述。本文将接着BRDF继续介绍一些更复杂的shading model，以及UE4和Disney的实现。

本文主要介绍理论部分，下篇将介绍UE4、Filament以及Disney对shading model的详细实现。

理论

主要介绍对各向异性、次表面、清漆和布料几种物理现象的shading model。

Anisotropy

上图左：各向同性金属球；右：各向异性金属球。图片来源：Filament文档。

将各向同性的shading model拓展到各向异性，需要解决两个问题，

1. 需要引入哪些新的参数来描述各向异性？
2. 需要在原有BRDF模型的基础上修改公式的哪些部分来描述各向异性？

对于第一个问题，既然是“各向”异性，必然要分别定义各个方向上的差异的属性，这里的属性一般是粗糙度。方向一般在切平面上选取两个轴，即tangent方向和bitagent方向，他们与normal方向共通构成一组正交基底（TBN坐标系）。对于需要精细控制切平面方向的模型，比如头发，会单独提供切向贴图，指明切向方向，类似于法向贴图。其他情况一般由系统（shader）确定切向方向。所以总结下来，需要额外增加一个参数，将原有的粗糙度替换为tangent和bitangent方向的粗糙度（记为${\alpha }_t,{\alpha }_b$或${\alpha }_x,{\alpha }_y$）。

实际操作中，一般是增加一个anisotropy各向异性参数，然后用粗糙度$\alpha$（$\alpha =roughness\ast roughness$）和anisotropy计算tangent/bitangent两个方向上的粗糙度，具体的映射方法有很多。

• Neubelt and Pettine

$$\begin{array}{rl}{\alpha }_x& =\alpha \\ {\alpha }_y& =lerp(0,\alpha ,1-anisotropy)\end{array}$$

• Burley

$$\begin{array}{rl}{\alpha }_x& =\frac{\alpha }{\sqrt{1-0.9\times anisotropy}}\\ {\alpha }_y& =\alpha \sqrt{1-0.9\times anisotropy}\end{array}$$

• Kulla

$$\begin{array}{rl}{\alpha }_x& =\alpha \times (1+anisotropy)\\ {\alpha }_y& =\alpha \times (1-anisotropy)\end{array}$$

对于第二个问题，BRDF公式取决于三个部分，法向分布函数$D$，遮挡项$G$，和表示微观BRDF的菲涅尔项$F$。其中$F$肯定是不会变的，因为各向异性是宏观的现象，不会影响微观表面的性质。至于$D$和$G$，由于$G$都是从$D$推导出来的，因此我们先看一下$D$需要作出怎样的修改。

形状不变性

法向分布函数$D$的形状不变性指的是这样一种特性，改变$D$的粗糙度系数相当于对微表面进行拉伸，而不改变微表面的形状。

要说明白形状不变性，我们引入一个新的函数$P^{22}(x_{\tilde{m}},y_{\tilde{m}})$，描述微表面的法向量$\mathbf{m}$的2D斜率分布。其中，$\mathbf{m}=(x_m,y_m,z_m)$并且

$$\begin{array}{rl}(x_{\tilde{m}},y_{\tilde{m}})& =(-\frac{x_m}{z_m},-\frac{y_m}{z_m})=-\tan {\theta }_m(\cos {\varphi }_m,\sin {\varphi }_m)\\ \mathbf{m}& =\frac{(-x_{\tilde{m}},-y_{\tilde{m}},1)}{\sqrt{x_{\tilde{m}}^2+y_{\tilde{m}}^2+1}}\end{array}$$

$\mathbf{m}$的球坐标系坐标为$(1,\theta ,\varphi )$，如下图所示。

既然是分布函数，$P^{22}(x_{\tilde{m}},y_{\tilde{m}})$满足：

$${\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }{P}^{22}(x_{\tilde{m}},y_{\tilde{m}})d{x}_{\tilde{m}}d{y}_{\tilde{m}}=1$$

并且$P^{22}(x_{\tilde{m}},y_{\tilde{m}})$与$D$的关系为：

$$D(\mathbf{m})=\frac{{\chi }^{+}(\mathbf{n},\mathbf{m})}{{\cos }^4{\theta }_m}{P}^{22}(x_{\tilde{m}},y_{\tilde{m}})=\frac{{\chi }^{+}(\mathbf{n},\mathbf{m})}{(\mathbf{n}\cdot \mathbf{m}{)}^4}{P}^{22}(x_{\tilde{m}},y_{\tilde{m}})$$

考虑到$P^{22}(x_{\tilde{m}},y_{\tilde{m}})$同时也是粗糙度$\alpha$的函数，因此可以写为$P^{22}(x_{\tilde{m}},y_{\tilde{m}},\alpha )$。

满足形状不变性的法向分布函数，它的$P^{22}(x_{\tilde{m}},y_{\tilde{m}},\alpha )$符合这种形式：

$$\begin{array}{rl}{P}^{22}(x_{\tilde{m}},y_{\tilde{m}},\alpha )& =\frac{1}{{\alpha }^2}f\left(\frac{\sqrt{x_{\tilde{m}}^2+y_{\tilde{m}}^2}}{\alpha }\right)=\frac{1}{{\alpha }^2}f\left(\frac{\tan {\theta }_m}{\alpha }\right)\\ & =\frac{1}{{\alpha }^2}f\left(\frac{\sqrt{1-(\mathbf{n}\cdot \mathbf{m}{)}^2}}{\alpha (\mathbf{n}\cdot \mathbf{m})}\right)\\ D(\mathbf{m})& =\frac{{\chi }^{+}(\mathbf{n},\mathbf{m})}{(\mathbf{n}\cdot \mathbf{m}{)}^4}{P}^{22}(x_{\tilde{m}},y_{\tilde{m}})=\frac{{\chi }^{+}(\mathbf{n},\mathbf{m})}{{\alpha }^2(\mathbf{n}\cdot \mathbf{m}{)}^4}f\left(\frac{\tan {\theta }_m}{\alpha }\right)\\ & =\frac{{\chi }^{+}(\mathbf{n},\mathbf{m})}{{\alpha }^2(\mathbf{n}\cdot \mathbf{m}{)}^4}f\left(\frac{\sqrt{1-(\mathbf{n}\cdot \mathbf{m}{)}^2}}{\alpha (\mathbf{n}\cdot \mathbf{m})}\right)\end{array}$$

对于任意大于零的实数$\lambda$，有：

$$P^{22}(x_{\tilde{m}},y_{\tilde{m}},\alpha )=\frac{1}{{\lambda }^2}{P}^{22}(\frac{x_{\tilde{m}}}{\lambda },\frac{y_{\tilde{m}}}{\lambda },\frac{\alpha }{\lambda })$$

也就是说，改变粗糙度仅相当于对$P^{22}(x_{\tilde{m}},y_{\tilde{m}},\alpha )$进行拉伸。

前面提到的三种法向分布函数：Beckmann、BlinnPhong和GGX，只有BlinnPhong不满足形状不变性。推广的GTR公式也同样不满足形状不变性。

各向异性D和G

对于具有形状不变性性质的法向分布函数$D$，它的形式变为：

$$D(\mathbf{m})=\frac{{\chi }^{+}(\mathbf{n},\mathbf{m})}{{\alpha }_x{\alpha }_y(\mathbf{n}\cdot \mathbf{m}{)}^4}f\left(\frac{\sqrt{\frac{(\mathbf{t}\cdot \mathbf{m}{)}^2}{{\alpha }_x^2}+\frac{(\mathbf{b}\cdot \mathbf{m}{)}^2}{{\alpha }_y^2}}}{(\mathbf{n}\cdot \mathbf{m})}\right)$$

相应的遮挡项$G$，只需要改变$\Lambda (\mathbf{s})$中的变量$c$：

$$c=\frac{\mathbf{n}\cdot \mathbf{s}}{\sqrt{{\alpha }_x^2(\mathbf{t}\cdot \mathbf{s}{)}^2+{\alpha }_y^2(\mathbf{b}\cdot \mathbf{s}{)}^2}}$$

这里的$\mathbf{s}$可以是入射光线$\mathbf{l}$，也可以是出射光线$\mathbf{v}$。

有了一般的推广形式，我们将其应用到具有形状不变性的Beckmann和GGX模型上，有：

$$\begin{array}{rl}{D}_{Beckmann}(\mathbf{h},{\alpha }_x,{\alpha }_y)& =\frac{{\chi }^{+}(\mathbf{n},\mathbf{h})}{\pi {\alpha }_x{\alpha }_y(\mathbf{n}\cdot \mathbf{h}{)}^4}\exp \left(-\frac{\frac{(\mathbf{t}\cdot \mathbf{h}{)}^2}{{\alpha }_x^2}+\frac{(\mathbf{b}\cdot \mathbf{h}{)}^2}{{\alpha }_y^2}}{(\mathbf{n}\cdot \mathbf{h}{)}^2}\right)\\ D_{GGX}(\mathbf{h},{\alpha }_x,{\alpha }_y)& =\frac{{\chi }^{+}(\mathbf{n},\mathbf{h})}{\pi {\alpha }_x{\alpha }_y{\left(\frac{(\mathbf{t}\cdot \mathbf{h}{)}^2}{{\alpha }_x^2}+\frac{(\mathbf{b}\cdot \mathbf{h}{)}^2}{{\alpha }_y^2}+(\mathbf{n}\cdot \mathbf{h}{)}^2\right)}^2}\end{array}$$

Clear Coat

Clear Coat模型是多层材质的一种最简单的形式。Clear Coat的中文一般翻译为“清漆”，但是这个中文词汇跟英文词汇一样让人摸不着头脑。简单来说，Clear Coat指的是材质表面的一层薄薄的透明薄膜，比如木制家具表面刷的漆，汽车的喷漆等等。Clear Coat Layer一般是透明的，光线会透过Clear Coat Layer照射到Base Layer，为了模拟Clear Coat材质，一般将这两层材质按照下图的方式建模：

Clear Coat Layer会产生Specular反射，剩余的光线照射到Base Layer产生Specular和Diffuse反射。因此，相比于基础BRDF模型，Clear Coat模型多一个由Clear Coat Layer产生的specular项。

$$f=C_{diffuse}\cdot f_d+C_{specular}\cdot f_s+C_{clearcoat}\cdot f_{sc}$$

三个系数$C_x$是为了保证能量守恒，不同的引擎会选取不同的计算方法。clear coat layer的specular项$f_{sc}$的计算公式与$f_s$类似，也是DGF三项乘积的形式。

如果要表示一个Clear Coat模型，需要两个额外的参数，一个是ClearCoat系数，另一个是ClearCoatRoughness，指明Clear Coat Layer的粗糙度。

Subsurface

Subsurface模型也是一种多层材质，它可以很复杂（如用于绘制皮肤的SSSS模型），也可以很简单（如Disney在Principled BRDF中的实现）。

Subsurface模型描述的是名为次表面散射（subsurface scattering）的物理现象。次表面散射现象指的是光线进入材质内部，经过多次反射、吸收，最终折射出材质表面的现象。这个物理现象在BRDF模型中我们是用Diffuse分量来描述的，但是无论是Lambert模型还是Burley模型，都无法表现材质对光线的吸收作用，或者说，无法表现出“通透”效果。例如玉或者皮肤在强光照射下呈现的半透现象，或者普通光照下的通透感。下图[2]展示了采用Subsurface scattering模型的皮肤（右上）和未采用（右下）的对比。

但是完整实现subsurface scattering太复杂了，也没有必要。因此一般的subsurface模型都是通过添加一层额外的diffuse layer来近似表现次表面散射现象。这一点刚好跟Clear Coat模型相反——Clear Coat模型添加了额外的specular层，而subsurface添加了额外的diffuse层。

$$f=f_d+f_s+C_{subsurface}\cdot f_{sc}$$

Cloth

布料模型和其他几个shading model很不一样，这是因为布料的组成形式是纤维相互堆叠，纤维之间还存在空隙，不符合微表面理论。相应的，Cloth model的公式也跟基于微表面理论shading model不太一样。模拟Cloth的BRDF模型大致可以分为三类：

• 基于观察的empirical models
• 基于微表面理论的模型
• 微圆柱体模型（micro-cylinder model）

第一种基于观察的模型很大程度上靠艺术家调整参数来模拟布料效果，涉及到的参数比较多。第三种微圆柱体模型其实类似于头发模型——都是用圆柱体对材质进行建模，效果相对精准，但是也更繁琐。这里主要介绍第二种基于微表面理论的模型，它是通过选取合适的法向分布函数和遮挡函数来模拟布料效果。

常用的布料模型中的法向分布函数是Ashikhmin[3]提出的，用逆高斯分布来拟合法向分布，经过后续的补充，最终的形态为：

$$\begin{array}{rl}{D}_{Ashikhmin}(\mathbf{m})& =\frac{{\chi }^{+}(\mathbf{n}\cdot \mathbf{m})}{\pi (1+k_{amp}{\alpha }^2)}\left(1+\frac{k_{amp}\exp \left(\frac{(\mathbf{n}\cdot \mathbf{m}{)}^2}{{\alpha }^2\left((\mathbf{n}\cdot \mathbf{m}{)}^2-1\right)}\right)}{{\left(1-(\mathbf{n}\cdot \mathbf{m}{)}^2\right)}^2}\right)\\ G_{Cloth}(\mathbf{l},\mathbf{v})& =\frac{1}{4(\mathbf{n}\cdot \mathbf{v}+\mathbf{n}\cdot \mathbf{l}-(\mathbf{n}\cdot \mathbf{v})(\mathbf{n}\cdot \mathbf{l}))}\\ f(\mathbf{l},\mathbf{v},\mathbf{h})& =(1-F(\mathbf{h},\mathbf{l}))\frac{{\rho }_{ss}}{\pi }+F(\mathbf{h},\mathbf{l})D(\mathbf{h})G(\mathbf{l},\mathbf{v})\end{array}$$

其中，$\alpha$控制逆高斯函数的宽度，$k_{amp}$控制幅度。$f(\mathbf{l},\mathbf{v},\mathbf{h})$已经包含了diffuse和specular两部分。

float D_Ashikhmin(float roughness, float NoH) {
    // Ashikhmin 2007, "Distribution-based BRDFs"
	float a2 = roughness * roughness;
	float cos2h = NoH * NoH;
	float sin2h = max(1.0 - cos2h, 0.0078125); // 2^(-14/2), so sin2h^2 > 0 in fp16
	float sin4h = sin2h * sin2h;
	float cot2 = -cos2h / (a2 * sin2h);
	return 1.0 / (PI * (4.0 * a2 + 1.0) * sin4h) * (4.0 * exp(cot2) + sin4h);
}

另一种布料模型是Estevez和Kulla提出的Charlie模型，Charlie模型采用的是正弦函数的幂来拟合，而不是逆高斯函数。它的实现更简单、外观更柔和。

$$D_{Charlie}=\frac{(2+\frac{1}{\alpha }){\sin }^{\frac{1}{\alpha }}\theta }{2\pi }$$

实现：

float D_Charlie(float roughness, float NoH) {
    // Estevez and Kulla 2017, "Production Friendly Microfacet Sheen BRDF"
    float invAlpha  = 1.0 / roughness;
    float cos2h = NoH * NoH;
    float sin2h = max(1.0 - cos2h, 0.0078125); // 2^(-14/2), so sin2h^2 > 0 in fp16
    return (2.0 + invAlpha) * pow(sin2h, invAlpha * 0.5) / (2.0 * PI);
}