光线追踪渲染实战（五）：低差异序列与重要性采样，加速收敛！

项目代码仓库：
GitHub：https://github.com/AKGWSB/EzRT
gitee：https://gitee.com/AKGWSB/EzRT

前言

使用低差异序列与重要性采样来加速光线追踪的收敛！下图是仅 1000 spp 的结果，且找不到任何噪点。而要知道以往博客的示例图样都是 4000 spp 起步的：

来看看我们究竟干了什么事情：

1. 低差异序列降噪：

2. 重要性采样（diffuse，specular）降噪：

3. 对 hdr 贴图重要性采样（从左到右依次是原图，样本，结果）：

4. 多重（chong）重（zhong）要性采样：

---

在上一篇博客 微平面理论与迪士尼 BRDF，严格遵循物理！ 中，我们学习了了微平面理论，并且通过给出的公式实现了完整的迪士尼原则的 BRDF，实现了基于物理的渲染

但是现阶段我们是通过半球均匀采样来计算光线追踪的。尽管对于漫反射来说，均匀采样能够得到很好的结果，但是对于 BRDF，尤其是粗糙度非常低的时候，均匀采样很难达到收敛，因为我们浪费了很多采样：

如图，对于上图的 10 个采样，真正有用的只有 3 个，效率很低。在这篇博客中我们将通过一些数学方法来强化采样策略：

1. 从样本上：用空间上更加均匀的 低差异序列 代替伪随机序列
2. 从分布上：用概率论中的 重要性采样 策略来加速渲染方程的收敛过程

不多说了，开始把

1. 低差异序列介绍

渲染方程是一个困难积分：

通常情况下，我们使用蒙特卡洛采样，来对困难积分的值进行估计，而对于一个函数的估计，需要这么几个步骤：

1. 选取积分域内的任意样本 $x$
2. 计算选取该样本的概率 $pdf(x)$
3. 计算样本 $x$ 对应的函数值 $f(x)$
4. 本次采样的贡献为 $\frac{f(x)}{pdf(x)}$

这里有一个结论，就是如果我们使用的样本，在样本空间内的分布越 均匀 ，那么积分的估计值就会越准确。先入为主的想，原因是因为均匀的样本 不会放过被积函数的任何一个角落 ，所以估计的值会更加准确

以 $y=x$ 为例，选取 5 个样本，假设 $pdf(x)$ 恒为 $\frac{1}{1.0-0.0}$，可以看到均匀的选取 0.2，0.4，0.6，0.8，1.0 得到的估计值，要准于不均匀的选取 0.1，0.7，0.8，0.9，1.0 的采样：

选取均匀的样本，对于采样的效果至关重要。于是我们引入差异（discrepancy）来衡量一组样本 $X$ 的均匀程度，差异越低的样本越均匀。对于一个拥有 $N$ 个样本和 $V$ 单位体积的样本空间的的点集 $X$，我们任取一体积为 $v$ 的子空间，其中包含 $n$ 个点，那么有：

$$Disc(X)=\left|\frac{n/N}{v/V}-1.0\right|$$

比如一杯水 100 ml，有 3000 个水分子，你一勺子打了 10 ml，那么你理应获得 300 个水分子，这说明水的分布是绝对均匀的，没有差异的

更加直观的，下面是二维平面上，使用较为均匀的低差异序列，和普通的随机数序列生成的点。可以直观的看到，在样本数量同为 256 个的时候，低差异序列对于空间的填充更加均匀：

可以看到低差异序列均匀的覆盖了样本区域，而对于伪随机数来说，存在很多的孔洞

注：
蒙特卡洛积分是无偏的，但是积分域的某些角落，在采样次数较少的时候难以达到，采样次数上去了还是可以采到的
而对于低差异序列，非常少量的采样就可以均匀照顾到区域，因此低差异序列在采样次数低的时候能够更快收敛

2. sobol 序列及其实现

关于 sobol sequence 的原理我看了 3 遍相也没懂，毕竟是数学的东西，都是怪物搞出来的。。。好在最后是勉强搞懂了那些参数和符号的意义，能够拼凑出代码实现了

这里 是讲 sobol 原理的论文地址，而 这里 是一篇不错的笔记，对代码实现有帮助。最后 这里 则是计算生成矩阵的生成多项式参数

原文比较复杂，下面仅挑一些重要的稍微概括一下，然后是代码实现

2.1. 生成矩阵与递推式

sobol 序列是一个递推序列，序列的第 i 个数，取决于序列的第 i-1 个数。此外 sobol 序列需要一个生成矩阵 $v$，它是一个 二进制 矩阵，通常以 一维数组 的形式出现：

$$v=\left[\begin{array}{cccc}0& 1& 0& 0\\ 0& 0& 1& 1\\ 1& 0& 0& 1\\ 0& 1& 1& 0\end{array}\right]=\left[\begin{array}{c}4\\ 3\\ 9\\ 6\end{array}\right]=\left[\begin{array}{c}{\vec{v}}_1\\ {\vec{v}}_2\\ {\vec{v}}_3\\ {\vec{v}}_4\end{array}\right]$$

此外 sobol 的表达式是基于按位异或的操作（bit xor），其中 $x_{i,j}$ 表示第 i 个 sobol 数的第 j 个 bit，而 $i_j$ 表示整数 i 的第 j 个 bit，$v_{i,j}$ 则表示生成矩阵的第 i 行第 j 列。第 i 个 sobol 数的第 j 个 bit 的表达式如下：

$$x_{i,j}=i_1{v}_{1,j}\oplus i_2{v}_{2,j}\oplus i_3{v}_{3,j}\cdots$$

对于每个 $x_i$ 的每一 bit 都要和所有的 $v_i$ 按位异或显得有些麻烦，但是存在一种特殊情况：如果 i 和 i-1 的二进制表示 只有 1 bit 的差异，比如 $5=(101{)}_2$ 和 $7=(111{)}_2$，那么通过 $x_{i-1}$ 异或上那一个差异的 bit 就可以得到 $x_i$，相当于我们得到了递推式：

$$x_{i,j}=x_{i-1,j}\oplus v_{c_{i-1},j}$$

其中 $c_{i-1}$ 是 i 和 i-1 的 发生差异 的那一 bit。对于一个整数从 0 变成 i，需要改变 k 次，k 为 i 的二进制表示中 1 的个数，下面以 0 变成 11 （二进制 1011）为例：

1. 改变第 1 个bit，由 0000 变为 0001，$c_0=0$
2. 改变第 2 个bit，由 0000 变为 0011，$c_1=1$
3. 改变第 4 个bit，由 0000 变为 1011，$c_2=3$

于是对于递推式有：

$$x_{i,j}=x_{i-1,j}\oplus v_{c_{i-1},j}=(x_{i-2,j}\oplus v_{c_{i-2},j})\oplus v_{c_{i-1},j}=\cdots$$

此外规定 $x_{0,j}=0$，而 0 异或任何数字都不影响结果，于是最终变成了：

$$x_{i,j}=v_{c_{i-1},j}\oplus v_{c_{}}$$