本文是 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) 称为候选人分布（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 = \int_A rho(x) Le(x) G(x) V(x)方程关于所有发光（不包含反光）表面积分，考虑其对当前交点 y （在式子中略去）的贡献，由此得到 y 点的直接光照（reflected radiance due to direct lighting）。其中：
- rho 是 BSDF（双向散射分布函数，是双向反射分布函数 BRDF 和双向透射分布函数 BTDF 之和，可以处理带有投射的表面如玻璃和空气的交接面）
- Le 是（下个交点 x 的）发光出射照度
- 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) 指导采样。
通过 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(n**k M) 个样本的信息，近似地等价于样本就是 p_hat 采样得到的，且伴随的 pdf 也正确。
时间复用的方法为，将当前池子与所有之前的池子中的一个（通过 Motion Vector 判断？）合并。
ReSTIR 算法的局限性和扩展
据实现过 ReSTIR 的同学告知，ReSTIR 实际使用中容易出现的问题是像素之间有相干性：一个显著更亮的光源样本可能会被相邻若干像素的样本池全部采用，因此产生亮暗条带的伪像（Artifact）。
一个可能的解决方法是，在确定每个像素的样本池后，将池中元素应用若干步 Metropolis-Hasting 算法，来减弱相邻像素样本的相干性（这种情况下 shadePixel 函数该如何改写？）
原论文提出的计算直接光照（Direct Illumination）的 ReSTIR 算法被称为 ReSTIR DI，后续它被扩展为支持全局光照（Global Illumination），称为 ReSTIR GI，这个情形下池子中存储的变量更多。后续 GI 被进一步改进，成为 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] 区间。
Leave a Reply