从零开始学BDPM

前置知识

渲染方程对路径积分

我们将基本的渲染方程递归形式展开：

可以写成无限求和的形式，其中：

Given above equation and a particular length , all that we need to do to compute a Monte Carlo estimate of the radiance arriving at due to paths of length is to sample a set of vertices with an appropriate sampling density in the scene to generate a path and then to evaluate an estimate of using those vertices. Whether we generate those vertices by starting a path from the camera, starting from the light, starting from both ends, or starting from a point in the middle is a detail that only affects how the weights for the Monte Carlo estimates are computed.

Russian Roulette and Splitting

RR

Russian roulette addresses the problem of samples that are expensive to evaluate but make a small contribution to the final result. Russian roulette makes it possible to skip tracing rays when the integrand’s value is very low but not necessarily 0, while still computing the correct value on average.

可以验证，RR保持期望不变，即$E[F^{\prime}]=E[F]$；但总是增大方差，除非$c=F$（不是$E[F]$）。选择合适的$q$可以得到较小的方差。

Splitting

Splitting addresses this problem by formalizing the approach of taking multiple samples in some of the dimensions of integration for each sample taken in other dimensions.

例如，为了计算直接光照，如下式：

使用Monte Carlo估计。每个像素进行$N$次采样，每次像素采样在与表面的首个交点都只对灯光采样一次。在直接光照比较复杂的情况下，$N$要比较大，但较大的$N$对于减少像素值的方差收益较低。因此，对于每次像素采样，都进行$M$次对灯光的采样。相比原来$M=1$，在对灯光的采样次数$N\times M$相同时，发射光线的总数（或路径总长度）$N+N\times M$却减小了。我们对积分的估计为：

Multiple Importance Sampling

在使用Monte Carlo估计积分时，如果采样分布$p(x)$与被积函数$f(x)$越相似，那么相同采样次数时结果越可能准确，从而提高采样效率。特别地，如果$p(x)=\frac{f(x)}{\int f(x) \,\mathrm{d}x}$，那么方差就是$0$。实践中，尽管我们无法获知$f(x)$，但至少获知其中的某一项，并由它构建采样分布，也是极好的。

但假如我们获得了被积函数$F$中的几项（例如BRDF$f_r$和直接光照$L_d$），但又不好按这几项的乘积来采样，那又如何呢？如果只根据其中一项的分布采样，这一项却可能与它们乘积的分布相去甚远。答案就是多重重要性采样（MIS）。它根据每一项不同的采样分布，分别作一定次数的采样，并使用一个权函数$w(x)$重新计算每次采样（Monte Carlo估计）的权重。以两种采样分布的情况为例，使用MIS的Monte Carlo估计是：

权值函数通常会考虑到样本$x$在每一种采样分布$p_s$下的概率。常用的banlance heuristic计算权函数的策略如下，注意加入这个权函数后估计依然是无偏的：

考虑这种情况：假设$p_f(X)$和$p_g(X)$分别是$f(X)$和$g(X)$合理的采样分布。如果我们根据采样概率$p_f(X)$采样到了X，当$p_f(X)$很小而不巧$g(X)$却很大，$f(X)\cdot g(X)$就会很大，从而引入很大的方差。而引入权函数后，由于$p_g(X)$会比较大，因此这次并不合理的采样的权重$w_f(X)$会很小——大胜利！

Path Tracing

我们的目标是计算无限条不同长度路径的贡献之和（见Integral over Paths）：

我们可以使用RR对上式的期望作得一个无偏估计。RR并没有解决无穷求和的问题——对于任何一个有限的步数，我们都有可能达到它。

在给定路径长度时，为了采样长度$i$的路径，我们首先想到（pbrt让我们想到）的方式是在所有物体的表面均匀采样$n$个点。对于任何一个点，它采样的概率密度函数都是$p_i=\frac{1}{\sum_j A_j}$。那么使用Monte Carlo估计这条路径的贡献就是:

注意这个式子实际上采样了$\mathrm{p_2}$到$\mathrm{p_n}$这$n-1$个点，这里PBRT假定$\mathrm{p_0}$和$\mathrm{p_1}$是已经被传感器位置和camera ray确定了。

但这种采样方式存在许多问题。注意分子上的几何项连乘，路径中任意相邻的两点间互不可见就会导致路径的贡献为0，另外也无法根据BRDF进行重要性采样。这导致很多采样都是低效的，增大了估计值的方差。因此一般采用递增的路径采样方式，每次在路径最后的点采样方向$\omega_i$。下式给出了固定路径长度$n$时使用路径递增采样方法的贡献计算，注意对于路径的终点$\mathrm{p_i}$我们还是在光源的表面采样：

注意上式中分母采样概率的变化。由于我们从表面采样改为立体角采用采样，因此需要应用下面的转换:（这样说反而更复杂了，直接从展开立体角形式的Monte Carlo估计渲染方程的递归积分角度理解还更容易一些）

Path-Space Measurement Equation

The measurement equation describes the value of an abstract measurement that is found by integrating over some set of rays carrying radiance. For example, when computing the value of a pixel in the image, we want to integrate over rays starting in the neighborhood of the pixel, with contributions weighted by the image reconstruction filter.

$W_e$就是$\mathrm{p_0}$和$\mathrm{p_1}$间的重要性(importance)。例如，对于针孔相机(不考虑DoF)而言，像素值度量的重要性函数表示传感器对场景中一点$x$的方向$w$光线的相应，它通常包括了delta分布，即只考虑传感器上$\mathrm{p_0}$射出相机的那一个方向：

这样带入LTE的路径积分形式，发现了$L_e$与$W_e$具有对称性：

这使得我们可以想象以下过程，即镜头发射假想的光线，击中光源后由$L_e$作为度量方程计算贡献（path tracing），传播importance；也可以想象逆过程，即光源向外发出射线（或粒子），击中传感器后由$W_e$计算贡献，这就是particle tracing——

Particle Tracing

Particle tracing算法会产生$N$个样本$(\mathrm{p_j}, \mathrm{\omega_j}, \beta_j)$，表示其位于顶点$\mathrm{p_j}$，有入射方向$\mathrm{\omega_j}$和权重$\beta_j$。权重包括了路径的throughput和采样概率。

我们的目标是使用这些粒子样本计算某个权重函数$W_e(\mathrm{p_j},\mathrm{\omega_j})$下的某个度量方程(measurement equation)，因此我们希望相同权函数下这些粒子权重的期望与该度量方程的值相同：

那么如何生成这些粒子，如何计算它们的权重以满足我们的目标呢？考察我们的计算目标（上式右侧），用对它使用Monte Carlo积分估计可得：

更进一步，我们在路径空间使用Monte Carlo估计上式中的$E\left[\frac{L\left(\mathrm{p}_{\mathrm{i}, 1} \rightarrow \mathrm{p}_{\mathrm{i}, 0}\right)}{p(\mathrm{p_{i,1}})}\right]$，其中$q_{i,j}$是第$j$步终止路径的概率。那么把$E$展开就是：（展开过程回顾上面path tracing）

粒子权重$\beta_{i,j}$的计算已经就写出来了。

说实话，PBRT这一节属实有些抽象。做完作业之后回头看看我写的，发现没有什么逻辑。还是直接看下一节比较好……

Photon Mapping

传统的光子映射分为两步：

(1) 从光源发射若干携带光通量$\Phi_p$的粒子，称为光子。与PT中的光线类似，光子在场景的表面间不断弹射，每经过一次碰撞光子携带的能量都会与表面的BRDF等相互作用。我们将光子击中表面时的位置、能量、入射方向等信息保存下来，就得到了光子图。

(2)从相机向场景发射光线(camera ray)，为了得到表面着色点的出射辐射度，我们可以利用第一步得到的光子图中的信息进行近似，如下式。其中$\kappa$是密度估计函数，例如$\frac{1}{\pi r^2}$。注意这个近似(密度估计)方式是有偏的，因为我们使用了交点$x$附近一定区域内的光子信息，而不是恰好处于该点的。

在原版的光子映射算法中还使用了一些技巧来提升收敛速度、降低噪声等，例如final gather，它是在camera ray与场景的交点处再采样并发射若干光线，通过Monte Carlo方法计算间接的入射辐射度。

Bidirectional Photon Mapping

双向光子映射主要是(1)将光子映射算法与路径追踪相结合，并且(2)在估计光子贡献时设计了一个多重重要性采样的权重函数。

基于Monte Carlo的路径追踪通过采样多个路径样本，并使用下式估计像素的值：

使用在上一节说明的、利用光子图估计辐射度的方法替换得到下式。其中T是camera path的throughput，就是BRDF值$f_r$、BRDF采样概率和Russian Roulette采样概率（如果使用RR）在每个交点$\mathrm{p}$的连乘：

上式还包含了一个MIS权重函数$w\left(X_{s(p), t}\right)$，其计算方式如下，即考虑在当前长度为$n$的路径$x$上，所有不同的光源子路径长度$s$和相机子路径长度$t$下的采样概率：

路径概率的计算方式如下：

由此我们需要分别计算相机子路径的采样概率$p^E_t$和光源子路径的采样概率$p^L_s$，它们的计算方式基本相同。以相机子路径为例：

概率$p\left(x_{i-1}, x_{i-1} \rightarrow x_{i}\right)$就是在$x_{i-1}$发射相应光线的采样概率，包含了BRDF采样概率和Roussian Roulette概率。而通过将相机子路径路径上顶点的标号扩展到完整的路径上：

就可以得到在$t<i<=n$时，从相机出发采样$i$段路径的概率。对于从光源出发路径的采样概率同理：

Implementation

BDPM与PM类似，算法依然可以分为两步，因此我们分别探讨每一步实现需要在mistuba上作出的修改。

First Pass

BDPM的第一步仍然是构建光子图。mitsuba已经实现了PM算法，并且实现了其所依赖的光子图类PhotonMap。与PM不同的是，BDPM需要记录路径上每个顶点的历史信息，并保存在光子的数据结构Photon中：

具体地，为了得到$p_i^L, \ 1<i<=s$, 需要在光子中保存路径从光源开始、到每个顶点的采样概率$p_i^L$。同时，为了得到$p_{t+1}^E, \ i\in \langle 2, s\rangle$，需要得到每个顶点沿反方向采样的概率，也就是$p\left(x_{t+i-1}^{\prime}, x_{t+i-2}^{\prime} \rightarrow x_{t+i-1}^{\prime} \rightarrow x_{t+i}^{\prime}\right)$。为了在每个顶点得到这两个概率，在光子中添加矢量数组用于记录。

Second Pass

第二步的总体框架是基于Path Tracing，最大的不同是在每一次发射光线与场景表面的交点上利用光子图中的信息进行一次辐射度估计，基本上完全参考了论文中Pass 2的伪代码进行实现。具体实现代码位于BDPMIntegrator类中的Li()函数，它采用递归的形式，每一层都进行一次辐射度估计，计算本次光子的贡献，之后采样并发射下一条光线进入新的递归。函数大致流程如下：

(1) 将当前光线与场景求交，如果找到交点，检查交点是否位于光源上并计入其贡献；如果没有交点，检查是否存在环境贴图并计入其贡献。

(2) 在交点处使用光子图估计出射辐射度，计算交点附近光子的贡献，并与相机子路径上的throughput相乘。如果交点BRDF是一个Delta函数，那么直接进入下一步，因为光子正好从该方向入射的概率为0，他们的贡献也就是0。

(3) 使用Russian Roulette判断光线是否被吸收。如果没有，则根据BRDF重要性采样，在交点处生成新的光线并进入下一次递归。

与第一步类似，在每次递归中，都需要记录当前相机子路径的采样概率，以及反方向采样当前BRDF的概率密度函数值。

MIS Radiance Estimate

我们基于BDPM::estimateRadiance函数实现新的辐射度估计函数estimateRadianceBDPM，并加入论文中提出的多重重要性采样权重函数，权重的计算过程同样包含在该函数中。

函数利用mitsuba已经实现的N近邻搜索方法，找到交点一定半径内的光子，随后计算每个光子的贡献，即光子能量$\times$BRDF函数$\times$权重函数。为了计算权重函数，需要利用在光子数据和相机子路径中保存的整条路径上的概率信息。