计算机图形学基础（8）——光线追踪

发表于 Waline：

之前我们介绍了在游戏领域中广泛使用的实时渲染技术——光栅化，本文我们来介绍一下在特效领域中广泛使用的离线渲染技术——光线追踪。

概述

那么有了光栅化渲染技术，为什么还要用光线追踪呢？根本原因在于光栅化的着色局部性，使得它无法解决很多全局效果，比如：

  • 软阴影（Soft shadows）
  • 光泽反射（Glossy reflection）
  • 间接光照（Indirect illumination）

软阴影的产生是由于光源的大小和距离。一个面积较大的光源照射一个物体，可以从多个角度投射光线到物体上。当光线从不同角度照射时，阴影边缘的重叠区域会因光线部分遮挡而出现阴影的渐变效果。

光泽反射是一种既不完全镜面反射也不完全漫反射的光反射现象。光泽度高、粗糙度低的表面产生的反射相对清晰，接近镜面反射；而光泽度低、粗糙度高的表面，反射图像就会更加模糊，趋向漫反射。

间接光照是一个场景中光线经过一次或多次反射后照射到物体表面的光照效果。在真实世界中，一个房间只通过一扇窗投过光照，但是房间里没有一个地方是完全黑色的，这就是光的多次反射造成的。

底层依据

光线追踪依赖以下几个基本的底层依据，分别是：

  • 光线沿直线传播：在微观角度，光是沿着波形传播；在宏观角度，光是沿着直线传播。
  • 光线相互不碰撞：当光线的传播路径发生交叠时，我们认为光线的传播互不干扰。
  • 光线具有可逆性：真实世界中，视觉成像是因为物体发射或反射光线，进入视网膜；在图形学中，屏幕成像是因为从相机向像素建立了反向的光线传播路径，采集路径上的所有着色信息。

基本原理

光线追踪的基本原理非常简单，其采用了 针孔相机模型（Pinhole Camera Model），分为 视线生成 和 像素着色 两个阶段。

对于视线生成，以相机为起点，以像素为锚点，利用光的可逆性，建立一条视觉射线，如下所示。此时，我们需要找到视线与空间中的物体所产生的最近的交点。

对于像素着色，当我们找到视线与物体之间的最近交点时，随即建立交点与光源的连线，判断交点是否在阴影之中，并根据结果来计算着色，写回像素结果。

光线追踪

上述基本原理只考虑了光源的直接影响，没有考虑间接影响。在实际情况中，光线会经过空间中物体的多次弹射，汇聚至交点上，从而对着色产生间接影响。

对此，经典光线追踪，也称 惠特式光线追踪（Whitted-Style Ray Tracing）或 递归光线追踪（Recursive Ray Tracing），引入了间接光照的处理。

在视线生成之后，我们找到最近的交点，然后根据光线的可逆性，找到反射光线（Reflected Ray）与折射光线（Refracted Ray）所产生的交点，依次类推，如下图所示。

之后，我们将所有交点分别与光源进行连线，判断这些交点是否在阴影之中，并根据结果计算着色，写回像素结果。当然，在反射和折射过程中会存在能量折损，因此，在计算着色时也会考虑光线折损的影响。

技术细节

光线追踪的整体原理非常简单，下面，我们来考虑其中所涉及到的一些技术细节。

交点判定原理

光线追踪中最重要的技术是交点判定，其中包含光线表示和几何表示两部分。

对于光线表示，我们只需要使用原点、方向向量即可进行表示，如下所示。

\[\begin{aligned} \vec{r}(t) = \vec{o} + t\vec{d} \end{aligned}\]

对于几何表示，在 《计算机图形学基础（6）——几何》 中，我们提到隐式几何表示和显式几何表示两种。

隐式几何的交点

隐式几何使用数学关系式表示，因此我们可以结合光线表示来求解方程组，即可计算得到交点，如下所示。其中，\(f\) 为隐式几何的数学关系式，将光线表示作为参数代入，求解 \(t\) 即可得到交点。

\[\begin{aligned} f(\vec{o} + t\vec{d}) = 0 \end{aligned}\]

显式几何的交点

显式几何直接或间接（通过参数映射的方式）定义点、线、面等元素集合，基本上最终都会转换成多边形网格来表示。这里，我们介绍平面和三角形两种情况下的交点判定技术。

平面交点判定

对于任意一个平面，我们可以使用一个向量和一个点来定义。其中，向量是法向量，平面上任意一个点与 \(\vec{p'}\) 所构成的向量都与法向量垂直，两者的点积为 0，因此得到如下所示的平面定义。

\[\begin{aligned} (\vec{p} - \vec{p'}) \cdot \vec{N} = 0 \end{aligned}\]

此时，我们再结合平面定义和光线表示得到如下所示的方程组，进而解得 \(t\) 的值，即交点。

\[\begin{aligned} (\vec{o} + t \vec{d} - \vec{p'}) \cdot \vec{N} = 0 \end{aligned}\]

三角形交点判定

对于三角形交点，我们可以通过 Moller Trumbore Algorithm 快速求解。算法利用重心坐标定义三角形内的一个点，以此得到如下关系式。

\[\begin{aligned} \vec{o} + t \vec{d} = (1 - b_1 - b_2)\vec{P_0} + b_1\vec{P_1} + b_2\vec{P_2} \end{aligned}\]

根据重心坐标的特性，当 \(1- b_1 - b_2 = 0\) 时，点在三角形内，因此关系式中存在三个变量 \(t\)、\(b_1\)、\(b_2\)。由于关系式中的向量都是三维空间中的向量，因此每个向量都由三个坐标值构成，可以进一步得到三个线性方程组，从而求解各个变量。

交点判定加速

在实际应用中，三维空间可能包含非常多的物体，每个物体可能包含非常多的三角形。假如，我们要判断光线与一个物体的交点，那么需要遍历物体的每一个三角形，如下所示。很显然，这是一个非常耗时的操作，尤其是场景中物体非常多的情况下。对此，图形学中也有一些针对交点判定加速的优化方法。

包围盒

包围盒（Bounding Box）的原理是使用一个简单的几何体包围一个复杂的物体，如果光线与包围盒没有交点，那么必然与内部物体没有交点，从而达到交点判定加速的目的。

轴对齐包围盒

包围盒通常是长方体。对于长方体，我们将它理解成 三对不同的平行面形成的交集，在实际应用中，我们会使用一种特殊的长方体，即三对面各自与坐标系的轴对齐，因此称为 轴对齐包围盒（Axis-Aligned Bounding Box，AABB）。

交点判定

很显然，接下来的问题就是如何判断光线是否与轴对齐包围盒相交。这里，我们首先考虑 2D 情况下的交点判定，然后进一步延伸至 3D 情况下的交点判定。

对于 2D 的情况，如下所示，分为三个步骤：

  • 考虑光线与 \(x_0\)、\(x_1\) 两个面的交点，根据方程组可以计算出一组 \(t_{min}\) 和 \(t_{max}\)。
  • 考虑光线与 \(y_0\)、\(y_1\) 两个面的交点，根据方程组又可以计算出一组 \(t_{min}\) 和 \(t_{max}\)。
  • 对于两组 \(t_{min}\) 和 \(t_{max}\)，我们必须保证 \(t_{min}\) 大于 0，否则表示交点在光源的背后，没有实际意义。两组 \(t_{min}\) 和 \(t_{max}\) 各自构成一个线段，我们只需要计算两者的交集，即可得到光线在长方形内的传播路径，而传播路径的两个端点正是光线与长方形的两个交点。

对于 3D 的情况，我们只需要额外延伸一步即可。

  • 考虑光线与 \(z_0\)、\(z_1\) 两个面的交点，根据方程组又可以计算出一组 \(t_{min}\) 和 \(t_{max}\)。
  • 求解三组 \(t_{min}\) 和 \(t_{max}\) 各自构成的线段，求线段的交集即可得到光线在包围盒内的传播路径，传播路径的两个端点正是光线与包围盒的两个交点。

最终求解得到的 \(t_{enter}\) 和 \(t_{exit}\) 的值可能存在一下几种情况：

  • 当 \(t_{exit}\) < 0 时，表示包围盒在光源的背后，没有交点。
  • 当 \(t_{enter}\) < 0 且 \(t_{exit}\) >= 0 时，表示光线在包围盒的内部。
  • 当且仅当 \(t_{enter}\) < \(t_{exit}\)，\(t_{exit}\) >= 0 时，光线和包围盒存在交点。

统一网格

统一网格（Uniform grids）是包围盒技术的延伸，它对包围盒内部的空间进行预处理。

  • 将包围盒拆分成统一的立体网格
  • 将与物体发生重叠的网格进行标识

在判断光线与包围盒内部物体的交点时，会先遍历所有网格。当网格与光线有交点时，会继续判断网格内的物体是否与光线有交点。

关于网格的划分，不同的划分方法，加速效果存在差异。

  • 当只有一个网格时，没有加速效果。
  • 当网格划分非常密集时，计算交点时遍历网格的开销会增大，甚至会出现性能降低的情况。

相对而言，一个加速效果良好的网格划分经验是：网格的数量等于物体的数量乘以一个常数（经验值是 27）。

然而，统一网格划分并不是万能的，它只适用于一些物体分布均匀的场景。如下所示的场景中，在餐桌上方存在大量没有物体的空间，此时使用统一网格划分的效果并不好。

空间划分

空间划分（Spatial partitions）主要是为了解决空间中物体分布不均匀的问题，其主要有 Oct-Tree、KD-Tree、BSP-Tree 等几种技术。

Oct-Tree，又称八叉树，其基本原理是把立方体分割成八等份，递归进行分割，直到格子为空或者物体足够少。

KD-Tree，与 Oct-Tree 类似，区别在于每次只将格子一分为二，总是沿着有个轴进行分割。相比于 Oct-Tree，其节点数量的复杂度不会随着维度而指数型增长。

BSP-Tree，对空间一分为二，每次选择一个方向进行分割。相比于 KD-Tree，其切割的方向并不一定与轴平行。

KD-Tree

这里我们着重介绍一下 KD-Tree 的工作原理。

首先，KD-Tree 会对包围盒空间进行预处理，对包围盒空间进行二分，同时构建二叉树，如下所示。其中，非叶子节点表示二分之前的整体空间，它不会存储物体；叶子节点表示存储物体的真实空间。

在进行交点判定时，本质上就是对二叉树进行先序遍历。对于上图中的例子，我们判定光线是否与整体的空间 A 有交点。

显然，上述例子中光线与空间 A 存在交点，那么它会节点 A 的两个叶子节点 1 和 B 分别进行交点判定，如下所示。如果交点存在则进一步判断内部物体是否与光线有交点。

下图所示，当与空间 B 存在交点时，则进一步遍历空间 B 的子节点。

以此类推，空间 C 与光线仍然存在交点，则继续遍历二叉树，直到判定光线与物体之间存在交点或不存在交点。

直观而言，KD-Tree 确实能够有效加速交点判定。但是 KD-Tree 也存在一些无法解决的问题，比如：

  • 存在无法判定的例子，比如：三角形反包围了包围盒。
  • 同一物体被划分在多个包围盒中，会有重复判断的问题。

物体划分

针对空间划分无法解决的问题，图形学采用物体划分来解决。下面，我们来介绍物体划分的典型技术——BVH。

BVH

BVH（Bounding Volume Hierarchy）的整体原理与 KD-Tree 类似，唯一的区别在于对于空间划分的方式不同。KD-Tree 是对空间进行划分，BVH 是对物体进行划分。

如下所示，空间中的物体彼此之间非常靠近。如果我们采用 KD-Tree 的划分方式，那么一个物体可能会包含在多个空间中，进而存在重复判定的情况。

对此，BVH 对物体进行划分，并重新计算包围盒，形成根节点的子节点。以此类推，最终在所有叶子节点存储物体列表。

交点渲染原理

当我们获取到光线直射、反射至物体表面的交点之后，我们可以进一步计算这些交点的着色，从而完成渲染。

光线追踪的问题

关于经典光线追踪，即 Whitted-Style Ray Tracing，对光线作出了如下假设。

  • 光线只进行镜面反射和折射
  • 光线在漫反射面停止弹射

对于第一种情况，经典光线追踪只能渲染镜面反射，无法渲染磨砂反射，如下所示。

对于第二种情况，对于漫反射物体，光线转播至表面时会停止弹射，如下所示。

既然经典光线追踪是存在问题的，那么我们该如何对交点进行渲染呢？答案是渲染方程。

渲染方程定义

关于渲染，我们在 上一篇文章 中介绍了辐射度量学，文中我们提到了渲染方程，如下所示。

\[\begin{aligned} L_r(p, {\omega}_r) = L_e(p, {\omega}_r) + \int_{H^2}L_i(p, {\omega}_i) f_r(p, {\omega}_i, {\omega}_r) (n \cdot {\omega}_i) d{\omega}_i \end{aligned}\]

渲染方程展开

渲染方程经过一系列变换可以简化成如下表达式。

\[\begin{aligned} l(u) = e(u) + \int l(v) K(u, v) dv \end{aligned}\]

我们将 \(K\) 作为反射操作符，可以进一步将简化为如下表达式，\(L\) 是一个递归项。

\[\begin{aligned} L = E + KL \end{aligned}\]

我们的目的是求解 \(L\)，前面说 \(K\) 是一个反射操作符，其实我们也可以将其理解为反射次数，反射次数越多，展开项越多，比如：\(K^2\) 表示光在空间中弹射两次。

\[\begin{aligned} L = & E + KL \\ IL - KL = & E \\ (I - K)L = & E \\ L = & (I - K)^{-1}E \\ L = & (I + K + K^2 + K^3 + ...)E \\ L = & E + KE + K^{2}E + K^{3}E + ... \\ \end{aligned}\]

我们可以非常容易地理解展开后的渲染方程中的各个项，如下所示。全局光照由 光源、直接光照、间接光照 组合而成，其中光栅化只包含了光源和直接光照两部分，难以实现间接光照；光线追踪则包含了全局光照。

下图所示，展示了全局光照包含不同数量的光照项的渲染对比效果。

求解理论基础

关于如何求解渲染方程，我们首先介绍两部分与之相关的内容。

  • 概率论基础
  • 蒙特卡洛积分

概率论基础

这里我们主要介绍概率、概率分布函数、期望等几个概念。

上图所示是一个的连续的概率分布曲线 \(X \sim p(x)\)，横坐标表示目标值 \(x\)，纵坐标表示目标值为 \(x\) 的概率 \(p(x)\)，即 概率分布函数（Probability Distribution Function，PDF），那么我们可以得到如下两个关系式：

  • 所有取值的概率之和为 \(1\)
  • 期望为概率与目标值乘积的积分 \(E[X]\)
\[\begin{aligned} \int p(x) dx = & 1; && p(x) \geq 0 \\ E[X] = & \int x p(x) dx \end{aligned}\]

蒙特卡洛积分

假设我们希望求解如下所示的不定积分，那么该如何求解？答案是蒙特卡洛积分。

蒙特卡洛积分（Monte Carlo Integration）的基本思想是：在积分域内不断采样 \(f(x)\)，不断地与 \(ab\) 构成一个个长方形，然后对所有的长方形的面积之和求平均值。由此我们可以得到离散形式的蒙特卡洛方程，如下所示。

\[\begin{aligned} F_N = \frac{1}{N} \sum^{1}_{i=1} \frac{f(x_i)}{p(x_i)} \end{aligned}\]

我们将其进一步表示为更加通用的连续形式的蒙特卡洛方程，如下所示。

\[\begin{aligned} \int_{a}^{b} f(x) dx = \frac{1}{N} \sum_{i=1}^{N} \frac{f(x_i)}{p(x_i)} \end{aligned}\]

对于蒙特卡洛积分，当 \(N\) 越大，求解的结果越精准。

渲染方程求解

在了解了概率论基础和蒙特卡洛积分之后，我们来正式求解渲染方程。

求解过程分析

假设，我们只考虑一个像素点的着色过程，如下图所示。其中，\({\omega}_o\) 为观测方向，即从着色点到观测点的方向；\({\omega}_i\) 为光线入射方向，即从光源到着色点的方向。

同时，我们先忽略着色点本身的发光项，那么可以得到一个简化的渲染方程，如下所示。

\[\begin{aligned} L_o(p, {\omega}_o) = \int_{H^2} L_i(p, {\omega}_i) f_r(p, {\omega}_i, {\omega}_o) (n \cdot {\omega}_i) d{\omega}_i \end{aligned}\]

我们将渲染等价为在半球面上进行采样。因此，任意一点的采样概率 \(p({\omega}_i)\) 为 \(1/2\pi\)，即概率分布函数（PDF）。对应，采样值 \(f({\omega}_i)\) 为 \(L_i(p, {\omega}_i) f_r(p, {\omega}_i, {\omega}_o) (n \cdot {\omega}_i)\)。由此可以用蒙特卡洛积分来求解渲染方程，如下所示。

\[\begin{aligned} L_o(p, {\omega}_o) = & \int_{H^2}L_i(p, {\omega}_i) f_r(p, {\omega}_i, {\omega}_o) (n \cdot {\omega}_i) d{\omega}_i \\ \approx & \frac{1}{N} \sum_{i=1}^{N} \frac{f({\omega}_i)}{p({\omega}_i)} \\ \approx & \frac{1}{N} \sum_{i=1}^{N} \frac{L_i(p, {\omega}_i) f_r(p, {\omega}_i, {\omega}_o) (n \cdot {\omega}_i) }{p({\omega}_i)} \end{aligned}\]

根据上述的渲染方程的求解关系式，我们可以实现对应的伪代码，如下所示。

上述，我们只处理了光的一次弹射，而全局光照还包括二次弹射、三次弹射...我们该如何处理呢？

上图所示，我们希望计算 Q 点反射到 P 点的光线。对此，我们可以转换一下思路：Q 点反射到 P 点的光线，不就等于从 P 点观测 Q 点，Q 点的着色吗？于是全局光照就转换成了一个递归着色的过程。如下所示，我们加入了处理全局光照的递归着色。

至此，着色算法的整体框架已经实现了，但是还存在两个问题：

  • 分治爆炸问题
  • 无限递归问题

分治爆炸问题

在上述算法实现中，对于每一个点的着色，我们都会采样各个方向，而每一个方向所经过的反射点又会采样各个方向。如果我们对于每个点采样 100 个方向，那么在下一次递归中将会采样 10000 个方向，依次类推，产生分治爆炸问题。

那么如何解决这个问题？很显然，只有当采样数量为 1 时，递归时才不会出现分治爆炸问题。在计算机图形学中，对于每个着色点，只通过一条射线进行光线追踪的方式被称为 路径追踪（Path Tracing）。

如下所示，为解决了分支爆炸问题的着色算法。我们不再采样各个方向，而是随机采样一个方向，进而计算着色。然而，我们知道在蒙特卡洛积分中，采样数越小，准确性越低。算法中将采样数降至 1，很显然，会出现很大的噪声。那么这又该如何解决呢？

对此，有一种方法另辟蹊径：每个像素生成多条光线，进行多次路径追踪，并求解平均值。

如下所示，我们定义了一个 ray_generation 方法，它以相机位置和像素点作为参数，内部对这个像素点进行多次路径追踪，即调用 shade 方法，最终求解平均值。

无限递归问题

对于无限递归问题，很显然，我们必须要找到一个停止递归的策略。对此，图形学借鉴了 俄罗斯轮盘赌（Russion Roulette，RR）的思想。

俄罗斯轮盘赌的基本原理是设定一个概率值 \(0 < P < 1\)。

  • 射出光线的概率为 \(P\)，由此返回的着色结果为 \(L_0/P\)
  • 不射出光线的概率为 \(1-P\)，由此返回的着色结果为 \(0\)

使用这种方式，我们仍然可以期望得到 \(L_0\)，如下所示为俄罗斯轮盘赌的期望公式.

\[\begin{aligned} E = P * (L_0 / P) + (1 - P) * 0 = L_0 \end{aligned}\]

于是，我们可以进一步优化着色算法，如下所示。

光源采样问题

经过上述改进，着色算法是正确了，但是它并不高效。如下所示，像素采样必须达到一定阈值才能得到相对高质量的结果。

这里的根本原因在于均匀采样，其中只有极少数方向能够打到光源，大多数方向则浪费掉了。

对此，我们可以考虑直接对光源进行采样。考虑到我们在求解蒙特卡洛积分时是在立体角进行采样，因此我们将空间中的面光源投影到球面上，如下所示。

我们得到立体角和光源微分面积的关系如下所示。

\[\begin{aligned} d\omega = \frac{dA cos{\theta}'}{||x' - x||^2} \end{aligned}\]

然后，我们微分立体角代入渲染方程，如下所示。

\[\begin{aligned} L_o(p, {\omega}_o) = & \int_{H^2}L_i(p, {\omega}_i) f_r(p, {\omega}_i, {\omega}_o) cos\theta d{\omega}_i \\ = & \int_{H^2} L_i(p, {\omega}_i) f_r(p, {\omega}_i, {\omega}_o) \frac{cos\theta cos{\theta}'}{||x' - x||^2} dA \end{aligned}\]

现在，我们认为着色结果来源于两部分：

  • 光源的直接作用：直接采样光源，无需俄罗斯轮盘赌
  • 光源的间接作用：反射、漫反射，需要俄罗斯轮盘赌

由此，我们结合两部分得到如下所示的优化算法。当然，对于光源的直接作用，我们需要判断光源与着色点之间是否存在遮挡，这一点一定要注意。

总结

本文我们主要介绍了光线追踪的原理。首先，我们我们介绍了光线追踪的基本原理，其主要包括：视线生成和像素着色两个阶段。

然后，我们介绍了经典光线追踪技术，也称为 Whitted-Style 光线追踪，其引入了间接光照的处理。

在光线追踪中，主要设计两个关键技术，分别是交点判定和交点渲染。对于交点判定，我们介绍了包围盒原理、统一网格、空间划分、物体划分等技术；对于交点渲染，我们在经典光线追踪的基础上引入了漫反射的处理，整体围绕渲染方程进行展开，使用蒙特卡洛积分进行求解。

在实现渲染方程算法的过程中，我们遇到了分治爆炸、无限递归、光源采样等问题，对此我们也依次进行了处理。最终实现了一个完整的着色算法。

参考

  1. 《GAMES 101》