本文讨论双向路径追踪算法的理论与实现。
LTE的路径和形式
对场景中表面上的某点$x$和某方向$\Phi$，用$\mathrm{Cast}(x, \Phi)$表示从$x$点沿方向$\vec\Phi$的射线与场景的首个交点。此时，光线传播方程（light transport equation，LTE）可以写作：
\[L(x \to \Theta) = L_e(x \to \Theta) + \int_{\mathcal S^2}f_s(\Phi \to x \to \Theta)L(\mathrm{Cast}(x, \Phi) \to -\Phi)\cos\langle N_x, \Phi\rangle d\omega_\Phi\]现在我们稍微扩张一下符号的含义：
\[\begin{aligned} L(x \to y) &= L(x \to \vec e_{xy}) \\ f_s(x \to y \to z) &= f_s(\vec e_{yx} \to y \to \vec e_{yz}) \end{aligned}\]再结合立体角微元和投影面积微元间的关系，便可以将LTE的积分域改写到场景表面$\mathcal M$上：
\[L(x' \to x) = L_e(x' \to x) + \int_{\mathcal M}f_s(x'' \to x' \to x)L(x'' \to x)G(x'' \leftrightarrow x')dA_{x''}\]其中：
\[G(x'' \leftrightarrow x') = V(x'' \leftrightarrow x')\frac{\cos\langle N_{x'}, \vec e_{x'x''}\rangle\cos\langle N_{x''}, \vec e_{x''x'}\rangle}{|x'' - x'|^2}\]$V(x’’ \leftrightarrow x’)$在$x’‘$和$x’$间没有遮挡物时为1，否则为0。现在我们不断把新LTE的左侧代换到它的右侧：
\[\begin{aligned} L(x_1 \to x_0) &= L_e(x_1 \to x_0) \\ & + \int_{\mathcal M}L_e(x_2 \to x_1)f_s(x_2 \to x_1 \to x_0)G(x_2 \leftrightarrow x_1)dA_{x_2} \\ & + \int_{\mathcal M}\int_{\mathcal M}L_e(x_3 \to x_2)f_s(x_3 \to x_2 \to x_1)G(x_3 \leftrightarrow x_2) \\ &~~~~~~~~~~~~~~~~\times f_s(x_2 \to x_1 \to x_0)G(x_2\leftrightarrow x_1)dA_{x_3}dA_{x_2} \\ & + \ldots \end{aligned}\]就得到了这样一个和式：
\[L(x_1 \to x_0) = \sum_{n = 1}^\infty P(\overline x_n)\]其中$\overline x_n$指长度为$n$的路径$x_n \to x_{n-1} \to \cdots \to x_1 \to x_0$，$P(\overline x_n)$是光源经$n$段传播过程后为最终结果作出的贡献：
\[\begin{aligned} P(\overline x_n) = \underbrace{\int_{\mathcal M}\cdots\int_{\mathcal M}}_{n-1}L_e(x_n \to x_{n-1})\left(\prod_{i=1}^{n-1}f_s(x_{i+1} \to x_i \to x_{i-1})G(x_{i+1}\leftrightarrow x_i)\right)dA_{x_n}\cdots dA_{x_2} \end{aligned}\]这就是LTE的路径和形式，它的含义是：$L(x_1 \to x_0)$可以拆成以下部分的和：
- 场景中的自发光经$0$次散射后从$x_1$点射向$x_0$的辐射亮度（其实就是$x_1$点的自发光）
- 场景中的自发光经$1$次散射后从$x_1$点射向$x_0$的辐射亮度
- 场景中的自发光经$2$次散射后从$x_1$点射向$x_0$的辐射亮度
- ……
- 场景中的自发光经$n$次散射后从$x_1$点射向$x_0$的辐射亮度
Measurement Equation
在之前的讨论中，我们总是说：只要对任意$x$和$\Phi$能求解$L(x \leftarrow \Phi)$，就算是大功告成了。但是仔细一想：把图像上的每个像素放到场景中都是一个矩形，该像素的颜色和函数$L(x \leftarrow \Phi)$之间的联系似乎也不是那么清晰？
事实上，设像素$j$在场景中对应的矩形区域是$\mathcal M_j$，$j$的颜色是由穿过$\mathcal M_j$的radiance flux决定的，也就是得在$\mathcal M_j$上积分。而且不是随便从哪个方向穿过$\mathcal M_j$的radiance都能被计入该flux中，只有特定方向的才能抵达摄像机镜头。以最基本的透视摄像机为例，设视点位置为$x_e$，对像素$j$上的每个点$x$，只有沿着方向$x \to x_e$穿过$x$的radiance才能计入其中。我们简单粗暴地取这些radiance的均值作为$j$的“颜色”：
\[I_j = \left. \int_{\mathcal M_j} L(x \leftarrow \vec e_{x_ex})d\omega_\Phi dA_x \middle / \int_{\mathcal M_j}dA_x \right.\]在这个例子中，每个点$x$仅对应唯一的方向$x_e \to x$。在更为一般的摄像机模型中（如需要考虑镜头半径、焦距等），每个$x$可能对应了一个方向的集合，记作立体角函数$\Omega_x$。$\Omega_x$中每个方向对$x$的“颜色”贡献可能是不一样的，我们用一个分布函数来表达，称为像素$j$的局部重要性函数$W_j$，此时$I_j$可以被写作：
\[I_j = \left. \int_{\mathcal M_j}\int_{\Omega_x}W_j(x \to \vec\omega)L(x \leftarrow \vec\omega)\cos\langle N_x, \vec\omega\rangle d\omega dA_x \middle / \int_{\mathcal M_j}\int_{\Omega_x}W_j(x \to \vec\omega)\cos\langle N_x, \vec\omega\rangle d\omega dA_x \right.\]具体到基本透视摄像机的例子里，根据：
\[\begin{aligned} &\left. \int_{\mathcal M_j} L(x \leftarrow \vec e_{x_ex})d\omega_\Phi dA_x \middle / \int_{\mathcal M_j}dA_x \right. = I_j = \\ &\left. \int_{\mathcal M_j}\int_{\Omega_x}W_j(x \to \vec\omega)L(x \leftarrow \vec\omega)\cos\langle N_x, \vec\omega\rangle d\omega dA_x \middle / \int_{\mathcal M_j}\int_{\Omega_x}W_j(x \to \vec\omega)\cos\langle N_x, \vec\omega\rangle d\omega dA_x \right. \end{aligned}\]可以解得：
\[W_j(x \to \vec\omega) = \frac{\delta(\vec\omega - \vec e_{x_ex})}{\cos\langle N_x, \vec\omega\rangle}\]我们不妨将像素$j$的局部重要性函数延拓到整个$\mathcal M$和整个$\mathcal S^2$，把$M_j$以外的$x$对应的函数值规定为0即可。此外，$I_j$中的分母对每个像素而言都是一个常量，不妨用它归一化一下好了，于是我们得到了这玩意儿：
\[W_e^{(j)}(x \to \vec\omega) = \begin{cases}\begin{aligned} &\left.W_j(x \to \vec\omega) \middle/ \left(\int_{\mathcal M_j}\int_{\Omega_x}W_j(x \to \vec\omega)\cos\langle N_x, \vec\omega\rangle d\omega dA_x\right)\right., &x \in \mathcal M_j, \vec\omega \in \Omega_x \\ &0, &\text{otherwise} \end{aligned}\end{cases}\]$W_e^{(j)}$就是所谓的重要性函数了，用$W_e^{(j)}$表示的$I_j$是这样的：
\[I_j = \int_{\mathcal M}\int_{\mathcal S^2}W_e^{(j)}(x \to \vec\omega)L(x \leftarrow \vec\omega)\cos\langle N_x, \vec\omega\rangle d\omega dA_x\]用立体角-面积转换把$\mathcal S^2$换成$\mathcal M$，记原来的$x$为$x_0$，用$x_1$表示换出来的$\mathcal M$上的点，就得到了所谓的Measurement Equation：
\[I_j = \int_{\mathcal M}\int_{\mathcal M}W_e^{(j)}(x_0 \to x_1)L(x_1 \to x_0)G(x_0 \leftrightarrow x_1) dA_{x_1}dA_{x_0}\]现在把LTE的路径和形式代入上式：
\[\begin{aligned} I_j &= \int_{\mathcal M}\int_{\mathcal M}W_e^{(j)}(x_0 \to x_1)L(x_1 \to x_0)G(x_0 \leftrightarrow x_1) dA_{x_1}dA_{x_0} \\ &= \int_{\mathcal M}\int_{\mathcal M}W_e^{(j)}(x_0 \to x_1)\left(\sum_{n = 1}^\infty P(\overline x_n)\right)G(x_0 \leftrightarrow x_1) dA_{x_1}dA_{x_0} \\ &= \sum_{n=1}^\infty \int_{\mathcal M}\int_{\mathcal M}W_e^{(j)}(x_0 \to x_1)P(\overline x_n)G(x_0 \leftrightarrow x_1) dA_{x_1}dA_{x_0} \\ &= \sum_{n=1}^\infty \underbrace{\int_{\mathcal M}\cdots\int_{\mathcal M}}_{n+1}W_e^{(j)}(x_0 \to x_1)T(\overline x_n)G(x_0 \leftrightarrow x_1)L_e(x_n \to x_{n-1}) dA_{x_n}\cdots dA_{x_0} \end{aligned}\]其中：
\[T(\overline x_n) = \prod_{i=1}^{n-1}f_s(x_{i+1} \to x_i \to x_{i-1})G(x_{i+1}\leftrightarrow x_i)\]路径采样
从LTE的路径和形式看来，我们可以在路径空间上直接采样路径，以对$I_j$进行估值。如果我们从摄像机出发一路采样出一条路径，就相当于是在使用路径追踪算法；如果我们从光源出发采样进入摄像机镜头的路径，就相当于是在使用“light tracing”。双向路径追踪则是两者的结合：从摄像机出发采样一条子路径（subpath），同时也从光源出发采样一条子路径，最后将这两条子路径连接起来。想法很简单，难点在于如何计算最后得到的路径的采样概率密度，以及一堆奇奇怪怪的corner cases。
首先描述一下从摄像机出发得到子路径的过程：我们让摄像机发射一条射线，一路按照与场景中其他物体的交点处的BSDF进行采样，且使用轮盘赌策略来结束采样过程。现用$x_0, x_1, \ldots$来表示该子路径上的点（$x_0$是镜头上的点），它们附带了一系列概率密度：
\[p_\mathrm{eye}(x_0), p(\vec e_{x_0x_1}\mid x_0), p(\vec e_{x_1x_2}\mid x_1), \cdots\]$p_\mathrm{eye}(x_0)$和$p(\vec e_{x_0x_1}\mid x_0)$是由摄像机给出的，后续的条件概率则是由每一次BSDF采样所使用的概率密度函数所决定的。注意条件概率中需要把轮盘赌的成功概率也算进去。类似地，对从光源出发的子路径$y_0, y_1, \ldots$，也有相似的计算过程。
最后，我们把$x_0x_1\cdots x_s$和$y_ty_{t-1}\cdots y_0$连接起来构成一条完整路径$\overline p_{s, t}$：
\[x_0x_1\cdots x_{s-1}x_sy_ty_{t-1}\cdots y_1y_0\]它的throughput值是：
\[\begin{aligned} T(\overline p_{s, t}) = &f_s(x_2 \to x_1 \to x_0) f_s(x_3 \to x_2 \to x_1) \cdots f_s(x_s \to x_{s-1} \to x_{s-2}) \\ \times &f_s(y_t \to x_s \to x_{s-1})G(x_s, y_t)f_s(y_{t-1} \to y_t \to x_s) \\ \times &f_s(y_{t-2} \to y_{t-1} \to y_t)\cdots f_s(y_1 \to y_2 \to y_3)f_s(y_0 \to y_1 \to y_2) \\ \times &W_e(x_0 \to \vec e_{x_0x_1})\cos\langle N_{x_0}, \vec e_{x_0x_1}\rangle L_e(y_0 \to \vec e_{y_0y_1}) \end{aligned}\]整条路径的采样概率则是上面那一大堆概率密度之积：
\[p(\overline p_{s, t}) = p_\mathrm{eye}(x_0)p(\vec e_{x_0x_1}\mid x_0)\cdots p(\vec e_{x_{s-1}x_s}\mid x_s)p(\vec e_{y_{t-1}y_t}\mid y_{t-1})\cdots p(\vec e_{y_0y_1} \mid y_0)p_\mathrm{light}(y_0)\]路径重用
在刚刚已经推出的方案中，我们每次只采样一条路径，将该路径的throughput除以其概率密度函数值，就可以拿来对$I_j$估值了。Veach很天才地想到：我们能不能把摄像机子路径的每个前缀和光源子路径的每个后缀都连接起来，从而得到一大堆子路径呢？这样做唯一的问题就是这一组路径间存在相关性，但只要我们重复这样的过程许多次，得到$N$组路径，那么某条长度为$n$的路径和这里面$N-1$个组中长度为$n$的路径都是无关的。随着$N$的增大，我们仍然可以保证估值无偏地收敛到正确的结果。并且当组数增大时，组内路径的相关性在全体路径中会显得愈发微不足道，即correlation也会收敛到0，不会对结果产生视觉可见的影响。
在使用了路径重用技术后，我们注意到，在每一轮采样（生成两条子路径并连接它们得到一堆路径，这样的过程称为一轮采样）中，一条长度为$k$的路径有$k + 1$个顶点，也就对应了$k+1$种采样方式。比如一条长度为2的路径有以下三种采样方式：
- 长度为0的摄像机路径和长度为2的光源路径连接起来
- 长度为1的摄像机路径和长度为1的光源路径连接起来
- 长度为2的摄像机路径和长度为0的光源路径连接起来
毫无疑问，这三种采样方式各有其擅长之处。譬如，1其实是light tracing的路径采样方式，擅长计算焦散等现象；3就是native path tracing；2是path tracing的直接照明项。我们能不能把这些采样方式的优点综合起来，而不是粗暴地求和呢？说到这里，一个名词已经呼之欲出了——多重重要性采样。
使用多重重要性采样的要点，在于给定按某种采样策略A得到的采样结果$a$时，如何计算$a$在B、C、D等其他采样策略下对应的概率密度函数值。在这个情景下这个问题其实很好办——沿着路径走一遍，把各种策略下的概率密度都算出来就行了。
几种特殊路径
s != 0, t = 0：采用MIS策略，而不是老老实实地连接什么子路径。连接子路径在这里相当于在以前的路径追踪中使用光源采样来计算直接照明。
s = 0, t != 0：这里连接子路径得到完整路径时，谁也不知道完整路径会对应到图像上的哪个像素去。因此要额外把这部分信息传递出去。这里会对系统的接口设计造成一点影响。
摄像机子路径击中了光源：完全不需要光源子路径的一种特殊情况。
光源子路径击中了镜头：要求镜头本身成为场景的一部分，不打算实现这个。
实现
（这篇文章说得有点不清不楚的，一是因为BDPT本来就很难讲，我建议直接去看Veach的PhD thesis，网上的博客之类的直接忽略；二是因为我这段时间花了些精力在毕业设计上，没太多心思去写博客。）
综上所述，我把BDPT算法分为以下两个步骤：
- 从摄像机和光源分别发射子路径，在子路径长度达到一定阈值后引入轮盘赌策略，记录下每个点相对于前一个点的条件概率。
- 对摄像机子路径的每个前缀和光源子路径的每个后缀，连接它们得到一条完整的路径，计算其throughput和MIS技术下的概率，加到最终结果上。
（施工中……）