图形学 - 光和材质

这部分的内容主要是围绕光源以及光的反射进行讨论,现实中光的反射方式有两种: 镜面反射(specular reflection)和漫反射(diffuse reflection),光如何反射是取决于反射面(也就是它照射的物体表面)的,如果反射面光滑,那么就是镜面反射;相反,反射面粗糙,那么就是漫反射.

可以看到在镜面反射中(左边),入射光(incoming ray)和反射光(reflected ray)的角度是一样的,观察者(viewer)实际上只能看到一个完美的反射光(因为反射角要等于入射角,符合这个条件的只有一个),这样的反射光一被叫做镜面反射高光(specular highlights).

当然现实中可不是只能看到一条完美的反射光,而是一束反射光(specular reflection cone),就像这样,

处于完美反射光角度上的所看到的反射光是最亮,越是偏离这个角度,看到的反射光越暗.反射光束的角度是由物体表面的光泽度(shininess)决定的,

这图上面的球体的光泽度从左到右体依次递增,白色区域就是镜面反射高光,光泽度越高,光束的角度越小,高光区域越小.

仔细回想一下自己是否在日常生活中见到过弯曲的表面上有一个白色光的圆,其实就是上面这图的实例,比如光照下的矿泉水瓶,红酒杯等等.

这种光叫做 镜面反射光(specular light).

而在理想的漫反射(右边)中,入射光被分散成多个不同方向并且能量大小相同的反射光,观察者可以在任何方位看到反射光,如果光是以平行射线的形式照射到表面,光源位置在无限远处就属于这种情况,比如室外的太阳光,那么表面就是被均匀照明(参考上面的图: Diffuse Reflection),

假设光是以非平行射线的形式照射到表面,比如物体附近的台灯(向各个方向照射的光源)或者光照射到表面弯曲的物体(曲面)都属于这种情况,表面并不是被均匀照明的,表面某个点上的照明程度取决于射线和表面的角度.

这种光叫做 漫反射光(diffuse light).

与弹球可以多次弹跳一样,光也可以被多次反射,在实际中,当光照射到物体表面时,部分光被吸收,部分光会被反射(除非物体完全不反光),然后这些被反射的光照射到别的物体上(比如,生物的眼睛,因此生物才能看到这个物体),一直循环这个过程直到再也没有多余的光能够被反射.

在这个过程中,被反射的光虽然也能照明其它物体,但并非最直接(也就是初始,没有经过"弹跳")的光源,这种间接光叫做 环境反射光(ambient light).

为了更好理解,举个现实例子:

一个房间里面打开唯一的灯,灯位于房间天花板的一个角落,这个时候大部份光都是可以确定直接来源,也就是灯;但是房间里面有一样物体的背面是没有直接被灯照到,但仍然可以看到它的背后,

这是因为光是会像弹球那样弹跳的,比如照射到墙壁弹跳一次,弹跳到另外一个东西上,不知道经过多少次弹跳最后弹到物体的背面,最后弹跳到眼睛上,就这样被看到了.

还有一个例子就是月亮,它的光就是从太阳借过来的,即使不开灯也能有月亮的光,这就是环境反射光的含义.

光在传播过程中被物体吸收从而减弱,这叫做衰减(attenuation),在现实中,光单纯地在空气中穿梭也会衰减,随着传播距离越远,光的强度越弱,想不衰减地则只有在真空中穿梭.

所以总结下来,按照反射的情况对光进行划分有三种光: 镜面反射光,漫反射光和环境光,它们三者的颜色之和就是光的颜色.

现实中每个物体对同样的光都有不同的反应,比如金属品和木制品,通常金属品的表面更加光滑,反光效果更好,

而木制品一般都比较粗糙,反光效果比较差. 甚至有些物体还会自己发光,比如汽车前灯,荧光棒等等,

这种自己发出的光叫做发射光(emissive), 物体对光的作出的反应取决于材质(material),可以理解为物体表面的特点,

比如颜色,还有在上面的 Phong lighting model 的例子中提到过的光泽度 \(\mathrm{shininess}\), 这些都是材质的属性.

在开始介绍光照模型前先思考一个问题, "光照计算" 具体实在计算什么呢?

光照计算就是为了得出物体在观察者眼中的颜色, 光照计算和阴影是没有关系的, 阴影的内容可以单独作为一个主题了.

图形学第一个提出的真实光照模型是 Phong 光照模型(Phong lighting model),

虽然它现在已经不是主流的照明方案, 但其核心思想依然是不少光照模型的基础,

比如 OpenGL 固定流水线(fixed function pipeline)中所使用的 Blinn-Phong shading model.

Phong 光照模型认为反射光分为三个部分: 环境反射光, 镜面反射光和漫反射光.

为了简单介绍这个模型是怎么计算的,假设现有一个发出 \(\mathrm{light\_color}\) 颜色的光源照射一个物体 \(o\) 上,

该物体颜色为 \(\mathrm{object\_color}\),现在要求计算出物体反射的光 \(\mathrm{result\_color}\).

其实总体思路就是要计算出光能够产生出多少 \(\mathrm{ambient}\), \(\mathrm{specular}\) 以及 \(\mathrm{diffuse}\),然后这三项经过物体的反射后得到的颜色分别是多少,

最后的总和就是物体的颜色 \(\mathrm{result\_color}\) 了.

有一点值得注意, 光照的计算可以发生在世界空间上, 也可以发生在视点空间上, 重点在与产于计算的参数必须位于同一个空间上, 否则计算就是有问题的.

不过, 这片笔记的的光照计算是假设在世界空间上的, 所以被照射的顶点, 光源位置, 法线向量以及相机位置全都是这个空间上的.

当光照射到物体表面时,物体材质的颜色就意味着是对光的反射率.

根据之前提到过的环境反射光,镜面反射光以及漫反射,材质的颜色分为环境反射色(ambient color),镜面反射色(specular color),以及漫反射色(diffuse color),此外还有一个没有对应的发射色(emissive color).

所以 OpenGL 的光照模型把光分为 4 种,从物体表面到眼睛中的颜色就是计算出4种光的混合结果:

\(\mathrm{result\_color} = \mathrm{emissive} + \mathrm{ambient} + \mathrm{diffuse} + \mathrm{specular}\).

除 \(\mathrm{emissive}\) 以外,每一项是光和材质共同的计算结果,这很好理解:因为光照射到物体,物体把部分光反射到眼睛上,反射了多少光就得看照射的光有多少以及材质的反射程度,而自发光是没有经过反射直接到眼睛的.

严格来说, \(\mathrm{emssive}\) 不是 Blinn-Phong shading model 里面定义的.

OpenGL 允许有多个光源,并且拥有一个 全局环境反射光(global ambient light),它不属于任何一个光源,也就是就算所有光源都关闭了,全局环境反射光还在,光源的环境反射光以及全局环境反射光共同构成整个场景的环境反射光.

OpenGL 提供 glLight*(light, pname, value) 一类的函数来对光的参数 pname, 以下是对参数 pname 的说明:

属性	默认值	描述
GL_AMBIENT	(0.0, 0.0, 0.0, 1.0)	形式如(x, y, z, w),环境反射光的 RGBA 强度
GL_DIFFUSE	(1.0, 1.0, 1.0, 1.0)	形式如(x, y, z, w),漫反射光的 RGBA 强度
GL_SPECULAR	(1.0, 1.0, 1.0, 1.0)	形式如(x, y, z, w),漫反射光的 RGBA 强度
GL_POSITION	(0.0, 0.0, 1.0, 0.0)	形式如(x, y, z, w),如果 w 为0,那么光源位置就是在无限远处,因此认为光线之间是平行的,(x, y, z) 表示 指向光源的方向,这种光源叫做 定向光源(directional light),而默认值的意思就是光源往 \(z\) 轴的负方向,也就是 (0, 0, -1) 发出光线;若 w 不为0,那么光源位置位于场景附近,位于 (x, y, z),以光源为中心光向各个反向出发,这种光源叫做 位置光源(positional light), 偶尔也被叫做 点光源(point light). 另外,把这种光源屏蔽相当的一部分就能够达到聚光灯的效果
GL_SPOT_DIRECTION	(0.0, 0.0, -1.0)	形式如(x, y, z), 聚光灯的方向,聚光灯本质就是位置光源,所以一般来说 GL_POSITION 的 w 不应该为 0,不过 OpenGL 并没有限制这么做,但是采用定向光加聚光灯设置这种组合不一定是你想要的结果
GL_SPOT_EXPONENT	0.0	[0,128] 之间的整数或者浮点数,聚光灯的强度分布,如果为 0,那么聚光灯照射范围内的所有光的强度都一样,如果为正数,光线离照射中心越远,光的强度越弱,如果该值越大,那么衰减程度就越大,最后光的中心区域和外围区域差别越明显
GL_SPOT_CUTOFF	180.0	, [0,90] 之间的整数或者浮点数以及特殊值 180.0,聚光灯光束角度的 \(\frac{1}{2}\),当为 180.0 的时候就表示没有屏蔽位置光源.
GL_CONSTANT_ATTENUATION	1.0	恒定衰减因素(factor),下文用 \(k_{c}\) 表示
GL_LINEAR_ATTENUATION	0.0	线性衰减因素,下文用 \(k_{l}\) 表示
GL_QUADRATIC_ATTENUATION	0.0	二次衰减因素,下文用 \(k_{q}\) 表示

看到上面这些属性基本上能对 OpenGL 里面的光有一个形象了,可以看出并没有直接设置光的颜色,

取而代之的是用了 \(\mathrm{GL\_AMBIENT}\), \(\mathrm{GL\_DIFFUSE}\) 和 \(\mathrm{GL\_SPECULAR}\) 来设置光的颜色,

这样在后续就不需要把光的每个部分求出来. 接下来就是一些其他的补充了.

首先是 光的衰减,对于定向光源来说,光的衰减计算是被禁掉的,因为光是随传播距离增加而减弱,而定向光源是位于无限远的,所以这是不能可能无限衰减的.

但是对于位置光源来说可以启用衰减计算,这里的材质发射光和全局环境反射光是不会衰减的,所以衰减的光就只有环境反射光,漫反射光以及镜面反射光.

另外, OpenGL 提供 glLightModel*(pname, param) 调整计算规则, 接下来是参数 pname 的说明:

直接从 glMaterial*(face, pname, param) 切入主题,下面是它能够设置的材质属性参数 pname 说明:

参数 face 用来指定被设定对象是正面还是背面,还是两面都要设置, 值分别为: GL_FRONT, GL_BACK 和 GL_FRONT_AND_BACK.

这里的颜色模式用的是 RGBA 模式.

先介绍一下颜色的运算,假设有两个颜色: \(A = (R1, G1, B1, A1)\) 和 \(B = (R2, G2, B2, A2)\):

颜色加法: \(A + B = (R1+R2, G1+G2, B1+B2, A)\),加法表示颜色之间混合;

颜色乘法: \(A \times B = (R1R2, G1G2, B1B2, A)\),乘法表示颜色经过缩放,物体反射入射光就是这么一个例子;

其中 \(A\) 等于顶点的漫反射色 \(\mathrm{GL\_DIFFUSE}_{material}\) 的 alpha 值,所以计算的时候可以不用理会 \(A\) 的计算,下文也是这么做.

当光照射一个顶点时,顶点颜色的计算过程如下:

\(\mathrm{color}_{vector} = \mathrm{Material\_Emission} + \mathrm{Scaled\_Global\_Ambient\_Light} + \mathrm{Contributions\_From\_Light\_Sources}\).

\(\mathrm{Material\_Emission}\) 就是材质的 \(\mathrm{GL\_EMISSION}\) 属性.

\(\mathrm{Scaled\_Global\_Ambient\_Light}\) 是指被材质反射过的全局环境光:

\(\mathrm{Scaled\_Global\_Ambient\_Light} = \mathrm{GL\_AMBIENT}_{material} \times \mathrm{GL\_LIGHT\_MODEL\_AMBIENT}\).

\(\mathrm{Contributions\_From\_Light\_Sources}\) 则是比较复杂,它指所有光源的总和.

一个光源由环境反射光,漫反射光以及镜面反射光共同构成(三项的和),另外还有两个因素需要考虑:

光有没有经过衰减,以及光源是否为聚光灯.

假设现有 n-1 个光源,把所有光源的构成全部加在一起,总的计算过程如下:

\(\mathrm{contribution} = \mathrm{attenuation\_factor} \times \mathrm{spotlight\_effect} \times (\mathrm{ambient} + \mathrm{diffuse} + \mathrm{specular})\)

\(\mathrm{Contributions\_From\_Light\_Sources} = \displaystyle\sum_{i=0}^{n-1}\mathrm{contribution}_{i}\).

\(\mathrm{attenuation\_factor}\) 就是在 Phong lighting model 提到过的系数 \(k\), 这里只是三项的系数都一样.

接下来分别对 \(\mathrm{contribution}\) 的每一项计算进行深入了解.

• \(\mathrm{attenuation\_factor}\)

  光线的衰减因素.

  如果光源是位置光源,那么计算如下: \(\mathrm{attenuation\_factor} = \frac{1}{k_{c} + k_{l}d + k_{q}d^{2}}\),

  其中,

  \(d\) 是指 \(\mathrm{GL\_POSITION} = (X, Y, Z, 1)\) 的光源到顶点 \((x_{v}, y_{v}, z_{v}, 1)\) 的距离: \(d = \sqrt{(X - x_{v})^{2} + (Y - y_{v})^{2} + (Z - z_{v})^{2} + (1 - 1)^{2}}\),

  \(k_{c} = \mathrm{GL\_CONSTANT\_ATTENUATION}\),

  \(k_{l} = \mathrm{GL\_LINEAR\_ATTENUATION}\)

  \(k_{q} = \mathrm{GL\_QUADRATIC\_ATTENUATION}\).

  如果光源是定向光源,那么结果直接固定为1: \(\mathrm{attenuation\_factor} = 1\).

• \(\mathrm{spotlight\_effect}\)

  聚光灯效果.

  分三种情况进行计算:

  当光源不是聚光灯,也就是说 \(\mathrm{GL\_SPOT\_CUTOFF} = 180.0\) 时, \(\mathrm{spotlight\_effect} = 1\).

  当光源是聚光灯,在这个前提下要分两种情况讨论: 顶点是否位于聚光灯的照射范围内,

  判定是哪种情况其实很简单,做法就是计算出 从光源到顶点的直线 和 聚光灯的方向 \(\mathrm{GL\_SPOT\_DIRECTION}\) 之间的夹角是不是小于 \(\mathrm{GL\_SPOT\_CUTOFF}\).(如果不理解的话回去看 \(\mathrm{GL\_SPOT\_CUTOFF}\) 参数的配图).

  假设有顶点 \(v = (x_{v}, y_{v}, z_{v})\), \(\mathrm{GL\_POSITION} = (x_{l}, y_{l}, z_{l}, 1)\) 的光源,

  先计算出光源到顶点的单位向量 \(\vec{n} = \frac{v - \mathrm{GL\_POSITION}}{|v - \mathrm{GL\_POSITION}|} = \frac{(x_{v} - x_{l}, y_{v} - y_{l}, z_{v} - z_{l})}{\sqrt{(x_{v} - x_{l})^{2} + (y_{v} - y_{l})^{2} + (z_{v} - z_{l})^{2}}}\),

  然后通过 \(\vec{n} \cdot \mathrm{GL\_SPOT\_DIRECTION} = |\vec{n}||\mathrm{GL\_SPOT\_DIRECTION}|\cos\theta\) 来求出夹角的 \(cosine\) 值.

  如果 \(\mathrm{GL\_SPOT\_DIRECTION}\) 是单位向量的话, \(\vec{n} \cdot \mathrm{GL\_SPOT\_DIRECTION} = \cos\theta\).

  (实际上方向也的确是用单位向量表示,后面就不赘述了,凡是方向统一暗示单位向量).

  因为 \(\mathrm{GL\_SPOT\_DIRECTION}\) 的取值可能是 \([0.0, 90.0]\) 以及 \(180.0\),不考虑 \(180.0\) 的话,

  整个计算结果应该是 \(max(\vec{n} \cdot \mathrm{GL\_SPOT\_DIRECTION}, 0)\), 如果这个结果比 \(\mathrm{GL\_SPOT\_CUTOFF}\) 的 \(cosine\) 值要大,

  根据 \(cosine\) 函数的性质可以得出顶点在聚光灯的照射范围内, 否则顶点就在照射范围内.

  如果 顶点不在聚光灯照射范围内,那么 \(\mathrm{spotlight\_effect} = 0\);

  如果 顶点在聚光灯照射范围内,那么 \(\mathrm{spotlight\_effect} = max(\vec{n} \cdot \mathrm{GL\_SPOT\_DIRECTION}, 0)^{\mathrm{GL\_SPOT\_EXPONENT}}\).

• \(\mathrm{ambient}\)

  就只是简单的环境光经过物体反射得到的环境反射色:

  \(\mathrm{ambient} = \mathrm{GL\_AMBIENT}_{\mathrm{light}} \times \mathrm{GL\_AMBIENT}_{\mathrm{material}}\).

• \(\mathrm{diffuse}\)

  漫反射光照射到物体表面得到的漫反射色.

  这需要计算光是否直接照射到顶点,这个计算过程还是很好理解的,只要 指向光源方向 和 顶点的法线,就可以断定光是直接照射到顶点上,这是亮度最大,越偏离这个方向亮度就越小.

  整个漫反射色的计算过程很好理解的,先假设目前有 \(\mathrm{GL\_POSITION} = (x_{l}, y_{l}, z_{l}, w)\) 的光源,以及被照射的顶点 \(v = (x_{v}, y_{v}, z_{v})\),顶点的法线为 \(N = (x_{n}, y_{n}, z_{n})\).

  先计算出指向光源的方向 \(L\).

  如果是定向光源(\(w = 0\)),那么 \(L = (x_{l}, y_{l}, z_{l})\);

  如果是位置光源(\(w \neq 0\)),那么指向光源的方向就是 从顶点到光源的方向,那么 \(L = \frac{\mathrm{GL\_POSITION} - v}{|\mathrm{GL\_POSITION} - v|} = \frac{(x_{l} - x_{v}, y_{l} - y_{v}, z_{l} - z_{v})}{|(x_{l} - x_{v}, y_{l} - y_{v}, z_{l} - z_{v}|}\).

  一旦计算出 \(L\) 以后,就进而可以得出 \(L\) 和 \(N\) 夹角 \(\theta\) 的 \(cosine\) 值了: \(\cos\theta = \frac{L \cdot N}{|L||N|}\), 这个值也表示了 \(L\) 和 \(N\) 不一致的时顶点接收到入射光的强度系数.

  如果 \(L \cdot N < 0\),也就是说光源位于物体表面的"错误"面,也就是光没有照射到顶点 \(v\) 上,可以通过 \(max(L \cdot N, 0)\) 来纠正错误为 0.

  最后整个计算过程就是 \(\mathrm{diffuse} = max(L \cdot N, 0) \times \mathrm{GL\_DIFFUSE}_{\mathrm{light}} \times \mathrm{GL\_DIFFUSE}_{\mathrm{material}}\).

  整个计算思路和 Phong lighting model 的那个例子一模一样,只是这里考虑了材质以及光的类型这两个因素.

• \(\mathrm{specular}\)

  镜面反射光照射到物体表面得到的镜面反射色.

  这一个步计算就是 Blinn-Phong shading model 和 Phong lighting model 的差别了,之所以有这个新的模型那必然是 Phong lighting model 存在一些限制的.

  先来假设一些变量来更好的做说明,设有和 \(\mathrm{diffuse}\) 里面一样的光源,顶点以及顶点的法线,光源同样分定向和位置光源两种情况进行讨论, \(L\) 和 \(N\) 的具体计算过程就不赘述了,

  另外多设一个反射向量 \(R\), 以及相机位置 \(v_{\mathrm{camera}} = (x_{\mathrm{camera}}, y_{\mathrm{camera}}, z_{\mathrm{camera}})\), 还需要一个从顶点 \((x_{v}, y_{v}, z_{v})\) 到 viewer 的方向 \(V\),

  最后再提醒一次,这些向量都是单位向量,否则不能保证计算结果的正确性.

  不过 \(V\) 需要分情况进行讨论:

  如果 \(\mathrm{GL\_LIGHT\_MODEL\_LOCAL} = \mathrm{GL\_TRUE}\), 那 \(V = v_{\mathrm{camera}} - v = (x_{\mathrm{camera}} - x_{v}, y_{\mathrm{camera}} - y_{v}, z_{\mathrm{camera}} - z_{v})\); 否则, \(V = (0, 0, 1)\).

  Phong lighting model 存在的限制就是 \(R\) 和 \(V\) 之间的角度不能超过 90 度,对于 \(\mathrm{diffuse}\) 来说, \(L\) 和 \(N\) 是不能超过 90 度的,

  但是对于 \(\mathrm{specular}\) 来说, \(R\) 和 \(V\) 是可以超过 90 度的,

  

  Figure 4: 镜面反射2(图片来源于网络)

  根据 Phong lighting model 的计算方法,如果 \(R\) 和 \(V\) 的夹角超过 90 度,那么 \(\mathrm{specular}\) 的结果就是 0,也就是说在 90 度的位置附近会看到一条明显的光暗分界线,这导致看起来不那么真实.

  James F. Blinn 在 1977 年提出了一个 \(\mathrm{specular}\) 计算方法: 不依赖反射向量,利用半角向量 \(H\) (halfway vector)进行角度的计算,半角向量就是 \(L\) 和 \(V\) 夹角的 \(\frac{1}{2}\) 的那个方向,

  

  Figure 5: 半角向量(图片来源于网络)

  只要把 \(L\) 或者 \(V\) 平移一下就可以发现 \(H = \frac{L + V}{|L + V|}\), 然后把原本计算 \(R \cdot V\) 改为计算 \(H \cdot N\).

  最终的计算结果和 Phong lighting model 的 \(\mathrm{specualr}\) 会稍微不一样,但依然是合理的,当 \(V\) 正好等于 \(R\),半角向量也正好等于 \(N\),这对应了 \(V\) 和 \(R\) 一致时亮度最高,也符合了 \(V\) 越时靠近 \(R\),亮度越高.

  所以 \(\mathrm{specular} = max(H \cdot N, 0)^{\mathrm{GL\_SHININESS}} \times \mathrm{GL\_SPECULAR}_{\mathrm{light}} \times \mathrm{GL\_SPECULAR}_{\mathrm{material}}\).

在 3D 建模中, 宏表面是实际存在的, 由顶点构成, 宏表面的法线由顶点法线决定;

微表面只在理论层面上存在, 由多个微平面(microfacets)构成, 这些微平面的法线是杂乱无章的, 因此, 只有一部分微平面才能反射光线.

在图形学中, 对于几何体外观的建模, 总会假设一定的建模尺度和观察尺度:

• 宏观尺度 (Macroscale): 通过三角形网格(mesh)对几何体外观进行建模, 由顶点法线(Vertex Normal)给每个顶点提供法线信息.
• 中尺度 (Mesoscale): 通过纹理(贴图)对几何外观进行建模, 由法线贴图(Normal Map)给每个像素提供法线信息.
• 微观尺度 (Microscale): 通过 BRDF 为几何外观进行建模, BRDF 会使用一个微表面的法线分布函数 NDF (Normal Distribution Function), 配合粗糙度贴图(Roughness Map)统计出一个宏表面上有多少微平面可以反射光线, 给每个亚像素(subpixel)提供法线.

亚像素

在数字成像领域，由于物理上已经无法在相邻的物理像素之间增加更多实际的感光单元，因此, 在软件上通过插值算法引入虚拟像素以在图像中提高测量的精度, 这些虚拟像素就是亚像素.

比如在下面的图中, 1 个方点就是 1 个物理像素, 1 个圆点就是 1 个亚像素.

软件上的 1 个像素是如下图那般由 4 个物理像素(1,2,3,4)构成的.

Figure 10: 亚像素

目前图形学中大多数微观尺度的结构都是高度场(heightfields),

每平面坐标都会有一个值来反应它的高度, 这些值的集合就一个高度场.

微表面就是一个高度场, 而高度场结构的微表面有一个特征: 微平面的法线 \(m\) 不能向下.

接下来就是分部地解析反射率方程 (The reflectance equation):

\(L_{o}(p, \omega_{o}) = \int_{\Omega}f_{r}(p, \omega_{i}, \omega_{o})L_{i}(p, \omega_{i})|w_{i} \cdot n|d\omega_{i}\)

能量守恒(Energy Conservation): 出射光的能量不能超过入射光的能量, 排除自发光的表面光.

为了遵守能量守恒, 需要对镜面反射光(specular light)和漫反射光(diffuse)进行明确的区分,

光在击中表面时会分成折射光(refraction)和反射光(reflection)两部分:

反射光是没有被表面吸收并且被直接反射的光, 也就是熟知的镜面反射光;

折射光是进入表面并且被吸收的光, 就是我们熟知道的漫反射光.

所谓的被吸收其实就是光在进入表面后, 与材质的粒子发生多次碰撞后全部转化成热能,

然而有一部分光会在全部转化成热能之前以随机的方向逃离出表面,

这一部分光被称为散射光(scattered light), 在逃出表面后成为漫反射光的一部分.

不过在目前介绍的 PBR 模型中会假设折射光被全部吸收, 忽视散射光在离开表面后的远处消失的效果,

不过有一项叫做次表面散射(subsurface scattering)的着色技术可以实现这些效果, 但伴随而来的是性能代价.

因为遵守能量守恒, 所以满足该关系: \(k_{\mathrm{specular}} + k_{\mathrm{diffuse}} = 1.0\),

其中 \(k_{\mathrm{specular}}\) 和 \(k_{\mathrm{diffuse}}\) 分别是镜面反射光占总能量的比例, 以及漫反射光占总能量的比例.

金属表面(metallic surface)只有反射光, 没有折射光, 换句话就是金属表面只显示镜面反射色, 不显示漫反射色;

非金属表面(non-metallic surface), 或者叫电解质(dielectric)

BRDF 是一个以入射光方向 \(\omega_{i}\), 出射光方向 \(\omega_{o}\), 宏表面法线 \(n\) 以及表面粗糙度 \(\alpha\) 为参数的函数,

计算出每个光线 \(\omega_{i}\) 在特定材质的不透明表面光上产生多少反射光.

如果一个表面和镜子一样光滑, 那么 BRDF 函数就会返回 \(0.0\) 来表示入射光的角度和出射光的角度是一样的,

BRDF 并不是一个具体的函数, 可以有很多种选择, 最早的 BRDF 是由 Fred Nicodemus 在 1965 年提出: \(f_r(\omega_i, \omega_o) = \frac{d L_o (\omega_o)}{d E_i(\omega_i)} = \frac{d L_o(\omega_o)}{L_i(\omega_i) \cos\theta_i d\omega_i}\).

这个公式之所以定义为辐射 \(L_o\) 和辐射度 \(E_i\) 之比, 而不是 \(L_o\) 和 \(L_i\) 或者 \(E_o\) 与 \(E_i\) 之比,

是因为在入射时需要考虑照射面积, 在反射时需要考虑每立体角的辐射通量, 反射角度以及投影面积.

目前业界主流的 PBR 渲染管线所使用的 BRDF 是 Cook-Torrance BRDF: \(f_{r} = k_{\mathrm{diffuse}}f_{\mathrm{lambert}} + k_{\mathrm{specular}}f_{\mathrm{cook-torrance}}\).

\(f_{\mathrm{lambert}} = \frac{c}{\pi}\), \(c\) 是表面颜色(surface color) 或者反射率(albedo), 就像 \(\mathrm{GL\_DIFFUSE}_{\mathrm{material}}\) 一样决定了反射多少漫反射光.

\(f_{\mathrm{cook-torrance}} = \frac{DFG}{4(\omega_{o} \cdot n)(\omega_{i} \cdot n)}\) 计算反射光, \(f_{\mathrm{cook-torrance}}\) 的定义源于最早的 BRDF, 其中 \(D\) 是法线分布函数, \(G\) 是几何函数, \(F\) 是菲涅耳方程.

接下来分别学习这 3 个部分, 并在最后介绍 \(f_{\mathrm{cook-torrance}}\) 是怎么推导出来的.

另外, 为了方便码字, 在解析 BRDF 这部分的内容里面会采用 \(v\) 代表 \(\omega_{o}\), \(l\) 代表 \(\omega_{i}\).

还有这部分的参数比较多, 时间久了可能会忘记这些参数的含义, 为了方便之后快速查询, 这里给出它们的汇总:

\(l\) 是入射光方向;

\(v\) 是出射光方向;

\(m\) 是微平面的法线向量, 使用 \(v\) 和 \(l\) 的半角向量 \(h\) 作为微平面的法线向量, 也就是 \(m = h\);

\(n\) 是宏表面的法线向量;

\(\alpha\) 是微平面的粗糙度

经过前面的推导后, 我们对反射率方程有了新的认识:

\(L_o(p, \omega_o) = \int_{\Omega} (k_{\mathrm{diffuse}} \frac{c}{\pi} + k_{\mathrm{specular}} \frac{F(\omega_i, h, F_0) D(m) G(\omega_o, \omega_i, h)}{4(\omega_o \cdot n)(\omega_i \cdot n)}) L_i(p, \omega_i) n \cdot \omega_i d\omega_i\), 其中 \(h = \frac{\omega_i + \omega_o}{|\omega_i + \omega_o|}\).

可以进一步变换, 把方程分成漫反射辐射度(diffuse irradiance)和镜面反射辐射度(specular irradiance):

\(\begin{equation*} \begin{aligned} L_o(p, \omega_o) &= \int_{\Omega} k_{\mathrm{diffuse}} \frac{c}{\pi} L_i(p, \omega_i) n \cdot \omega_i d\omega_i + \int_{\Omega} k_{\mathrm{specular}} \frac{F(\omega_i, h, F_0) D(m) G(\omega_o, \omega_i, h)}{4(\omega_o \cdot n)(\omega_i \cdot n)} L_i(p, \omega_i) n \cdot \omega_i d\omega_i \\ &= k_{\mathrm{diffuse}} \frac{c}{\pi} \int_{\Omega} L_i(p, \omega_i) n \cdot \omega_i d\omega_i + k_{\mathrm{specular}} \int_{\Omega} \frac{F(\omega_i, h, F_0) D(m) G(\omega_o, \omega_i, h)}{4(\omega_o \cdot n)(\omega_i \cdot n)} L_i(p, \omega_i) n \cdot \omega_i d\omega_i \end{aligned} \end{equation*}\)

在了解 PBR 的理论以后, 就可以描述美术人员如何为 PBR 方程创作制作表面的物理参数, 每个参数都可以通过纹理贴图来定义或建模, 并且纹理贴图能够以逐片元的精度(per-fragment)控制表面的物理属性.

以下是 PBR 管线中常用的贴图类型:

• 反射贴图 (Albedo texture)

  反射贴图为每个纹素(texel)指定表面的颜色, 如果纹素为金属材质则对应其基础反射率.

  作用基本上与中的漫反射贴图 (diffuse texture) 一致, 漫反射贴图通常包含了细微阴影或暗部凹陷, 这正是反射贴图需要避免的部分, 只能包含表面颜色或者金属的基础反射率.

  反射贴图是 \(f_{\mathrm{lambert}} = \frac{c}{\pi}\) 中的 \(c\) 的来源.

• 法线贴图 (Normal map)

  法线贴图逐片元地指定宏表面的法线, 是 \(n\) 的来源, 用来模拟表面的凹凸效果.

• 金属贴图 (Metallic map)

  金属贴图定义纹素是否金属, 从金属贴图取得的值被称为金属度,

  通常为灰度值(grayscale value), 或直接使用非黑即白的二值化形式.

  如果金属度不为 1, 那么材质就被认为是非金属材质, 通常渲染引擎会采用一个表示非金属材质的值作为默认值.

  金属度影响的是 PBR 菲涅耳反射 \(F\) 中的镜面反射率 \(F_0\) 的计算, 以及漫反射系数 \(k_{\mathrm{diffuse}}\) 的计算.

  严格来说, \(F_0\) 是由反射率 \(c\) 和金属度 \(m\) 共同影响的: \(F_0 = c_0 + (c - c_0) \times m\).

  其中 \(c_0\) 为默认反射率, 这个值通常表示是 \((0.03, 0.03, 0.03)\), \(c\) 和 \(c_0\) 均为向量.

  这条公式实际上是近似计算, 在原本的计算中, 需要讨论材质是金属才是非金属,

  如果是金属材质, 那么认为表面不存在漫反射, 镜面反射为最大, 也就是直接以 \(c\) 作为镜面反射率: \(F_0 = c\);

  如果是非金属材质, 那么 \(F_0 = \frac{(n_1 - n_2)^2}{(n_1 + n_2)^2}\), \(n_1\) 和 \(n_2\) 分别是两种介质的折射率(index of refraction),

  通常 \(n_1\) 是空气的折射率 1, \(n_2\) 为非金属表面的折射率, \(F_0\) 意味着从空气射入到表面中时的镜面反射率.

  要找出特定非金属材质的金属度, 先在 IOR LIST 找到该材质的折射率 \(n_2\) 算出 \(F_0\), 最后计算出金属度 \(m = \frac{F_0 - s}{c - s}\).

  \(k_{\mathrm{diffuse}}\) 也存在类似关系, 如果金属度为 1, 那么就认为漫反射为 0: \(k_{\mathrm{diffuse}} = c + ((0, 0, 0) - c) \times m\).

• 粗糙度贴图 (Roughness map)

  粗糙度贴图以逐纹素的方式定义表面的粗糙度. 粗糙贴图是 BRDF 中 \(D\) 和 \(G\) 项粗糙度 \(\alpha\) 的来源.

  表面越粗糙, 反射越宽泛且模糊, 而越光滑, 反射越集中且清晰, 简单点说就是影响高亮光圈的大小以及亮度.

  部分引擎使用光滑度贴图 (Smoothness map) 而非粗糙度贴图, 可以通过 \(1.0 - \alpha\) 转换成粗糙度.

• 环境光遮蔽贴图 (Ambient occlusion map, 简称 AO map)

  环境光遮蔽贴图为表面以及其周围几何体设置额外的阴影因子 \(ao\) , 给反射贴图补充阴影信息, 模拟几何遮挡导致的环境光衰减.

  注意, 是只影响环境光部分, 不影响直接光照, 大概为 \(ao \times ambient\), 其中 \(ambient\) 是环境光.

该实现的 Fragment Shader 来源于 LearnOpenGL, 但作了一些调整, 使得可用于 WebGL.

以下 Shader 作为场景中每个物体的材质.

Vertex Shader:

#version 100
#if defined(USE_TEXTURE)
attribute vec3 tangent;
attribute vec3 aBitangent;
#endif
attribute vec3 normal;
attribute vec3 position;
attribute vec2 uv;

varying vec2 vUV;
varying vec3 vPos;

#if defined(USE_TEXTURE)
varying mat4 TBN;
#else
varying vec3 Normal;
#endif

uniform mat4 modelMatrix;
uniform mat4 viewMatrix;
uniform mat4 projectionMatrix;

void main() {
  #if defined(USE_TEXTURE)
  vec3 T = normalize(vec3(modelMatrix * vec4(tangent, 0.0)));
  vec3 N = normalize(vec3(modelMatrix * vec4(normal, 0.0)));
  vec3 B = normalize(cross(N, T));
  TBN = mat3(T, B, N);
  #else
  Normal = normal;
  #endif

  vUV = uv;
  vPos = modelMatrix * vec4(position, 1.0);
  gl_Position = projectionMatrix * viewMatrix * vPos;
}

Fragment Shader:

#version 100
varying vec2 vUV;
varying vec3 vPos;

#if defined(USE_TEXTURE)
varying mat4 TBN;
// material textures
uniform sampler2D albedoMap;
uniform sampler2D normalMap;
uniform sampler2D metallicMap;
uniform sampler2D roughnessMap;
uniform sampler2D aoMap;

#else
// material parameters
varying vec3 Normal;
uniform vec3 albedo;
uniform float metallic;
uniform float roughness;
uniform float ao;
#endif

// lights
uniform vec3 lightPositions[4];
uniform vec3 lightColors[4];

uniform vec3 cameraPosition;

const float PI = 3.14159265359;

float DistributionGGX(vec3 N, vec3 H, float roughness) {
  float a      = roughness*roughness;
  float a2     = a*a;
  float NdotH  = max(dot(N, H), 0.0);
  float NdotH2 = NdotH*NdotH;

  float num   = a2;
  float denom = (NdotH2 * (a2 - 1.0) + 1.0);
  denom = PI * denom * denom;

  return num / denom;
}

float GeometrySchlickGGX(float NdotV, float roughness) {
  float r = (roughness + 1.0);
  float k = (r*r) / 8.0;

  float num   = NdotV;
  float denom = NdotV * (1.0 - k) + k;

  return num / denom;
}

float GeometrySmith(vec3 N, vec3 V, vec3 L, float roughness) {
  float NdotV = max(dot(N, V), 0.0);
  float NdotL = max(dot(N, L), 0.0);
  float ggx2  = GeometrySchlickGGX(NdotV, roughness);
  float ggx1  = GeometrySchlickGGX(NdotL, roughness);

  return ggx1 * ggx2;
}

vec3 fresnelSchlick(float cosTheta, vec3 F0) {
  return F0 + (1.0 - F0) * pow(clamp(1.0 - cosTheta, 0.0, 1.0), 5.0);
}

#if defined(USE_TEXTURE)
vec3 getNormalFromNormalMap() {
  // obtain normal from normal map in range [0,1]
  vec3 normal = texture(normalMap, vUV).rgb;
  // transform normal vector back to range [-1,1]
  normal = normal * 2.0 - 1.0;
  return normalize(TBN * normal);
}
#endif

void main()
{
  #if defined(USE_TEXTURE)
  vec3 albedo     = pow(texture(albedoMap, vUV).rgb, 2.2); // with Crt Gamma
  vec3 N          = getNormalFromNormalMap();
  float metallic  = texture(metallicMap, vUV).r;
  float roughness = texture(roughnessMap, vUV).r;
  float ao        = texture(aoMap, vUV).r;
  #else
  vec3 N = normalize(Normal);
  #endif

  vec3 V = normalize(cameraPosition - vPos);

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

  // reflectance equation
  vec3 Lo = vec3(0.0);
  for(int i = 0; i < 4; ++i)
    {
      vec3 L = normalize(lightPositions[i] - vPos);
      vec3 H = normalize(V + L);

      // calculate per-light radiance
      float dist        = length(lightPositions[i] - vPos);
      float attenuation = 1.0 / (dist * dist);
      vec3 radiance     = lightColors[i] * attenuation;

      // cook-torrance brdf
      float NDF = DistributionGGX(N, H, roughness);
      float G   = GeometrySmith(N, V, L, roughness);
      vec3 F    = fresnelSchlick(max(dot(H, V), 0.0), F0);

      vec3 kS = F;
      vec3 kD = vec3(1.0) - kS;
      kD *= 1.0 - metallic;

      vec3 numerator    = NDF * G * F;
      float denominator = 4.0 * max(dot(N, V), 0.0) * max(dot(N, L), 0.0) + 0.0001;
      vec3 specular     = numerator / denominator;

      // add to outgoing radiance Lo
      float NdotL = max(dot(N, L), 0.0);
      Lo += (kD * albedo / PI + specular) * radiance * NdotL;
    }

  vec3 ambient = vec3(0.03) * albedo * ao;
  vec3 color = ambient + Lo;

  color = color / (color + vec3(1.0));
  color = pow(color, vec3(1.0/2.2)); // Gamma Correction

  gl_FragColor = vec4(color, 1.0);
}

从 Camera Plane 到 Shadow Func 之间的距离: Sample Pos, 表示像素从阴影贴图上对应的深度值 \(d_{c}\);

红色虚实线表示在 fragment shader 中实时计算的深度值 \(P_{normal}.z\).

图像 Shadow Func 反映了一个区域上若干像素的实时深度值在不同位置上的变化, 以及阴影贴图深度值在不同位置上都是一致的,

这是阴影贴图是对连续对象进行采样(离散化)的结果, 对连续对象进行离散化后, 离散结果与连续对象相比必然会有精度差异的,

无法避免的, 提升采样精度也不可能解决这个问题, 这就是自遮挡产生的原因:

同一个区域的像素有一部分的实时深度值 \(P_{normal}.z\) 大于贴图上的深度值 \(d_{c}\).

因此, 可以立马想到一个解决自遮挡问题的方法: 让实时深度值 \(P_{normal}.z\) 减去一个偏移值 \(b\) 后再与贴图深度值 \(d_{c}\) 进行对比,

又或者让贴图深度值加上偏移值 \(b\) 后与实时深度值进行对比, 其中 \(b\) 是一个很小的值, 这个方法叫做 shadow bias.

但是, 这个偏移值 \(b\) 应该如何决定呢? 从图中可以发现, 偏移值大小和光线与 Geometry 之间的夹角 \(\theta \in [0, \pi]\) 有关.

当 \(\theta \to \frac{\pi}{2}\) 时, 自遮挡现象就呈消失趋势, 因为实时深度值和贴图深度值大小是一致的, 因此, \(b\) 接近于 0;

当 \(\theta \to 0\) 或 \(\theta \to \pi\), 自遮挡现象开始出现, \(|b|\) 越大.

这符合 \(cosine\) 函数的描述, 所以 \(b = \cos\theta\).

但通常只能拿到光线方向 \(l\) 和 Geometry 的法线 \(n\), 这个时候 \(b = 1 - \cos\theta_{n}\), 其中 \(\theta_{n} \in [0, \frac{\pi}{2}]\) 是 \(l\) 和 \(n\) 的夹角, \(\cos\theta_{n} = l \cdot n\).

不过 \(b \in [0, 1]\), 不满足很小的要求, 所以这里对它重新映射: \([0.0005, 0.005]\).

float ShadowCalculation(vec4 fragPosLightSpace, vec3 normal, vec3 lightDir) {
  // parameter 'normal' and 'lightDir' need to be normalized beforehand
  // perform perspective divide, be NDC
  vec3 fragPosLightNDC = fragPosLightSpace.xyz / fragPosLightSpace.w;
  // transform to [0,1] range, be UV
  vec3 fragPosLightNormal = fragPosLightNDC * 0.5 + 0.5;
  float closestDepth = texture2D(shadowMap, fragPosLightNormal.xy).r;
  float currentFragDepth = fragPosLightNormal.z;
  float bias = max(0.005 * (1.0 - dot(normal, lightDir)), 0.0005);
  float shadow = currentFragDepth - bias > closestDepth ? 1.0: 0.0;
  return shadow;
}

void main() {
  // ...
  float shadow = ShadowCalculation(fs_in.fragPosLightSpace, normal, lightDir);
  // ...
}

需要注意的是, 偏移值不能太小, 否则无法完全解决自遮挡问题.

虽然 Shadow Bias 解决了自遮蔽, 但也引入了一个新问题: 阴影发生了偏移, 这个问题被称为 Peter Panning.

由于阴影贴图的分辨率是固定且有限的, 这导致阴影有锯齿.

通常第一反应就是通过增加阴影贴图的分辨率进行解决, 又或者让光源的视椎体尽可能贴近场景.

然而这些都不是好方法, 锯齿可以通过反走样处理(anti-aliasing)进行解决, 对阴影进行反走样被成为 PCF.

PCF 有很多种实现, 最简单的实现为 Box Filtering, 也被成为盒装模糊.

以下是经过改进的实现:

float ShadowCalculation(vec4 fragPosLightSpace, vec3 normal, vec3 lightDir) {
  // parameter 'normal' and 'lightDir' need to be normalized beforehand
  // perform perspective divide, be NDC
  vec3 fragPosLightNDC = fragPosLightSpace.xyz / fragPosLightSpace.w;
  // transform to [0,1] range, be UV
  vec3 fragPosLightNormal = fragPosLightNDC * 0.5 + 0.5;
  float currentFragDepth = fragPosLightNormal.z;
  float bias = max(0.005 * (1.0 - dot(normal, lightDir)), 0.0005);
  float shadow = 0.0;
  // vec2 texelSize = 1.0 / textureSize(shadowMap, 0); // use this if textureSize is available
  vec2 texelSize = 1.0 / vec2(shadowMap_size); // shadowMap_size is the size of shadow map
  for (int x = -1; x <= 1; x++) {
    for (int y = -1; y <= 1; y++) {
      float closestDepth = texture2D(shadowMap, fragPosLightNormal.xy + vec2(x, y) * texelSize).r;
      shadow += currentFragDepth - bias > closestDepth ? 1.0: 0.0;
    }
  }
  shadow /= 9.0;
  return shadow;
}

事实上有些图形库不使用 CubeMap 作为阴影贴图, 比如 three.js 就是使用 2D 贴图作为阴影贴图.

它同样进行 6 次渲染, 只不过把每次渲染结果分别按照不同位置储存在同一张 2D 贴图里面,

读取的时候根据光照计算出 UV 坐标在进行读取. 就像是立方体映射的张开图那样.

但 three.js 没有采用这套布局, 该布局需要 \(4 \times 3\) 个 \(2 \times 2\) 规格的方块储存 6 张贴图的内容,

每个方块的纹理范围是 \([-1, 1]^2\) 的范围, 整个贴图的有效利用率只有 \(\frac{1}{2}\).

为了节省空间, 从 three.js 的 shader 函数 cubeToUV (r179) 的描述中可以发现它的布局是这样的,

以 \(-Z\) 作为中心, 只需要用到 \(4 \times 2\) 个 \(2 \times 2\) 规格的方块, 还能把有效利用率提高到 \(\frac{3}{4}\).

这是 three.js 的实现:

vec2 cubeToUV( vec3 v, float texelSizeY ) {
  // v 是射线方向, 一般经过单位化
  // texelSizeY 是纹数尺寸

  // Number of texels to avoid at the edge of each square

  vec3 absV = abs( v );

  // Intersect unit cube

  float scaleToCube = 1.0 / max( absV.x, max( absV.y, absV.z ) );
  absV *= scaleToCube;

  // Apply scale to avoid seams
  // 划分过的边界纹理坐标是不等于整数的, 比如
  //      |
  //      |     方格 A
  //      |
  //    0.99    <--- 决策边界
  //      |
  //      |     安全边界, 不做处理
  //      |
  //      1(0)  <--- 理论边界, 自身方格的结束边界, 相邻方格的起始边界
  //      |
  //      |
  //      |     方格 B
  //      |
  //    0.01
  // 需要创建一个安全边界来防止采样错误

  // 这里的安全边界宽度为 2 个纹素大小, NEAREST 滤波只需要一个纹素
  // two texels less per square (one texel will do for NEAREST)
  v *= scaleToCube * ( 1.0 - 2.0 * texelSizeY );

  // Unwrap

  // space: -1 ... 1 range for each square
  //
  // #X##               dim    := ( 4 , 2 )
  //  # #               center := ( 1 , 1 )

  vec2 planar = v.xy;

  float almostATexel = 1.5 * texelSizeY;
  float almostOne = 1.0 - almostATexel;

  if ( absV.z >= almostOne ) {

    if ( v.z > 0.0 )
      planar.x = 4.0 - v.x;

  } else if ( absV.x >= almostOne ) {

    float signX = sign( v.x );
    planar.x = v.z * signX + 2.0 * signX;

  } else if ( absV.y >= almostOne ) {

    float signY = sign( v.y );
    planar.x = v.x + 2.0 * signY + 2.0;
    planar.y = v.z * signY - 2.0;

  }

  // Transform to UV space: [0, 1]

  // scale := 0.5 / dim => 0.5 / (4, 2) = (0.125, 0.25)
  // translate := ( center + 0.5 ) / dim => ((1, 1) + 0.5) / (4, 2) = (0.375, 0.75)
  return vec2( 0.125, 0.25 ) * planar + vec2( 0.375, 0.75 );

}

官方给出不使用 CubeMap 理由在这个 2017 年的 PR 中可以找到.

以极角 \(\theta\) 和方位角 \(\phi\) 作为参数进行参数化可以把极坐标转换到斜率空间上.

做法是先把极坐标转换成笛卡尔坐标: \(\begin{equation*}\left\{ \begin{aligned} x &= \sin\theta \cos\phi \\ y &= \sin\theta \sin\phi \\ z &= \cos\theta \end{aligned} \right.\end{equation*}\).

把笛卡尔坐标转换成极坐标: \(\begin{equation*}\left\{ \begin{aligned} \theta &= \arccos(z) \\ \phi &= \arctan(\frac{y}{x}) \end{aligned} \right.\end{equation*}\).

然后把笛卡尔坐标转换到斜率空间上: \(\begin{equation*}\left\{ \begin{aligned} \widetilde{x} &= -\tan\theta \cos\phi \\ \widetilde{y} &= -\tan\theta \sin\phi \end{aligned} \right.\end{equation*}\).

从斜率空间转换会极坐标上稍微需要推导一下,

根据 40, \(|\widetilde{n}| = \sqrt{\widetilde{x}^{2} + \widetilde{y}^{2}} = \sqrt{\frac{x^{2} + y^{2}}{z^{2}}} = \frac{\sqrt{x^{2} + y^{2}}}{z} = \frac{\sqrt{\sin^2\theta (\cos^2\theta + \sin^2\theta)}}{\cos\theta} = \tan\theta\) 是 \(n\) 在 \(xy\) 平面上的投影长度, 所以 \(\theta = \arctan(\sqrt{\widetilde{x}^{2} + \widetilde{y}^{2}})\).

根据笛卡尔坐标到斜率空间上的转换可以得出 \(\tan\phi = \frac{-\widetilde{y}}{-\widetilde{x}} = \frac{-\frac{y}{z}}{\frac{x}{z}} = \frac{y}{x}\).

最终, 得到从斜率空间转换会极坐标上的关系: \(\begin{equation*}\left\{ \begin{aligned} \theta &= \arctan(\sqrt{\widetilde{x}^{2} + \widetilde{y}^{2}}) \\ \phi &= \arctan(\frac{-\widetilde{y}}{-\widetilde{x}}) \end{aligned} \right.\end{equation*}\).

1. 微表面是一个高度场, 法线被限制在球体的上半部分中. 而我们的斜率空间是通过参数化球体上半部分(\(z \gt 0\))得到的, 两者非常契合.

2. 通过调整微平面的法线向量 \(n\) 来改变微平面的方向是限制的: 不能让微平面面向球体上半部分以外的范围.

   而斜率空间是没有边界的, 因为斜率空间是通过参数化球体上半部分得到的, 所以不管怎么调整 \(\widetilde{n}\) 都不需要担心出现越界问题.

3. 很多 针对表面的线性变换都能在斜率空间上找到对应的变换, 而且两者的对应关系十分简单.

   斜率空间的这些变换通常被参数化 BRDF 模型用来实现粗糙度, 各向异性.

   • 微表面在水平方向上缩放 \((\alpha_{x}, \alpha_{y})\) 倍, 微平面 \(n\) 在斜率空间上的 \(\widetilde{n}\) 就缩放 \((\frac{1}{\alpha_{x}}, \frac{1}{\alpha_{y}})\) 倍.

     这里需要透过法线和平面垂直的关系才能看出来, 假设 \(p_{1} = (x_{1}, y_{1}, z_{1})\) 和 \(p_{0} = (x_{0}, y_{0}, z_{0})\) 是平面上的两个点,

     这两个点构成的直线 \(l = p_{1} - p_{0} = (x_{1} - x_{0}, y_{1} - y_{0}, z_{1} - z_{0}) = (x_{l}, y_{l}, z_{l})\) 与法线向量 \(n = (x, y, z)\) 垂直:

     \(n \cdot l = (x, y, z) \cdot (x_{l}, y_{l}, z_{l}) = 0\).

     让平面在水平方向进行缩放会导致直线 \(l\) 在 \(x\) 和 \(y\) 分量上缩放得到: \(l^{'} = (\alpha_{x} x_{l}, \alpha_{y} y_{l}, z_{l})\),

     为了保证 \(n \cdot l^{'} = 0\) 成立, \(n\) 需要做如下变换 \((\frac{x}{\alpha_{x}}, \frac{y}{\alpha_{y}}, z)\).

     最后把变换后的 \(n\) 转换到斜率空间上得到 \(\widetilde{n} = (-\frac{x}{z \alpha_{x}}, -\frac{y}{z \alpha_{y}}) = (\frac{\widetilde{x}}{\alpha_{x}}, \frac{\widetilde{y}}{\alpha_{y}})\).

     接下来其它的变换都是遵守类似的规律: 微表面发生变换, 从而影响微平面的法线变量 \(n\), 进而影响 \(\widetilde{n}\).

     因此, 后面会一直使用法线 \(n\) 和直线 \(l = p_{1} - p_{0}\) 作为讨论例子, 不再重复声明, 并且进一步简化描述.

   • 沿着表面的 \(z+\) 方向缩放 \(\alpha\) 倍, 微平面 \(n\) 在斜率空间上的 \(\widetilde{n}\) 就缩放 \(\alpha\) 倍.

     对微平面的 \(z\) 分量进行缩放, 得到变换后的直线 \(l^{'} = (x_{l}, y_{l}, \alpha z_{l})\).

     为了保证 \(n \cdot l^{'} = 0\), \(n\) 需要变换为 \((x, y, \frac{z}{\alpha})\), 最终 \(\widetilde{n} = (-\frac{\alpha x}{z}, -\frac{\alpha y}{z}) = (-\alpha \widetilde{x}, -\alpha \widetilde{y})\).

   • 围绕 \(z\) 轴对表面旋转 \(\theta\) 度, 斜率空间的斜率也是旋转 \(\theta\) 角.

     先假设表面的旋转矩阵为 \(M_{\mathrm{r}} = \left( \begin{array}{c} \cos\theta & -\sin\theta & 0 \\ \sin\theta & \cos\theta & 0 \\ 0 & 0 & 1 \end{array} \right)\), 表面的变换的逆矩阵是其面法线的变换矩阵.

     因为 \(M_{r}\) 旋转矩阵(旋转矩阵是正交矩阵), 所以它的转置矩阵就是它的逆矩阵.

     由此可得, 法线 \(n\) 的变换矩阵为 \(M_{r}^{-1} = M_{r}^{T} = \left( \begin{array}{c} \cos\theta & \sin\theta & 0 \\ -\sin\theta & \cos\theta & 0 \\ 0 & 0 & 1 \end{array} \right)\).

     所以, \(n\) 经过变换后为 \(M_{r}^{-1} n = \left( \begin{array}{c} \cos\theta & \sin\theta & 0 \\ -\sin\theta & \cos\theta & 0 \\ 0 & 0 & 1 \end{array} \right) \left( \begin{array}{c} x \\ y \\ z \end{array} \right) = \left( \begin{array}{c} x \cos\theta + y \sin\theta \\ -x\sin\theta + y \cos\theta \\ z \end{array} \right)\).

     最终 \(\widetilde{n} = (-\frac{x \cos\theta + y \sin\theta}{z}, -\frac{-x\sin\theta + y \cos\theta}{z}) = \left( \begin{array}{c} \cos\theta & -\sin\theta \\ \sin\theta & \cos\theta \end{array} \right) \left( \begin{array}{c} \widetilde{x} = -\frac{x}{z} \\ \widetilde{y} = -\frac{y}{z} \end{array} \right)\),

     这个矩阵 \(\left( \begin{array}{c} \cos\theta & -\sin\theta \\ \sin\theta & \cos\theta \end{array} \right)\) 就是往逆时针方向旋转 \(\theta\) 度.

     其实前面的两个变换也可以通过矩阵进行分析, 而且从矩阵角度分析的效率更高, 有兴趣的可以自行尝试.

   • 围绕 \(y\) 轴对表面旋转, 或者围绕 \(x\) 轴对表面旋转, 在斜率空间上就是透视投影.
   • 从水平方向上对微表面进行剪切变换 \(\left( \begin{array}{c} 1 & k_{2} \\ k_{1} & 1 \end{array} \right)\)

在 BRDF 模型中, NDF 就是先在斜率空间上定义的: \(D(\widetilde{n})\), 最后再把 \(D(\widetilde{n})\) 变换到其它空间上, 比如笛卡尔坐标系 \(D(n_{\mathrm{cartesian}})\) 或极坐标系上: \(D(n_{\mathrm{polar}})\).

因为斜率空间是一个无边界的 2D 平面, 所以可以在上面引入一元或二元概率密度函数, 比如说引入正态分布函数. 有时候斜率空间上的概率密度函数被称为 \(P_{22}(\widetilde{n})\) 函数.

为了把 \(P_{22}(\widetilde{n})\) 变换到其它空间上得到概率密度函数 \(P(n)\), 就得找出 \(P_{22}(\widetilde{n})\) 和 \(P(n)\) 之间的关系.

那么该如何找出这个变换呢? 可以尝试从概率密度函数的一些特性入手.

概率密度函数的值就是 \(\frac{符合特定条件的样本数量}{\text{总样本数量}}\), 对概率密度函数在特定范围内进行积分, 得到是 \(\frac{\text{特定范围内的样本数量}}{\text{总样本数量}}\).

比如以下例子,

这个比例也等于 \(\frac{\text{阴影面积}}{\text{总面积}}\), 为 \(\int_{60}^{70} f(x) dx\).

切换空间其实就和切换单位一样, 切换单位后, 函数在定义域上发生改变, 并且函数的定义也要适应定义域的改变而发生改变.

假设这个例子的函数 \(y = f(x)\) 的定义域为 \([50, 80]\), 在把单位从英寸切换成厘米后, 得到定义域为 \([50 * 2.64, 80 * 2.64] = [132, 211.2]\) 的函数 \(g(x^{'})\).

不管是什么单位, 特定身高范围内的女性占比是不会变的, 也就是说 \(\int_{60}^{70} f(x) dx = \int_{60 * 2.64}^{70 * 2.64} g(x^{'}) dx^{'}\).

你也可以这么想, 在一元函数的积分中, \(f(x)dx\) 是一个矩形的面积, 当 \(f(x)\) 转换成 \(g(x^{'})\) 后, 这个矩形的面积和 \(g(x^{'})dx^{'}\) 是一样的,

这两个矩形只是形状发生了改变而已, 面积大小没变.

在二元函数的积分中, \(f(x)dx\) (\(x\) 是二维向量) 就是体积, 也是等同于 \(g(x^{'})dx^{'}\), 就不赘述了.

因此, 我们可以得出这一个关系: \(f(x) dx = g(x^{'}) dx^{'}\).

也就是说, 在 \(x \mapsto x^{'}\) 前提下, 找出 \(f(x)\) 和 \(g(x^{'})\) 之间的关系就可以得到想要的变换了,

比如从 \(f(x)\) 变换到 \(g(x^{'})\) 上: \(f(x) dx = g(x^{'}) dx^{'} \Rightarrow g(x^{'}) = \frac{dx}{dx^{'}} f(x)\), 其中 \(\frac{dx}{dx^{'}}\) 就是从 \(f(x)\) 到 \(g(x^{'})\) 的变换.

同理, 回到 \(P_{22}\) 和 \(P\) 之间的转换, 可以得到这个变换关系 \(\begin{equation*} P_{22}(\widetilde{n}) d\widetilde{n} = P(n) dn \Rightarrow \left\{ \begin{aligned} & P(n) = \frac{d\widetilde{n}}{dn} P_{22}(\widetilde{n}) \\ & P_{22}(\widetilde{n}) = \frac{dn}{d\widetilde{n}}P(n) \end{aligned} \right. \end{equation*}\).

那么现在的关键是如何找出 \(\frac{d\widetilde{n}}{dn}\) 或 \(\frac{dn}{d\widetilde{n}}\) 呢?

\(\frac{d\widetilde{n}}{dn}\) 和 \(\frac{dn}{d\widetilde{n}}\) 刚好又是雅可比行列式(the determinant of Jacobian Matrix), 所以问题就变成了: 如何解出雅可比行列式.

这里先直接给出答案:

从斜率空间到极坐标系的雅可比行列式是 \(\frac{d\widetilde{n}}{dn} = \frac{1}{\cos^{4}\theta}\);

从极坐标系到斜率空间的雅可比行列式是 \(\frac{dn}{d\widetilde{n}} = \cos^{4}\theta\).

在得到雅可比行列式后, 就可以在斜率空间上引入想要的分布函数, 改变其参数并对函数乘以 \(\frac{d\widetilde{n}}{dn}\) 得到在论文中看到的分布函数.

至于雅可比矩阵的计算方法, 就本人所知而言有两种:

首先是 参考文章, 但该解法涉及了微分形式(differential forms)的形式论和几何代数, 本人没这方面的了解;

其次是根据变换关系直接求解, 这里主要介绍这一种方法, 具体思路如下:

首先, 没法直接从 \((x, y, z)\) 和 \((\widetilde{x} = -\frac{x}{z}, \widetilde{y} = -\frac{y}{z})\) 的关系中找到关于从 \(f\) 变换到 \(\widetilde{f}\) 的雅可比矩阵,

因为 \(f\) 和 \(\widetilde{f}\) 的定义是不确定的, 没法直接知道它们之间的关联, 不过它们都可以与极坐标之间进行变换,

所以可以分别找出极坐标到笛卡尔空间的雅可比矩阵 \(J_{f}\) 以到斜率空间的雅可比矩阵 \(J_{\widetilde{f}}\),

最后间接地得出从 \(f\) 到 \(\widetilde{f}\) 之间的雅可比矩阵行列式 \(\frac{|J_{\widetilde{f}}|}{|J_{f}|}\).

需要注意的是, 根据 NDF的定义, 还要分别找出从极坐标到 \(xy\) 平面的投影面积 \(dxdy\) 以及到 \(\widetilde{x} \widetilde{y}\) 平面的投影面积 \(d\widetilde{x} d\widetilde{y}\),

因此, 在计算 \(J_{f}\) 时不需要考虑 \(z\) 轴的变换, 最终得到一个 \(2 \times 2\) 的矩阵;

最后 \(\frac{d\widetilde{x}d\widetilde{y}}{dxdy}\) 就是我们想要的结果. 接下来开始计算.

计算从极坐标系到笛卡尔空间的 \(xy\) 平面投影的面积 \(dxdy\):

\(f(\theta, \phi) = (\sin\theta \cos\phi, \sin\theta \sin\phi, \cos\theta)\)

\(J_{f} = \left( \begin{array}{c} \frac{\partial x}{\partial \theta} & \frac{\partial x}{\partial \phi} \\ \frac{\partial y}{\partial \theta} & \frac{\partial y}{\partial \phi} \end{array} \right) = \left( \begin{array}{c} \cos\theta \cos\phi & -\sin\theta\sin\phi \\ \cos\theta \sin\phi & \sin\theta \cos\phi \end{array} \right)\)

\(|J_{f}| = \frac{\partial x}{\partial \theta} \frac{\partial y}{\partial \phi} - \frac{\partial x}{\partial \phi} \frac{\partial y}{\partial \theta} = \sin\theta\cos\theta \cos^2\phi + \sin\theta \cos\theta \sin^2\phi = \sin\theta \cos\theta\)

\(dx dy = |J_{f}| d\theta d\phi = \sin\theta \cos\theta d\theta d\phi\)

计算从极坐标系到斜率空间投影的面积 \(d\widetilde{x}d\widetilde{y}\):

\(\widetilde{f}(\theta, \phi) = (-\tan\theta \cos\phi, -\tan\theta \sin\phi)\)

\(J_{\widetilde{f}} = \left( \begin{array}{c} \frac{\partial \widetilde{x}}{\partial \theta} & \frac{\partial \widetilde{x}}{\partial \phi} \\ \frac{\partial \widetilde{y}}{\partial \theta} & \frac{\partial \widetilde{y}}{\partial \phi} \end{array} \right) = \left( \begin{array}{c} -\frac{1}{\cos^2 \theta} \cos\phi & \tan\theta \sin\phi \\ -\frac{1}{\cos^2\theta} \sin\phi & -\tan\theta \cos\phi \end{array} \right)\)

\(|J_{\widetilde{f}}| = \frac{\partial \widetilde{x}}{\partial \theta} \frac{\partial \widetilde{y}}{\partial \phi} - \frac{\partial \widetilde{x}}{\partial \phi} \frac{\partial \widetilde{y}}{\partial \theta} = \frac{1}{\cos^2 \theta} \cos^2\phi \tan\theta + \frac{1}{\cos^2 \theta} \sin^2 \phi \tan\theta = \frac{\tan\theta}{\cos^2\theta}\)

\(d\widetilde{x} d\widetilde{y} = |J_{\widetilde{f}}| d\theta d\phi = \frac{\tan\theta}{\cos^2\theta} d\theta d\phi\)

最后计算最终结果:

\(\frac{d\widetilde{x} d\widetilde{y}}{dx dy} = \frac{|J_{f}|}{|J_{\widetilde{f}}|} = \frac{\tan\theta}{\sin\theta\cos^3 \theta} = \frac{1}{\cos^4\theta}\)

比如把二元正态分布函数引入斜率空间得到: \(D(\widetilde{m}, \sigma) = \frac{1}{2 \pi \sigma^{2}} \times e^{-\frac{|\widetilde{m}|^{2}}{2 \sigma^{2}}}\),

让 \(D(\widetilde{m}, \sigma)\) 乘以 \(\frac{1}{\cos^{4}\theta}\) 就可以得到 NDF: \(D(m = (\theta, \phi), \sigma) = \frac{1}{\cos^{4}\theta}D(\widetilde{m}, \sigma) = \frac{1}{2 \pi \sigma^{2} \cos^{4}\theta} \times e^{-\frac{\tan^{2}\theta}{2 \sigma^{2}}}\).

这个 NDF 就是各向同性的 Beckmann Distribution.

\(\sigma\) 就是粗糙度, 但也有文献上用 \(\alpha\) 来替代, 看到它们其中任意一个要意识到实际上是一个东西.

你可能在其它文献中看到的 NDF 这么写: \(D(m, \alpha) = \frac{1}{\pi \alpha^{2} \cos^{4}\theta} \times e^{-\frac{\tan^{2}\theta}{\alpha^{2}}}\), 这是因为令 \(\alpha = \sqrt{2}\sigma\).

另外, 斜率空间本身就限制了 \(m\) 不可能超出上半球范围, 但是笛卡尔坐标或者极坐标并没有该限制, 所以我们需要注意 \(\theta\) 的范围:

\(\theta \in [-\frac{\pi}{2}, \frac{\pi}{2}]\), 范围外的统一不进行计算.

所以 NDF 的完整定义是: \(D(m = (\theta, \phi), \sigma) = \frac{\chi^{+}(\cos\theta)}{2 \pi \sigma^{2} \cos^{4}\theta} \times e^{-\frac{\tan^{2}\theta}{2 \sigma^{2}}}\), 其中 \(\begin{equation*}\chi^{+}(x) = \begin{cases} 1, & x \gt 0 \\ 0, & x \le 0 \end{cases}\end{equation*}\).

通过观察可以看到 \(D(m, \alpha)\) 的定义和 \(m\) 的 \(\phi\) 没有任何关系, 而 \(\theta\) 又刚好是微平面面法线 \(m\) 和宏表面法线 \(n\) 的夹角, 所以 \(D(m, \alpha)\) 又有另外两种写法:

• \(D(\theta, \alpha) = \frac{\chi^{+}(\cos\theta)}{\pi \alpha^{2} \cos^{4}\theta} \times e^{-\frac{\tan^{2}\theta}{\alpha^{2}}}\)
• \(D(m = (x_{0}, y_{0}, z_{0}), n = (x_{1}, y_{1}, z_{1}), \alpha) = \frac{\chi^{+}(m \cdot n)}{\pi \alpha^{2} (m \cdot n)^{4}} \times e^{\frac{(m \cdot n)^2 - 1}{\alpha^{2} (m \cdot n)^{2}}}\)

把二元T分布函数引入斜率空间得到: \(D(\widetilde{m}, \sigma) = \frac{\sigma^{2}}{\pi (\sigma^{2} + |\widetilde{m}|^{2})^{2}}\),

同样乘以 \(\frac{1}{\cos^{4}\theta}\) 得到 NDF: \(D(m = (\theta, \phi), \sigma) = \frac{\sigma^{2}}{\pi \cos^{4}\theta (\sigma^{2} + \tan^{2}\theta)^{2}}\).

再限制一下 \(\theta\) 的范围: \(D(m, \sigma) = \frac{\sigma^{2}\chi^{+}(\cos\theta)}{\pi \cos^{4}\theta (\sigma^{2} + \tan^{2}\theta)^{2}}\).

这个 NDF 就是各向同性的 GGX (Trowbridge-Reitz) Distribution.

另外, 有些文档会用下面的形式来表示 GGX Distribution:

\(\begin{equation*} \begin{aligned} D(m = (\theta, \phi), \sigma) &= \frac{\sigma^{2}\chi^{+}(\cos\theta)}{\pi \cos^{4}\theta (\sigma^{2} + \tan^{2}\theta)^{2}} \\ &= \frac{\sigma^{2}\chi^{+}(\cos\theta)}{\pi ((\cos^{2}\theta) \sigma^{2} + \sin^{2}\theta)^{2}} \\ &= \frac{\sigma^{2}\chi^{+}(\cos\theta)}{\pi ((\cos^{2}\theta) \sigma^{2} + (1 - \cos^{2}\theta))^{2}} \\ &= \frac{\sigma^{2}\chi^{+}(\cos\theta)}{\pi ((\cos^{2}\theta)(\sigma^{2} - 1) + 1)^{2}} \end{aligned} \end{equation*}\)

\(D(m, \sigma)\) 同样和 \(m\) 的 \(\phi\) 没有任何关系, 所以类似前面的 Beckmann Distribution 一样可以改写成另外两种表示:

• \(D(\theta, \sigma) = \frac{\sigma^{2}\chi^{+}(\cos\theta)}{\pi \cos^{4}\theta (\sigma^{2} + \tan^{2}\theta)^{2}}\)
• \(D(m = (x_{0}, y_{0}, z_{0}), n = (x_{1}, y_{1}, z_{1}), \sigma) = \frac{\sigma^{2} \chi^{+}(m \cdot n)}{\pi (m \cdot n)^{4} (\sigma^{2} - \frac{(m \cdot n)^{2} - 1}{(m \cdot n)^{2}})^{2}} = \frac{\sigma^{2}\chi^{+}(m \cdot n)}{\pi ((m \cdot n)^{2}(\sigma^{2} - 1) + 1)^{2}}\)