采样技巧和 ReSTIR 渲染算法

2024-09-16

本文是 ReSTIR DI 算法的一些前置知识和论文阅读笔记

关于采样技巧

蒙特卡洛积分和重要性采样

蒙特卡洛积分是通过采样来估算积分的方法：

$$ \int_x f(x)dx = \frac{1}{N} \Sigma_{i=1}^N f(x_i)/p(x_i), x_i \sim p(x_i) $$

一个性质是，当样本的分布 p 接近 f 的形状时，估计值的方差就小；当 p 经过放缩等于 f 时，方差变为0，即仅用一个样本得到的估计值就是真实值。

以 $ f(x) = x, x \in [1,3]$ 为例，若选取 p(x) = x/4，则任意 Xi, f(Xi)/p(Xi) = Xi / (Xi / 4) = 4，为 f(x) 在 [1,3] 上积分的真实值。

重要性采样 (Importance Sampling) 就是努力让 p 接近 f 来减少估计值的方差的采样方法。

多重重要性采样 (MIS)

多重重要性采样 (Multiple Importance Sampling) 就是采用来自若干不同分布 p1(x), ..., pn(x) 的样本组合。注意，它在数学上不等价于使用 p1(x), ..., pn(x) 的混合分布采样多次得到样本（见附录1）。

通用地，假设每个分布采样 ni 个，则估计值

$$ F = \Sigma_{i=1}^n \frac{1}{n_i} \Sigma_{j=1}^{n_i} w_i(X_{ij}) f(X_{ij}) / P_i(X_{ij}) $$

当对于任意 X 满足 $\Sigma_i w_i(X) = 1$，该估计值是无偏的。证明：不妨设 $n_i=1$，免去下标 j，

$$ \begin{align*} E[F] &= \Sigma_{i=1}^n \int_X w_i(X) \frac{f(X)}{P_i(X)} P_i(X) dX \\ &= \Sigma_{i=1}^n \int_X w_i(X) f(X) dX \\ &= \int_X (\Sigma_{i=1}^n w_i(X)) f(X) dX \\ &= \int_X f(X) dX \end{align*} $$

一个常见的 $w_i$ 取值是 Balance Heuristic:

$$ w_i(X) = \frac{p_i(X)} {\Sigma_k p_k(X)} $$

重要性重采样 (RIS)

重要性重采样 (Importance Resampling 或 Resampled Importance Sampling, RIS) 的方法允许你从一个可以是未被归一化的分布 $\hat{p} (x)$ 中采样。它的方法为：

  1. 从一个方便采样的分布 p(x) 开始，采样得到 M 个样本 x1, x2,..., xM

  2. 给每个样本分配权重 $w(x_i) = \hat{p} (x_i) / p(x_i)$

  3. 依据归一化后的权重，随机采样得到其中一个样本，记为 y。当 M 足够大时，y 的分布就足够接近 $\hat{p}(x)$ （为此，p(x) 需要满足一定条件，见附录2）

上述步骤中, p(x) 称为候选分布（candidate distribution），而 $\hat{p} (x)$ 称为目标分布（target distribution）。

进一步，它可以被用来估计积分：设估计量 F 为

$$ F = \frac{f(y)} {\hat{p} (y)} * \frac{\Sigma w(x_i)}{M} $$

F 的期望为

$$ \begin{align*} E[F] &= E_{x_1,...,x_M} E_y F \\ &= E_{x_1,...,x_M} \Sigma_{i=1}^M p(Y=X_i) * \frac{f(x_i)} {\hat{p} (x_i)} * \frac{\Sigma_j w(x_j)}{M} \\ &= E_{x_1,...,x_M} \Sigma_{i=1}^M \frac{w(x_i)} {\Sigma_j w(x_j)} * \frac{f(x_i)} {\hat{p} (x_i)} * \frac{\Sigma_j w(x_j)}{M} \\ &= E_{x_1,...,x_M} \Sigma_{i=1}^M \frac{w(x_i) f(x_i)}{M \hat{p} (x_i)} \\ &= E_{x_1,...,x_M} \Sigma_{i=1}^M \frac{f(x_i)}{M p(x_i)} \\ &= \frac{1}{M} \Sigma_{i=1}^M E_{x_i} \frac{f(x_i)}{p(x_i)} \\ &= \int_x f(x) \end{align*} $$

加权池采样 (WRS)

在 x1,...,xM 中加权采样（权重函数为 w(xi)）N 个样本，是否可以只用 O(N) 的附加存储空间（即 x1,...,xM 以流的形式输入算法）？

以 N=1 的情形为例，可以的确 "流式" 处理：仅保留一个样本的池子，并维护一个已经处理过的元素的权重和，当要处理元素 xi 时，以

$$ P_i = w(x_i) / \Sigma_{j=1}^i w(x_j) $$

的概率将元素池中的样本替换为 xi，最终处理完所有元素后，池子中的样本就等价于加权采样的结果。这个算法叫做加权池采样 (Weighted Reservior Sampling, WRS)

ReSTIR 渲染算法

参考[1]，[2]，[3]

论文：
"Spatiotemporal reservoir resampling for real-time ray tracing with dynamic direct lighting"
利用时空域加权池采样做动态直接光照下的实时光线追踪

受到 RIS 和 WRS 算法的启发，论文提出了流式算法 ReSTIR（REservior-based Spatial Temporal Importance Resampling，基于池子的时空重要性重采样），解决了 在百万光源的场景高效地计算出直接光照 (direct lighting) 的问题。

求解直接光照的渲染方程

ReSTIR 要求解的渲染方程如下：

$$ L(y,w) = \int_{x\in A} \rho(y, yx \leftrightarrow w) L_e(x\rightarrow y) G(x \leftrightarrow y) V(x \leftrightarrow y) $$

方程关于所有发光（不包含反光）表面积分，发光面元的坐标为 x。我们考虑坐标 x 的发光面元 - 物体 y - 相机 这条光路，由此计算得到的是 y 点的直接光照（reflected radiance due to direct lighting）。其中：

文章希望根据被积函数 $f(x) = \rho(x) L_e(x) G(x) V(x)$ 来指导 RIS 的采样，这里使用的是 RGB 的加权和即“亮度” 作为指导采样的单值函数。然而，由于判断可见性需要追踪一条光线，开销太大，所以实际使用刨除 V(x) 的 $\hat{p} (x) = \rho(x) L_e(x) G(x)$ 指导采样。

至此，思考两个问题：
Q: $L_e(x\leftrightarrow y)$ 和 $G(x\leftrightarrow y)$ 的计算为什么不需要追踪一条光线？
A：首先，x 是从所有表面中采样得到的，而非先在 y 点表面的半球上选择一个方向再去追踪光线求交得到的（即与一般的渲染流程相反）。然后，Le 和 G 只与向量 xy （其方向和长度）有关，与路径中发生了什么无关，所以不用追踪光线。
Q: 在这个 RIS 采样中，候选分布是谁？目标分布是谁？
A: 候选分布是一个方便采样且能支撑目标分布的分布，比如说可以是按照发光光源面积的均匀分布（随口说的，如有错误请指出）；目标分布是 $\hat{p} (x)$，它可能亦被称作“光源分布”。

通过 WRS 做流式 RIS

使用 RIS 可以将多次采样获得的信息 “浓缩” 到一个样本点及其配套的池子权重之中，而这个过程因为 WRS 而可以流式（streaming）地进行（代码复制自 paper 的 Algorithm 3，有修改）：

class Reservior
  y, w_sum, M = 0, 0, 0
  def update(xi, wi)
    w_sum, M = w_sum + wi, M + 1
    if rand() < wi / w_sum
      y = xi

def RIS(q)
  Reservior r
  for i in 1..M
    generate xi ~ p
    r.update(xi, p_hat(xi)/p(xi))
  return r
def shadePixel(Reservior r, q)
  return f(r.y)/p_hat(r.y) (1/M r.w_sum) # Only trace a shadow ray here, since f(r.y)/p_hat(r.y) is V(r.y)
def render()
  for pixel q in image
    image[q] = shadePixel(RIS(q), q)

注：以上方法在 Ray tracing Gems 的章节中有所涉及

合并 RIS 池子来达到时空复用

注意到临近像素或临近帧的 $\hat{p}$ 函数很可能接近，因此我们可以将这些 RIS 池子合并到本像素的池子中（权重和相加，按权重和抽取一个代表样本）

本算法额外地想要在将 q' 的池子应用于 q 时，进行一个比例为 $ \hat{p}_{q} (r.y) / \hat{p}_{q'} (r.y) $ 的权重调整：令每个池子中额外存储一个新变量 W： r.W = r.w_sum / p_hat(r.y) / r.M，这样调整后的权重和就不需要用到 q' 了，传给 update 函数的第二个参数（权重调整后的 w_sum）可以写成和 q' 无关的形式：

r.w_sum * p_hat(q, r.y) / p_hat(q', r.y)
= ( r.W * p_hat(q', r.y) * r.M ) * p_hat(q, r.y) / p_hat(q', r.y)
= r.W * p_hat(q, r.y) * r.M

整理后的算法为 paper 中的 Algorithm 4：

def combineReservoirs(q, r1..rk)
  Reservior s
  for r in r1..rk
    s.update(r.y, p_hat(q, r.y) * r.W * r.M)
  s.M = r1.M + r2.M + ... + rk.M
  s.W = s.w_sum / p_hat(s.y) / s.M
  return s

每轮迭代中，让每个像素的池子更新为合并了 k 个周围的像素的池子。如果进行 n 轮迭代，花费的计算代价仅仅为 O(nk+M), 但每个像素因此在最好情况下获得了 $ O(k^n M)$ 个样本的信息，近似地等价于样本就是 p_hat 采样得到的，且伴随的 pdf 也正确。

时间复用的方法为，将当前池子与所有之前的池子中的一个（通过 Motion Vector 判断？）合并。

ReSTIR 算法的局限性和扩展

据实现过 ReSTIR 的同学告知，ReSTIR 实际使用中容易出现的问题是像素之间有相干性：一个显著更亮的光源样本可能会被相邻若干像素的样本池全部采用，因此产生亮暗条带的伪像（Artifact）。

一个可能的解决方法是，在确定每个像素的样本池后，将池中元素应用若干步 Metropolis-Hasting 算法，来减弱相邻像素样本的相干性（这种情况下 shadePixel 函数该如何改写？）

原论文提出的计算直接光照（Direct Illumination）的 ReSTIR 算法被称为 ReSTIR DI，后续它被扩展为支持全局光照（Global Illumination），称为 ReSTIR GI，这个情形下池子中存储的变量更多。再之后，ReSTIR 的思想被借鉴到一个复用 path 的算法，即 ReSTIR PT (Path-Tracing) 中，可以参考这篇博文。

附录

附录1：将多个分布的样本组合不等于从混合分布中采样多次

比如说，p1(x) = U[0, 1], p2(x) = U[2, 3]。下面两种方法：

  1. 从 p1(x) 采样 x1，从 p2(x) 采样 x2，并打乱顺序

  2. 从 ( p1(x)+p2(x) ) / 2 中采样 x1, x2

得到的 x1, x2 联合分布并不相同。前者必然得到一个 [0, 1] 间的样本，一个 [2, 3] 间的样本；后者有概率两个样本全部在 [0, 1] 区间。

附录2：候选分布应满足的条件

候选分布的 support 需要覆盖目标分布的 support。support 可以理解为“函数非零区域” 或者说 “样本可能出现的区域”。