AdaPT - Monte Carlo Path Tracer II

Posted on Edited on In Views: Symbols count in article: 9.2k Reading time ≈ 8 mins.

Ada Path Tracer II

​ 本文为 Ada Path Tracer I 的续文。在上一篇博客中，我只是说明了 Monte Carlo 积分在path tracing中的应用以及path tracing的理论、实现流程。并没有深入讨论 PDF 如何计算、采样如何进行、PDF的选择如何影响结果以及如何使得 path tracer 具有更高的效率。笔者在本文中对上文未讨论清楚的问题进行分析。分析过程中我果真也发现自己之前的某些理解是不对的，实现时出了问题，例如在：Lambertian BSDF的实现、MIS 对无面积光源的处理、MIS weight计算等等方面均有瑕疵。

​ 本部分学习的结束标志着【基本 path tracer 实现】的完成（V1.0）。除了对BSDF进行小的修改之外（表面反射定义过于简单），下一步我将主要围绕两部分进行深入研究：

  • Bidirectional path tracing。双向 tracer 对处理 directional / collimated 光源十分有用，并且在caustic建模上有着更好的效果（虽然不是最好的）
  • Volumetric path tracing。主要研究匀质介质散射吸收。
Figure 1. AdaPT 2min增量式渲染（CPU，GPU时间除以8-10）：点光源/顶部面光源与球面光源

II. PDF Confusion

2.1 讨论对象分析

​ 有如下问题需要明确：Path tracer中 PDF 对结果的影响是什么？选取不正确的PDF是否会导致结果成为有偏估计？PDF针对分子中哪一项起作用？

​ 已知MC积分公式： \[ \begin{equation}\label{mc_est} \frac 1N \sum_{i=1}^N \frac{f(X_i)}{p(X_i)} \end{equation} \] ​ 求期望可以得到： \[ \begin{equation}\label{mc} \frac 1N \sum_{i=1}^N E\left(\frac{f(X_i)}{p(X_i)}\right) \end{equation} \] ​ 而上式内部的期望为： \[ \begin{equation}\label{exp} E\left(\frac{f(X_i)}{p(X_i)}\right) = \int_\Omega \frac{f(x)}{p(x)}p(x)dx \end{equation} \] ​ 注意，根据随机变量函数的期望公式： \[ E(f(X)) = \int_\Omega f(x)p(x)dx \] ​ 可知公式\(\eqref{exp}\)求解的期望是随机变量\(x\)的函数\(f(x)/p(x)\)的期望。其中的\(p(x)\)应当仍然是\(x\)的分布。那么，在Monte Carlo积分中，或者在某个 vertex 处，\(x\) 代表了什么意义？我们举个例子来说明：在 direct component estimation （next event部分），我们采样M条shadow rays。则PDF是：\(p(x_l|e_i)p(e_i)\)，注意此处的PDF是定义在面积微元上的（differential area），而一般来说，我们需要的是“某个方向的PDF”（也即立体角定义）。首先我们需要记住立体角的物理概念：

The point from which the object is viewed is called the apex of the solid angle, and the object is said to subtend its solid angle at that point.

​ 接下来可以进行推导，见下图：

Figure 2. 微分立体角定义下的PDF与微分面积定义下的PDF（非光源采样）
Figure 3. 微分立体角定义下的PDF与微分面积定义下的PDF（光源采样 - 实际上与上一情况类似）

注意我们在direct component estimation 中，我们的实际估计目标是： \[ \int_\Omega f_r(x, w_i, w_o)L_d(x, w_i)\cos\theta dw_i\rightarrow \frac{1}{N}\sum_{i=1}^N\frac{f_r(x, w_i, w_0)L_d(x, w_i)\cos\theta}{p} \] ​ 也即我们实际应该获知各个方向上（微分立体角）的入射情况。而直接采样光源（特别是面光源）得到的PDF并不是定义在“方向上”的，而是在面积上的（面光源上某个微分面积与击中位置张成立体角），故应当进行转化。下面给出上式中分母\(p\)的形式。由于\(p\)实际上由面积定义的PDF转化而来，并且注意我们使用的并不是投影立体角 --- 我们需要的就是某一方向的PDF。\(p\) 自然应该是 \(p(w_i)\) 的形式，个人认为这与随机变量函数的期望计算有一定关系：随机变量是什么就用什么PDF（此处随机变量是立体角而不是投影立体角），并且在MIS中，我们需要与BSDF sampling使用同一域下的PDF（BSDF sampling使用立体角而非投影立体角）。非投影立体角下，\(p(w_e)\)应该是： \[ \begin{equation}\label{convert} p(x_l|e_i)p(e_i) \frac{\Vert x_l - x\Vert^2}{\cos(w_o, n_e)} \end{equation} \] ​ 故在笔者的实现内仅仅使用了光源的normal，而没有使用击中点的normal。

​ 故综合以上的例子来看，公式\(\eqref{mc}\)中的PDF到底是谁的PDF？就目前笔者的理解而言，笔者认为是【方向】。Path tracing 中的积分问题几乎都可以归结为在方向上的积分（而且有时必须归结为在方向上的问题，如MIS中统一在立体角方向域上讨论）。我们不妨进一步讨论BSDF常用的两个方法中，返回的PDF：

  • sample 方法：由BSDF采样获得反射、透射方向以及对应的PDF
  • pdf 方法：直接衡量某个方向的PDF，用以MIS

​ 既然以上两种方法都是在衡量某个方向的PDF，那么为什么诸如 Lambertian 这样的BSDF，以上两个函数返回的PDF却与入射方向有关（cosine weighted）？按理来说，Lambertian BSDF 向半球内的任意方向反射的可能都应相等。

至于为什么Lambertian PDF不是 \(1/2\pi\) 而是 \(\cos\theta /\pi\)，原因貌似是 importance sampling... 我们大可以使用\(1/2\pi\)对应的 uniform hemisphere sampling，但效率不如 cosine weighted hemisphere sampling 高。

​ 不过，这样的PDF是否会导致有偏估计？这就设计到 importance sampling 里的一些知识了。不过在此之前，请记住：

Lambertian uniform PDF 为\(1/2\pi\)是针对\(w\)（立体角）而言，如果变换到球坐标系下，应当为\(p(\theta, \phi) = \sin\theta /\pi\)。其他同理，故一定要 明确PDF在哪个空间下讨论。

It is important to be clear which PDF is being evaluated—for example, for a direction on the hemisphere, we have already seen these densities expressed differently in terms of solid angle and in terms of \((\theta, \phi)\). For hemispheres (and all other directional sampling), these functions return values with respect to solid angle. For the hemisphere, the solid angle PDF is a constant \(p(w) = 1/2\pi\). [1]

此外，在sampling阶段，PDF的选择是任意的。但一定需要保证：PDF 的计算结果与 采样方法 是一致的，这一点将在importance sampling中详谈。

2.2 Importance Sampling

As we will see in Section 13.10, it is often useful to sample from a distribution that has a shape similar to that of the integrand being estimated. For example, because the scattering equation weights the product of the BSDF and the incident radiance with a cosine term, it is useful to have a method that generates directions that are more likely to be close to the top of the hemisphere, where the cosine term has a large value, than the bottom, where the cosine term is small. [1]

​ Importance sampling 其实是老朋友了，还记得在特龙的《概率机器人》中，particle filter里就有“目标分布”和“建议部分”之说，并且作者称：建议分布越接近目标分布，采样效果越好，收敛越快。由于都是 Monte Carlo 方法（PF定位是 Monte Carlo 定位，Path Tracing 是 Monte Carlo 积分），故Path Tracing中如果可以获得类似公式\(\eqref{mc_est}\)中的分子项的PDF，则可以提升采样效率。

​ 不过需要记住的是：importance sampling 的思想内核是“贡献（contribution）大的部分多采样，贡献小的部分少采样”。这完全不同于“值（value）大的地方多采样，值小的地方少采样”，我们举与PBR book同样的例子：scattering equation： \[ \int_\Omega f_r(x, w_i, w_o)L_d(x, w_i)\cos\theta dw_i \] ​ 上式中真正产生贡献的是\(f_r\times \cos\theta\)，在此乘积项大的位置应当多采样，小的位置少采样（毕竟小的位置贡献小，极端一点的例子就是镜面反射）。采样时，我们既可以选择与\(f_r\)相似的 PDF，也可以选择与\(f_r\times \cos\theta\) 相似的PDF，这是因为，在公式\(\eqref{exp}\)中，分母人为除的PDF将会与当前sample出现的概率抵消。也即：

不考虑效率、认为样本数量无穷时，PDF可以是任意的，只要采样的方法与计算的PDF之间一致即可。（Dirac-delta分布除外）

​ 也即，我们如果使用\(1/2\pi\)作为PDF，就必须要使用 uniform hemisphere sampling，否则“抵消”将无法进行。在 Monte Carlo 积分估计式\(\eqref{mc_est}\) 中，这种抵消是隐式的（所以我花了很长时间才明白）！【样本\(X_i\)】按PDF \(p(x)\)定义的概率 \(P(X)\) 生成了出来，样本 \(x_i\) 或者与此类似的样本出现的频率（可能性）其实已经包含了PDF 信息\(p(x)\)，试想一个3面是4，余下三面为1、2、3的骰子，我们对其采样（投掷），得到的结果\(f(X_i)\) 确实看似不包含骰子的信息，但1、2、3、4在多次采样中出现的次数却包含PDF信息。

​ 故，有偏估计？有偏估计只会出现在采样方式与PDF不符的情况中，至少对于sample方法来说是这样的。sample方法获得的PDF并不一定与BSDF的形状相关，还可能与其他附加项有关（cosine term）。那么这里就有一个新的问题：pdf方法返回的PDF是否可以与其他项有关？因为我个人的想法是这样的：pdf项只会使用在MIS中，MIS需要联系两种不同的域，如BSDF sampling (solid angle)和light sampling (differential area)，而 light sampling PDF 在经过公式\(\eqref{convert}\)等号右边的分数项变换后将得到立体角定义下的 PDF，此 PDF 是否只应当与 BSDF 项相关？毕竟按照直觉理解，MIS就是为了针对【光源】以及【BSDF】都尝试计算某个方向入射是否可能，加上 cosine term 就不单单是在衡量BSDF而是在衡量 interactions 的整体贡献了。此问题的答案我们将在MIS中进行分析。

2.3 MIS

2.3.1 基本概念

​ Multiple importance sampling 是一种“信息融合”方法。它在 path tracing 中的一个最简单用法就是融合【BSDF sampling】以及【Emitter sampling】，融合的要求是：

Expressing all densities with respect to the same measure. [2]

​ 下面给出 path tracing 中的两种最常用的 ray direction 采样方法：

​ BSDF采样的意思就是：根据BSDF进行采样，采样更倾向于在BSDF值大的位置进行。这样做的好处非常显然：假设BSDF本身是一个窄带分布，在衡量直接光照时如果只是在光源上取点，这些点很有可能落在BSDF峰值区域外。这将会导致我们在低贡献区域浪费大量 samples。则可知，在BSDF分布较窄、光源较大时，BSDF采样结果较好。

​ Emitter采样是现阶段 AdaPT估计直接光照时的采样方法。由于对于漫反射表面而言，入射光线经过反射后将分布在整个半球上，当光源较小时（特别是点光源时），通过BSDF采样可能很难击中光源。那么我们也可以直接在光源上选点，使用所选点以及当前击中点构建shadow rays。这种方法在BSDF分布较宽、光源较小时的采样结果较好。

Figure 4. MIS结合两种采样方法有助于提升渲染质量

​ 复杂场景自然可能存在大/小光源以及光滑/粗糙/透射/散射的BSDF，则我们最好是找到一种方法，结合两种采样方法得到结果。这里有两种方案：

  • 简单方案：对于采样过程，我们只使用 emitter sampling。但在结果计算阶段，原始 Monte Carlo 积分将用\(p(x)\)除\(f(x)\)，而此方案中，我们对结果进行加权，不计算\(f(x)/p(x)\)，而是计算\(w(x)f(x)\)。\(w(x)\)称为 MIS 权重。
  • 复杂方案：采样过程同时实现 BSDF sampling 以及 emitter sampling，对于两种方法采样的结果，分别计算权重，最后的结果为：\(w_1(x)f_{\text{BSDF}}(x)+w_2(x)f_{\text{emitter}}(x)\)。仅从理论上看，这样的方法应该可以有更好的收敛性，因为上一方案仅仅在BSDF/emitter PDF结果冲突时进行重加权，并不帮助获得更好的样本。例如我们不管怎么在光滑表面上采样，落在specular lobe中的样本总是不够多的，虽然最后的权重低，不会有太大错误，但无助于加速收敛。本方案确实也相对复杂，一条shadow rays需要两次 ray intersection 操作。不过我们可以通过随机策略选择的方式，随机使用其中一种策略以降低计算量。

​ 如果每次从多种采样方法的每种方法中采样至少1个样本，则对应方法被成为“multi-sample estimator”[2]，形式如下。其中，\(n_i\geq 1\) 为第i种采样方法的样本数，\(w_i(·)\)为对应方法的权重函数。 \[ \begin{equation}\label{multi} F = \sum_{i=1}^n\frac {1}{n_i}\sum_{j=1}^{n_i}w_i(X_{i,j})\frac{f(X_{i,j})}{p_i(X_{i,j})} \end{equation} \] ​ 以上方法在满足一定条件下一定可以根据每一种无偏采样估计方法构建新的无偏估计方法，基本条件为：(1) \(\sum_i w_i(x) = 1, \text{ if } f(x)\neq 0\)， (2) \(w_i(x) = 0, \text{ if } p_i(x) = 0\)。具体的证明以及推导见 [2]

​ 如果每次随机选择某种采样方法，采样一个样本，则对应方法被称为“one-sample estimator”[2]。AdaPT与mitsuba(2/3)采用的则是这种方法，但区别是 mitsuba 使用的是 power heuristic 而 AdaPT 采用的是 balance heuristic（之后会说）。具体形式如下： \[ \begin{equation}\label{mono} F = \frac{w_I(X_I) f(X_I)}{c_I p_I(X_I)} \end{equation} \] ​ \(I\) 为随机变量，表示随机选择多种采样策略中的一种，\(c_I\) 为选择策略\(I = i\) 的概率。。AdaPT与mitsuba(2/3)采用的方法为 one-sample 的极端情况：只用一种方法进行sampling，但对sampling结果进行MIS 加权。相当于emitter、BSDF sampling同时存在，记为策略I与II，而选择策略I的概率为1，选择策略2的概率为0。

2.3.2 Heuristics

​ 权重函数的选择对于结果的方差有很大的影响。最极端的两个例子：

  • 某一策略的权重为1，其他策略为0。相当于固定某一种策略进行采样，这当然会引起一些问题，见图4
  • 策略权重一致，假设有N种策略，则每个策略的权重为常数 \(1/N\)。这种... 并不好，由于方差具有可加性[2]

​ 下面，笔者不加证明地给出一个较好的权重启发函数，其被称为“balance heuristic”: \[ \begin{equation} w_i(x)=\frac{n_ip_i(x)}{\sum_k n_kp_k(x)} \end{equation} \] ​ 也即：假设某个样本\(x\)原本来自策略\(i\)，我们不仅需要计算样本在此策略下被采样出的PDF，也需要计算在其他策略下被采样的PDF（所以pdf方法就有用了），最后的权重是本策略PDF在所有策略PDF之和中的比重。这样的函数在机器学习领域比较常见（可能会包装一下，比如softmax），所以很直观，其背后的intuition（话说，这个词翻译过来就是没那个味道）也很合理：如果某个策略拥有较大的PDF，而其他策略相对较小，则优先选用对应策略。

​ 注意，在one-sample 的极端情况中，有些策略存在但从来不进行采样，但这样的PDF也是需要计入的。例如AdaPT与mitsuba中，在计算直接光照时都使用 emitter sampling，但采样结果都会用BSDF pdf 计算BSDF PDF。最后的权重则为： \[ \begin{equation}\label{weights} w_{\text{AdaPT}} = \frac{P_{\text{emitter}}}{P_{\text{emitter}} + P_{\text{BSDF}}},w_{\text{mitsuba}} = \frac{P^2_{\text{emitter}}}{P^2_{\text{emitter}} + P^2_{\text{BSDF}}} \end{equation} \] ​ [2] 分析了 balance heuristic 与其他heuristic函数的方差关系，本文在此处引用其公式： \[ \begin{equation} V[F_{\text{balance}}] - V[F] \leq \left(\frac{1}{\min_i n_i} - \frac{1}{\sum_i n_i} \right)\mu^2 \end{equation} \] ​ 其中\(\mu = E[F]\)，也即需要估计的积分。一般来说，balance heuristic 已经是一个很好的heuristic函数了，其他heuristic函数相对 balance heuristic 而言具有方差下界（方差自然是越小越好，但就算是下界也差不了多少）。

​ 其他的 heuristic 函数就深入了，我只简单讲讲其中的power heuristic。power heuristic函数正如公式\(\eqref{weights}\)中 mitsuba 部分，它在原有的PDF上进行了幂操作。背后的intuition是：希望PDF变小的速度更快，以使得某一小PDF策略对权重产生的影响变得更小。具体的案例分析见[2]的【9.2.3.1 Low-variance problems: examples and analysis】。在这里，我不得不说，英语确实很重要，图形学有很多很棒的资料都没有中文翻译。那些教育局的领导拍拍脑袋就觉得高考取消英语有助于提高学生对于传统文化的认同与民族自信心，这种想法纯属多余，对理工类学生来说尤其如此。如果自信不是比出来的，而是闭门造车养出来的，那这种自信无疑非常脆弱。

​ 不过[2]中也分析了，在 one-sample 的情况下，balance heuristic 的表现将比 power heuristic 好（所以为什么mitsuba不用呢）。

2.3.3 上文问题的答案

sample方法获得的PDF并不一定与BSDF的形状相关，还可能与其他附加项有关（cosine term）。那么这里就有一个新的问题：pdf方法返回的PDF是否可以与其他项有关？

​ 答案是，主不在乎。对MIS而言，PDF只需要是对同一定义域下的某个随机变量而言即可。MIS权重只是一个系数，只要满足条件就不会影响无偏估计性。加cosine项一般意味着 BSDF PDF的峰值整体变窄了（\(f_r \times \cos\theta\) 只会比 \(\cos\theta\) 窄），从直觉上理解也即在此情况下，glossy reflection边缘（与入射与法向角度较大时）更倾向于使用 emitter sampling。这也是合理的，毕竟 glossy reflection 的占比与 diffusive reflection 占比在边缘处很接近。

Reference

[1] PBR-Book: 13.6 2D Sampling with Multidimensional Transformations

[2] Stanford graphics: Chapter 9 Multiple Importance Sampling