深入理解渲染方程

在图形学入门（三）：基础着色中，我们讨论了 Phong 反射模型，当时我们提到过 Phong 反射模型不是一个物理模型，而是一个经验模型，这意味着这个模型对光照效果的模拟是不准确的。即便在简单情况下它能近似出一些不错的效果，但随着场景的复杂度提升（例如复杂的光照、复杂的材质等），要想继续用 Phong 反射模型达到很强的真实感就变得越来越困难。例如下面的这幅图 ¹ 中，士兵和长官的铠甲上都投影出了电梯里非常复杂的灯光，在后面的长官的铠甲上还能看到前面两个士兵的投影：

要想在计算机中模拟出这样真实的效果，我们首先要知道真实世界中的光线和物体是如何作用的，例如一束光线照射到物体表面，这束光线的量有多少，光线照射到物体表面后发生了什么，最后光线反射的方向和量又有多少，我们用什么物理单位来衡量这些量等等。这些准确定义的信息是计算出正确的结果的前提，只有先理解这些内容，我们才能更好地模拟光与物体相互作用的过程。

接下来我们会讨论关于辐射度量学（Radiometry）² 的内容，这是一个研究各种电磁辐射强弱的学科。在本文后面我们将会看到，基于辐射度量学构建的渲染方程（Rendering Equation）³ 可以准确描述光如何在场景中流动，并在理论上给出了一个完美的结果。而各种各样的渲染技术，就是这个理想结果的一个近似。事实上，只要涉足高质量的实时渲染，渲染方程几乎是绕不开的前置知识。深入理解渲染方程包含的物理意义可以使我们更好地学习高质量真实感渲染的进阶内容。

在正式开始讨论前，我们首先明确要讨论的对象，也就是光的性质。我们知道，光其实有波粒二象性 ⁴，但我们后面的讨论会基于几何光学 ⁵，既不考虑其波动性，也不考虑其粒子性，仍然认为光沿直线传播。

辐射度量学 #

辐射通量（Radiant Flux） #

让我们从最简单的物理量开始讨论，辐射能（Radiant Energy）⁶ $Q$ 是指电磁辐射具有的能量，它的单位是焦耳 $\mathrm{J}$。能量这个概念我们很熟悉，但是它在我们的讨论中几乎用不上，我们更关心的是单位时间的能量，也就是功率（Power）。在辐射度量学中，我们一般称其为辐射通量（Radiant Flux）⁷ $\Phi$，单位是瓦特 $\mathrm{W}$。

$$ \Phi \equiv \frac{\mathrm{d} Q}{\mathrm{d} t} $$

我们更关心单位时间的能量而不是总能量是很符合常理的，例如一个电灯打开的时间越长肯定辐射出的光的总能量越多，但我们关心的主要还是这个灯有多亮，这也就是单位时间的光的能量。

辐射通量是辐射度量学中最基本的物理量。一个表面的辐射通量可以用「通过」该表面的坡印廷向量（Poynting Vector）⁸ $S$ 表示（这也是为什么要叫「通量」）。所谓坡印廷向量其实就是有向的能量流，它的方向为电磁能的传递方向，大小为能量密度（单位面积的能量传输速率）。在我们的讨论里，我们认为它是光线就可以了。设所求表面的表面法线为 $\hat{n}$，表面面积为 $A$，$\hat{n}$ 和 $S$ 的夹角为 $\theta$，那么表面 $\sigma$ 的辐射通量可以表示为：

$$ \Phi = \int_{\sigma} S \cdot \hat{n} \ \mathrm{d} A = \int_{\sigma} \left| S \right| \cos{} \theta \ \mathrm{d} A $$

有没有觉得这个说法听起来很熟悉？没错，这正是我们在图形学入门（三）：基础着色中计算漫反射分量时提到的内容。当时我们通过下图 ⁹ 直观了解入射光和表面法线存在 $\theta$ 夹角时，漫反射能量的减少情况：

通过刚刚的讨论，我们就可以明确其中的物理意义了，即光作为一个有向的能量流流过一个表面的时候表面通过的能量并不是光全部的能量，而和角度有关。如果光的能量流比较抽象，我们可以将其想象为水流，如果一块平面平行于水流方向，那么它不会受到水流的力的作用，随着角度倾斜，受到的力会增大，直到垂直于水流方向，此时受到的力最大。

除了辐射通量之外，我们常关注的物理量还有如下三个：

• 辐射强度（Radiant Intensity）：用于描述光源发出光的情况
• 辐照度（Irradiance）：用于描述物体表面接收到的光的情况
• 辐射率（Radiance）：用于描述一条光线

下面我们来分别讨论它们。

辐射强度（Radiant Intensity） #

辐射强度 $I$ 用于描述光源发出光的情况，它的定义是光源在每单位立体角（Solid Angle）的辐射通量 ¹⁰，单位是瓦特每球面度（Steradian）$\mathrm{W} / \mathrm{sr}$。

$$ I(\omega) \equiv \frac{\mathrm{d} \Phi}{\mathrm{d} \omega} $$

这里提到了「立体角」这个概念，这是一个物体对特定点的三维空间的角度，如下图 ¹¹ 所示：

这描述看起来就非常难理解，怎么这个角度是「一个物体相对于一个点」的角度？这和我们在二维空间中认知的角度定义很不一样，在二维空间中角度似乎和别的东西没关系，我们不会说角度是一个形状相对于一个点的角度。而且，如果这角度是个圆锥之类的形状可能还比较好理解，上图中这个物体居然具有不规则的形状，那它张开的角度到底应该怎么理解和计算呢？

理解这个概念的关键在于，我们不要把它理解为是类似二维空间那样张开的角度。而应该将计量立体角的球面度作为二维空间中弧度的延伸。在二维空间中我们是如何定义弧度的呢？是用弧长除以半径。对于半径为 $r$ 的圆上的一段圆弧 $l$ 而言，其对应的弧度 $\theta$ 表示为：

$$ \theta = \frac{l}{r} $$

一个圆的总弧度是 $2\pi$，即周长与半径的比例。对于单位圆来说，由于 $r = 1$，因此在单位圆上，我们有 $\theta = l$，我们可以说：「平面角是单位圆上的一段弧长」。如果你已经回忆起来了上面说的这个知识，那么应该可以很容易将这个概念推广到三维空间中，得出这样的定义：「立体角是单位球面上的一块面积」。对于上面的定义公式，我们也可以将「弧长」变为「面积」，将「半径」变为「半径的平方」，如下图 ¹² 所示：

那么，我们就可以得出：

$$ \Omega = \frac{A}{r^2} $$

类似于圆的总弧度，一个球的总球面度是 $4 \pi$，这也就是球的表面积 $4 \pi r^2$ 和 $r^2$ 的比值。另外可以一提的是，球面度的英文 Steradian 是希腊语「立体」（Stereos）和弧度（Radian）的组合，从字面上也可以了解到它们之间的关联。

回到辐射强度的定义公式 $I(\omega) \equiv \mathrm{d} \Phi / \mathrm{d} \omega$。从这个定义中可以知道，为了求辐射强度，我们需要求球面度的微分，那么这个微分又应该怎么求呢？根据立体角的定义，我们知道，要想求球面度的微分，就需要求球面上面积的微分。参考下图 ¹²，我们可以用球坐标 $(r,\ \theta,\ \varphi)$ 表示三维空间中的一个点 ¹³，其中，$r$ 为点到原点的径向距离（Radial Distance），$\theta$ 为点与原点连线和 $z$ 轴形成的极角（Polar Angle），$\varphi$ 为点与原点连线在 $xy$ 平面的投影与 $x$ 轴形成的方位角（Azimuthal Angle）。那么球面上面积的微分就相当于让 $\theta$ 和 $\varphi$ 都产生极小的变化（$\mathrm{d} \theta$ 和 $\mathrm{d} \varphi$），此时这个变化会在球面上形成一个极小的矩形区域。因为它在纵向上角度的变化是 $\mathrm{d} \theta$，根据弧长和弧度的关系我们知道，它的高为 $r \ \mathrm{d} \theta$，又由于这是一个单位球，因此 $r = 1$，也就是这个矩形区域实际的高度就是 $\mathrm{d} \theta$。同理，它的宽为 $r \sin{} \theta \ \mathrm{d} \varphi$。这里的宽度之所以会多一个 $r \sin{} \theta$ 是因为随着球面上的点从「极点」移动到「赤道」的过程中，球的横切面的圆的半径是不同的，这使得每张开单位大小的 $\varphi$ 角度，它在球面上对应的长度不同，在「极点」时最小，在「赤道」时最大，而这个切面的半径长度就是 $r \sin{} \theta$：

因此我们知道球面上面积的微分为：

$$ \begin{align} \mathrm{d} A &= (r \ \mathrm{d} \theta) \ (r \sin{} \theta \ \mathrm{d} \varphi) \newline &= r^2 \sin{} \theta \ \mathrm{d} \theta \ \mathrm{d} \varphi \end{align} $$

因此，球面度的微分就为：

$$ \mathrm{d} \omega = \frac{\mathrm{d} A}{r^2} = \sin{} \theta \ \mathrm{d} \theta \ \mathrm{d} \varphi $$

在辐射度量学中，我们也会用 $\omega$ 表示三维空间中的一个方向，此时这个 $\omega$ 就表示由 $\theta$ 和 $\varphi$ 确定的方向上的单位向量。

辐照度（Irradiance） #

辐照度 $E$ 用于描述物体表面接收到的光的情况，它的定义是光在入射于物体表面时，每单位面积的辐射通量 ¹⁴，单位是瓦特每平方米 $\mathrm{W} / \mathrm{m^{2}}$。

$$ E(\mathrm{x}) \equiv \frac{\mathrm{d} \Phi (\mathrm{x})}{\mathrm{d} A} $$

其中，$\mathrm{x}$ 表示物体表面上的一个点。从辐照度的定义我们也可以看出，在针对物体表面的一点计算辐照度的时候，也需要像计算辐射通量的时候一样，乘上光照方向和物体表面法线之间的夹角的余弦值。

辐照度常被称为「强度」（Intensity），但为了和前面我们讨论的「辐射强度」（Radiation Intensity）区分开来，因此起了「Irradiance」这个名字 ¹⁵。

使用辐照度可以很容易解释我们日常生活中一些常见的现象。例如，我们可以观察到灯泡的亮度会随着距离的增加而逐渐衰减，那么这个衰减的东西是什么呢？是能量消失了么？其实这里衰减的量就是辐照度。对于一个点光源来说，我们可以以其为球心作球面，那么，只要功率不变，那么每个球面上的总辐射通量就是不变的，从下图中我们可以很明显地看出，当辐射通量相同的时候，随着照亮的距离增大，由于球面的表面积增大，因此单位面积的辐射通量（即辐照度）就减少了，这就使得我们感受到灯的亮度降低。由于球面的面积和 $r^2$ 成正比，因此这个亮度的变化就与 $r^2$ 成反比：

辐射率（Radiance） #

辐射率 $L$ 用于描述光源为非点光源时，光源单位面积的辐射强度，它的定义是在指定方向上的单位立体角和垂直此方向的单位面积上的辐射通量 ¹⁶，单位是瓦特每球面度每平方米 $\mathrm{W} / (\mathrm{sr} \cdot \mathrm{m^2})$。我们看到，这里描述的是一个极小的面积上对某一个方向一个极小的夹角的辐射，这意味着事实上我们可以认为它描述的就是一条「光线」。我们可以从下图 ¹⁷ 中直观了解辐射率描述的内容：

对于点 $\mathrm{x}$ 和辐射方向 $\omega$ 而言，若表面法线方向和 $\omega$ 夹角为 $\theta$，那么我们定义辐射率 $L(\mathrm{x}, \ \omega)$ 如下：

$$ L(\mathrm{x}, \ \omega) \equiv \frac{\mathrm{d}^2 \Phi (\mathrm{x}, \ \omega)}{\mathrm{d}\omega \ \mathrm{d}A \cos{}\theta} $$

在上面的定义公式中，因为需要同时计算单位面积和单位立体角的辐射通量，因此我们要对辐射通量 $\Phi$ 进行两次微分操作，一次是面积 $A$，一次是立体角 $\omega$。另外，由于辐射率是定义在「垂直 $\omega$ 方向」上的辐射通量，因此我们需要对面积乘上 $\cos{}\theta$。

虽然辐射率的定义公式看起来就比之前定义公式的要复杂很多，但我们也不难发现它和前面我们讨论的辐射强度和辐照度的关联，我们前面提到：

• 辐照度是每单位面积的辐射通量
• 辐射强度是每单位立体角的辐射通量

因此我们也可以认为：

• 辐射率是每单位立体角的辐照度
• 辐射率是每单位面积的辐射强度

虽然准确来说，辐照度是指入射表面的辐射通量，如果我们想讨论发射的辐射通量，我们应该使用的物理概念是辐射出射度（Radiant Exitance）¹⁸，但如果我们不考虑方向性的问题，只关心的是这里的量有多少的情况下，我们可以直接认为「辐射率是每单位立体角的辐照度」，并将该「入射辐射率」$L_i$ 表示为：

$$ L_i(\mathrm{x},\ \omega) = \frac{\mathrm{d}E(\mathrm{x})}{\mathrm{d}\omega \cos{}\theta} $$

如下图所示，辐照度 $E$ 指一个单位面积上收到的辐射通量，那么如果我们考虑这个单位面积收到的来自某一个方向上的辐射通量，我们就可以得到这个方向上的入射辐射率：

类似地，我们可以说「辐射率是每单位面积的辐射强度」：

$$ L_e(\mathrm{x},\ \omega) = \frac{\mathrm{d}I(\mathrm{x},\ \omega)}{\mathrm{d}A \cos{}\theta} $$

如下图所示，辐射强度 $I$ 是每单位立体角的辐射通量。那么如果我们只考虑光源某个单位面积上对这个立体角方向上发射的辐射通量，我们就可以得到这个方向上的出射辐射率：

这样一来，我们就可以将几个物理概念联系在一起，更好地理解其中的意义。

双向反射分布函数（BRDF） #

有了上面的基本物理量的定义之后，我们就可以讨论如何更准确地来描述反射这件事了。前面提到，在 Phong 反射模型中，我们对反射进行了简单的近似处理，但这个近似不具有物理意义。这种不精确的近似手段导致了 Phong 反射模型无法达到更真实的渲染效果，尤其是对于复杂的表面材质和复杂的光照情况。那么，我们应该如何精确地描述光的反射呢？

对于反射这个过程，我们既可以将其理解为光线照射到物体表面后发生弹射，也可以将其理解为是物体表面吸收了某个方向来的入射光的能量后，向某个方向辐射出这个能量的过程。根据前面讨论的内容，我们知道我们可以使用辐照度 $E$ 来描述物体表面某个单位面积上从四面八方接收到的光，那么从 $\omega_i$ 这个方向接收到的光就可以用 $\mathrm{d}E(\omega_i)$ 表示，类似地，由于这个单位面积上接收到光后，可能会向各个方向发射光，那么由于入射 $\mathrm{d}E$ 而在 $\omega_r$ 方向发射出的反射光我们可以用 $\mathrm{d}L_r(\omega_r)$ 来描述。如果我们以这种形式理解反射，那么我们就可以将反射和前面讨论的辐照度 $E$ 以及辐射率 $L$ 的物理概念联系在一起了，如下图 ¹⁷ 所示：

我们前面提到，可以认为「辐射率是每单位立体角的辐照度」。那么根据前面讨论的公式 $L_i(\mathrm{x},\ \omega) = \mathrm{d}E(\mathrm{x}) / (\mathrm{d}\omega \cos{}\theta)$，我们可以得知图中位置 $\mathrm{x}$ 的 $\mathrm{d}E(\omega_i)$ 可以表示为：

$$ \mathrm{d}E(\mathrm{x},\ \omega_i) = L_i(\mathrm{x},\ \omega_i) \cos{}\theta_i \ \mathrm{d}\omega_i $$

而这样一个「如何在某个 $\omega_r$ 方向上分配 $\omega_i$ 方向来的入射光」方式，我们可以用一个函数 $f_r(\omega_i \rightarrow \omega_r)$ 来描述，这就是「双向反射分布函数」（Bidirectional Reflectance Distribution Function, BRDF）¹⁹。

BRDF 描述了光在一个不透明物体表面上反射时的性质，它本质上是一个四元函数，参数中的入射方向 $\omega_i$ 和反射方向 $\omega_r$ 本身是用 $\theta$ 和 $\varphi$ 表示的。除此之外，我们用 $\hat{n}$ 表示表面的法线。如下图 ¹⁷ 所示：

BRDF 的意义是在 $\omega_r$ 方向上反射光线的辐射率和同一点上从 $\omega_i$ 方向上入射的辐照度的比值，它的单位是 $1 / \mathrm{sr}$，即：

$$ f_r(\omega_i \rightarrow \omega_r) = \frac{\mathrm{d}L_r(\omega_r)}{\mathrm{d}E_i(\omega_i)} = \frac{\mathrm{d}L_r(\omega_r)}{L_i(\omega_i) \cos{}\theta_i \ \mathrm{d}\omega_i} $$

正是因为 BRDF 描述了光如何在一个表面上被反射，因此 BRDF 本身就决定了一个物体表面的材质属性，不同材质的物体表面对光的反射情况不同，就可以用不同的 BRDF 来描述。

渲染方程（Rendering Equation） #

有了 BRDF 后，我们就知道了反射光与入射光的关系，这样一来，我们终于可以讨论渲染方程的定义了。首先，我们可以将 BRDF 的公式稍作变换，得到平面上任意一点 $\mathrm{x}$ 上，$\mathrm{d}L_r$ 和 $L_i$ 的关系：

$$ \mathrm{d}L_r(\mathrm{x},\ \omega_r) = f_r(\mathrm{x},\ \omega_i \rightarrow \omega_r) \ L_i(\mathrm{x},\ \omega_i) \cos{}\theta_i $$

那么对于该位置 $\mathrm{x}$ 而言，我们可以对其法线方向所在的单位半球面 $H^2$ 上的入射光的能量进行积分（也就是将各个 $\omega_i$ 方向上的入射光的能量累加起来），就可以得到 $\omega_r$ 方向上反射光的能量：

$$ L_r(\mathrm{p},\ \omega_r) = \int_{H^2} f_r(\mathrm{x},\ \omega_i \rightarrow \omega_r) \ L_i(\mathrm{x},\ \omega_i) \cos{}\theta_i \ \mathrm{d}\omega_i $$

我们之所以只需要在半球上进行积分，是因为我们认为被求表面的是不透明的，其背面的光不对这个表面的反射产生任何贡献。上式已经告诉了我们某个方向上的因为反射而发出光的情况，我们还可以在此基础上加入自发光项 $L_e$ 以覆盖物体本身会发光的情况，让其变得更加通用，这就是渲染方程：

$$ L_o(\mathrm{x},\ \omega_o) = L_e(\mathrm{x},\ \omega_o) + \int_{H^2} f_r(\mathrm{x},\ \omega_i \rightarrow \omega_o) \ L_i(\mathrm{x},\ \omega_i) \cos{}\theta_i \ \mathrm{d}\omega_i $$

由于我们知道了物体表面的法线，因此我们会将 $\cos{}\theta_i$ 替换为 $\hat{n} \cdot \omega_i$，得到：

$$ L_o(\mathrm{x},\ \omega_o) = L_e(\mathrm{x},\ \omega_o) + \int_{H^2} f_r(\mathrm{x},\ \omega_i \rightarrow \omega_o) \ L_i(\mathrm{x},\ \omega_i) \ (\hat{n} \cdot \omega_i) \ \mathrm{d}\omega_i $$

渲染方程的物理基础是能量守恒定律，在一个特定的位置和方向，出射光 $L_o$ 是该点本身发射光 $L_e$ 与该点的反射光之和。渲染方程是渲染领域中的一个核心理论概念，在各个场景下如何去解渲染方程是高质量真实感渲染的核心挑战。由于积分在计算机中非常难算，有时甚至没有解析解，因此在实际应用中，渲染方程主要是用两类数值解法来求解。一类基于有限元法（Finite Element Method）²⁰，例如辐射度算法（Radiosity）²¹。另一类基于蒙特卡洛方法（Monte Carlo Method）²²，如路径追踪（Path Tracing）²³、光子映射（Photon Mapping）²⁴ 等。这些具体的解法超出了本文的讨论范围，但相信经过了前面的讨论，我们可以更容易去学习和理解这些进阶的知识。

参考资料 #