路径追踪是光线追踪的方法之一，可以做到几乎 100%的真实（照片级真实感）

1 渲染方程与全局光照

[[02 PBR理论#渲染方程和反射方程|渲染方程如下：]]

$$L_o(p,\omega_o) = L_{e}(p,\omega_{o})+\int\limits_{\Omega} f_r(p,\omega_i,\omega_o) L_i(p,\omega_i) n \cdot \omega_i d\omega_i$$

直接光

点光源的反射方程：
点光源只有一个方向有入射光，所以不用积分，多个点光源直接累加即可。

面光源的反射方程：
面光源是点光源的集合，因此对面光源所在立体角进行积分可以得出面光源的反射方程

间接光

光源不一定都是直接光，还有反射来的间接光。将其他物体反射过来的光当成光源，得到渲染方程

• $w_i$ 本来规定从点 $x$ 指向光源，是为了方便光照计算。这里不使用这个规定。将 $w_i$ 规定为视为从反射点指向点 $x$，我们将这个反射来的光也参与光照计算，那么在计算的时候需要将其取负号。
• 只需要计算反射光，其余参数均已知

我们可以将方程简写：
$$
I(u)=e(u)+\int I(v)\boxed{K(u,v)dv}
$$
写成算子形式：
$$
L=E+KL
$$

• 在一个场景内能量守恒
• $L$ ：反射能量
• $E$ ：直接光辐射的能量
• $K$：反射符，可以将光源反射能量 $L$ 再次反射
• $KL$：光源反射能量 $L$ 再次反射

解出： $L=E+KL→L=(I-K)^{-1}E$

通过二项式定理得到下面的式子
$$
L=(1+K+K^2+K^3+\ldots)E
$$
$$
L=E+KE+K^2E+K^3E+\ldots
$$

• $L=E$ 为直接光辐射的能量
• $L=E+KE$ 为直接光辐射的能量+第一次反射的能力 （第一次反射为直接光照反射，之后的就是间接光反射）
• $L=KE+K²E$ 为直接光辐射的能量+一次反射的能量+二次反射的能量
• 以此类推，所有弹射加起来就是全局光照

2 概率论基础知识

$X$：随机变量。 表示潜在值的分布
$X\sim p (x)$：概率密度函数。 描述随机过程选择值的相对概率

随机变量的期望值： 从随机分布中反复抽取样本所得到的平均值。

概率密度函数（PDF）： 一个随机变量 $X$，它可以取一组连续值中的任意一个，其中某个特定值的相对概率由连续概率密度函数 $p (X)$ 给出。


求函数的期望：

3 蒙特卡洛积分 Monte Carlo Integration

有些积分没有解析式，那么怎么求解呢？
用数值方法计算一个复杂定积分的近似值

蒙特卡洛积分是一种数值积分技巧，允许我们使用有限数量的随机样本估计任何积分。此外，蒙特卡洛保证收敛到正确的解决方案——采样越多，效果越好。

原理：对函数的积分域多次采样求均值，作为积分的近似值
每次采样的积分域为宽为 $b-a$，高为 $f (x_i)$ 的长方形。

先看下图：

这幅图表明了什么意思呢？我们知道，计算 $[a,b]$ 内的定积分就是求曲线 $f(x)$、直线 $x=a,x=b$ 以及 $x$ 轴围成的形状的面积，因此，如果我们在曲线上随机地选取 $N$ 个点，计算如图所示的粉红色长方形面积之和，再求个平均，其实就得到了定积分的近似值。点的数量取得越多，这个平均值就越逼近定积分的真实值。

数学表达如下：
$$
\frac1N[(b-a)\times f(X_1)+(b-a)\times f(X_2)+\cdots+(b-a)\times f(X_N)]=\frac{b-a}N\sum_{i=1}^Nf(X_i)
$$

结果即蒙特卡洛积分公式：
$$
F^N≈\frac1N\sum_{i=1}^N\frac{f(X_i)}{p(Xi)}
$$
这个式子没有积分符号 $∫$，但是它却叫做 “积分” 公式，这是因为这个式子求的是积分的近似值——当 N 越大的时候，计算出的值就越接近定积分的真实值。

$p(Xi)$ 是概率密度函数，它表示 $Xi$ 这个点，在某个分布下取 $Xi$ 这个值的概率。从图中很明显看到，我们取的 $X_i$ 在 $[a, b]$ 范围内，即取 $Xi$ 这个值的概率为 $\frac{1}{b-a}$。所以 $p(Xi)=\frac{1}{b-a}$

4 路径追踪 Path Tracing

Witted-Style 光线追踪的局限

Witted-Style 中不正确的部分：

• 无法表现 Glossy 的材质（Glossy 材质上的反射方向不稳定，不符合入射角等于出射角的关系）
• 漫反射后仍然会多次反射，需要引入全局光照

路径追踪解决了这些问题

红色墙面的漫反射，反射到左边的柱子，使其变红。（color bleeding）

蒙特卡洛积分解渲染方程

$$
L_o(p,\omega_o)=L_e(p,\omega_o)+\int_{\Omega^+}L_i(p,\omega_i)f_r(p,\omega_i,\omega_o)(n\cdot\omega_i)\mathrm{d}\omega_i
$$
计算全局光照：

1. 求解半球上的积分
2. 递归执行

使用蒙特卡洛方法解积分，即将某个渲染点 $p$ 的 $\omega_i$ 方向（作为随机变量）进行多次采样找到反射光源的角度得到蒙特卡洛积分：

我们想要计算该积分$$
L_o(p,\omega_o)=\int_{\Omega^+}L_i(p,\omega_i)f_r(p,\omega_i,\omega_o)(n\cdot\omega_i)\mathrm{d}\omega_i
$$

1. 根据蒙特卡洛积分：$\displaystyle\int_{a}^{b}f(x)\mathrm{d}x\approx\frac{1}{N}\sum_{k=1}^{N}\frac{f(X_{k})}{p(X_{k})}\quad X_{k}\sim p(x)$
2. 找到 $f(x)$：$\displaystyle L_i(p,\omega_i)f_r(p,\omega_i,\omega_o)(n\cdot\omega_i)$
3. 找到 $p(x)$ (概率密度函数)：$p(\omega_{i})=\frac{1}{2\pi}$ 。因为半球的立体角为 $2\pi sr$，我们取的 $X_i$ 在 $[0, 2\pi]$ 范围内，即取 $Xi$ 这个值的概率为 $\frac{1}{2\pi}$。所以 $p(Xi)=\frac{1}{2\pi}$
4. 带入原积分式：（注意该式中 $p$ 是着色点，$p(\omega_{i})$ 是概率密度函数）$$
   \begin{aligned}
   L_{o}(p,\omega_{o})& =\int_{\Omega^{+}}L_{i}(p,\omega_{i})f_{r}(p,\omega_{i},\omega_{o})(n\cdot\omega_{i})\mathrm{d}\omega_{i} \
   &\approx\frac1N\sum_{i=1}^N\frac{L_i(p,\omega_i)f_r(p,\omega_i,\omega_o)(n\cdot\omega_i)}{p(\omega_i)}\
   &\approx\frac1N\sum_{i=1}^N{2\pi L_i(p,\omega_i)f_r(p,\omega_i,\omega_o)(n\cdot\omega_i))}
   \end{aligned}

$$

全局光照算法

直接光照的算法

1
2
3
4
5
6
7
8

shade(p, wo)
    Randomly choose N directions wi~pdf //随机选择N个wi方向
    Lo = 0.0 //初始化L0
    For each wi //对每个wi方向
        Trace a ray r(p, wi) //从p点往wi发出一条光线
        If ray r hit the light //如果光线击中光源,说明该点对光源是可见的
            Lo += (1 / N) * L_i * f_r * cosine / pdf(wi) //计算着色，累加L0
    Return Lo

全局光照的算法

如图，如果使用直接光照算法，$p$ 点反射到 $Q$ 上的光线因为没有 hit 光源，所以是不考虑的。引入间接光照，我们要考虑弹射到光源的光线。

思路：将 $Q$ 点视为光源！这样就相当于我们在 $p$ 点计算 $Q$ 点的直接光照

1
2
3
4
5
6
7
8
9
10

shade(p, wo)
    Randomly choose N directions wi~pdf //随机选择N个wi方向
    Lo = 0.0 //初始化L0
    For each wi //对每个wi方向
        Trace a ray r(p, wi) //从p点往wi发出一条光线
        If ray r hit the light //如果光线击中光源,说明该点对光源是可见的
            Lo += (1 / N) * L_i * f_r * cosine / pdf(wi) //计算着色，累加L0
        Else If ray r hit an object at q //如果光纤击中其他物体q
            Lo += (1 / N) * shade(q, -wi) * f_r * cosine / pdf(wi) //⭐
    Return Lo

⭐解释第 9 行：
$L_i(p,\omega_i)$ 是入射光辐射率，这里将 $L_i$ 替换为 shade(q, -wi)，shade(q, -wi) 递归执行该算法结果返回 $Q$ 点的 $L_o$ ，这就很合理了，将 $Q$ 点反射的光作为 $P$ 点的入射光，就好像 $Q$ 点是光源一样！

传入 shade 函数的第二个参数为什么是 $-w_i$？
渲染方程定义 $w_i$ 从着色点指向光源，这里的 $w_i$ 是是相对于 $P$ 点来说的，即从 $P$ 点指向光源的方向为 $w_i$。我们知道 $shade(p, wo)$ 函数第二个参数传入的是出射方向，所以对于被当作光源的 $Q$ 点来说，该点的出射方向应该是 $Q$ 指向指向 $P$，方向和 $w_i$ 是相反的，所以要加负号。

算法优化->路径追踪

实际上，上述算法中一根光线打到物体后会反射很多个光线到同一个物体，以此类推产生指数爆炸。

如果蒙特卡洛积分采样次数为 1 则不会出现指数爆炸的现象。(1 的指数幂还是 1)，由此引出路径追踪算法：
只采样一次的光线追踪被称之为路径追踪（即一条光线形成了一条路径），上面采样 N 次的被称之为分布式光线追踪（Distributed Ray Tracing）

1
2
3
4
5
6
7

shade(p, wo)
    Randomly choose ONE direction wi~pdf(w) //随机选择1个wi方向
    Trace a ray r(p, wi)
    If ray r hit the light
        Return L_i * f_r * cosine / pdf(wi)
    Else If ray r hit an object at q
        Return shade(q, -wi) * f_r * cosine / pdf(wi

只采样一次会产生噪声，我们可以对每一个像素生成生成多个路径进行路径追踪，然后平均他们的结果：

对于每个像素发射 N 条光线（采样）做以下的算法，将每个射出去的采样接收到能量的点做蒙特卡洛积分得到平均值

1
2
3
4
5
6
7
8

ray_generation(camPos, pixel) 
    Uniformly choose N sample positions within the pixel //在该像素内选择N个采样点
    pixel_radiance = 0.0 //初始化反射率
    For each sample in the pixel //对于该像素内的每个采样点
        Shoot a ray r(camPos, cam_to_sample) //从相机位置发射光线
        If ray r hit the scene at p //如果击中场景中的物体
            pixel_radiance += 1 / N * shade(p, sample_to_cam) //执行上面提到的路径追踪算法，做蒙特卡洛积分得到平均值
    Return pixel_radiance

解决死循环->俄罗斯轮盘赌算法

上面的算法会无限递归，解决办法：俄罗斯轮盘赌（RR）算法

• 手动指定一个结束的概率 $P$，范围 0~1
• 每次调用 shade 函数的时候做一个随机数判断
• 若大于随机数则 return 0，反之继续执行 shade 函数并返回 $\frac{L_0}{P}$
• 这样算下来数学期望不会变，$E (L_o)=P∗(L_o/P)+(1−P)∗0=L_o$
• 且路径追踪总会停下来

该算法会产生噪声，但结果误差很小

1
2
3
4
5
6
7
8
9
10
11

shade(p, wo)
    Manually specify a probability P_RR  //指定结束概率
    Randomly select ksi in a uniform dist. in [0, 1] //取随机数
    If (ksi > P_RR) return 0.0;
    
    Randomly choose ONE direction wi~pdf(w) //随机选择1个wi方向
    Trace a ray r(p, wi)
    If ray r hit the light
        Return L_i * f_r * cosine / pdf(wi) / P_RR //结果除以P_RR
    Else If ray r hit an object at q
        Return shade(q, -wi) * f_r * cosine / pdf(wi) / P_RR //结果除以P_RR

提高效率->直接采样光源

前面介绍的算法是从着色点向每个 $w_i$ 方向发出光线去寻找光源。如果我们在着色点对半球均匀发射光线，会造成浪费。 比如下图中最左边光源比较大，只需要 5 个光线就能找到光源。最右边的例子光源比较小，需要 50000 个光线才能找到光源，大多数的光线被浪费了

在着色点采样寻找光源效率太低了，我们可以直接采样光源！直接计算光源对着色点的贡献
假设在面积为 A 的面光源上均匀采样：$pdf = 1/A$
原来的渲染方程是对着色点半球的积分，蒙特卡洛积分要求在什么地方采样就在什么地方积分。我们对光源采样，就要在光源上积分，我们就要对渲染方程进行改写。

1. 我们需要将 $d\omega$ 换成换成 $dA$，根据他们的位置关系和立体角定义，可以得
   $$
   d\omega=\frac{dA\cos\theta^{\prime}}{|x^{\prime}-x|^2}
   $$
2. 重写渲染方程，将渲染方程写成对光源的积分 $$
   \begin{aligned}
   L_o(x,\omega_o)& =\int_{\Omega^+}L_i(x,\omega_i)f_r(x,\omega_i,\omega_o)\cos\theta\mathrm{d}\omega_i \
   &=\int_AL_i(x,\omega_i)f_r(x,\omega_i,\omega_o)\frac{\cos\theta\cos\theta^{\prime}}{|x^{\prime}-x|^2}\mathrm{d}A
   \end{aligned}
   $$

这样就可以进行蒙特卡洛积分了

1. 根据蒙特卡洛积分：$\displaystyle\int_{a}^{b}f(x)\mathrm{d}x\approx\frac{1}{N}\sum_{k=1}^{N}\frac{f(X_{k})}{p(X_{k})}\quad X_{k}\sim p(x)$
2. 找到 $f(x)$：$L_i(x,\omega_i)f_r(x,\omega_i,\omega_o)\frac{\cos\theta\cos\theta^{\prime}}{|x^{\prime}-x|^2}$
3. 找到 $p(x)$ (概率密度函数)：$p(\omega_{i})=\frac{1}{A}$ 。
4. 带入原积分式：（注意该式中 $p$ 是着色点，$p(\omega_{i})$ 是概率密度函数）$$
   \begin{aligned}
   L_{o}(p,\omega_{o})& =\int_AL_i(x,\omega_i)f_r(x,\omega_i,\omega_o)\frac{\cos\theta\cos\theta^{\prime}}{|x^{\prime}-x|^2}\mathrm{d}A \
   &\approx\frac1N\sum_{i=1}^N\frac{L_i(x,\omega_i)f_r(x,\omega_i,\omega_o)\frac{\cos\theta\cos\theta^{\prime}}{|x^{\prime}-x|^2}}{p(\omega_i)}
   \end{aligned}
   $$

路径追踪算法总结

最终路径追踪分为两部分

• 直接光 (光源的直接光照)：直接对光源采样进行计算，不需要俄罗斯轮盘赌
• 间接光（光源照射其他物体，从其他物体反射来的光照）：还是用原来对着色点采样的方法，使用俄罗斯轮盘赌

如果光源中有物体挡住则不能渲染光源的直接光照，通过一个 if 来解决

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

shade(p, wo)
    # 直接光
    L_dir = 0.0
    Uniformly sample the light at x’ (pdf_light = 1 / A) //光源处均匀采样
    Shoot a ray from p to x’ //从着色点向光源发出光线
    If the ray is not blocked in the middle //如果路径中没有物体遮挡
        L_dir = L_i * f_r * cos θ * cos θ’ / |x’ - p|^2 / pdf_light //计算着色
    # 间接光
    L_indir = 0.0
    Test Russian Roulette with probability P_RR
    Uniformly sample the hemisphere toward wi (pdf_hemi = 1 / 2pi)
    Trace a ray r(p, wi)
    If ray r hit a non-emitting object at q //如果光纤击中场景中的物体
        L_indir = shade(q, -wi) * f_r * cos θ / pdf_hemi / P_RR
    Return L_dir + L_indir