本文是 ReSTIR DI 算法的一些前置知识和论文阅读笔记
更新策略：直接修改本文内容
关于采样技巧
蒙特卡洛积分和重要性采样
蒙特卡洛积分是通过采样来估算积分的方法：
\int_x f(x)dx = 1/N \Sigma_{i=1}^N f(x_i)/p(x_i), x_i~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 (1/ni) \Sigma_{j=1}^ni w_i(Xij) f(Xij) / P_i(Xij)该估计值是无偏的，当对于任意 X 满足 \Sigma_i w_i(X) = 1:
E[F]
= \Sigma_{i=1}^n \int_X w_i(X) f(X) / P_i(X) P_i(X) dX
= \Sigma_{i=1}^n \int_X w_i(X) f(X) dX
= \int_X f(X) dX一个常见的 w_i 取值（Balance Heuristic）是
w_i(X) = p_i(X) / \Sigma_k p_k(X)重要性重采样 (RIS)
重要性重采样 (Importance Resampling 或 Resampled Importance Sampling, RIS) 的方法允许你从一个可以是未被归一化的分布 p_hat(x) 中采样。它的方法为：
- 从一个方便采样的分布 p(x) 开始，采样得到 M 个样本 x1, x2,…, xM
- 给每个样本分配权重 w(xi) = p_hat(xi) / p(xi)
- 依据归一化后的权重，随机采样得到其中一个样本，记为 y。当 M 足够大时，y 的分布就足够接近 p_hat(x) （为此，p(x) 需要满足一定条件，见附录2）
上述步骤中, p(x) 称为候选分布（candidate distribution），而 p_hat(x) 称为目标分布（target distribution）。
进一步，它可以被用来估计积分：
F = f(y) / p_hat(y) * (1/M) \Sigma w(x_i)
E[F]
= E_{x1,...,xM} E_y F
= E_{x1,...,xM} \Sigma_{i=1}^M f(xi) / p_hat(xi) * p(Y=Xi) * (1/M) \Sigma w(x_i)
= E_{x1,...,xM} (1/M) \Sigma_{i=1}^M f(xi) / p_hat(xi) * w_(xi) / \Sigma w(x_i) * \Sigma w(x_i)
= E_{x1,...,xM} (1/M) \Sigma_{i=1}^M f(xi) / p_hat(xi) * w_(xi)
= E_{x1,...,xM} (1/M) \Sigma_{i=1}^M f(xi) / p(xi)
= \int_x f(x)加权池采样 (WRS)
在 x1,…,xM 中加权采样（权重函数为 w(xi)）N 个样本，是否可以只用 O(N) 的附加存储空间（即 x1,…,xM 以流的形式输入算法）？
以 N=1 的情形为例，可以的确 “流式” 处理：仅保留一个样本的池子，并维护一个已经处理过的元素的权重和，当要处理元素 xi 时，以
Pi = w(xi) / Sigma_{j=1}^i w(xj)的概率将元素池中的样本替换为 xi，最终处理完所有元素后，池子中的样本就等价于加权采样的结果。这个算法叫做加权池采样 (Weighted Reservior Sampling, WRS)
ReSTIR 渲染算法
论文：
“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_A rho(y, yx<->w) Le(x->y) G(x<->y) V(x<->y)方程关于所有发光（不包含反光）表面积分，发光面元的坐标为 x。我们考虑坐标 x 的发光面元 – 物体 y – 相机 这条光路，由此计算得到的是 y 点的直接光照（reflected radiance due to direct lighting）。其中：
- rho 是物体 y 点的 BSDF（双向散射分布函数，是双向反射分布函数 BRDF 和双向透射分布函数 BTDF 之和，可以处理带有透射的表面如玻璃和空气的交接面）。传入的方向是发光面元到 y 的入射方向，以及 y 到相机的方向 w
- Le 是面元的发光出射照度
- G 是 x 和 y 的平方反比距离乘上余弦项
- V 是 x 和 y 间的可见性
文章希望根据 f(x) = rho(x) Le(x) G(x) V(x) 来指导 RIS 的采样，这里使用的是 RGB 的加权和即“亮度” 作为指导采样的单值函数。然而，由于判断可见性需要追踪一条光线，开销太大，所以实际使用刨除 V(x) 的 p_hat(x) = rho(x) Le(x) G(x) 指导采样。
至此，思考两个问题：
Q: Le(x<->y) 和 G(x<->y) 的计算为什么不需要追踪一条光线？
A：首先，x 是从所有表面中采样得到的，而非先在 y 点表面的半球上选择一个方向再去追踪光线求交得到的（即与一般的渲染流程相反）。然后，Le 和 G 只与向量 xy （其方向和长度）有关，与路径中发生了什么无关，所以不用追踪光线。
Q: 在这个 RIS 采样中，候选分布是谁？目标分布是谁？
A: 候选分布是一个方便采样且能支撑目标分布的分布，比如说可以是按照发光光源面积的均匀分布（随口说的，如有错误请指出）；目标分布是 p_hat(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 = xidef render()
for pixel q in image
image[q] = shadePixel(RIS(q), q)
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)注：以上方法在 Ray tracing Gems 的章节中有所涉及
时空复用就是 RIS 池子合并
注意到临近像素或临近帧的 p_hat 函数很可能接近，因此我们可以将这些 RIS 池子合并到本池子中（权重和相加，按权重和抽取一个代表样本）
本算法额外地想要在将 q’ 的池子应用于 q 时，进行一个比例为 p_hat(q, r.y) / p_hat(q’, r.y) 的权重调整：令每个池子中存储 r.W = r.w_sum / p_hat(r.y) / r.M，这样调整后的权重和就不需要用到 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，这个情形下池子中存储的变量更多。后续 GI 被进一步改进，成为 再之后，ReSTIR 的思想被借鉴到一个复用 path 的算法，即 ReSTIR PT (Path-Tracing) 中，可以参考这篇博文。
附录
附录1：将多个分布的样本组合不等于从混合分布中采样多次
比如说，p1(x) = U[0, 1], p2(x) = U[2, 3]。下面两种方法：
- 从 p1(x) 采样 x1，从 p2(x) 采样 x2，并打乱顺序
- 从 ( p1(x)+p2(x) ) / 2 中采样 x1, x2
得到的 x1, x2 联合分布并不相同。前者必然得到一个 [0, 1] 间的样本，一个 [2, 3] 间的样本；后者有概率两个样本全部在 [0, 1] 区间。
附录2：候选分布应满足的条件
候选分布的 support 需要覆盖目标分布的 support。support 可以理解为“函数非零区域” 或者说 “样本可能出现的区域”。
Leave a Reply