Physical Based Rendering入门
基本概念(Fundamentals)
在图形学的早期发展中,业界大佬们一直在研究光照模型,各种光照模型层出不穷,但是一直苦于无法将其整合成一个统一的公式,直到1986年,James Kajiya等人发表了一篇称为《The Rendering Equation》的论文,论文标题短小却直击本质,这篇论文提出了大名鼎鼎的渲染方程(The Rendering Equation),为现代图形学奠定了基础。
光的反射模型
先简单复习一下光的反射模型,我们考虑一下一束光打到光滑表面上的情况,此时光线会反射到一个单一方向上(镜面反射),但是实际上并不存在绝对光滑的物体,真实的情况是,除了镜面反射的那部分光线,微观上物体表面会更加粗糙,会有一部分光线会被反射到各个方向上(漫反射),甚至还有有一部分进入物体,和物体的分子之间相互作用,一部分可能被物体吸收,一部分可能会穿过物体(半透明物体),另一部分可能经过多次碰撞后又重新反射出表面(次表面散射)。
一般而言,会有大部分的光线集中在靠近镜面反射方向的地方,这部分我们称之为高光反射(Specular),高光也有资料里面称为Spread,被微观粗糙的部分反射到各个方向上,这部分我们称之为漫反射(Diffuse),如图1所示。
渲染方程(The Rendering Equation)的定义
渲染方程表达了什么呢?从结论上来看非常简单,渲染方程计算了某个点在某个出射光方向上的光的总能量是多少,如图2所示。
有了渲染方程,我们就可以计算出视角方向上收到的所有能量,而颜色正是能量的反馈,因此我们就能计算出正确的颜色值。
渲染方程将所有光照模型统一成一个公式,即:
\[ L_o(\textbf x,\omega_o)=L_e(\textbf x,\omega_o)+L_r(\textbf x,\omega_o) \]
其中\(L_o(\textbf x,\omega_o)\)表示在目标点\(\textbf x\)处出射方向\(\omega_o\)上接受到的能量,\(L_e(\textbf x,\omega_o)\)表示目标点\(\textbf x\)自身在出射方向\(\omega_o\)上发射的能量(自发光),\(L_r(\textbf x,\omega_o)\)表示在目标点\(\textbf x\)处所有在出射方向\(\omega_o\)上接受到的能量(其他光源的能量)。
\(L_r(\textbf x,\omega_o)\)是这个方程最关键的部分,它的具体定义为:
\[ L_r(\textbf x,\omega_o)=\int_{\omega_i\in\Omega}f_r(\textbf x,\omega_i,\omega_o)L_i(\textbf x,\omega_i)(\omega_i\cdot\textbf n)d\omega_i \]
这个等式也有其自己的名字,称为反射方程(The Reflectance Equation),也是这篇文章要着重讲解的部分,为了理解反射方程,我们需要先了解一些基本概念。
立体角(Solid Angle)
在讲立体角之前,我们需要先复习一下弧度(radian)的概念,我们定义单位圆上弧度为\(\alpha\)对应的弧长也为\(\alpha\),那么可以推导出单位弧度为\(1/2\pi\),假设其对应半径\(r\)的圆上的弧长为\(s\),那么可以得出\(\alpha/s=2\pi/2\pi r\),从而推导出\(\alpha=s/r\),如图3所示。
同理我们将这个概念扩展到三维空间下,三维空间下的立体角称为steradian(简称sr),我们定义单位球上立体角为\(\omega\)对应的面积也为\(\omega\),那么可以推导出单位立体角为\(1/4\pi\),假设其对应半径\(r\)的球上的面积为\(A\),那么可以得出\(\omega/A=4\pi/4\pi r^2\),从而推导出\(\omega=A/r^2\),如图4所示。
立体角的优点不仅在于其形成了球面,假想光的传播形成了一个个球面,而根据能量守恒定律单位面积上的能量会随着距离的增加而减小(与\(r^2\)成反比),但是立体角上的能量却是恒定的,如图5所示。
有了立体角之后,我们就不需要考虑光沿直线传播的衰减问题,而通过对立体角微分来表达光的传播方向,因此接下来要介绍微分立体角。
根据之前立体角的讲解可以得知\(\omega=A/r^2\),那么对两边微分就可以得到\(d\omega=dA/r^2\),那么接下来的问题就转到了如何求\(dA\)上。
首先我们需要复习一下球面坐标系(Spherical Coordinate System),球面坐标系由极角\(\theta\)(polar angle)、方位角\(\varphi\)(azimuthal angle)和半径\(r\)构成,也就是\((r,\theta,\varphi)\),如图6所示。
在渲染中,我们需要将法线\(\vec n\)对齐\(z\)轴,经过目标点且垂直于法线的平面对齐\(xy\)平面,球面上某一点与法线的夹角即为极角\(\theta\),将目标点投影到\(xy\)平面上后与\(x\)轴的夹角为方位角\(\phi\)。
我们可以用球面坐标系的两个角度来表示微分立体角,然后通过极角来确定方向。假设某个微分立体角的极角为\(d\theta\)、方位角为\(d\varphi\),球面坐标为\((r,\theta,\varphi)\)。
在微分情况下,我们可以将\(dA\)近似成一个矩形,根据图7可以看到其两条边分别为\(rd\theta\)和\(r\sin\theta d\varphi\),因此可以推导出\(dA=r^2\sin\theta d\theta d\varphi\)。
这里也可以简单验证一下,我们对整个球面的\(dA\)进行积分,得到的应该是\(4\pi r^2\),如下所示:
\[ \int_{\varphi=0}^{2\pi}\int_{\theta=0}^\pi r^2\sin\theta d\theta d\varphi=2\pi r^2(-\cos\theta)\mid_0^\pi=4\pi r^2 \]
证明公式没有问题后,我们又知道\(d\omega=dA/r^2\),代入可以得到\(d\omega=\sin\theta d\theta d\varphi\)。
辐射度量学(Radiometry)
此外由于渲染方程得到的是能量,我们还需要了解光的能量和颜色之间的关系,才能更好地理解它,为此我们需要引入一个概念:辐射度量学(Radiometry)。
Radiometry衡量的是一种电磁辐射,辐射是通过波传播的,电磁波拥有不同的波长,其波长可以小到100nm以下,也可以大到上万千米以上,然而人眼所能感知到的波长范围仅在400nm-700nm之间,如图8所示。
Radiometry使用的量化单位具体如下:
Radiant Flux
表示单位时间内的辐射能量,用\(\phi\)表示,单位是\(W(Watt)\),之前以及下文提到的能量都是指这个单位。
\[ \phi=\frac{dQ}{dt} \]
Irradiance
表示单位面积上的能量,用\(E\)表示,单位是\(W/m^2\)。
\[ E=\frac{d\phi}{dA} \]
Radiant Intensity
表示单位立体角上的能量,用\(I\)表示,单位是\(W/sr\)。
\[ I=\frac{d\phi}{d\omega} \]
Radiance
表示单位面积单位立体角上的能量,用\(L\)表示,单位是\(W/(m^2sr)\)。
\[ L=\frac{d^2\phi}{dAd\omega}=\frac{dE}{d\omega} \]
在渲染方程中,需要关注的是Irradiance和Radiance这两个概念,根据渲染方程的定义,我们首先要知道目标点上总共的能量,也就是Irradiance,然后再计算这些能量分布到\(\omega_o\)方向上的能量,也就是Radiance,我们所得到的颜色信息来自Radiance。
那么渲染过程就可以这样理解:眼睛和屏幕定义了一个点和一组方向(例如,通过每个像素的光线),并且在每个方向上在眼睛处进行评估,最终得到的叠加后的Radiance就是我们绘制的颜色值。
或许看过相关资料的朋友会看到有些时候Radiance也有叫Luminance的,各种单位也有另外一套叫法,这其实并没有什么问题,之前提到过光是由许多不同波长的波构成的,但是人眼能感知到的波长是有限的,所谓的Luminance其实也属于另一个度量学称作Photometry,它指的是人眼所能感知到的部分,但是实际上这两种度量学的概念是相同的,并且它们之间的转换是一种常量转换,因此我们在渲染中可以假设能量只包含人眼感知到的那部分,从而简化模型,否则就还要考虑波长等因素了,因此这两种度量方式都是正确的。
下表列出了Radiometry和Photometry对应的量以及单位。
Radiometric Quantity: Units | Photometric Quantity: Units |
---|---|
radiant flux: \(W(Watt)\) | luminous flux: \(lm(lumen)\) |
irradiance: \(W/m^2\) | illuminance: \(lx(lux)\) |
radiant intensity: \(W/sr\) | luminous intensity: \(cd(candela)\) |
radiance: \(W/(m^2sr)\) | luminance: \(nit=cd/m^2\) |
反射方程(The Reflectance Equation)
了解了上述基础概念之后,我们重新回到反射方程,也就是这部分:
\[ L_r(\textbf x,\omega_o)=\int_{\omega_i\in\Omega}f_r(\textbf x,\omega_i,\omega_o)L_i(\textbf x,\omega_i)(\omega_i\cdot\textbf n)d\omega_i \]
等式中的积分式是指对以目标点为球心,法线为\(z\)轴的上半球面进行积分,积分的内容是半球中的每个方向反射到\(\omega_o\)方向上的能量。
其中\(L_i(\textbf x,\omega_i)\)表示某个入射方向\(\omega_i\)在目标点\(\textbf x\)处贡献的能量,\(f_r(\textbf x,\omega_i,\omega_o)\)表示BRDF,\(\omega_i\cdot\textbf n\)则是由Lambert’s Cosine Law推导而来,我们先来了解一下Lambert’s Cosine Law。
Lambert’s Cosine Law
Lambert提出单位面积上接收到的能量和入射方向与法线夹角的余弦值成正比,通过图9可以更直观地理解它。
从数学模型上来看的话,我们看图10:
由于\(d\omega\)方向上Radiance实际上是基于\(dA'\)的,假设\(d\omega_i\)方向上的微元面积为\(dA_i\),那么\(dA\)处的Irradiance实际上应该为:
\[ E=\frac{\sum_i^\infty L_i(\omega_i)d\omega_idA_i}{dA} \]
而\(dA_i\)由于白色虚线上的轴被缩放为\(\cos\theta_i\),因此\(dA_i=\cos\theta_idA\),由此代入上式可以得到:
\[ E=\sum_i^\infty L_i(\omega_i)\cos\theta_id\omega_i \]
而\(\cos\theta_i\)就是\(\omega_i\cdot\textbf n\),使用积分式来表达的话,最终可以概括为点\(\textbf x\)处接收到的整个半球方向的Irradiance为:
\[ E(\textbf x)=\int_{\omega_i\in\Omega}L_i(\textbf x,\omega_i)(\omega_i\cdot\textbf n)d\omega_i \]
推导到这一步之后,对于反射方程就只差最后一个部分了,那就是计算出点\(\textbf x\)处接受到的所有能量中有多少贡献给了目标方向\(\omega_o\)上,这里就引出了BRDF的概念。
The BRDF
BRDF全称Bidirectional Reflectance Distribution Function,即双向反射分布方程。BRDF描述了每个入射方向\(\omega_i\)的能量贡献了多少比例给出射方向\(\omega_o\)上,离散来看大致就是这样的:
\[ f_r(\textbf x,\omega_o,\omega_i)=\begin{cases} a_0 & \omega_i=\omega_0 \\ a_1 & \omega_i=\omega_1 \\ \vdots & \vdots \\ a_n &\omega_i=\omega_n \\ \end{cases} \]
上述函数中一定有\(0\leq a_i\leq1\),这里的\(a_i\)一般可以理解为反射率,即有多少能量反射了出去,至此我们通过BRDF对每个\(\omega_i\)方向定义了一个反射率,那么将其加入被积式子中,得到的就是最终\(\omega_o\)方向上的Radiance了,BRDF的数学定义为:
\[ f_r(\textbf x,\omega_o,\omega_i)=\frac{dL_o(\textbf x,\omega_o)}{dE_i(\textbf x,\omega_i)}=\frac{dL_o(\textbf x,\omega_o)}{L_i(\textbf x,\omega_i)(\omega_i\cdot\textbf n)d\omega_i} \]
BRDF描述的是材质的一种反射特性,根据材质的性质还可以将BRDF分为各向同性(Isotropic)的BRDF和各向异性(Anisotropic)的BRDF,这一部分会在微表面部分讲解。
严格来说,BRDF会受到物理性质的约束,使其遵循两点规律,即Helmholtz互易性和能量守恒。
Helmholtz互易性表示BRDF的入射方向和出射方向可以交换且函数的返回值相同,即:
\[ f_r(\textbf x,\omega_o,\omega_i)=f_r(\textbf x,\omega_i,\omega_o) \]
换个说法,我们可以考虑半球范围\(\Omega\)的入射能量如何分布到某个出射方向\(\omega_o\)上,也可以考虑某个入射方向\(\omega_i\)上的能量如何分布到半球范围\(\Omega\)中,而它们的分布方式是相同的。然而实际上在实时渲染中,这个性质一般并没有完全遵守而只是近似,这样就已经能模拟出比较真实的结果了。
能量守恒意味着出射能量不能大于入射能量(不包括自发光),离线渲染算法,如路径追踪,需要准确的能量守恒以确保收敛性,然而对于实时渲染,近似的能量守恒很重要,但不需要完全精确,如果使用违反能量守恒的BRDF渲染表面,它将显得不够真实。
结合BRDF我们可以计算一些其他的属性,比如半球-方向反射率(Hemispherical-Directional Reflectance),其测量了半球范围内的能量经过目标点后反射到某个方向上的能量损失,其数学定义为:
\[ R(\omega_o)=\int_{\omega_i\in\Omega}f_r(\textbf x,\omega_i,\omega_o)(\omega_i\cdot\textbf n)d\omega_i \]
同样的也有方向-半球反射率(Directional-Hemispherical Reflectance),其测量了某个方向经过目标点后反射到半球范围内的能量损失,其数学定义为:
\[ R(\omega_i)=\int_{\omega_o\in\Omega}f_r(\textbf x,\omega_i,\omega_o)(\omega_o\cdot\textbf n)d\omega_o \]
很显然,如果BRDF是互易的,那么这两种反射率相同:\(R(\omega_i)=R(\omega_o)\)。
我们也可以使用球面坐标系的极角和方位角来表示反射方程,我们前面已经推导过\(d\omega=\sin\theta d\theta d\phi\),代入反射方程可以得到:
\[ \int_{\varphi_i=0}^{2\pi}\int_{\theta_i=0}^{\pi/2} f_r(\textbf x,\theta_o,\varphi_o,\theta_i,\varphi_i)L_i(\textbf x,\theta_i,\varphi_i)\cos\theta_i\sin\theta_id\theta_id\varphi_i \]
这样整个模型就变成了考虑极角和方位角的情形,之后推导其他公式为了简洁会将\(\textbf x\)省略。
我们可以再结合图11理解一下BRDF的含义。
我们接下来来讲一个最简单的BRDF:Lambertian BRDF。
Lambertian BRDF是一种最简单的光照模型,它假设光线击中表面后向所有方向均匀反射,所以BRDF是个常数,我们定义目标点处的半球-方向反射率为\(\rho_{ss}\),则有:
\[ \begin{align*}&\rho_{ss}=\int_{\varphi_i=0}^{2\pi}\int_{\theta_i=0}^{\pi/2} f_r\cos\theta_i\sin\theta_id\theta_id\varphi_i\\\Rightarrow &\rho_{ss}=2\pi f_r\int_{\theta_i=0}^{\pi/2}\sin\theta_id(\sin\theta_i)\\\Rightarrow &\rho_{ss}=2\pi f_r(\frac12\sin^2\theta)|_0^{\pi/2}\\\Rightarrow &f_r=\frac{\rho_{ss}}{\pi}\end{align*} \]
这样就推导出了Lambertian BRDF为:
\[ f_{lambertian}=\frac{\rho_{ss}}{\pi} \]
微表面理论(Microfacet Theory)
由于微观结构的尺寸远小于像素,因此无法直接建模,我们需要对其进行统计建模,而微表面理论(Microfacet Theory)则是基于微观几何的一种模型,即我们把目标点视为由许许多多个微表面构成的,每个微表面拥有自己的法线,并各自独立地反射光线,微表面模型如图12所示:
一般我们用\(\textbf m\)表示微表面法线,入射角和反射角用\(\theta_m\)表示。
对于大多数表面,微几何表面的法线分布是连续的,并在宏观表面法线处具有强烈的峰值,这个分布的紧密程度是由表面的粗糙度(roughness)决定的。
粗糙度(Roughness)
宏观上的粗糙度比较容易理解,带来的视觉效果具有很明显的凹凸感,随着观测尺度越来越小,小到比一个像素还小的时候,这个时候粗糙度会带来怎样的视觉效果呢?可以想象出许多光线在这个像素内受粗糙度的影响而反射到各个方向,宏观和微观的光线反射方式如图13所示。
而随着微观粗糙度的增加,反射的光线会变得越分散,也就会导致反射出去的环境信息更模糊更暗,如图14所示:
假想有一个小而亮的光源,设置一个粗糙度的值,并将衡量它的尺度逐渐减小,得到的高光光斑如图15所示:
粗糙度正是微表面法线分布的一种体现,这也体现了研究微观法线分布的重要性。
Masking, Shadowing, Interreflecting
在微观几何中,除了要关注法线分布,还需要考虑到其他效应。微表面的分布形式会产生一些遮挡关系,如图16所示:
正常来说,我们对微表面建模时这三点都是需要考虑的,不过由于interreflecting比较复杂,很多BRDF中只计算了shadowing和masking这两部分,尽管这可能使最终的结果变暗。
如果微观几何的高度和表面法线之间存在相关性,那么shadowing和masking会有效地影响法线分布,例如想象一个表面,其中凸起的部分非常平滑(比如被风化),而较低的部分仍然非常粗糙,这个时候在斜视角度下,表面的较低部分往往会被遮挡从而导致表面变得更加平滑,如图17所示:
对于所有表面类型,随着入射角度\(\theta_i\)的增加,表面可见的不规则性(Irregularity)会减小,这些效应与菲涅尔效应(以后的分享中会讲)相结合,使表面在观察和照明角度接近90度时呈高度反射状态。
法线分布函数(Normal Distribution Function)
微表面模型的一个重要属性是微表面法线\(\textbf m\)的统计分布,称为法线分布函数(Normal Distribution Function, NDF),用\(D(\textbf m)\)表示。
虽然\(D(\textbf m)\)表示的是法线\(\textbf m\)的分布比例,但为了能兼容到宏观的BRDF中,需要做一些规定:所有\(D(\textbf m)\)在宏观表面上的投影面积之和必须是\(dA\),也就是说\(D(\textbf m)\)需要满足:
\[ \int_{\textbf m\in\Theta}D(\textbf m)(\textbf n\cdot\textbf m)d\textbf mdA=dA\\ \Rightarrow \int_{\textbf m\in\Theta}D(\textbf m)(\textbf n\cdot\textbf m)d\textbf m=1 \]
根据这个条件,我们可以将\(D(\textbf m)\)理解为\(\textbf m\)对应的微表面的面积与宏观表面面积的比值,如图18所示:
同理我们可以推导出:对于任意方向\(\textbf v\),微表面和宏观表面在垂直于\(\textbf v\)的平面上的投影是相等的,即:
\[ \int_{\textbf m\in\Theta}D(\textbf m)(\textbf v\cdot\textbf m)d\textbf m=\textbf v\cdot\textbf n \]
我们可以结合图19来理解:
NDF一般以粗糙度作为参数来构建,不同NDF模型中的粗糙度参数的含义一般都不太一样,而且并不那么线性,一般需要通过二次拟合来达到均匀感知的效果。
NDF根据其形状可以分为各向同性的(Isotropic)和各向异性的(Anisotropic),各向同性的NDF具有旋转对称性,也就是外观只和极角有关而与方位角无关,大部分的物体表面都是各向同性的,而各向异性的NDF一般具有偏好的某个法线方向,使得我们在不同的方位角位置观察物体呈现出来的颜色不一样。
各向同性的NDF只需要一个粗糙度参数就能定义,各向异性的NDF则需要切线方向和副切线方向上的两个粗糙度参数定义,具体的NDF会在以后的分享中阐述。
几何函数(Geometry Function)
在渲染时我们只关心可见的微表面,即多个重叠微表面中最靠近相机的那部分,为了计算宏观表面上的实际生效面积,我们可以定义一个几何函数\(G_1(\textbf m,\textbf v)\)(Masking Function),它表达的是法线为\(\textbf m\)的微表面的可见部分在视线\(\textbf v\)上的投影面积占总投影面积的比例,我们将其加入球面积分中,得到的即为宏观表面在垂直于\(\textbf v\)的平面上的投影面积比例:
\[ \int_{\textbf m\in\Theta}G_1(\textbf m,\textbf v)D(\textbf m)(\textbf v\cdot\textbf m)^+d\textbf m=\textbf v\cdot\textbf n \]
由于微表面的背面我们是看不到的,因此\(\textbf v\cdot\textbf m\)需要将小于零的部分裁减掉,在式子中表示为\((\textbf v\cdot\textbf m)^+\),\(G_1(\textbf m,\textbf v)D(\textbf m)\)这一部分我们将其称为可见法线分布(distribution of visible normals),可以结合图20理解:
同理我们可以推导出\(G_1(\textbf m,\textbf l)\)(Shadowing Function),它表达的是法线为\(\textbf m\)的微表面的可见部分在视线\(\textbf v\)上的投影面积占总投影面积的比例:
\[ \int_{\textbf m\in\Theta}G_1(\textbf m,\textbf l)D(\textbf m)(\textbf l\cdot\textbf m)^+d\textbf m=\textbf l\cdot\textbf n \]
实际上我们需要同时考虑masking和shadowing的情况,我们用\(G_2\)函数表示它,假设masking和shadowing是相互独立的,那么我们可以简单地通过累乘来计算出\(G_2\):
\[ G_2(\textbf m,\textbf l,\textbf v)=G_1(\textbf m,\textbf l)G_1(\textbf m,\textbf v) \]
但实际情况往往是masking和shadowing存在重叠部分,此时\(G_2\)表达的就是同时考虑masking和shadowing的可见比例。
由于masking和shadowing的可见法线分布的积分各自不可能超过\(\textbf n\cdot\textbf l\)和\(\textbf n\cdot\textbf v\),因此有:
\[ \int_{\textbf m\in\Omega}G_2(\textbf m,\textbf l,\textbf v)D(\textbf m)(\textbf m\cdot\textbf l)^+(\textbf m\cdot\textbf v)^+d\textbf m\leq|\textbf n\cdot\textbf l||\textbf n\cdot\textbf v| \]
基于微表面的BRDF
最后我们基于微观几何定义一个微表面的BRDF:\(f_\mu(\textbf m,\textbf l,\textbf v)\),这个BRDF需要满足:
\[ f_r(\textbf l,\textbf v)=\int_{\textbf m\in\Omega}f_\mu(\textbf m,\textbf l,\textbf v)d\textbf m \]
这个时候我们再考虑可见法线分布的比例进一步过滤\(\textbf m\),根据给定的NDF \(D(\textbf m)\)和\(G_2(\textbf m,\textbf l,\textbf v)\),推导得到的宏表面BRDF为:
\[ f_r(\textbf l,\textbf v)=\int_{\textbf m\in\Omega}f_\mu(\textbf m,\textbf l,\textbf v)G_2(\textbf m,\textbf l,\textbf v)D(\textbf m)\frac{(\textbf m\cdot\textbf l)^+}{|\textbf n\cdot\textbf l|}\frac{(\textbf m\cdot\textbf v)^+}{|\textbf n\cdot\textbf v|}d\textbf m \]
之前提到过可见法线分布的最大值为\(|\textbf n\cdot\textbf l||\textbf n\cdot\textbf v|\),因此我们这里计算比例需要除以这个值。这个式子就是基于微观几何推导出的BRDF公式了。
总结(Summary)
通过了解辐射度量学和微分立体角的概念,我们了解了宏观概念下的渲染方程,并引出了BRDF的概念,这样我们就可以将各种光照模型的焦点转移到各自的BRDF中,然后我们又通过微观几何的概念建立了微表面理论,推导出了基于微表面理论的BRDF,这个时候又将焦点转移到了NDF和几何函数上面,渲染方程的伟大之处便在于此,我们将整个光照模型抽象成了多个子模型,让我们能专注在建立子模型来达到自己想要的渲染效果。