Metatron Dev. V: 抽样优化

Path Guiding
#

路径引导基于已有样本估计场景贡献, 通过MIS选择BSDF或PG, 结果无偏.

ReSTIR PG
#

局部路径引导拟合的理想分布为后缀贡献:

\begin{equation} p(\omega_i|\mathbf{p}, \omega_o) \propto f(\mathbf{p}, \omega_o, \omega_i)L_i(\mathbf{p}, \omega_i)|\cos\theta| \end{equation}

ReSTIR PT将路径贡献作为目标分布 $\hat{p}(\mathbf{x})=f(\mathbf{x})$ , 记 $\mathcal{A}$ 为场景表面, $W_e$ 为传感器响应, $h_x$ 为像素 $x$ 的重建滤波器, $L_e(\mathbf{p}_n)$ 为末端顶点自发光, 路径贡献函数为:

\begin{equation} f(\mathbf{x}) = W_e(\mathbf{p}_0 \to \mathbf{p}_1)G(\mathbf{p}_0 \leftrightarrow \mathbf{p}_1) \prod_{i=1}^{n-1} f(\mathbf{p}_{i+1} \to \mathbf{p}_i \to \mathbf{p}_{i-1})G(\mathbf{p}_i \leftrightarrow \mathbf{p}_{i+1}) L_e(\mathbf{p}_n) \end{equation}

定义逐顶点的传输因子 $T$ 并简化上式:

\begin{equation} T(\mathbf{p}_i) = \begin{cases} f(\mathbf{p}_{i+1} \to \mathbf{p}_i \to \mathbf{p}_{i-1})G(\mathbf{p}_i \leftrightarrow \mathbf{p}_{i+1}), & i > 1\\ W_e(\mathbf{p}_0 \to \mathbf{p}_1)G(\mathbf{p}_0 \leftrightarrow \mathbf{p}_1)f(\mathbf{p}_2 \to \mathbf{p}_1 \to \mathbf{p}_0)G(\mathbf{p}_1 \leftrightarrow \mathbf{p}_2), & i = 1 \end{cases} \end{equation}

\begin{equation} f([\mathbf{p}_0, \dots, \mathbf{p}_n]) = \prod_{i=1}^{n-1}T(\mathbf{p}_i)\ L_e(\mathbf{p}_n) \end{equation}

记像素 $x$ 的路径空间为 $\Omega_x$ :

\begin{equation} C_x = \frac{1}{\int_{\Omega_x} f(\mathbf{x})\mathrm{d}\mathbf{x}} \end{equation}

全体像素构成对整个路径空间 $\Omega=\bigcup_{n=2}^\infty \mathcal{A}^n$ 的分层抽样, 其密度为各像素分布的混合:

\begin{equation} p(\mathbf{x}) = \Phi(\mathbf{p}_0, \mathbf{p}_1)f(\mathbf{x}), \quad \Phi(\mathbf{p}_0, \mathbf{p}_1) = \frac{1}{N}\sum_{x=1}^N C_x[h_x(\mathbf{p}_0 \to \mathbf{p}_1) > 0] \end{equation}

记 $\mathrm{d}\mathbf{p}_{a:b}=\mathrm{d}\mathbf{p}_a \cdots \mathrm{d}\mathbf{p}_b$ , 约定 $a > b$ 时积分退化为被积函数本身:

\begin{equation} p(\mathbf{p}_{i+1}|\mathbf{p}_0, \dots, \mathbf{p}_i) = \frac{p(\mathbf{p}_0, \dots, \mathbf{p}_{i+1})}{p(\mathbf{p}_0, \dots, \mathbf{p}_i)} = \frac{\sum_{n=i+1}^\infty \int_{\mathcal{A}^{n-i-1}} p(\mathbf{x})\mathrm{d}\mathbf{p}_{i+2:n}} {\sum_{n=i+1}^\infty \int_{\mathcal{A}^{n-i}} p(\mathbf{x})\mathrm{d}\mathbf{p}_{i+1:n}} \end{equation}

$\Phi$ 只依赖 $\mathbf{p}_0$ 与 $\mathbf{p}_1$ , 条件 $i \geq 1$ 已固定二者, 同时 $\prod_{t=1}^{i-1}T(\mathbf{p}_t)$ 与积分变量无关, 一并约去:

\begin{equation} \begin{aligned} p(\mathbf{p}_{i+1}|\mathbf{p}_0, \dots, \mathbf{p}_i) &= \frac{\sum_{n=i+1}^\infty \int_{\mathcal{A}^{n-i-1}} \prod_{t=1}^{n-1}T(\mathbf{p}_t)L_e(\mathbf{p}_n)\mathrm{d}\mathbf{p}_{i+2:n}} {\sum_{n=i+1}^\infty \int_{\mathcal{A}^{n-i}} \prod_{t=1}^{n-1}T(\mathbf{p}_t)L_e(\mathbf{p}_n)\mathrm{d}\mathbf{p}_{i+1:n}}\\ &= \frac{T(\mathbf{p}_i)\sum_{n=i+1}^\infty \int_{\mathcal{A}^{n-i-1}} \prod_{t=i+1}^{n-1}T(\mathbf{p}_t)L_e(\mathbf{p}_n)\mathrm{d}\mathbf{p}_{i+2:n}} {\sum_{n=i+1}^\infty \int_{\mathcal{A}^{n-i}} \prod_{t=i}^{n-1}T(\mathbf{p}_t)L_e(\mathbf{p}_n)\mathrm{d}\mathbf{p}_{i+1:n}} \end{aligned} \end{equation}

分子中 $T(\mathbf{p}_i)$ 之后的求和即入射辐亮度:

\begin{equation} L_i(\mathbf{p}_{i+1} \to \mathbf{p}_i) = \sum_{n=i+1}^\infty \int_{\mathcal{A}^{n-i-1}} \prod_{t=i+1}^{n-1}T(\mathbf{p}_t)L_e(\mathbf{p}_n)\mathrm{d}\mathbf{p}_{i+2:n} \end{equation}

分母相比入射辐亮度少了自发光:

\begin{equation} L_i(\mathbf{p}_i \to \mathbf{p}_{i-1}) - L_e(\mathbf{p}_i \to \mathbf{p}_{i-1}) = \sum_{n=i+1}^\infty \int_{\mathcal{A}^{n-i}} \prod_{t=i}^{n-1}T(\mathbf{p}_t)L_e(\mathbf{p}_n)\mathrm{d}\mathbf{p}_{i+1:n}\\ \end{equation}

分母只依赖 $\mathbf{p}_{i-1}$ 与 $\mathbf{p}_i$ , 相对 $\mathbf{p}_{i+1}$ 为常数, 因此 $\mathbf{p}_0, \dots, \mathbf{p}_{i-2}$ 对条件分布没有影响:

\begin{equation} p(\mathbf{p}_{i+1}|\mathbf{p}_{i-1}, \mathbf{p}_i) \propto f(\mathbf{p}_{i+1} \to \mathbf{p}_i \to \mathbf{p}_{i-1})G(\mathbf{p}_i \leftrightarrow \mathbf{p}_{i+1})L_i(\mathbf{p}_{i+1} \to \mathbf{p}_i) \end{equation}

转立体角测度结果如下, 即服从路径贡献分布的路径, 局部弹射方向的条件分布为理想的局部引导分布. ReSTIR样本只是无偏权重, 初始样本质量差时相关性较强.

\begin{equation} p(\omega_i|\mathbf{p}_i, \omega_o) \propto f(\mathbf{p}_i, \omega_o, \omega_i)L_i(\mathbf{p}_i, \omega_i)|\cos\theta| \end{equation}

$p(\omega_i|\mathbf{p}, \omega_o)$ 需要拟合7D模型因此实时不可行. 划分场景为空间网格, 在格内拟合求平均以消去 $\mathbf{p}$ , 降维到4D的 $p(\omega_i|\omega_o)$ . 漫反射BSDF为常数, 光滑镜面更依赖BSDF抽样, 因此通过积分消去 $\omega_o$ 后, 只有粗糙镜面效果较差, 可进一步降维到2D:

\begin{equation} p(\omega_i) = \int_{\mathcal{H}^2} p(\omega_i|\omega_o)p(\omega_o)\mathrm{d}\omega_o \propto L_i(\omega_i)|\cos\theta| \int_{\mathcal{H}^2} f(\omega_o, \omega_i)p(\omega_o)\mathrm{d}\omega_o \end{equation}

$p(\omega_i)$ 以vMF拟合, $\mu_k \in \mathcal{S}^2$ 为单位平均方向, $\kappa_k \geq 0$ 为集中度, $\sum_k \pi_k = 1$ :

\begin{equation} \mathcal{V}(\omega;\Theta) = \sum_{k=1}^K \pi_k v(\omega;\mu_k, \kappa_k) = \sum_{k=1}^K \frac{\kappa}{4\pi\sinh\kappa}e^{\kappa\mu \cdot \omega} \end{equation}

展开 $\sinh$ 可得vMF随着 $\kappa$ 增长而收窄集中于 $\mu$ 的波瓣:

\begin{equation} v(\omega;\mu, \kappa) = \frac{\kappa}{2\pi(1 - e^{-2\kappa})}e^{-\kappa(1 - \mu \cdot \omega)} \approx \frac{\kappa}{2\pi}e^{-\kappa(1 - \mu \cdot \omega)} \end{equation}

以EM迭代求解, E步固定参数, 计算责任:

\begin{equation} \gamma_k(\omega_n) = \frac{\pi_k v(\omega_n;\mu_k, \kappa_k)}{\sum_{j=1}^K \pi_j v(\omega_n;\mu_j, \kappa_j)} \end{equation}

M步固定责任, 更新参数, $\epsilon = 0.01$ 防止分量权重归零:

\begin{equation} \begin{aligned} w_k = \sum_{n=1}^N \gamma_k(\omega_n), \quad r_k = \sum_{n=1}^N \gamma_k(\omega_n)\omega_n\\ \pi_k = \frac{w_k + \epsilon}{\sum_{j=1}^K (w_j + \epsilon)}, \quad \mu_k = \frac{r_k}{\|r_k\|}, \quad \bar{R}_k = \frac{\|r_k\|}{w_k} \end{aligned} \end{equation}

方向越集中 $\bar{R}_k \in [0, 1]$ 越接近1, 由其反解 $\kappa_k$ :

\begin{equation} \kappa_k \approx \frac{\bar{R}_k(3 - \bar{R}_k^2)}{1 - \bar{R}_k^2} \end{equation}

ReSTIR PT无偏权重包含足够信息, 不复用历史分布, 记录最终路径所有顶点以更新引导.

VXPG
#

将第二次弹射顶点的NEE采样结果注入体素, 估计漫反射辐照度:

\begin{equation} E(v_i) = \frac{1}{N}\sum_{\mathbf{x}_2 \in v_i} L_l(\mathbf{x}_2) \end{equation}

以 $32 \times 32$ 图块中心像素为矩心, 单次SLIC聚类计算超像素, $\mathbf{p}$ , $\mathbf{n}$ 为世界空间位置和法线, $\mathbf{u}$ 为像素坐标, $w_u$ 为超参数, 距离函数如下:

\begin{equation} \text{dist}_p(x, y) = |\mathbf{p}_x - \mathbf{p}_y|^2 + w_u |\mathbf{u}_x - \mathbf{u}_y|^2 + \begin{cases} \begin{aligned} &0,&\ \text{if}\ \mathbf{n_x}\cdot\mathbf{n_y} > 0.1\\ &1000000,&\ \text{otherwise} \end{aligned} \end{cases} \end{equation}

在屏幕空间分层抽样128条路径, 所有体素向所有样本路径的 $\mathbf{x}_2$ 发射光线以验证可见性, 得128位位域 $R$ , 执行K-means聚类. $\oplus$ 为异或, $w_E$ 为超参数, 距离函数如下:

\begin{equation} \text{dist}_v(x, y) = \text{countbits}(R_x \oplus R_y) + w_E |E_x - E_y| \end{equation}

聚类后依据超像素与超体素组成的元组分桶, 从桶中抽取32个路径样本统计平均路径通量:

\begin{equation} \bar{T}(p', v') = \frac{1}{N}\sum_{i = 1}^N f(\mathbf{p}_2 \rightarrow \mathbf{p}_1 \rightarrow \mathbf{p}_0) G(\mathbf{p}_1 \leftrightarrow \mathbf{p}_2) \end{equation}

依据平均路径通量重要性抽样得超像素与超体素, 再抽样桶中的体素. 光栅化执行体素化, 得体素中所有三角形的包围盒, $A$ 为包围盒最大面的面积, 抽样权重如下:

\begin{equation} \phi(v_i) = E(v_i)A(v_i) \end{equation}

得到体素后, 使用球面三角形抽样确定与包围盒可见面的相交点, 发射光线求交.

RCPG
#

基于辐射度级联摆放探针, 使用八面体纹理存储辐射度, 转为立体角的Jacobian如下:

\begin{equation} |J| = \left(1 + 2\min(1 - |u| - |v|, 0) - 2|u| - 2|v| + 2|uv| + 2u^2 + 2v^2\right)^{-\frac{3}{2}} \end{equation}

贴近物体的探针命中距离短, 需要的采样率低, 反之采样率高. 不同LOD的探针对应不同的距离区间, 基于区间远端与探针体素内嵌球形成的锥体, 可估计角采样Nyquist频率. 令 $w$ 为体素边长, $d$ 为与体素中心的距离, $c$ 为超参数, 采样频率估计如下:

\begin{equation} \theta_{\min} = 2\arcsin\frac{w}{2d_{\max}}, \quad f_{\max} = \frac{1}{\theta_{\min}} < \frac{f_s}{2}, \quad d_{\max} = \frac{cw}{2\sin\frac{1}{f_s}} \end{equation}

所有探针纹素每帧更新, 只追踪距离区间, 命中确定二值透明. 下采样父级探针, 做透明度混合填充区间外辐射度. 探针通过指数混合累积历史, 命中时查询历史探针, 存储无限弹射辐射度. 每帧对场景执行体素化, 剔除空体素对应的探针.

遍历表面对应的LOD 0探针的纹素, 作为解析面光通过LTC计算贡献, 执行功率重要性抽样, 逐级查询父级纹素LTC并抽样. 令 $\mathbf{v}$ 为几何顶点, 面光辐照度如下:

\begin{equation} E = \frac{1}{2\pi}\sum_{i=1}^m \arccos(\mathbf{v}_i \cdot \mathbf{v}_j) \frac{\mathbf{v}_i \times \mathbf{v}_j}{|\mathbf{v}_i \times \mathbf{v}_j|} \cdot \mathbf{n}, \quad j = (i + 1) \bmod m \end{equation}

使用三参数LTC以减少计算量:

\begin{equation} M^{-1} = \begin{pmatrix} a & 0 & b\\ 0 & 1 & 0\\ 0 & 0 & c \end{pmatrix}, \quad a, b, c \in [0, 1] \end{equation}

MCPG
#

基于混合vMF估计分布, 链状态为单个波瓣, 存于多分辨率哈希网格, 同时有均匀静态网格保证LOD边界的状态交换. 基于三线性权重抽取 $N_\mathrm{mc}$ 个候选顶点, 以亮度为权重执行RIS.

$f_\mathrm{mc}$ 为当前样本入射亮度估计, 路径尾部查询辐照度缓存, 接受概率如下, 分母可加速预热:

\begin{equation} p_\mathrm{accept} = \min\left(\frac{f_\mathrm{mc}}{\mathrm{sum}_\mathrm{mc}/N_\mathrm{mc}}, 1\right) \end{equation}

以最大似然估计更新状态, 混合因子 $\alpha = \max(\frac{1}{N}, \alpha_\mathrm{min})$ . 更新后抽样写回位置, 读写两侧的随机访问使状态扩散. 转移步骤不满足细致平衡, 稳态分布不保证收敛到目标分布.

NASGPG
#

归一化各向异性球面高斯(NASG)基于正交坐标系 $[\mathbf{x}, \mathbf{y}, \mathbf{z}]$ 定义, $\mathbf{z}$ 为波瓣轴, $\lambda$ 为锐度, $a$ 为各向异性, $a = 0$ 时退化为球面高斯. 记 $u = \frac{\mathbf{v} \cdot \mathbf{z} + 1}{2}$ , $e = \frac{a(\mathbf{v} \cdot \mathbf{x})^2}{1 - (\mathbf{v} \cdot \mathbf{z})^2}$ , 形式如下:

\begin{equation} G(\mathbf{v};[\mathbf{x}, \mathbf{y}, \mathbf{z}], \lambda, a) = \begin{cases} \begin{aligned} &\exp\left(2\lambda u^{1 + e} - 2\lambda\right)u^e, &\ \mathbf{v} \neq \pm\mathbf{z}\\ &1, &\ \mathbf{v} = \mathbf{z}\\ &0, &\ \mathbf{v} = -\mathbf{z} \end{aligned} \end{cases} \end{equation}

$u^e$ 为球坐标换元的Jacobian, 因此NASG有闭式积分, 可归一化:

\begin{equation} K = \int_{\mathcal{S}^2} G(\mathbf{v};[\mathbf{x}, \mathbf{y}, \mathbf{z}], \lambda, a)\mathrm{d}\omega = \frac{2\pi(1 - e^{-2\lambda})}{\lambda\sqrt{1 + a}} \end{equation}

正交坐标系以欧拉角 $\theta, \phi, \tau$ 参数化, $\mathbf{y} = \mathbf{z} \times \mathbf{x}$ , 单个NASG分量只需 $\cos\theta$ , $\sin\phi$ , $\cos\phi$ , $\sin\tau$ , $\cos\tau$ , $\lambda$ , $a$ 七个标量表示:

\begin{equation} \mathbf{z} = \begin{pmatrix} \cos\phi\sin\theta\\ \sin\phi\sin\theta\\ \cos\theta \end{pmatrix}, \quad \mathbf{x} = \begin{pmatrix} \cos\theta\cos\phi\cos\tau - \sin\phi\sin\tau\\ \cos\theta\sin\phi\cos\tau + \cos\phi\sin\tau\\ -\sin\theta\cos\tau \end{pmatrix} \end{equation}

神经网络为4层128宽无偏置MLP, 输入为 $\mathbf{p}$ , $\omega_o$ , $\mathbf{n}$ , $\mathbf{p}$ 归一化后应用one-blob编码, 即分为 $k$ 个等宽区间, 每个区间对应 $\sigma = \frac{1}{k}$ 的高斯核, 执行积分:

\begin{equation} \text{ob}(x)_i = \int_{\frac{i - 1}{k}}^{\frac{i}{k}} \frac{1}{\sqrt{2\pi}\sigma}e^{-\frac{(t - x)^2}{2\sigma^2}}\mathrm{d}t, \quad i = 1, \dots, k \end{equation}

输出为 $8N + 1$ 维, N为NASG分量数量, 8为NASG7个标量与其权重 $A$ , 1为MIS抽样概率 $c$ . 引导分布为归一化分量的混合, 与BSDF按 $c$ 做MIS:

\begin{equation} \begin{aligned} \hat{q}(\omega_i) &= cq(\omega_i) + (1 - c)p_f(\omega_i)\\ &= c\sum_{k=1}^N A_k\frac{G_k(\omega_i)}{K_k} + (1 - c)p_f(\omega_i) \end{aligned} \end{equation}

以KL散度衡量 $q$ 与目标分布 $p(\omega_i) = Cf(\mathbf{p}, \omega_o, \omega_i)L_i(\mathbf{p}, \omega_i)|\cos\theta|$ 的差异, 令 $\gamma$ 为NASG参数, $p$ 与 $\gamma$ 无关, 梯度如下:

\begin{equation} \begin{aligned} \nabla_\gamma D_{KL}(p\|q) &= \nabla_\gamma\int_{\mathcal{S}^2} p(\omega_i)(\log p(\omega_i) - \log q(\omega_i;\gamma))\mathrm{d}\omega_i\\ &= -\int_{\mathcal{S}^2} p(\omega_i)\nabla_\gamma\log q(\omega_i;\gamma)\mathrm{d}\omega_i \end{aligned} \end{equation}

基于 $\hat{q}$ 得单样本估计. 归一化常数 $C$ 未知, 但只对梯度整体缩放, Adam归一化矩后可约去:

\begin{equation} \nabla_\gamma D_{KL}(p\|q) \approx -\frac{p(\omega_i)}{\hat{q}(\omega_i)}\nabla_\gamma\log q(\omega_i;\gamma) \end{equation}

$D_{KL}(p\|q)$ 与 $c$ 无关, $D_{KL}(p\|\hat{q})$ 可对 $c$ 求导, 但只优化 $\hat{q}$ 时 $c$ 会倾向0, 因为初始 $q$ 质量较差. 因此以 $D_{KL}(p\|q)$ 为主项保证 $q$ 持续更新, 混入 $D_{KL}(p\|\hat{q})$ 以学习 $c$ , $e = 0.2$ :

\begin{equation} \text{loss} = eD_{KL}(p\|\hat{q}) + (1 - e)D_{KL}(p\|q) \end{equation}

NPMPG
#

使用神经网络隐式表示场景, 以 $\Phi$ 为可训练参数, 连续地将位置映射为vMF混合参数:

\begin{equation} \text{NPM}(\mathbf{x}\mid\Phi) = \hat{\Theta}(\mathbf{x}), \quad \mathcal{V}(\omega_i\mid\hat{\Theta}(\mathbf{x})) \propto L_i(\mathbf{x}, \omega_i) \end{equation}

使用可训练多分辨率空间编码处理高频, 定义 $L$ 级均匀LOD网格, 体素存 $F$ 维特征, 查询时拼接各层的三线性插值结果得到 $G(\mathbf{x})$ :

\begin{equation} G(\mathbf{x}\mid\Phi_E) = \bigoplus_{l=1}^L \text{trilinear}(\mathbf{x}, V_l[\mathbf{x}]) \end{equation}

MLP为3层64宽ReLU, 输出 $K$ 个vMF. $\kappa', \lambda', \theta', \phi'$ 为原始输出, $(\theta, \phi)$ 为 $\mu$ 的归一化球坐标:

\begin{equation} \begin{aligned} \kappa_i = \exp(\kappa_i'), \quad \lambda_i = \frac{\exp(\lambda_i')}{\sum_{j=1}^K\exp(\lambda_j')}\\ \theta_i = \frac{1}{1 + \exp(-\theta_i')}, \quad \phi_i = \frac{1}{1 + \exp(-\phi_i')} \end{aligned} \end{equation}

以KL散度为目标做小批量随机梯度下降. 令 $\mathcal{D} \propto L_i$ 为目标分布, 路径顶点更新附近体素, $\tilde{p}$ 为BSDF与引导分布组合的实际抽样分布, 梯度估计如下:

\begin{equation} \begin{aligned} \nabla_\Theta D_{KL}(\mathcal{D}\|\mathcal{V};\Theta) &= \nabla_\Theta\int_\Omega \mathcal{D}(\omega)\log\frac{\mathcal{D}(\omega)}{\mathcal{V}(\omega\mid\hat{\Theta})}\mathrm{d}\omega\\ &\approx \nabla_\Theta\frac{1}{N}\sum_{j=1}^N\frac{\mathcal{D}(\omega_j)}{\tilde{p}(\omega_j\mid\hat{\Theta})}\log\frac{\mathcal{D}(\omega_j)}{\mathcal{V}(\omega_j\mid\hat{\Theta})}\\ &= -\frac{1}{N}\sum_{j=1}^N\frac{\mathcal{D}(\omega_j)\nabla_\Theta\mathcal{V}(\omega_j\mid\hat{\Theta})}{\tilde{p}(\omega_j\mid\hat{\Theta})\mathcal{V}(\omega_j\mid\hat{\Theta})} \end{aligned} \end{equation}

学习完整被积函数时额外输入 $\omega_o$ , 目标分布改为 $f_s L_i\cos\theta_i$ , 余弦项以固定vMF波瓣近似. $\mathbf{n}$ 与粗糙度 $r$ 作为辅助特征输入, $\omega_o$ 与 $\mathbf{n}$ 使用球谐编码.

Radiance Cache
#

辐射度缓存在探针中存储辐射度, 命中后直接查询缓存, 因此有偏.

ORCA
#

根据BSDF抽样概率和粗糙度决定舍弃概率, 根据预算归一化以避免光滑场景光线超支:

\begin{equation} s'_i=\frac{b s_i}{\sum_{i=1}^N s_i} \end{equation}

稀疏光线完整追踪, 根据第二次弹射顶点信息计算hash, 更新最细LOD体素的累积辐亮度, LOD与相机距离相关. 下采样以更新LOD体素, 其余光线单次弹射并查询缓存, 模拟重连接. 只用本帧数据, 不做时域累积.

SHARC
#

依据世界空间顶点和LOD计算hash, 因此跨帧hash一致, 体素辐亮度逐帧累积.

NRC
#

单个MLP缓存散射辐亮度 $L_s(\mathbf{x}, \omega)$ , 路径足迹足够大时误差被模糊, 令 $p$ 为BSDF抽样概率, $\theta_1$ 为主顶点处视线与法线夹角, $c = 0.01$ , 路径足迹定义如下:

\begin{equation} \begin{aligned} &a(\mathbf{x}_1 \cdots \mathbf{x}_n) = \left(\sum_{i=2}^n\sqrt{\frac{\|\mathbf{x}_{i-1} - \mathbf{x}_i\|^2}{p(\omega_i\mid\mathbf{x}_{i-1}, \omega)|\cos\theta_i|}}\right)^2\\ &a_0 = \frac{\|\mathbf{x}_0 - \mathbf{x}_1\|^2}{4\pi\cos\theta_1}, \quad a > ca_0 \end{aligned} \end{equation}

屏幕分块后每块抽样一条路径更新缓存, 高学习率与每帧多步导致闪烁, 因此推理时使用权重的指数移动平均, $\alpha = 0.99$ , $\eta_t$ 修正初期偏差, 不反馈到训练:

\begin{equation} \bar{W}_t = \frac{1 - \alpha}{\eta_t}W_t + \alpha\eta_{t-1}\bar{W}_{t-1}, \quad \eta_t = 1 - \alpha^t \end{equation}

MLP为7层64宽无偏置, 输出RGB. $\omega$ , $\mathbf{n}$ 转球坐标, 与 $1 - e^{-r}$ 一同做4区间one-blob编码, 漫反射与镜面反射率 $\alpha$ , $\beta$ 直接输入. 位置微小变化引起辐亮度剧变, 改用频率编码:

\begin{equation} \text{freq}(x) = \left(\sin(2^0\pi x), \sin(2^1\pi x), \dots, \sin(2^{11}\pi x)\right) \end{equation}

MLP输出乘 $\alpha + \beta$ 得到近似出射辐射度. 由于 $L_s$ 为无偏估计量, 使用相对L2损失保证梯度无偏, $\text{sg}$ 为停止梯度, 损失函数如下:

\begin{equation} \mathcal{L}_2(L_s, \hat{L}_s) = \frac{(L_s - \hat{L}_s)^2}{\text{sg}(\hat{L}_s)^2 + \epsilon} \end{equation}

Markov Chain
#

Metropolis-Hastings
#

令 $\mathbf{x}_k$ 为顶点数为 $k$ 的路径, $p_i$ 为 $X_i$ 的密度, 转移函数 $K(\mathbf{x} \to \mathbf{y})$ 为 $\mathbf{x}$ 转移至 $\mathbf{y}$ 的概率密度, 满足 $\sum_{k=1}^\infty\int_{\mathbf{y}_k} K(\mathbf{x} \to \mathbf{y}_k)\mathrm{d}\mathbf{y}_k = 1$ , Markov链的演化如下:

\begin{equation} p_i(\mathbf{x}) = \sum_{k=1}^\infty\int_{\mathbf{y}_k} K(\mathbf{y}_k \to \mathbf{x})p_{i-1}(\mathbf{y}_k)\mathrm{d}\mathbf{y}_k \end{equation}

稳态分布(stationary distribution) $p_\infty$ 为不动点即 $p_i = p_{i-1}$ . Metropolis-Hastings以提议密度 $T(\mathbf{x} \to \mathbf{y})$ 生成候选, 以接受概率 $a(\mathbf{x} \to \mathbf{y})$ 决定是否接受:

\begin{equation} p_i(\mathbf{x}) = p_{i-1}(\mathbf{x})(1 - \sum_{k=1}^\infty\int_{\mathbf{y}_k} T(\mathbf{x} \to \mathbf{y}_k)a(\mathbf{x} \to \mathbf{y}_k)\mathrm{d}\mathbf{y}_k) + \sum_{k=1}^\infty\int_{\mathbf{y}_k} p_{i-1}(\mathbf{y}_k)T(\mathbf{y}_k \to \mathbf{x})a(\mathbf{y}_k \to \mathbf{x})\mathrm{d}\mathbf{y}_k \end{equation}

给定目标分布 $\hat{p}$ , 细致平衡(detailed balance)条件如下:

\begin{equation} \hat{p}(\mathbf{x})T(\mathbf{x} \to \mathbf{y})a(\mathbf{x} \to \mathbf{y}) = \hat{p}(\mathbf{y})T(\mathbf{y} \to \mathbf{x})a(\mathbf{y} \to \mathbf{x}) \end{equation}

归一化常数 $\|\hat{p}\| = \sum_{k=1}^\infty\int_{\mathbf{x}_k} \hat{p}(\mathbf{x}_k)\,\mathrm{d}\mathbf{x}_k$ , 设 $p_{i-1} = \frac{\hat{p}}{\|\hat{p}\|}$ , 代入演化得 $\frac{\hat{p}}{\|\hat{p}\|}$ 为不动点:

\begin{equation} \begin{aligned} p_i(\mathbf{x}) &= p_{i-1}(\mathbf{x}) - \frac{1}{\|\hat{p}\|}\sum_{k=1}^\infty\int_{\mathbf{y}_k}( \hat{p}(\mathbf{x})T(\mathbf{x} \to \mathbf{y}_k)a(\mathbf{x} \to \mathbf{y}_k) - \hat{p}(\mathbf{y}_k)T(\mathbf{y}_k \to \mathbf{x})a(\mathbf{y}_k \to \mathbf{x}) )\mathrm{d}\mathbf{y}_k\\ &= p_{i-1}(\mathbf{x}) \end{aligned} \end{equation}

按如下方式定义 $a$ , 使得较大的 $\hat{p}(\mathbf{x})T(\mathbf{x} \to \mathbf{y})$ 总是被接受:

\begin{equation} a(\mathbf{x} \to \mathbf{y}) = \min(1, \frac{\hat{p}(\mathbf{y})T(\mathbf{y} \to \mathbf{x})}{\hat{p}(\mathbf{x})T(\mathbf{x} \to \mathbf{y})}) \end{equation}

不动点不蕴含唯一性, 例如取 $\Omega = \{1, 2, 3\}$ , 转移矩阵行为起点列为终点:

K = \begin{bmatrix} \frac{1}{2} & \frac{1}{2} & 0\\ \frac{1}{2} & \frac{1}{2} & 0\\ 0 & 0 & 1 \end{bmatrix}

令 $p_0 = (\alpha, \beta, \gamma)$ , 一步后即为不动点, 可见不动点与初始状态相关:

p_1 = \left(\frac{\alpha + \beta}{2}, \frac{\alpha + \beta}{2}, \gamma\right)

从任意 $\mathbf{x}$ 出发均能在有限步内到达任何 $\hat{p} > 0$ 的区域则不动点唯一. 若满足细致平衡, 已知 $p_0 = \frac{\hat{p}}{\|\hat{p}\|}$ 不动点为 $\frac{\hat{p}}{\|\hat{p}\|}$ , 因此任意初值均收敛到它:

\begin{equation} \lim_{i \to \infty}p_i = \frac{\hat{p}}{\|\hat{p}\|} \end{equation}

Metropolis Light Transport
#

基于BDPT实现MLT, 令 $p_0$ 为BDPT的PDF, 初始权重 $W_0 = \frac{\hat{p}(X_0)}{p_0(X_0)}$ , 变异保持 $W_i = W_{i-1}$ , 记 $\rho_i(w, \mathbf{x})$ 为 $(W_i, X_i)$ 的联合密度, 加权均衡条件(weighted equilibrium condition)如下:

\begin{equation} \int_\mathbb{R} w\,\rho_i(w, \mathbf{x})\,\mathrm{d}w = \hat{p}(\mathbf{x}) \end{equation}

$i = 0$ 时 $W_0$ 由 $X_0$ 确定, $\rho_0(w, \mathbf{x}) = \delta(w - \frac{\hat{p}(\mathbf{x})}{p_0(\mathbf{x})})p_0(\mathbf{x})$ , 代入直接满足. 变异与 $W$ 无关, 故 $\rho_i$ 与 $p_i$ 的演化相同:

\begin{equation} \rho_i(w, \mathbf{x}) = \rho_{i-1}(w, \mathbf{x})(1 - \sum_{k=1}^\infty\int_{\mathbf{y}_k} T(\mathbf{x} \to \mathbf{y}_k)a(\mathbf{x} \to \mathbf{y}_k)\mathrm{d}\mathbf{y}_k) + \sum_{k=1}^\infty\int_{\mathbf{y}_k} \rho_{i-1}(w, \mathbf{y}_k)T(\mathbf{y}_k \to \mathbf{x})a(\mathbf{y}_k \to \mathbf{x})\mathrm{d}\mathbf{y}_k \end{equation}

两侧乘 $w$ 并对 $w$ 积分, 基于加权均衡条件与细致平衡可得:

\begin{equation} \begin{aligned} \int_\mathbb{R} w\,\rho_i(w, \mathbf{x})\,\mathrm{d}w &= \hat{p}(\mathbf{x}) - \sum_{k=1}^\infty\int_{\mathbf{y}_k}( \hat{p}(\mathbf{x})T(\mathbf{x} \to \mathbf{y}_k)a(\mathbf{x} \to \mathbf{y}_k) - \hat{p}(\mathbf{y}_k)T(\mathbf{y}_k \to \mathbf{x})a(\mathbf{y}_k \to \mathbf{x}) )\mathrm{d}\mathbf{y}_k\\ &= \hat{p}(\mathbf{x}) \end{aligned} \end{equation}

令 $h$ 为滤波器权重, $I_j$ 为像素 $j$ 的真值, 可验证无偏:

\begin{equation} \begin{aligned} E[W_i h_j(X_i)] &= \sum_{k=1}^\infty\int_{\mathbf{x}_k}\int_\mathbb{R} w\,h_j(\mathbf{x}_k)\rho_i(w, \mathbf{x}_k)\,\mathrm{d}w\,\mathrm{d}\mathbf{x}_k\\ &= \sum_{k=1}^\infty\int_{\mathbf{x}_k} h_j(\mathbf{x}_k)\hat{p}(\mathbf{x}_k)\,\mathrm{d}\mathbf{x}_k\\ &= I_j \end{aligned} \end{equation}

ReSTIR MCMC
#

添加MCMC后的无偏权重如下:

\begin{equation} W_i = \frac{\hat{p}(X_{i-1})}{\hat{p}(X_i)}W_{i-1} = \frac{\hat{p}(X_0)}{\hat{p}(X_i)}W_0 \end{equation}

归纳证明 $W_i$ 无偏, 基于GRIS已知 $W_0$ 无偏, 设 $W_{i-1}$ 无偏, 记候选 $Z_{i-1} \sim T(X_{i-1} \to \cdot)$ :

\begin{equation} \begin{aligned} E[f(X_i)W_i] =\ & E[(1 - a(X_{i-1} \to Z_{i-1}))f(X_{i-1})W_{i-1}] + E[a(X_{i-1} \to Z_{i-1})f(Z_{i-1})\frac{\hat{p}(X_{i-1})}{\hat{p}(Z_{i-1})}W_{i-1}]\\ =\ & E[f(X_{i-1})W_{i-1}] + E[f(Z_{i-1})a(X_{i-1} \to Z_{i-1})\frac{\hat{p}(X_{i-1})}{\hat{p}(Z_{i-1})}W_{i-1}] - E[f(X_{i-1})a(X_{i-1} \to Z_{i-1})W_{i-1}] \end{aligned} \end{equation}

后两项含候选 $Z_{i-1}$ , 消去它才能归纳. $W_{i-1}$ 仅依赖 $X_{i-1}$ , 可从条件期望中作为常数提出:

\begin{equation} \begin{aligned} E[g(X_{i-1}, Z_{i-1})W_{i-1}] &= E\big[E[g(X_{i-1}, Z_{i-1})W_{i-1}|\mathcal{F}_{i-1}]\big]\\ &= E\big[W_{i-1}E[g(X_{i-1}, Z_{i-1})|\mathcal{F}_{i-1}]\big]\\ &= E\left[\left(\int_\Omega g(X_{i-1}, \mathbf{z})T(X_{i-1} \to \mathbf{z})\mathrm{d}\mathbf{z}\right)W_{i-1}\right] \end{aligned} \end{equation}

此时仅与 $X_{i-1}$ 有关, 代入即得二重积分:

\begin{equation} \begin{aligned} E[f(X_i)W_i] &= \int_\Omega f(\mathbf{x})\mathrm{d}\mathbf{x}\\ &+ \int_\Omega\int_\Omega f(\mathbf{z})a(\mathbf{x} \to \mathbf{z})\frac{\hat{p}(\mathbf{x})}{\hat{p}(\mathbf{z})}T(\mathbf{x} \to \mathbf{z})\mathrm{d}\mathbf{z}\mathrm{d}\mathbf{x}\\ &- \int_\Omega\int_\Omega f(\mathbf{x})a(\mathbf{x} \to \mathbf{z})T(\mathbf{x} \to \mathbf{z})\mathrm{d}\mathbf{z}\mathrm{d}\mathbf{x} \end{aligned} \end{equation}

改写细致平衡, 可抵消后两项, 即估计无偏:

\begin{equation} T(\mathbf{x} \to \mathbf{z})a(\mathbf{x} \to \mathbf{z}) = \frac{\hat{p}(\mathbf{z})}{\hat{p}(\mathbf{x})}T(\mathbf{z} \to \mathbf{x})a(\mathbf{z} \to \mathbf{x}) \end{equation}

变异为对 $\mathbf{u}_{i-1}$ 做高斯扰动, 即 $j \neq i$ 时 $\mathbf{p}^y_j = \mathbf{p}^x_j$ , 只需变换 $\mathbf{u}_{\{i - 1, i, i + 1\}}$ . $\mathbf{p}^y_{\{i+1, i+2\}}$ 固定在 $\mathbf{p}^x_{\{i+1, i+2\}}$ , 以变异后的前缀 $\mathbf{p}^y_{\{0, \dots, i\}}$ 为条件可得条件密度:

\begin{equation} p(\mathbf{p}^y_{\{i+1, i+2\}}|\mathbf{p}^y_{\{0, \dots, i\}}) = \delta(\mathbf{p}^y_{i+1} - \mathbf{p}^x_{i+1})\delta(\mathbf{p}^y_{i+2} - \mathbf{p}^x_{i+2}) \end{equation}

$a$ 需要提议密度比值, $\mathbf{u}_{i-1}$ 的密度只与 $|\mathbf{u}^x_{i-1} - \mathbf{u}^y_{i-1}|$ 有关, 在比值中约去. $\mathbf{u}^y_{\{i, i+1\}}$ 的密度为上式的变换 $p(\mathbf{p}^y_{\{i+1, i+2\}}|\mathbf{p}^y_{\{0, \dots, i\}})\left|\frac{\partial\mathbf{p}^y_{\{i+1, i+2\}}}{\partial\mathbf{u}^y_{\{i, i+1\}}}\right|$ , Dirac delta函数同样约去, 此时可得:

\begin{equation} \begin{aligned} \frac{T(\mathbf{y} \to \mathbf{x})}{T(\mathbf{x} \to \mathbf{y})} &= \frac{\left|\frac{\partial\mathbf{p}^x_{\{i+1, i+2\}}}{\partial\mathbf{u}^x_{\{i, i+1\}}}\right|}{\left|\frac{\partial\mathbf{p}^y_{\{i+1, i+2\}}}{\partial\mathbf{u}^y_{\{i, i+1\}}}\right|}\\ &= \left|\frac{\partial\mathbf{u}^y_{\{i, i+1\}}}{\partial\omega^y_{\{i, i+1\}}}\right| \left|\frac{\partial\omega^y_{\{i, i+1\}}}{\partial\mathbf{p}^y_{\{i+1, i+2\}}}\right| \left|\frac{\partial\mathbf{p}^x_{\{i+1, i+2\}}}{\partial\omega^x_{\{i, i+1\}}}\right| \left|\frac{\partial\omega^x_{\{i, i+1\}}}{\partial\mathbf{u}^x_{\{i, i+1\}}}\right|\\ &= \frac{p_{y_i}(\omega^y_i)G(\mathbf{p}^y_i \to \mathbf{p}^x_{i+1})}{p_{x_i}(\omega^x_i)G(\mathbf{p}^x_i \to \mathbf{p}^x_{i+1})}\frac{p_{y_{i+1}}(\omega^y_{i+1})}{p_{x_{i+1}}(\omega^x_{i+1})}\\ &= J_{\mathbf{x} \to \mathbf{y}} \end{aligned} \end{equation}

重连接后 $\mathbf{u}_{i-1}$ 为隐式, 需要BSDF抽样可逆即由 $\omega^x_{i-1}$ 恢复 $\mathbf{u}^x_{i-1}$ . 若以对称核 $k(\omega^x_{i-1}, \omega^y_{i-1})$ 扰动立体角, 变换到 $\mathbf{u}_{i-1}$ 后为 $k(\omega^x_{i-1}, \omega^y_{i-1})\left|\frac{\partial\omega^y_{i-1}}{\partial\mathbf{u}^y_{i-1}}\right| = \frac{k(\omega^x_{i-1}, \omega^y_{i-1})}{p_{y_{i-1}}(\omega^y_{i-1})}$ , 此时比值变为:

\begin{equation} \frac{T(\mathbf{y} \to \mathbf{x})}{T(\mathbf{x} \to \mathbf{y})} = \frac{p_{y_{i-1}}(\omega^y_{i-1})}{p_{x_{i-1}}(\omega^x_{i-1})}J_{\mathbf{x} \to \mathbf{y}} \end{equation}

$\mathbf{p}_{i+2} = \mathbf{p}_e$ 时, $\mathbf{p}_{i+1} \to \mathbf{p}_e$ 确定, $\mathbf{u}_{i+1}$ 不存在, 去除 $\mathbf{u}_{i+1}$ 的变换:

\begin{equation} \frac{T(\mathbf{y} \to \mathbf{x})}{T(\mathbf{x} \to \mathbf{y})} = \frac{p_{y_{i-1}}(\omega^y_{i-1})}{p_{x_{i-1}}(\omega^x_{i-1})} \frac{p_{y_i}(\omega^y_i)G(\mathbf{p}^y_i \to \mathbf{p}^x_{i+1})}{p_{x_i}(\omega^x_i)G(\mathbf{p}^x_i \to \mathbf{p}^x_{i+1})} \end{equation}

$\mathbf{p}_{i+1} = \mathbf{p}_e$ 时, $\mathbf{p}_i \to \mathbf{p}_e$ 确定, $\mathbf{u}_i$ 不存在, 提议只包含 $\mathbf{u}_{i-1}$ 的变换:

\begin{equation} \frac{T(\mathbf{y} \to \mathbf{x})}{T(\mathbf{x} \to \mathbf{y})} = \frac{p_{y_{i-1}}(\omega^y_{i-1})}{p_{x_{i-1}}(\omega^x_{i-1})} \end{equation}

Correlated Sampling
#

Control Variates
#

控制变量引入与 $f$ 相关的辅助函数 $h$ , 积分 $H = \int_{\Omega_x} h(\mathbf{x})\mathrm{d}\mathbf{x}$ 已知, 可改写估计量:

\begin{equation} \langle I_x \rangle = \alpha H + \left(f(X) - \alpha h(X)\right)W_X \end{equation}

图像空间控制变量以相邻像素 $y$ 为辅助函数, 要求 $I_y$ 可低方差估计, 估计量如下:

\begin{equation} \langle I_x \rangle_{\leftarrow y} = \alpha\langle I_y \rangle + \langle I_x - \alpha I_y \rangle \end{equation}

记 $A = f(X)W_X$ , $B = f(Y)W_Y$ , 分别无偏估计 $I_x$ 与 $I_y$ , $A - \alpha B$ 的方差为:

\begin{equation} \mathrm{Var}[A - \alpha B] = \mathrm{Var}[A] + \alpha^2\mathrm{Var}[B] - 2\alpha\mathrm{Cov}[A, B] \end{equation}

对 $\alpha$ 求导得最优系数, 令 $A$ , $B$ 的相关系数为 $\rho = \frac{\mathrm{Cov}[A,B]}{\sqrt{\mathrm{Var}[A]\mathrm{Var}[B]}}$ , 代回得方差缩减与 $\rho$ 相关:

\begin{equation} \alpha^* = \frac{\mathrm{Cov}[A, B]}{\mathrm{Var}[B]}, \quad \mathrm{Var}[A - \alpha^* B] = (1 - \rho^2)\mathrm{Var}[A] \end{equation}

若 $X$ 与 $Y$ 独立则 $\mathrm{Cov}[A, B] = 0$ , 任意 $\alpha \neq 0$ 均使方差增加. 对于 $\alpha\langle I_y \rangle + (A - \alpha B)$ , 若 $\langle I_y \rangle$ 为 $B$ 则退化为 $A$ ; 若为独立样本 $B'$ , 方差为 $\mathrm{Var}[A] + \alpha^2(\mathrm{Var}[B] + \mathrm{Var}[B'])$ . 两种情况都不低于 $\mathrm{Var}[A]$ , 因此独立样本无法缩减方差, $X$ 与 $Y$ 必须相关.

令 $Y = T_{x \to y}(X)$ , 场景连续处 $f(\mathbf{x}) \approx f(T_{x \to y}(\mathbf{x}))$ , 此时 $\rho \to 1$ .

\begin{equation} \begin{aligned} I_x - \alpha I_y &= \int_{\Omega_x} f(\mathbf{x})\mathrm{d}\mathbf{x} - \alpha\int_{\Omega_y} f(\mathbf{y})\mathrm{d}\mathbf{y}\\ &= \int_{\Omega_x} \left(f(\mathbf{x}) - \alpha f(T_{x \to y}(\mathbf{x}))J_{\mathbf{x} \to \mathbf{y}}\right)\mathrm{d}\mathbf{x} - \alpha\int_{\Omega_y \setminus T_{x \to y}(\Omega_x)} f(\mathbf{y})\mathrm{d}\mathbf{y} \end{aligned} \end{equation}

不保证 $T_{x \to y}(\Omega_x) \supseteq \Omega_y$ , 需要MIS:

\begin{equation} \begin{aligned} \langle I_x - \alpha I_y \rangle &= m_x(X)\left(f(X) - \alpha f(T_{x \to y}(X))J_{X \to Y}\right)W_X\\ &+\ m_y(Y)\left(f(T_{y \to x}(Y))J_{Y \to X} - \alpha f(Y)\right)W_Y \end{aligned} \end{equation}

ReSTCV
#

$x$ 重投影到 $y$ , $X$ 为当前帧新样本, $\langle I_x \rangle_\mathrm{init} = f(X)W_X$ . 时域估计为:

\begin{equation} \langle I_x \rangle = \frac{M_y\langle I_x \rangle_{\leftarrow y} + M_\mathrm{init}\langle I_x \rangle_\mathrm{init}}{M_y + M_\mathrm{init}} \end{equation}

$\mathcal{N}(x)$ 为包含 $x$ 自身的空域像素集合, 空域估计为:

\begin{equation} \langle I_x \rangle = \frac{\sum_{y \in \mathcal{N}(x)} M_y\langle I_x \rangle_{\leftarrow y}}{\sum_{y \in \mathcal{N}(x)} M_y} \end{equation}

ReSTIR PT下 $\langle I_x \rangle_{\leftarrow y}$ 的参数可从蓄水池获取, GRIS估计量只根据无偏权重调整样本亮度, 控制变量估计量包含多个通道, 可有效降低复杂色彩光照或光谱渲染的方差.

Histogram Stratification
#

设 $X \sim p$ , $p$ 覆盖 $f$ 的支撑集, 贡献 $C = f(X)/p(X) = f(X)W_X$ 可积, $C$ 的CDF为:

\begin{equation} F_C(c) = P(C \leq c) = \int_{\{\mathbf{x} \in \Omega_x:\,p(\mathbf{x}) > 0,\,f(\mathbf{x})/p(\mathbf{x}) \leq c\}} p(\mathbf{x})\mathrm{d}\mathbf{x} \end{equation}

定义 $C$ 的广义逆用于逆变换法:

\begin{equation} F_C^{-1}(u) = \inf\{c \in \mathbb{R}: F_C(c) \geq u\}, \quad u \in (0, 1) \end{equation}

令 $U \sim \mathcal{U}(0, 1)$ , 此时可得:

\begin{equation} I_x = \int_{\Omega_x} f(\mathbf{x})\mathrm{d}\mathbf{x} = E[C] = \int_0^1 F_C^{-1}(u)\mathrm{d}u = E[F_C^{-1}(U)] \end{equation}

令有限候选 $\mathcal{S} = \{(X_i, W_i)\}_{i=1}^N$ , $C_i = f(X_i)W_i$ , 定义抽样函数:

\begin{equation} G_{\mathcal{S}}(u) = C_{1 + \lfloor Nu \rfloor}, \quad u \in [0, 1) \end{equation}

取独立于候选的均匀随机数 $U$ , 每个候选对应长度 $1/N$ 的区间:

\begin{equation} \begin{aligned} E[G_{\mathcal{S}}(U) \mid \mathcal{S}] &= \int_0^1 G_{\mathcal{S}}(u)\mathrm{d}u = \frac{1}{N}\sum_{i=1}^N C_i\\ E[G_{\mathcal{S}}(U)] &= E\left[E[G_{\mathcal{S}}(U) \mid \mathcal{S}]\right] = \frac{1}{N}\sum_{i=1}^N E[C_i] = I_x \end{aligned} \end{equation}

按 $C_i$ 从小到大排序得置换 $\pi$ , 将列表重排后, 抽样函数成为经验逆CDF, 不影响无偏性:

\begin{equation} G_{\mathcal{S}}(u) = C_{\pi(1 + \lfloor Nu \rfloor)} \end{equation}

取独立于候选的 $T \sim \mathcal{U}[0, 1/4]$ , 构造四个相关样本, $U_k$ 在对应的 $\frac{1}{4}$ 区间内均匀分布:

\begin{equation} U_1 = T, \quad U_2 = \frac{1}{2} - T, \quad U_3 = \frac{1}{2} + T, \quad U_4 = 1 - T \end{equation}

给定候选 $\mathcal{S}$ , 记 $Z_k = G_{\mathcal{S}}(U_k)$ , 平均估计量为 $\langle I_x\rangle_\mathrm{anti} = \frac{1}{4}\sum_{k=1}^4 Z_k$ . 由四层等宽可得:

\begin{equation} \begin{aligned} &E[Z_k\mid\mathcal{S}] = 4\int_{(k-1)/4}^{k/4}G_{\mathcal{S}}(u)\mathrm{d}u\\ &E\left[\langle I_x\rangle_\mathrm{anti}\middle|\mathcal{S}\right] = E\left[\frac{1}{4}\sum_{k=1}^4 Z_k\middle|\mathcal{S}\right] = \frac{1}{4}\sum_{k=1}^4 E[Z_k\mid\mathcal{S}] = \frac{1}{N}\sum_{i=1}^N C_i \end{aligned} \end{equation}

条件方差展开为:

\begin{equation} \mathrm{Var}[\langle I_x\rangle_\mathrm{anti}\mid\mathcal{S}] = \frac{1}{16}\left(\sum_{k=1}^4\mathrm{Var}[Z_k\mid\mathcal{S}] + 2\sum_{i<j}\mathrm{Cov}[Z_i,Z_j\mid\mathcal{S}]\right) \end{equation}

令 $z_k(t) = G_{\mathcal{S}}(U_k(t))$ , $L = \frac{1}{4}$ , $T$ 在 $[0,L]$ 上均匀分布, 由协方差定义得:

\begin{equation} \begin{aligned} \mathrm{Cov}[Z_i,Z_j\mid\mathcal{S}] &= \frac{1}{L}\int_0^L z_i(t)z_j(t)\mathrm{d}t - \frac{1}{L^2}\left(\int_0^L z_i(t)\mathrm{d}t\right)\left(\int_0^L z_j(t)\mathrm{d}t\right)\\ &= \frac{1}{L^2}\int_0^L\int_0^L\left(z_i(t)z_j(t)-z_i(t)z_j(s)\right)\mathrm{d}s\mathrm{d}t\\ &= \frac{1}{2L^2}\int_0^L\int_0^L\left(z_i(t)z_j(t)+z_i(s)z_j(s)-z_i(t)z_j(s)-z_i(s)z_j(t)\right)\mathrm{d}s\mathrm{d}t\\ &= \frac{1}{2L^2}\int_0^L\int_0^L \left(z_i(t)-z_i(s)\right)\left(z_j(t)-z_j(s)\right)\mathrm{d}s\mathrm{d}t \end{aligned} \end{equation}

排序后的 $G_{\mathcal{S}}$ 使 $z_1,z_3$ 单调不减, $z_2,z_4$ 单调不增, 因此:

\begin{equation} \begin{aligned} \mathrm{Cov}[Z_i,Z_j\mid\mathcal{S}] &\leq 0, &&(i,j)\in\{(1,2),(1,4),(2,3),(3,4)\}\\ \mathrm{Cov}[Z_i,Z_j\mid\mathcal{S}] &\geq 0, &&(i,j)\in\{(1,3),(2,4)\} \end{aligned} \end{equation}

与独立均匀样本比较, 四样本平均的条件方差为:

\begin{equation} \begin{aligned} &\mathrm{Var}[G_{\mathcal{S}}(U)\mid\mathcal{S}] = \frac{1}{4}\sum_{k=1}^4\left(\mathrm{Var}[Z_k\mid\mathcal{S}]+\left(E[Z_k\mid\mathcal{S}]-\frac{1}{N}\sum_{i=1}^N C_i\right)^2\right)\\ &\mathrm{Var}[\langle I_x\rangle_\mathrm{iid}\mid\mathcal{S}] = \frac{1}{16}\sum_{k=1}^4\left(\mathrm{Var}[Z_k\mid\mathcal{S}]+\left(E[Z_k\mid\mathcal{S}]-\frac{1}{N}\sum_{i=1}^N C_i\right)^2\right)\\ &\mathrm{Var}[\langle I_x\rangle_\mathrm{anti}\mid\mathcal{S}] -\mathrm{Var}[\langle I_x\rangle_\mathrm{iid}\mid\mathcal{S}] = \frac{1}{16}\left(2\sum_{i<j}\mathrm{Cov}[Z_i,Z_j\mid\mathcal{S}]-\sum_{k=1}^4\left(E[Z_k\mid\mathcal{S}]-\frac{1}{N}\sum_{i=1}^N C_i\right)^2\right) \end{aligned} \end{equation}

以上均固定候选集合, 但 $\mathcal{S}$ 本身也需要随机抽样. 令 $\langle I_x\rangle$ 表示任意一种四样本平均估计量, 基于全方差公式可估计方差, 候选抽样方式直接影响第二项, 得到候选后再执行RIS:

\begin{equation} \begin{aligned} \mathrm{Var}[\langle I_x\rangle] &=\mathrm{Var}\left[\frac{1}{N}\sum_{i=1}^N C_i\right]+E\left[\mathrm{Var}[\langle I_x\rangle\mid\mathcal{S}]\right] \end{aligned} \end{equation}

Russian Roulette
#

Zero Variance
#

令 $c_i = \sum_{j > i}\beta_j L_e(\mathbf{p}_j)$ , 由 $L_i = L_r + L_e$ 可得:

\begin{equation} E[c_i|\mathbf{p}_{0 \dots i}] = E\big[\beta_{i+1}L_i(\mathbf{p}_{i+1} \to \mathbf{p}_i)\big|\mathbf{p}_{0 \dots i}\big]\\ = \beta_i L_r(\mathbf{p}_i \to \mathbf{p}_{i-1}) \end{equation}

零方差(zero-variance)抽样密度 $p_\mathrm{zv}$ 即前文条件分布, 抽样密度正比于方向贡献:

\begin{equation} p_\mathrm{zv}(\mathbf{p}_{i+1}|\mathbf{p}_{i-1}, \mathbf{p}_i) = \frac{f(\mathbf{p}_{i+1} \to \mathbf{p}_i \to \mathbf{p}_{i-1})G(\mathbf{p}_i \leftrightarrow \mathbf{p}_{i+1})L_i(\mathbf{p}_{i+1} \to \mathbf{p}_i)}{L_r(\mathbf{p}_i \to \mathbf{p}_{i-1})} \end{equation}

以 $p_\mathrm{zv}$ 抽样得 $\mathbf{p}_{i+1}$ , 代入上式形式如下. 光源不再弹射, 故到达光源顶点时贡献为像素真值.

\begin{equation} \begin{aligned} E[c_i|\mathbf{p}_{0 \dots i}] &= \beta_i L_r(\mathbf{p}_i \to \mathbf{p}_{i-1})\\ &= \beta_i\frac{f(\mathbf{p}_{i+1} \to \mathbf{p}_i \to \mathbf{p}_{i-1})G(\mathbf{p}_i \leftrightarrow \mathbf{p}_{i+1})}{p_\mathrm{zv}(\mathbf{p}_{i+1}|\mathbf{p}_{i-1}, \mathbf{p}_i)}L_i(\mathbf{p}_{i+1} \to \mathbf{p}_i)\\ &= \beta_{i+1}L_i(\mathbf{p}_{i+1} \to \mathbf{p}_i)\\ &= \beta_{i+1}L_e(\mathbf{p}_{i+1}) + E[c_{i+1}|\mathbf{p}_{0 \dots i+1}] \end{aligned} \end{equation}

由前文可得 $fGL_i = L_r p_\mathrm{zv}$ , 任意抽样密度 $p$ 下方差如下, 可见 $p=p_{zv}$ 时方差最小.

\begin{equation} \begin{aligned} \mathrm{Var}\left[\frac{fGL_i}{p}\right] &= \int_\mathcal{A}\left(\frac{fGL_i}{p} - L_r\right)^2 p\,\mathrm{d}\mathbf{p}_{i+1}\\ &= L_r^2\int_\mathcal{A}\left(\frac{p_\mathrm{zv}}{p} - 1\right)^2 p\,\mathrm{d}\mathbf{p}_{i+1}\\ &= L_r^2\int_\mathcal{A}\frac{(p_\mathrm{zv} - p)^2}{p}\mathrm{d}\mathbf{p}_{i+1} \end{aligned} \end{equation}

$\langle L_i \rangle$ 为后缀估计, 由方差定义得 $E[\langle L_i \rangle^2|\mathbf{p}_{i+1}] = L_i^2(\mathbf{p}_{i+1} \to \mathbf{p}_i) + \mathrm{Var}[\langle L_i \rangle|\mathbf{p}_{i+1}]$ . 分裂为 $n = \lfloor q \rfloor$ 条后缀, 再以概率 $q - n$ 的Bernoulli试验决定最后一次分裂, 每条权重乘 $\frac{1}{q}$ , $q < 1$ 时退化为轮盘赌. 记路径数 $N = n + B$ , $B \sim \mathrm{Bernoulli}(q - n)$ , $E[N] = q$ , $E[N^2] = q^2 + (q - n)(n + 1 - q)$ . 估计量 $\langle L_r \rangle_q = \frac{fG}{p}\frac{1}{q}\sum_{k=1}^N\langle L_i \rangle_k$ , 各后缀独立且与 $N$ 独立, 此时可得:

\begin{equation} \begin{aligned} \mathrm{Var}[\langle L_r \rangle_q] &= E[\langle L_r \rangle_q^2] - E[\langle L_r \rangle_q]^2\\ &= E\big[E[\langle L_r \rangle_q^2|\mathbf{p}_{i+1}]\big] - L_r^2\\ &= E\left[\left(\frac{fG}{p}\right)^2\frac{E\left[\left(\sum_{k=1}^N\langle L_i \rangle_k\right)^2\middle|\mathbf{p}_{i+1}\right]}{q^2}\right] - L_r^2\\ &= E\left[\left(\frac{fG}{p}\right)^2\frac{E\left[N\mathrm{Var}[\langle L_i \rangle|\mathbf{p}_{i+1}] + N^2L_i^2\middle|\mathbf{p}_{i+1}\right]}{q^2}\right] - L_r^2\\ &= E\left[\left(\frac{fG}{p}\right)^2\frac{q\mathrm{Var}[\langle L_i \rangle|\mathbf{p}_{i+1}] + \left(q^2 + (q - n)(n + 1 - q)\right)L_i^2}{q^2}\right] - L_r^2\\ &= \mathrm{Var}\left[\frac{fGL_i}{p}\right] + E\left[\left(\frac{fG}{p}\right)^2\left(\frac{\mathrm{Var}[\langle L_i \rangle|\mathbf{p}_{i+1}]}{q} + \frac{(q - n)(n + 1 - q)}{q^2}L_i^2\right)\right] \end{aligned} \end{equation}

记 $\mathbf{p}_i$ 处相对方差 $r_i = \frac{\mathrm{Var}[\langle L_r \rangle_q]}{L_r^2}$ , 设 $\mathbf{p}_{i+1}$ 非光源, 由 $L_i(\mathbf{p}_{i+1} \to \mathbf{p}_i) = L_r(\mathbf{p}_{i+1} \to \mathbf{p}_i)$ 得后缀相对方差 $r_{i+1} = \frac{\mathrm{Var}[\langle L_i \rangle|\mathbf{p}_{i+1}]}{L_i^2}$ . 记 $w = \frac{p_\mathrm{zv}}{p}=\frac{fGL_i}{pL_r}$ , 代入方差可得:

\begin{equation} \mathrm{Var}[\langle L_r \rangle_q] = L_r^2\left(\int_\mathcal{A}\frac{(p_\mathrm{zv} - p)^2}{p}\mathrm{d}\mathbf{p}_{i+1} + E\left[w^2\left(\frac{r_{i+1}}{q} + \frac{(q - n)(n + 1 - q)}{q^2}\right)\right]\right) \end{equation}

第一项与 $q$ 无关; 取整项分情况讨论, 均不超过 $\frac{1}{4}$ 与 $q$ , 因此第二项不超过 $E\left[\frac{w^2(r_{i+1} + 1)}{q}\right]$ .

\begin{equation} (q - n)(n + 1 - q) = \begin{cases} q(1 - q) \leq \min\left(\frac{1}{4}, q\right), & q < 1\\ (q - n)(n + 1 - q) \leq \frac{1}{4} \leq q, & q \geq 1 \end{cases} \end{equation}

预算 $E[q] = 1$ 下应用Cauchy-Schwarz不等式, 等号成立当且仅当 $\frac{w\sqrt{r_{i+1} + 1}}{\sqrt{q}} \propto \sqrt{q}$ 即 $w\sqrt{r_{i+1} + 1} = cq$ , 由 $E[q] = 1$ 得 $c = E\left[w\sqrt{r_{i+1} + 1}\right]$ , 此时上界最小:

\begin{equation} \begin{aligned} E\left[w\sqrt{r_{i+1} + 1}\right]^2 &= E\left[\frac{w\sqrt{r_{i+1} + 1}}{\sqrt{q}}\sqrt{q}\right]^2\\ &\leq E\left[\frac{w^2(r_{i+1} + 1)}{q}\right]E[q]\\ &= E\left[\frac{w^2(r_{i+1} + 1)}{q}\right] \end{aligned} \end{equation}

此时可得最优的 $q$ :

\begin{equation} q^* = \frac{w\sqrt{r_{i+1} + 1}}{E\left[w\sqrt{r_{i+1} + 1}\right]} \end{equation}

ADRR
#

假设 $r_{i+1}$ 为常数, 此时可以约去 $\sqrt{r_{i+1} + 1}$ , 由 $E[w] = \int_\mathcal{A} p_\mathrm{zv}\mathrm{d}\mathbf{p}_{i+1} = 1$ 可得 $q^* = w$ . 令 $q$ 为期望贡献与像素真值之比, 此时可得最优的 $q$ :

\begin{equation} q = \frac{\beta_{i+1}L_r(\mathbf{p}_{i+1} \to \mathbf{p}_i)}{I} = \frac{fGL_i(\mathbf{p}_{i+1} \to \mathbf{p}_i)}{pL_r(\mathbf{p}_i \to \mathbf{p}_{i-1})} = \frac{p_\mathrm{zv}}{p} = w \end{equation}

记 $\chi^2 = \int_\mathcal{A}\frac{(p_\mathrm{zv} - p)^2}{p}\mathrm{d}\mathbf{p}_{i+1}$ , 代入 $q = w$ 得:

\begin{equation} \begin{aligned} \mathrm{Var}[\langle L_r \rangle_q] &= \mathrm{Var}[L_r q] + E\left[L_r^2 q^2\left(\frac{r_{i+1}}{q} + \frac{(q - n)(n + 1 - q)}{q^2}\right)\right]\\ &= L_r^2\left(\int_\mathcal{A}\frac{(p_\mathrm{zv} - p)^2}{p}\mathrm{d}\mathbf{p}_{i+1} + E\big[(q - n)(n + 1 - q)\big] + \int_\mathcal{A} p_\mathrm{zv}\,r_{i+1}\,\mathrm{d}\mathbf{p}_{i+1}\right)\\ &= L_r^2\left(\chi^2 + r_{i+1} + E\big[(q - n)(n + 1 - q)\big]\right) \end{aligned} \end{equation}

$q \equiv 1$ 时形式如下:

\begin{equation} \begin{aligned} \mathrm{Var}[\langle L_r \rangle_1] &= \mathrm{Var}\left[\frac{fGL_i}{p}\right] + E\left[\left(\frac{fG}{p}\right)^2\mathrm{Var}[\langle L_i \rangle|\mathbf{p}_{i+1}]\right]\\ &= L_r^2\chi^2 + E\left[\left(\frac{fGL_i}{p}\right)^2 r_{i+1}\right]\\ &= L_r^2\left(\chi^2 + r_{i+1}E[w^2]\right)\\ &= L_r^2\left(\chi^2 + r_{i+1}\int_\mathcal{A}\frac{p_\mathrm{zv}^2}{p}\mathrm{d}\mathbf{p}_{i+1}\right)\\ &= L_r^2\left(\chi^2 + r_{i+1}\int_\mathcal{A}\frac{(p_\mathrm{zv} - p)^2 + 2p_\mathrm{zv}p - p^2}{p}\mathrm{d}\mathbf{p}_{i+1}\right)\\ &= L_r^2\left(\chi^2 + r_{i+1}(1 + \chi^2)\right) \end{aligned} \end{equation}

记 $\{w\}$ 为 $w$ 的小数部分, 任意 $w \geq 0$ 有 $\{w\}(1 - \{w\}) \leq |w - 1|$ , 此时可得:

\begin{equation} \begin{aligned} E\big[(q - n)(n + 1 - q)\big] &= E\big[\{w\}(1 - \{w\})\big]\\ &\leq E\big[|w - 1|\big]\\ &\leq \sqrt{E\big[(w - 1)^2\big]}\\ &= \chi \end{aligned} \end{equation}

因此 $r_{i+1}\chi^2 \geq \min\left(\frac{1}{4}, \chi\right)$ 时轮盘赌更优, 后缀噪声越大, $p$ 偏离 $p_\mathrm{zv}$ 越远, 收益越大.

Variable Rate Shading
#

DOOM: The Dark Ages
#

基于可见性缓冲分发每个材质的绘制调用, 材质序号作为深度比较是否相等. 材质无空间连贯性, Hi-Z无法剔除, 性能与Forward+持平.

DOOM将多种相似材质打包为单个着色器, 为每个着色器分配隐式命令缓冲. 计算分发基于图块判断块内像素需要调用的着色器和材质, 以子组为单位原子填充各着色器的命令缓冲, 即根据子组大小对齐某种材质的像素数量.

对于不支持隐式命令缓冲的平台, 对每种着色器根据每个图块所请求的工作组大小分桶, 增加一个维度用于标记图块序号. 由于桶的数量已知, CPU上为每个桶分别分发命令即可.

依据亮度梯度确定着色率, 填充命令缓冲时只统计主像素并附加复制位. 若判断为四元组只渲染0号样本, 再添加更严格的判断, 以避免体积雾等效果的抖动失效.

Call of Duty: Modern Warfare
#

基于MSAA实现软件VRS, 颜色附件宽高减半并开启4xMSAA. 深度预通道以像素频率着色, 插值中心为0号样本, 基于梯度还原其余样本的重心坐标, 执行透明度测试以提升植被渲染效果并修改覆盖掩码.

着色率表示为2位, 由于0号样本总是着色, 只需存储水平, 垂直或其余三个像素是否被着色. 每帧轮换0号像素, 使用覆盖掩码更新着色率, 历史4帧着色率存储为8位, 对其求或以利用时域统计的着色频率, 例如4帧都只着色0号样本, 当前帧保持即可.

Forward+绘制逐样本着色, 着色率通过修改模板缓冲剔除低着色率区域. 计算着色器缺少硬件子组打包, 以图块为单位将存活样本按平面紧凑打包以降低分歧.

Call of Duty: Modern Warfare 4
#

基于A-SVGF时间梯度调整采样率. 将屏幕以 $s \times s$ 分层, 每层选一个历史帧表面样本 $G_{i-1}$ , 前向投影/重放到当前帧得 $\vec{G}_{i-1}$ . 梯度与归一化因子定义如下:

\begin{equation} \begin{aligned} \delta &= f_i(\vec{G}_{i-1}, \xi_{i-1}) - f_{i-1}(G_{i-1}, \xi_{i-1})\\ \Delta &= \max\big(f_i(\vec{G}_{i-1}, \xi_{i-1}), f_{i-1}(G_{i-1}, \xi_{i-1})\big) \end{aligned} \end{equation}

滤波并归一化得历史权重 $\lambda$ , 调整全局 $\alpha$ 为混合系数 $\alpha_i$ , 取 $3 \times 3$ 邻域最大值弥补重建误差:

\begin{equation} \lambda = \min\left(1, \frac{|\hat{\delta}|}{\hat{\Delta}}\right), \quad \alpha_i = (1 - \lambda)\alpha + \lambda \end{equation}

$N$ 为历史样本数, $N_t = \frac{1}{\alpha_i}$ 为可容忍的历史样本数, $N_\mathrm{max} = \frac{1}{\alpha}$ 为目标样本数, 此时可计算被保留的历史样本数 $N'$ :

\begin{equation} N' = \min\left(N, N_t - 1\right) \end{equation}

历史未满时超采样, 历史已满时估计方差并降采样. $\rho$ 为依据方差调整的分段常数, 决定降采样追踪概率. 请求样本数为连续实数, 小数部分随机抽样光线, 样本数如下:

\begin{equation} n = \begin{cases} \min(N_\mathrm{max} - N', n_\mathrm{max}), & N' < N_\mathrm{max} - 1\\ \rho(\sigma^2), & N' = N_\mathrm{max} - 1 \end{cases} \end{equation}

对 $n$ 分桶统计全屏光线, 档位为 $\frac{1}{64}, \frac{1}{63}, \cdots, \frac{1}{2}, 1, 2, \cdots, 16$ , 选桶流程如下:

\begin{equation} r_{b(n)} = \begin{cases} \lfloor n \rfloor, & n \geq 1\\ 1/\lceil 1/n \rceil, & n < 1 \end{cases} \end{equation}

超出预算时自顶向下合并, 最高非空桶并入下一档直至不超过预算, 反之则自底向上合并. 末次合并时整档调整粒度过粗, 对桶内像素统一降低/提升样本数以刚好满足预算.

Path Guiding #

ReSTIR PG #

VXPG #

RCPG #

MCPG #

NASGPG #

NPMPG #

Radiance Cache #

ORCA #

SHARC #

NRC #

Markov Chain #

Metropolis-Hastings #

Metropolis Light Transport #

ReSTIR MCMC #

Correlated Sampling #

Control Variates #

ReSTCV #

Histogram Stratification #

Russian Roulette #

Zero Variance #

ADRR #

Variable Rate Shading #

DOOM: The Dark Ages #

Call of Duty: Modern Warfare #

Call of Duty: Modern Warfare 4 #