线性独立路径_蒙特卡洛路径追踪（Monte Carlo Pathing Tracing）(下篇)

【文章来源于斯坦福大学课件http://graphics.stanford.edu/courses/cs348b-01/course29.hanrahan.pdf，以下为翻译，大量使用google翻译，如有错误还请指出。上一篇https://zhuanlan.zhihu.com/p/110076660】

三：随机游走和Markov链

为了进一步了解路径跟踪为什么起作用，让我们考虑一个更简单的问题：离散随机游走。考虑由位置和方向等连续变量组成的，具有n个状态的离散物理系统，而不是具有连续变量的物理系统。如上所述的路径跟踪就是随机游走的示例，从一个采样点到另一个采样点，而这些采样点采样自连续的概率分布。在离散随机游走里，我们从一个状态转移到另一个状态，而采样时从离散的概率分布中得来。

一个随机游走初始化可以被分为三步：

1. 设从状态i开始的概率为
2. 设从状态i到状态j的概率为
3. 设在状态i终止的概率为

离散的随机游走包括如下几步：

第一步：创建一个随机粒子，它在状态i的概率是

第二步：在状态i终止的概率是

，如果恰好终止的话，则让状态i加一分，然后返回第一步。

第三步：没有终止，那么根据概率分布随机一个新状态j，并让目前的状态等于i，返回第2步。

随机游走也称为马尔可夫链。 马尔可夫链是由随机过程生成的一系列状态。 当我们讨论Metropolis算法时，我们将返回到马尔可夫链。 马尔可夫链还提出了贝叶斯推理和学习理论（例如隐马尔可夫模型）。 睁大眼睛看马尔可夫链； 您将看到在计算机图形学中越来越多地使用了这些技术。

给出一组符合随机游走规则的粒子，那么最终目标就是计算出各个粒子可能在哪个状态上终止。为了解决这个问题，我们引入一个新随机变量

，也就是在n次游走后在状态i的概率。由于这次游走的方向和上次游走的方向是完全独立的，所以我们只需要将概率简单叠加上去就行了。

定义一个矩阵，其第i行第j列的元素为

，那么上面的计算步骤可以看成矩阵乘以向量：

而上式的迭代形式又可以进一步简化为矩阵等式

因为

假设矩阵是概率分布，则此过程将始终收敛。 这样做的根本原因是，概率总是小于1，因此概率的乘积很快趋于零。因此，假设矩阵是概率转移矩阵，则随机游走提供了一种求解线性方程组的方法。 请注意，当我们使用Neumann级数求解渲染方程时，矩阵的离散迭代与连续算子的迭代应用相似。

这种使用离散随机游动来求解矩阵方程的方法可以直接应用于光能传递问题。 在光能传递公式中，

其中

这个形式因子

就是光束从某一个表面

到达另一个表面

的概率。因此，这个形式因子F也可以被解释为概率矩阵。在这个等式里，

就是漫反射率。反射率必须是正数并小于1。而吸收率因此等于

。

这些发现导致了一种使用离散随机游动来求解矩阵辐射度方程的方法。在使用辐射度解决方案技术时，经常遇到的主要问题是计算形式因子F矩阵。由于环境的复杂性和进行精确的可见表面确定的困难，该过程很耗费资源且容易出错。形式因子矩阵也很大。例如，由一百万个表面元素组成的场景将需要一百万个平方矩阵。因此，在实践中，形式因子矩阵通常是即时计算的。假设某个粒子位于某个表面元素i上，则可以从一个随机方向发出一个输出粒子，其中从余弦加权分布（此处的余弦相对于表面元素法线）中选择随机方向。然后对粒子进行射线追踪，并找到表面元素j上的相交闭合点。该随机过程大致等效于从已知形状因子矩阵生成的过程。

有趣的是，证明随机游走提供了线性方程组的近似解，不会偏差很多。 尽管此证明有点形式化，但值得探讨其背后的数学细节。

第一步是在所有路径的空间上定义一个随机变量。 让我们将长度为k的路径表示为

该路径涉及从状态

到状态

的一系列转换，最终在状态ik中的k个转换之后结束。 不带下标的随机变量α是长度为1到无穷大的所有路径的集合。

下一步是为每个路径选择一个估计量W（α）。 然后，我们可以通过给定路径被采样的概率p（α）加权W来计算估计量的期望值。

在最后一行中，我们将相同长度的所有路径归为一组。 每个索引i的总和从1到n-系统中的状态数。 因此，存在

个长度为k的路径，并且路径当然可以具有无限的长度。 需要考虑的路径太多了！

以状态

结束的路径

的概率

是多少？ 假设上述离散随机游走过程，是粒子在状态i1中产生的概率，乘以进行必要转换以到达状态

的概率，乘以在状态

中终止的概率

通过这些仔细的定义，可以计算出期望值

回想一下，我们的估算器计算的是计数以状态j终止的粒子的数量。 在数学上，我们可以用一个增量函数来描述这个计数过程，

，此增量函数仅对以

结尾的粒子评分。 期望值就是

该总和可以被视为矩阵方程的j分量

这是线性方程组的理想解。

请注意，我们必须谨慎对待估算器。 如果我们不将粒子计数除以终止事件的可能性，则期望值将不会等于正确答案。 当选择蒙特卡洛技术时，为复杂的采样过程选择错误的估计器（即导致错误的期望值的估计器）是最常见的错误之一。 在您没有丰富经验之前，值得您说服自己的估计量是公正的。

该技术最初是由蒙特卡洛方法的创始人冯·诺伊曼和乌拉姆开发的。 他们使用的估算器通常称为吸收估算器，因为仅对被吸收的粒子进行计数。 Wasow开发的一个有趣的变体是计算通过状态j的所有粒子的数量（包括终止和过渡的粒子）。 这被称为碰撞估计器，因为它计算与表面碰撞的所有粒子。 这是一个有趣的练习，它表明碰撞估计器还提供了线性方程解的无偏估计。 推导碰撞估计器何时以及是否比吸收估计器更有效的条件也更具挑战性，但更有趣。

这种基本的证明技术很容易推广到连续分布，但表示法很混乱。 有关详细信息，请参见《 Spanier and Gelbard》一书中子传输[？]的第3章，这是最权威的资料，如果您想了解蒙特卡洛技术背后的理论来解决积分方程和运输问题。

四：伴随方程式（Adjoint Equations ）和重要性抽样

回想一下，像素的具体颜色，等于长度为n的路径上的总和

在这里我们切换了表示法并将术语写为

如上所述，在光和传感器的互换下，该方程是对称的。 用R代替

，并注意

在第二步中，我们从几何形状的对称性中注意到

由于相互原则，BRDF也是对称的

这些对称性意味着我们可能从光线或眼睛发出光线轨迹；这两种方法都将有相同的积分。

假设现在我们在点k处断开路径。 产生给k的光量为

以类似的方式，将传感器视为虚拟光源，我们可以计算出来自传感器的光量达到k

然后

请注意，使用符号LS和LR来指示来自光源的辐射“发射”与接收器。

具有传输方程的前向和后向（伴随）解的路径。 从源项S生成前向解，从接收项R生成后向解。对于在路径反转下传输方程不变的物理情况，前向方程和后向方程相同。

我们对此方程进行观察。 首先，可以将该方程视为两个辐射函数的内积。 如果我们认为辐射度是射线

的函数，则如果我们具有函数f（r）和g（r），则f和g的内积为

此处

是光线的估计值。自然的测量两个表面元素A和A'直接的光线是

同样的，考虑r被位置x和方向w参数化，也就是

其次，这种积分自然会导致一种对路径进行重要性采样的方法。 假设我们正在追踪光线并到达表面k。 为了计算传感器响应，我们需要将L与R进行积分。从这个意义上讲，R被认为是采样下一方向的重要函数，因为我们希望与R成比例的采样技术能够实现低方差。 但是R是逆向传输方程的解，如果我们要追踪来自传感器的射线，则可以计算出该反向方程。 R告诉我们来自传感器的光将达到这一点。 因此，后向解决方案为前向解决方案提供了重要的功能，反之亦然。 这是双向射线追踪之间的关键思想。

使用算子可以很容易地处理伴随方程。 用算子表示，一个积分方程就是

我们要估计由测量方程式给出的积分，它是两个函数的内积

f作为传感器的响应，K◦g作为渲染方程的解。 该方程式可能会重新排列

注意下面两个等式的不同：

和

一个积分位于第一个变量上方，另一个积分位于第二个变量上方。 当然，如果

则这两个积分是相同的，这时候

，并且算子是自伴的。

该表示法提供了一种简洁的方法来证明前向估计等于渲染方程的后向估计。 回忆

我们也可以在另一个方向上写一个对称方程

那么

即使算子不是自伴的，该结果仍然成立。 我们将把证明作为练习。

这个结果很漂亮，在实际计算时也很有用。 伴随方程在数学物理的所有领域中都有许多应用。 它们允许您执行的操作是创建输出敏感算法。 通常，当您求解方程式时，您到处都会找到答案。 想想光能传递； 使用有限元方法时，可以解决所有曲面上的光能传递问题。 这同样适用于光线追踪或经典的离散随机游动； 您可以解决粒子在任何状态下着陆的可能性。 但是，在许多问题中，您只想在几点找到解决方案。 在图像合成的情况下，我们只需要计算我们看到的或落在胶片上的辐射率。 仅在可以观察到其影响的情况下，才需要计算其他位置的辐射。

我们可以将解决方案子集的选择建模为响应函数的内积乘以辐射度。 如果我们只想观察解决方案的一小部分，则可以在不关心的位置将响应函数设为零。 现在考虑当所有表面都充当源并且只有胶片平面贡献了非零响应的情况。 从源头开始运行粒子跟踪算法将非常低效，因为粒子很少会终止在胶片平面上。 但是，以相反方向运行该算法非常有效，因为所有粒子都将在源上终止。 因此，每个粒子都提供有用的信息。 逆向思维会让算法更加有效。

仅解决问题的子集是蒙特卡罗技术的一大优势。 实际上，在算法开发的早期，蒙特卡洛技术就用于求解方程的线性系统。 事实证明，如果您只想求解一个变量，这个方法将非常有效。 但请注意：如果您想寻求完整的解决方案，那么像Gaussian elimination这样的更常规的技术会更加有效。