論文閱讀2020:ReSTIR(Spatiotemporal reservoir resampling for real-time ray tracing with dynamic direct lighting)（一):主體思想

問題引入

背景

渲染是計算機圖形學的最重要方向之一，是模擬場景/模型成像的過程。渲染在建築，游戲，模擬，電影電視特效，以及醫學可視化中都有廣泛的應用。目前基於物理渲染Ray Tracing仍然難以達到實時，學術界和工程界持續進行渲染性能研究和探索。ReSTIR這篇英偉達2020的文章,在圖形學領域重新發現Rediscover了加權蓄水池采樣（Weighted Resorvior sampling）方法，帶來采樣性能的巨大提升，而過去加權蓄水池主要用於金融、統計等領域。本篇文章將介紹蒙特卡洛光追和蓄水池采樣用於蒙特卡洛計算的基本原理。

基於蒙特卡洛的光線追蹤

渲染的基本問題是求解渲染方程，求出物體在光源的照射下，出射到人眼方向的輻射亮度(Radiance，\(L_0\))。

渲染方程表達如下，要計算任意一個表面位置\(x\),朝向為\(\omega_0\)的輻射亮度\(L_0(x,\omega_o)\),等於位置\(x\)處的自發光\(L_e\)，與半球面所有方向上三個乘積項乘積的積分的和，這三個乘積項分別是入射輻射亮度\(L_i\)，雙向散射分布函數BRSF \(f\)，以及入射方向與表面法線的夾角的余弦值。該渲染方程展開是一個無限積分，很難求出解析解，一般采用蒙特卡洛方法進行估計。

\[L_{o}(x,\omega_o)=L_{e}\left(x, \omega_{o}\right)+\int_{\Omega} L_{i}\left(x, \omega_{o}\right) f\left(\omega_{i}, \omega_{o}\right)\langle\cos \theta\rangle d \omega_{i} \]

蒙特卡洛積分的正確性推導

• 首先我們按照概率密度函數\(p(x)\)在區間\([a,b]\)上進行采樣得到樣本集\(\{X_1,X_2,...,X_N\}\)
• 再構造新的\(F_N\)函數：

\[F_{N}=\frac{1}{N} \sum_{i=1}^{N} \frac{f\left(X_{i}\right)}{p\left(X_{i}\right)} \]

• 計算\(F_N\)的數學期望如下：

\[E\left[F_{N}\right]=E\left[\frac{1}{N} \sum_{i=1}^{N} \frac{f\left(X_{i}\right)}{p\left(X_{i}\right)}\right]=\frac{1}{N} \sum_{i=1}^{N} E\left[\frac{f\left(X_{i}\right)}{p\left(X_{i}\right)}\right]\\=\frac{1}{N} \sum_{i=1}^{N} \int_{a}^{b} \frac{f\left(X_{i}\right)}{p\left(X_{i}\right)} p\left(X_{i}\right) d x=\frac{1}{N} \sum_{i=1}^{N} \int_{a}^{b} f(x) d x\\=\int_{a}^{b} f(x) d x \]

• 經過上述推導以及大數定理,可以得到蒙特卡洛積分方法如下：

\[\int_{D} f(x) d x=\lim _{N \rightarrow \infty} \frac{1}{N} \sum_{i=1}^{N} \frac{f\left(X_{i}\right)}{p\left(X_{i}\right)} \]

蒙特卡洛積分法的收斂性和收斂速度分析

求上述\(F_N\)的方差如下：

\[\sigma^{2}\left(F_{N}\right)=\sigma^{2}\left[\frac{1}{N} \sum_{i=1}^{N} \frac{f\left(X_{i}\right)}{p\left(X_{i}\right)}\right]=\frac{1}{N^{2}} \sum_{i}^{N} \sigma^{2}\left[\frac{f\left(X_{i}\right)}{p\left(X_{i}\right)}\right] \]

令：

\[Y=\frac{f(X)}{p(X)} \]

則：

\[\sigma^{2}\left(F_{N}\right)=\frac{1}{N^{2}} \sum_{i}^{N} \sigma^{2}\left[Y_{i}\right]=\frac{1}{N^{2}}\left(N \sigma^{2}[Y]\right)=\frac{1}{N} \sigma^{2}[Y] \]

則：

\[\sigma\left[F_{N}\right]=\frac{1}{\sqrt{N}} \sigma[Y] \]

• 由此可知當采樣個數N 趨向無窮大時, \(\sigma\left[F_{N}\right]\)趨向於0，這也證明蒙特卡洛積分法的收斂性
• 由方差公式可以知道蒙特卡洛積分法收斂速度和穩定性由\(Y\)決定，如果\(Y\)變化劇烈，就會造成方差較大，估計值的收斂速度越慢；若\(f(x)\)和\(p(x)\) 越接近(即若他們的形狀越接近)，則收斂更快更穩定。最理想的情況是\(p(x) = f(x)\)或成正比，則方差為0，快速收斂。

因此，選擇任意概率密度分布函數\(p(x)\)，當收斂次數足夠大時，蒙特卡洛積分都會收斂；若\(f(x)\)和\(p(x)\)形狀接近，成正比，則收斂更快更穩定。

利用蒙特卡洛積分計算渲染方程

利用蒙特卡洛積分計算渲染方程可寫作：

\[L_{o}=L_{e}\left(x, \omega_{o}\right)+\frac{1}{N} \sum_{j=1}^{N} \frac{L_{i}\left(x, \omega_{o}\right) f\left(\omega_{i}, \omega_{o}\right)\langle\cos \theta\rangle}{p d f\left(\omega_{i}\right)} \]

當\(pdf(\omega_i)\)的分布與被積函數，即三個乘積項\(L_{i}\left(x, \omega_{o}\right) f\left(\omega_{i}, \omega_{o}\right)\langle\cos \theta\rangle\)相似時，能有效減小方差。

待解決問題

本文試圖通過結合多種采樣方法，例如重采樣重要性采樣、蓄水池采樣等多種方法，達到渲染方程積分計算快速收斂，提升積分計算穩定性的目的。

方法

現有技術

重要性采樣

如上述對蒙特卡洛積分收斂方差和速度的分析，盡可能選擇概率密度分布函數與被積函數形狀類似，或成正比時，積分收斂更快。重要性采樣即一種選擇特殊概率密度分布函數\(p(x)\)進行采樣的策略。

RIS(Resampled Importance Sampling) 重采樣重要性采樣

Importance Resampling重要性重采樣

假設我們對 被積函數 \(f\) 有了一個比較好的估計 \(g\)，需要生成滿足分布 \(g\) 的 樣本，有四種方法，分別是（1）逆變換采樣（2）拒絕采樣（3）Metropolis 采樣（4）Importance resampling。這里主要介紹Importance resampling。

Importance Resamping是計算統計學中從復雜分布生成樣本的一種常用方法。它在文獻4，2015年Talbot的碩士論文中充分闡述。主要思想是，由於生成\(g\) 的分布很復雜，我們先生成M個樣本，再按照一定的權重從M中取N個樣本，可以產生和\(g\)近似的分布。它的主要步驟是：

• 由於生成\(g\) 的分布很復雜或者並不歸一化，我們選擇另一個容易的分布 \(p\)，從概率分布函數\(p\)中選取M個樣本candidate，得到\(X=\{x_1,x_2,...,x_M\}\)
• 在所有M個樣本中，對任意第\(j\)個樣本，均計算權重\(w_{j}=\frac{g\left(x_{j}\right)}{p\left(x_{j}\right)}\)
• 從M中生成N(N ≤ M)個樣本：\(\{y_1,y_2,...,y_N\}\)，樣本對應的權重\(w\)作為概率比例

這樣生成的N個樣本\(\{y_1,y_2,...,y_N\}\)和\(g\)的分布近似。

RIS(Resampled Importance Sampling)

套用Importance Resamping獲得N個樣本，來做蒙特卡洛積分的方法我們稱之為 RIS。

由上文，蒙特卡洛積分可寫作以下形式，需要的前提是\(g\)的樣本容易獲得且\(g\)的分布歸一化。

\[L \approx \hat{L}_{M C}=\frac{1}{N} \sum_{i=1}^{N} \frac{f\left(y_{i}\right)}{g\left(y_{i}\right)} \]

首先，我們通過Importance Resamping獲得關於分布\(g\)樣本，Importance Resamping在\(g\)的分布不歸一化時，也適用。Importance Resamping生成的樣本近似於分布$g \(,為了維持蒙特卡洛積分的無偏性，必須引入一個權重函數\)\omega$，這個權重函數是M個Candidate，以及N個樣本的函數，可寫作：

\[\hat{L}_{r i s}=\frac{1}{N} \sum_{i=1}^{N} \frac{f\left(y_{i}\right)}{\hat{q}\left(y_{i}\right)} \cdot w^{*}\left(X_{1}, \ldots, X_{M}, y_{1}, \ldots, y_{N}\right) \]

推導(見附錄)可得：

\[w^{*}\left(X_{1}, \ldots, X_{M}, y_{1}, \ldots, y_{N}\right)=\frac{1}{M} \sum_{j=1}^{M} w\left(x_{j}\right) \]

結合上述兩個方程可得：

\[\hat{L}_{r i s}=\frac{1}{N} \sum_{i=1}^{N}\left(\frac{f\left(y_{i}\right)}{\hat{q}\left(y_{i}\right)}\right) \cdot \frac{1}{M} \sum_{j=1}^{M} w\left(x_{j}\right) \]

滿足以下兩個條件，可保證無偏，

• 在\(f\)的任意非零位置， \(g\)和 \(p\)都是正的。
• M和N均大於0

光線追蹤中的應用

對於光線追蹤求蒙特卡洛積分，被積函數\(L_{i}\left(x, \omega_{o}\right) f\left(\omega_{i}, \omega_{o}\right)\langle\cos \theta\rangle\)，計算生成M個Candidate，及其對應的\(\omega\)，這里N取1，每次取一個樣本，通過RIS可獲得關於被積函數的無偏估計。

Reservoir sampling蓄水池采樣

蓄水池采樣是一類采樣算法，給定 M 個元素，一次性在其中隨機選擇 N個元素作為一個子集。通常，N是 定義為一個小的常數，M 通常非常大，對應的數據無法本地存儲，或者M個元素是一個數據流。

基本方法和推導

假設數據序列的規模為M，或數據大小未知時，當前數據流一共通過了\(M_{sofar}\)，需要采樣的數量的為 N。

• 首先構建一個可容納 N 個元素的數組，將序列的前 N個元素放入數組中
• 然后從第 N+1 個元素開始，以\(N/M_{so far}\)的概率來決定該元素是否被替換到數組中,數組中被替換的概率相同
• 當遍歷完所有元素之后，數組中剩下的N個元素即為所需采取的樣本

推導如下：

持續不斷的序列流過，對於蓄水池N個元素中第 i 個數（i≤N）:

• 在N步之前，被選中的概率為 1。
• 當走到第 N+1 步時，被第 N+1 個元素替換位置i處的概率為：N+1 個元素被選中的概率 * \(i\)號位被選中替換的概率，即為 \(\frac{N}{N+1} \times \frac{\mathbf{1}}{N}=\frac{1}{N+1}\)，則被保留的概率為

\[1-\frac{N}{N+1} \times \frac{\mathbf{1}}{N}=\frac{N}{N+1} \]

• 依次類推，不被 N+2 個元素替換的概率為:

\[1-\frac{N}{N+2} \times \frac{\mathbf{1}}{N}=\frac{N+1}{N+2} \]

• 則運行到第 n 步時，被保留的概率 = 被選中的概率 * 不被替換的概率，即：

\[1 \times \frac{N}{N+1} \times \frac{N+1}{N+2} \times \frac{N+2}{N+3} \times \ldots \times \frac{M-1}{M}=\frac{N}{M} \]

因此對於每個元素，被保留的概率都是\(N/M\).當M的數量不確定時，累計的數據流個數為\(M_{sofar}\), 蓄水池中每個元素的概率為\(N/M_{so far}\).

加權蓄水池采樣

加權蓄水池采樣是蓄水池采樣算法的拓展，第$i \(個元素有對應的概率\)w_i\(，對應的第\)i \(個元素被保留的概率為\)\sum_{0}^{{i-1}w_{j}/\sum_{0}}w_{j}$。當蓄水池大小為1，即N = 1，對應計算的偽代碼如下：

結合加權蓄水池采樣的RIS

結合加權蓄水池采樣的RIS采樣比較直觀，步驟如下：

• 每個像素、相鄰像素、以及時域可復用的像素生成樣本流
• 每次更新Reservoir中樣本、權重和以及總的樣本數
• 拿到更新后的Reservoir，根據RIS無偏估計，進行采樣

實驗結果

下面給出原始文獻的數據，將應用了蓄水池流式采樣的方法ReSTIR與2015年Talot的RIS同樣時間所得到的效果做比較。可以看出利用蓄水池采樣的方法，相比原有方法提升20%~30%。

參考文獻

[1]Talbot J F. Importance resampling for global illumination[J]. 2005.

[2]Bitterli, Benedikt, et al. "Spatiotemporal reservoir resampling for real-time ray tracing with dynamic direct lighting." ACM Transactions on Graphics (TOG) 39.4 (2020): 148-1.

[3]Chris Wyman."WEIGHTED RESERVOIR SAMPLING:RANDOMLY SAMPLING STREAMS".Ray Tracing Gems II

[4]Justin F. Talbot. 2005. Importance Resampling for Global Illumination. Masters Thesis.Brigham Young University. https://diglib.eg.org/handle/10.2312/EGWR.EGSR05.139-146

附錄

RIS采樣無偏的證明：

由

\[E_{X Y}(X Y)=E_{X}\left(E_{Y}(Y \mid X) \cdot X\right) \]

可得：

\[E\left(\hat{L}_{r i s}\right)=E\left[E\left(\frac{1}{N} \sum_{i=1}^{N} \frac{f\left(Y_{i}\right)}{g\left(Y_{i}\right)} \mid X_{1}, \ldots, X_{M}\right) \cdot \frac{1}{M} \sum_{j=1}^{M} w\left(X_{j}\right)\right] \]

\[E\left(\hat{L}_{r i s}\right)=E\left[E\left(\frac{f(Y)}{g(Y)} \mid X_{1}, . ., X_{M}\right) \cdot \frac{1}{M} \sum_{j=1}^{M} w\left(X_{j}\right)\right] \]

由於Y是從樣本里隨機選出，因此可以寫成下面的求和形式：

\[\begin{aligned} E\left(\hat{L}_{r i s}\right) &=E\left[\sum_{k=1}^{M}\left(\frac{f\left(X_{k}\right)}{g\left(X_{k}\right)} \cdot \frac{w\left(X_{k}\right)}{\sum_{j=1}^{M} w\left(X_{j}\right)}\right) \cdot \frac{1}{M} \sum_{j=1}^{M} w\left(X_{j}\right)\right] \\ &=E\left[\frac{1}{M} \sum_{k=1}^{M}\left(\frac{f\left(X_{k}\right)}{g\left(X_{k}\right)} \cdot w\left(X_{k}\right)\right)\right] \\ &=E\left(\frac{f(X)}{g(X)} \cdot w(X)\right) \end{aligned} \]