光线追踪
光线追踪可以解决解决光栅化中没有解决的一些问题：
- 全局的光照效果不好表示；
- 软阴影效果的产生；
- 光泽反射（Glossy reflection）；
- 间接光照，在漫反射场景中，有些光线在到达眼睛前反射不止一次。
光栅化是一种快速、近似的一种渲染方式，用于实时渲染；光线追踪是比较准确但是比较慢的渲染方式，用于离线渲染。
基本光线追踪方法
我们假设光线满足下面 3 点要求：
- 光线沿直线传播；
- 光线与光线直接「无碰撞」；
- 光线从光源射入人眼中，但在路径反转的情况下，物理性质是不变的（光路的可逆性）。
光的可逆性是光线追踪的重要思想（尽管它不符合实际物理学）。
- 通过每像素投射一条光线生成图像
- 通过向灯光发送光线来检查阴影
视眼睛（相机）为一个点，光源为点光源。对于每一个像素，我们从眼睛向像素画出一条光线，如果光线和物体有交点（这个交点要求是光线上的第一个交点，因为后面的交点都会被第一个交点遮挡），将这个交点和光源进行连线，如果可以连接到光源（即连线上无遮挡），那么这一点会被光源照亮，在这一点上要计算着色。
Whitted-Style 光线追踪
从眼睛沿着像素连接一条光线，我们认为光线在传播过程中会发生反射和折射，且假设反射是完美的镜面反射。将所有反射、折射点和光源连接起来，如果可以连接，则这一点会被光源照亮，那么这一点的着色要叠加到这个像素上。对于一个像素，可能对应多个被照亮的点，需要将这些点的着色全部叠加到这个像素上。
我们认为光线在传播的过程中发生反射和折射现象。假设反射是完美的镜面反射。当我们将所有的反射（折射）点和光源连接起来。如果这一点可以被光源照亮，那么我们认为这一点的着色应当叠加在这个像素上。对于光线我们认为存在能量的消逝，不会一直反射。
在上图中，通过了像素但还没通过反射、折射点的光线称为 primary ray，经过了反射、折射点的光线称为 secondary rays，从反射、折射点连接到光源的光线称为 shadow rays.
我们需要求出光线与物体表面的交点，才能确定反射点和折射点。
光线-表面求交
光线方程
光线方程通过其光源实在的位置和一个方向向量定义，如下所示：
光线方程为： \[ \mathbf{r}\left(t\right)=\mathbf{o}+t \mathbf{d} \quad 0 \leq t<\infty \]
光线与隐式表面求交
先来看光线如何和球面求交，已知光线方程为： \[ \mathbf{r}\left(t\right) = \mathbf{o} + t\mathbf{d} \quad 0 \leq t<\infty \] 球面方程为： \[ \mathbf{p}: \left(\mathbf{p - c}\right)^2 - R^2 = 0 \]
那么光线和球面的交点满足： \[ \left(\mathbf{o} + t\mathbf{d} - \mathbf{c}\right)^2 - R^2 = 0 \] 这是一个二次方程，可以写为： \[ \begin{align*} &at^2 + bt^2 + c = 0 \\ &a = \mathbf{d} \cdot \mathbf{d}\\ &b = 2 \left(\mathbf{o} - \mathbf{c}\right) \cdot \mathbf{d}\\ &c = \left(\mathbf{o}-\mathbf{c}\right) \cdot \left(\mathbf{o}-\mathbf{c}\right) - R^2 \end{align*} \] 这个方程的解为： \[ t = \frac{-b \pm \sqrt{b^2 -4ac}}{2a} \] 解出的 \(t\) 为正实数是光线与圆就有交点。
对于一般的隐式表面，其方程为： \[ \mathbf{p} : f\left(\mathbf{b}\right) = 0 \] 带入光线方程的： \[ f\left(\mathbf{o} + t\mathbf{d}\right) = 0 \] 只要求解出的 \(t\) 为正实数，则光线与表面有交点。
光线与显示三角形面求交
对于一个封闭的曲面，我们可以通过求光线与曲面的交点的个数来判断光源实在这个曲面的内部还是内部。如果光线与曲面的交点个数为奇数，则光源在曲面内部；如果光线与曲面的交点个数为偶数，则光源在曲面外部。
三角形一点处于一个平面上，光线与三角形求交分为两步：
- 求出光线与三角形所在平面的交点；
- 判断这个交点是否在三角形内部。
对于一个平面，我们可以使用平面的法向量和平面上的一点来表示这个平面的方程：
假设平面的法向量为 \(\mathbf{N}\)，平面上的一点为 \(\mathbf{p'}\)，这平面方程为： \[ \mathbf{p}: \left(\mathbf{p} - \mathbf{p}'\right) \cdot \mathbf{N} = 0 \] 将光线方程带入平面方程得： \[ \left(\mathbf{o} + t \mathbf{d} - \mathbf{p}\right)\cdot \mathbf{N} = 0 \] 解得： \[ t = \frac{\left(\mathbf{p}' - \mathbf{o}\right) \cdot \mathbf{N}}{\mathbf{d} \cdot \mathbf{N}} \] 然后检查 \(t\) 是否为正实数。
Möller Trumbore 算法
对于 \(\triangle P_0P_1P_2\)，在三角形平面上的点满足：
\[ \mathbf{P} = (1 - b_1 - b_2) \mathbf{P_0} + b_1\mathbf{P_1} + b_2\mathbf{P_2} \]
求交点得：
\[ \mathbf{O} + t\mathbf{D} = (1 - b_1 - b_2) \mathbf{P_0} + b_1\mathbf{P_1} + b_2\mathbf{P_2} \]
解得：
\[ \begin{bmatrix} t \\ b_1\\ b_2 \end{bmatrix} = \frac{1}{\mathbf{S_1} \cdot \mathbf{E_1}}\begin{bmatrix} \mathbf{S_2} \cdot \mathbf{E_2} \\ \mathbf{S_1} \cdot \mathbf{S} \\ \mathbf{S_2} \cdot \mathbf{D} \end{bmatrix} \]
其中： \[ \begin{align*} \mathbf{E_1} &= \mathbf{P_1} - \mathbf{P_0} \\ \mathbf{E_2} &= \mathbf{P_2} - \mathbf{P_0} \\ \mathbf{S} &= \mathbf{O} - \mathbf{P_0} \\ \mathbf{S_1} &= \mathbf{D} \times \mathbf{E_2} \\ \mathbf{S_2} &= \mathbf{S} \times \mathbf{E_1} \end{align*} \] 最后检查 \(b1 \ge 0\)、\(b2 \ge 0\)、\(1 - b_1 - b_2 \ge 0\) 是否成立，成立则说明解是合理的。
光线-表面求交加速
简单的光线表面求交算法每一像素需要对每一个三角形都进行测试，然后找出深度最小的那个点，这个算法需要计算 \(\text{\#pixels} \times \text{\#traingles}\) 次，是非常慢的。
包围盒（Bounding Volumes）
将物体使用一个简单的包围盒包起来，如果光线与包围盒都没有交点，则一定与这个物体表面没有交点，所以我们可以先测试包围盒，如果光线与包围盒有交点，再测试物体。
我们常用的包围盒是轴对齐包围盒（Axis-Aligned Bounding Box, AABB），包围盒的长宽高和坐标轴都是平行的。我们认为包围盒是 3 对无限大的平面。
先来看一下 2 维情况下光线与包围盒的求交：
我们求出光线在 x 平面上的两个交点以及 y 平面上的两个交点，那么最终光线在包围盒内部的部分是这些交点区间的交集。
可以认为：
- 光线进入到三个面中，才可以进入到包围盒；
- 光线离开任意一个面，就会离开包围盒。
在 3 维的情况下，我们对每对平面求一个 \(t_{\text{min}}\) 和 \(t_{\text{max}}\)，然后令 \(t_{\text{enter}} = \max{\left\{t_{\text{min}}\right\}}\)，\(t_{\text{exit}} = \min{\left\{t_{\text{max}}\right\}}\).
如果 \(t_{\text{exit}} \lt 0\)，则包围盒在光源的「背面」；如果 \(t_{\text{exit}} \ge 0\) 且 \(t_\text{enter} \lt 0\)，则光源在包围内。
综上所述，当 \(t_\text{enter} < t_\text{exit}\) 且 \(t_\text{exit} \gt 0\) 时，光源和包围盒有交点。
使用轴对齐包围盒可以非常快速地计算出 \(t\)，例如对于垂直于 \(x\) 轴的平面：
\(t\) 求出来为： \[ t = \frac{\mathbf{p}'_x - \mathbf{o}_x}{\mathbf{d}_x} \]
空间划分
为了可以加速光线和物体求交点，我们可以使用大小相同的网格将原本比较大的包围盒进行划分。我们进行网格划分分为以下几个步骤：
- 找到场景中的包围盒；
- 将包围盒划分成一个一个小格子；
- 存储哪些小格子中包含物体（我们认为物体都是非实心的面，只记录包含面的小格子，物体内部不包含面的小格子不计入）；
- 判断光线是否和格子相交，如果光线和某个格子相交并且这个格子中包含物体，那么我们要对存储在这个网格内的所有物体进行求交。
在这个里，我们有两个假设：
- 判断光线和格子是否相交是很快的；
- 我们可以使用类似于光栅化直线的方式来判断直线与那些网格是相交的。
一般来说，格子的划分不可以太稀疏也不可以太稠密，需要对格子的量进行控制。网格的数量 \(\text{\#cells} = C * \text{\#objs}\)，在 3 维情况下，\(C\approx 27\).
这种划分方式对于物体分布均匀的场景比较合适。对于物体分布稀疏的场景需要多次和格子进行相交判断，这种方法相对不合适。
还有其他空间划分方法，包括八叉树（Oct-Tree），KD 树（KD-Tree）以及 BSP 树（BSP-Tree）。
- 八叉树（Oct-Tree），将一个包围盒切成八块（图中是二维情况，只有 4 块）。对于每一个小格子我们会继续进行划分直到小格子中没有物体或者物体的数量比较少；
- KD 树（KD-Tree），每一次都进行一次水平划分或者竖直划分，将包围盒分成两部分。可以形成一个二叉树的存储结构。水平划分和竖直划分交替进行，保证划分的空间是均匀的；
- BSP 树（BSP-Tree），每一次选择一个方向进行一次划分，并不是沿着轴平行方向划分。、
三种划分方法，KD 树更常用并且使用起来比较方便，可以用二叉树来存储。在每一个二叉树节点中，我们都要储存以下信息：
- 如果是非叶子结点，需要存储划分轴，划分的位置以及孩子节点的指针；
- 如果是叶子结点，需要存储格子中包含的物体。
实际划分出的盒子均在叶子结点上。
我们通过类似于二分查找的方式计算交点：
- 如果光线和一个节点有交点，那么它和它的子节点也有交点；
- 如果光线和叶子结点有交点，那么它需要和格子内的所有物体求交点。
这种方法存在两个问题，首先，格子和一个三角形面是否相交的判断比较复杂。其次，一个物体可能会在多个不同的格子中，需要多次存储。因此我们会使用更常用的物体划分的方式。
物体划分
物体划分（Bounding Volume Hierarchy，BVH）的主要思想是对物体进行进行划分，并重新计算包围盒。BVH 中，每一个物体只属于一个包围盒。
BVH 的划分主要分为以下几步：
- 划分包围盒；
- 将物体组成的集合划分为两个子集合；
- 计算每个子集合的包围盒；
- 当叶子结点的三角形面数量足够少的时候，停止划分。
当我们进行划分的时候，也有不同的划分技巧：
- 沿着最长的轴划分为两半（让长轴变短，分割更加均匀）；
- 取中间的物体进行划分，可以保证两边三角形数量差不多（找到第 \(k\) 个物体的算法可以在 \(\Omicron(n)\) 的时间内解决，被称作快速选择算法）；
- 当包围盒中的物体数量小于一定数量的时候停止划分。
非叶子节点存储：
- 包围盒
- 孩子节点的指针。
叶子节点存储：
- 包围盒；
- 包围盒中的物体。
BVH 的遍历：
1 | Intersect(Ray ray, BVH node) { |
辐射度量学（Basic Radiometry）
Radiant Energy And Flux
Radiant energy：指的是电磁辐射的能量（The energy of electromagnetic radiation），单位为焦耳，符号记为： \[ Q\,\left[\mathrm{J}=\text { Joule }\right] \] Radiant flux（power）：指的是单位时间内的能量（The energy emitted, reflected, transmitted or received, per unit time.），单位是瓦特： \[ \Phi \equiv \frac{\mathrm{d} Q}{\mathrm{d} t}\,\left[\mathrm{W}=\mathrm{Watt}\right]\left[\operatorname{lm}=\text { lumen }\right]^* \]
Radiant flux 是对单位时间内流过传感器的光子数目的度量。
Radiant Intensity
Radiant Intensity 指的是光源在单位立体角上的功率（The power per unit solid angle）。数学定义为： \[ I(\omega) \equiv \frac{\mathrm{d} \Phi}{\mathrm{d} \omega} \,\left[\frac{\mathrm{W}}{\mathrm{sr}}\right]\left[\frac{\mathrm{lm}}{\mathrm{sr}}=\mathrm{cd}=\text { candela }\right] \]
立体角（solid angle）
在 2 维中，弧度制的定义如下： \[ \theta = \frac{l}{r} \] 单位为 \(rad\). 整圆的角度为 \(2\pi\, rad\).
在 3 维中，立体角的定义如下： \[ \omega = \frac{A}{r^2} \]
单位为 \(sr\)，整球对应的立体角为 \(4 \pi \, sr\).
接下来我们推出单位立体角（Differential Solid Angles）的公式：
\[ \begin{align*} \mathrm{d}A &= \left(r \mathrm{d} \theta \right)\left(r \sin {\theta}\mathrm{d}\phi\right) \\ &=r^2 \sin {\theta} \mathrm{d} \theta \mathrm{d} \phi \\ \mathrm{d}\omega &= \frac{\mathrm{d} A}{r^2} = \sin \theta \mathrm{d} \theta \mathrm{d}\phi \end{align*} \]
对整个球面的单位立体角进行积分： \[ \Omega=\int_{S^2} d \omega=\int_0^{2 \pi} \int_0^\pi \sin \theta d \theta d \phi=4 \pi \]
对 \(\mathrm{d} \Phi\) 积分得：
\[ \begin{align*} \Phi & = \int_{S^2} \mathrm{d} \Phi \\ & = \int_{S^2} I \mathrm{~d} \omega \\ & =4 \pi I \\ \end{align*} \]
所以： \[ I =\frac{\Phi}{4 \pi} \]
Irradiance
Irradiance：是指单位面积上所接收的功率（The power per unit area incident on a surface point.）。 \[ E(\mathbf{x}) \equiv \frac{\mathrm{d} \Phi(\mathbf{x})}{\mathrm{d} A \cos \theta} \, \left[\frac{W}{\mathrm{~m}^2}\right]\left[\frac{\operatorname{lm}}{\mathrm{m}^2}=\operatorname{lux}\right] \] 上式中，\(\theta\) 是光线与平面法线的夹角。
Lambert’s 余弦定理可以使用 Irradiance 来解释。
假设光线一功率 \(\Phi\) 均匀向外均匀发出，一个球面上的 Irradiance 和这个球的半径成平方反比。假设在半径为 1 的球面上的 Irradiance 为 \(E\)，则在半径为 \(r\) 的球面上的 Irradiance 为： \[ \begin{align*} E' &= \frac{\Phi}{4 \pi r^2} \\ &= \frac{E} {r^2} \\ \end{align*} \]
Radiance
Radiance 是描述光在环境中发布的基本场量。
- Radiance 是与光线相关的量
- 渲染就是在就算 Radiance.
Radiance 是单位立体角、单位面积上的功率（The power emitted, reflected, transmitted or received by a surface, per unit solid angle, per projected unit area.）。 \[ L\left(\mathrm{p}, \omega\right) \equiv \frac{\mathrm{d}^2 \Phi(\mathrm{p}, \omega)}{\mathrm{d} \omega \mathrm{d} A \cos \theta} \, \left[\frac{\mathrm{W}}{\mathrm{sr} \mathrm{m}^2}\right]\left[\frac{\mathrm{cd}}{\mathrm{m}^2}=\frac{\mathrm{lm}}{\mathrm{sr} \mathrm{m}^2}=\mathrm{nit}\right] \] 我们可以从两个方面来理解 Radiance。
- 从入射角度来说（Incident Radiance），我们认为 Radiance 是单位立体角下的 Irradiance：
\[ L\left(\mathrm{p}, \omega\right) = \frac{\mathrm{d} E\left(\mathrm{p}\right)}{\mathrm{d} \omega \cos \theta} \]
Incident radiance is the irradiance per unit solid angle arriving at the surface.
- 从出射角度来说（Exiting Radiance），我们认为 Radiance 是单位面积下的 Intensity：
\[ L\left(\mathrm{p}, \omega\right) = \frac{\mathrm{d} I\left(\mathrm{p}, \omega\right)}{\mathrm{d} A \cos \theta} \]
Exiting surface radiance is the intensity per unit projected area leaving the surface.
Irradiance vs. Radiance
Irradiance: total power received by area \(\mathrm{d}A\).
Radiance: power received by area \(\mathrm{d}A\) from "direction" \(\mathrm{d} \omega\).
那么 Irradiance 可以表示为 Radiance 在所有角度上的积分： \[ \begin{align*} \mathrm{d} E(\mathrm{p}, \omega) &= L_i\left(\mathrm{p},\omega\right) \cos \theta \mathrm{d} \omega \\ E\left(\mathrm{p}\right)& = \int_{H^2} L_i\left(\mathrm{p},\omega\right) \cos \theta \mathrm{d} \omega \end{align*} \] 我们这里只对上半求进行积分，下半球方向的光线对这一点没有任何贡献。
BRDF
BRDF 介绍
BRDF 的全称是：Bidirectional Reflectance Distribution Function. 翻译过来就是：双向反射分布函数。
BRDF 是一个函数，这个函数描述了一个点上的反射：将从 \({\omega}_i\) 方向来的 radiance 转换为 \(\mathrm{d} A\) 上接收的功率 \(E\)，任何功率 \(E\) 会变成其他任意方向 \({\omega}_{0}\) 方向的 radiance.
我们可以将反射分为两步：
- 光线入射某一点，得到了一部分能量；
- 这个点将得到的能量发射出去。
Differential irradiance incoming: \[ \mathrm{d}E\left(\omega_i\right) = L\left(\omega_i\right) \cos \theta_i \mathrm{d} \omega_i \] Differential radiance exiting (due to \(\mathrm{d} E\left(\omega_i\right)\)): \[ \mathrm{d}L_r\left(\omega_r\right) \] The BRDF represents how much light is reflected into each outgoing direction \(\omega_r\) from each incoming direction.
BRDF 的定义如下： \[ f_r\left(\omega_i \rightarrow \omega_r \right) = \frac{\mathrm{d} L_r\left(\omega_r\right)}{\mathrm{d} E_i\left(\omega_i\right)} = \frac{\mathrm{d}L_r\left(\omega_r\right)}{L\left(\omega_i\right) \cos \theta_i \mathrm{d} \omega_i} \,\left[\frac{\text{1}}{\text{sr}}\right] \]
反射方程（The Reflection Equation）
对于某一个出射方向，对于的能量应该是所有入射方向光线的叠加，于是乎我们可以得到反射方程：
\[ L_r\left(\mathrm{p}, \omega_r\right) = \int_{H^2}f_r\left(\mathrm{p}, \omega_i \rightarrow \omega_r \right)L_i\left(\mathrm{P}, \omega_i\right) \cos \theta_i \,\mathrm{d} \omega_i \]
反射 radiance 依赖于一个 入射 radiance，这个入射 radiance 又依赖于另外一个点出的反射 radiance，这是一个递归的过程。
渲染方程（The Rendering Equation）
反射方程为： \[ L_r\left(\mathrm{p}, \omega_r\right) = \int_{H^2}f_r\left(\mathrm{p}, \omega_i \rightarrow \omega_r \right)L_i\left(\mathrm{P}, \omega_i\right) \cos \theta_i \,\mathrm{d} \omega_i \] 有些物体本来就会发光，于是乎我们加上发光项，便得到了渲染方程： \[ L_o\left(\mathrm{p}, \omega_o\right) = L_e\left(\mathrm{p}, \omega_o\right) + \int_{\Omega^+}L_i\left(\mathrm{p}, \omega_i\right) f_r\left(\mathrm{p}, \omega_i, \omega_o\right)\left(n \cdot \omega_i \right)\mathrm{d} \omega_i \]
注意，所有的方向都是向外的。
我们可以认为 \(L_i\) 包含其他点光源、面光源以及其他物体二次反射光线，最终结果我们使用积分的方式进行叠加。渲染方程可以简写为： \[ l\left(u\right) = e\left(u\right) + \int l\left(v\right)K\left(u, v\right) \mathrm{d} v \] 我们可以通过矩阵再一次简化公式为： \[ L = E + KL \] 其中 \(K\) 是反射操作符，上述方程是一个离散化的简单矩阵方程，\(L\)、\(E\) 都是向量，\(K\) 是 light transport matrix.
求解得： \[ L = (I - K)^{-1}E \] 泰勒展开得： \[ \begin{align*} L &= (I + K + K^2 + K^3 + \cdots)E \\ &= E + KE + K^2 E + K^3E + \cdots \end{align*} \] 每一项分别代表的是：物体直接发出的光，光源经过一次反射的光，光源经过两次反射得到的间接光照……
着色是直接光照，对应 \(E\) 和 \(KE\)；\(KE\) 对应直接光照，\(K^2E\) 对应间接光照；\(KE + K^2E + K^3E + \cdots\) 对应全局光照。
蒙特卡罗积分（Monte Carlo Integration）
对于任意一个函数，无论其表达式复杂与否，我们都希望能够求出这个函数得积分值。在微积分中我们学习过黎曼积分，即：分隔、近似、求和、取极限，这对于表达式较为简单的函数很好操作，但对于表达式比较复杂的函数甚至表达式无法写出来的函数来说是不可行的。对于表达式比较复杂或者表达式无法写出来的函数我们可以采用数值方法求出其积分值，蒙特卡洛法就是这样的一种数值方法。
蒙特卡洛积分会进行多次随机采样，将采样点对应函数值除以采样点对应的概率，然后做一个平均，将这个结果作为积分结果。
对于下面这个积分： \[ \int_a^b f\left(x\right) \,\mathrm{d} x \] 取一个随机变量： \[ X_i \sim p\left(x\right) \]
则蒙特卡洛积分为： \[ F_N = \frac{1}{N} \sum^N_{i = 1}\frac{f\left(X_i\right)}{p\left(X_i\right)} \] 例如 \(X_i\) 满足下列分布： \[ X_i\sim p\left(x\right) = \frac{1}{b - a} \] 则蒙特卡洛积分为： \[ F_N = \frac{b- a}{N} \sum^N_{i = 1}f\left(X_i\right) \] 蒙特卡洛积分对于任何一种采样分布都是成立的。只要进行采样就可以得到对应的积分。使用蒙特卡洛积分要注意两点：
- 采样的次数越多，得到的结果越准确；
- 如果对 \(x\) 进行采样，则要对 \(x\) 进行积分。
路径追踪（Path Tracing）
Whitted-style 光线追踪有两个问题：
- 只做镜面反射（折射），这样不能很好地表示 glossy 材质的物体。
- 在漫反射材质直之间没有反射
所以 Whitted-style 的路径追踪是有问题的。但渲染方程是绝对正确的： \[ L_o\left(\mathrm{p}, \omega_o\right) = L_e\left(\mathrm{p}, \omega_o\right) + \int_{\Omega^+}L_i\left(\mathrm{p}, \omega_i\right) f_r\left(\mathrm{p}, \omega_i, \omega_o\right)\left(n \cdot \omega_i \right)\mathrm{d} \omega_i \] 我们现在就要求出这个积分。
积分求解
假设我们现在下面场景中渲染一个只有直接光照的像素：
我们先忽略渲染方程中的非积分项好了： \[ L_o\left(\mathrm{p}, \omega_o\right) = \int_{\Omega^+}L_i\left(\mathrm{p}, \omega_i\right) f_r\left(\mathrm{p}, \omega_i, \omega_o\right)\left(n \cdot \omega_i \right)\mathrm{d} \omega_i \] 我们使用蒙特卡洛积分求出这个积分值： \[ \int_a^b f\left(x\right) \, \mathrm{d}x \approx = \frac{1}{N}\sum^N_{k = 1}\frac{f\left(X_k\right)}{p\left(X_k\right)} \quad X_k \sim p\left(x\right) \] \(f\left(x \right)\) 是： \[ L_i\left(\mathrm{p}, \omega_i\right) f_r\left(\mathrm{p}, \omega_i, \omega_o\right)\left(n \cdot \omega_i \right) \] pdf 选取半球上均匀分布： \[ p\left(\omega_i\right) = \frac{1}{2\pi} \] 所以积分值求出来为： \[ \begin{align*} L_o\left(\mathrm{p}, \omega_o\right) &= \int_{\Omega^+}L_i\left(\mathrm{p}, \omega_i\right) f_r\left(\mathrm{p}, \omega_i, \omega_o\right)\left(n \cdot \omega_i \right)\,\mathrm{d} \omega_i \\ & \approx \frac{1}{N}\sum_{i = 1}^N \frac{L_i\left(\mathrm{p}, \omega_i\right) f_r\left(\mathrm{p}, \omega_i, \omega_o\right)\left(n \cdot \omega_i \right)}{p\left(\omega_i\right)} \end{align*} \] 这个算法就是一个适用于直接光照的正确着色算法，写出伪代码是：
1 | shade(p, wo) |
介绍全局光照
对于间接光照，对于下图中的 \(P\) 点，如果我们想要求出来自 \(Q\) 点的光线，我们可以看作我们从 \(P\) 点看向 \(Q\) 点来自光源的光线反射得到的结果。那么以上的伪代码我们可以改写为递归的形式：
伪代码为： 1
2
3
4
5
6
7
8
9
10shade(p, wo)
Randomly choose N directions wi~pdf
Lo = 0.0
For each wi
Trace a ray r(p, wi)
If ray r hit the light
Lo += (1 / N) * L_i * f_r * cosine / pdf(wi)
Else If ray r hit an object at q
Lo += (1 / n) * shade(q, -wi) * f_r * cosine / pdf(wi)
Return Lo
介绍路径追踪
指数爆炸问题
如果我们采样 \(N\) 次，那么在经过 \(r\) 次反射后，光线的数目可以达到 \(N^r\) 条，会使计算量大大增加。
当且仅当 \(N = 1\) 时，才不会出现指数爆炸的情况。也就是说，对每个着色点只追踪一条光线。这样做会产生噪声，no problem，只需在每个像素中追踪更多路径，然后做一个平均即可。
1 | shade(p, wo) |
则会就是路径追踪，即 \(N = 1\)，对于 \(N \neq 1\) 的称为分布式光线追踪。
对于路径追踪，使用如下的光线生成算法（与光线追踪中的光线投射很像）：
1 | ray_generation(camPos, pixel) |
这个算法还有一个问题，递归永远不会结束，这样复合现实情况，但在程序里这样肯定是不行的，我们必须选取一个合适的时机让算法停下来。
递归停止问题
引入俄罗斯轮盘赌（RR）来解决这个问题。对于，假设光线有 \(P\) 的概率能射出，这时候返回 Lo / P，那么会有 1 - P 的概率不能射出，这时候返回 0 即可。
使用这种方式的到的期望是：E = P * (Lo / P) + (1 - P) * 0 = Lo.
引入 RR 后 shade 的伪代码是：
1 | shade(p, wo) |
这个算法就是正确的路径追踪算法！
但很没用效率：
如上图所示，在 Low SPP 下效果很差。
对光源采样
对于同一个点来说，光源面积大，那么我们使用较少的光线就可以接触到光源，但是如果光源面积太小，我们必须使用较多的光线才可以和光源发生接触。那么对于小光源计算量会上升。
为了解决这个问题，我们直接在光源上采样，这样就不会有光线被浪费。
假设对光源均匀采样，则： \[ \text{pdf} = \frac{1}{A} \quad(\int \text{pdf} \, \mathrm{d} A = 1) \] 但是蒙特卡洛积分的采样必须在积分域上。此时积分域发生了变换，我们需要找到 \(\mathrm{\omega}\) 和 \(\mathrm{d} A\) 之间的关系。\(\mathrm{\omega}\) 是 \(\mathrm{d} A\) 在对应单位球上的投影，因此： \[ \mathrm{d} \omega = \frac{\mathrm{d} A \cos \theta'}{\left\|x' - x\right\| ^ 2} \] 如何就可以写出对光源的渲染方程啦： \[ \begin{align*} L_o\left(\mathrm{x}, \omega_o\right) &= \int_{\Omega^+}L_i\left(\mathrm{x}, \omega_i\right) f_r\left(\mathrm{x}, \omega_i, \omega_o\right)\cos \theta \,\mathrm{d} \omega_i \\ &= \int_A L_i\left(\mathrm{x}, \omega_i\right)f_r\left(\mathrm{x}, \omega_i, \omega_o\right)\frac{\cos \theta \cos \theta'}{\left\|x' - x\right\| ^ 2} \, \mathrm{d} A \end{align*} \] 以前，我们假设光是通过均匀半球采样「意外」射出的，现在我们考虑来自两个部分的辐射：
- 直接光照，不需要做 RR;
- 间接光照，需要做 RR.
对光源采样的代码如下：
1 | shade(p, wo) |
最后一件事，考虑遮挡：
1 | shade(p, wo) |
这就是最终的路径追踪算法！
最后放一张 Cornell box 的图：
左边的是使用路径追踪的渲染图，右边是照片，可以看到渲染图和照片非常像。