一个光源会在三维空间以及角度空间中随机产生光线分布。比方说,一个点光源的发出位置 origin 不变,而光线方向按照余弦来均匀分布。当进行渲染时,必须发出足够多个光线,才能精确描述这个点光源。
光线的随机产生通常使用随机数产生器,随机数产生器的目标是制造一系列互相无关的随机数,然后追踪大量的随机光线,在运算量较小的情况下获得尽可能高的渲染结果。
所有基于现代 CPU 的随机数生成算法都是伪随机的(quasi-random)。它们受限于一个周期。当超过周期后就会重复出现,而不再是相互无关的随机数。这个周期的最终限定是由电脑的位数来决定的,因此,没有一个内建的随机数是“真正”随机的。
Sobol 采样使用了不同的方式来采样。比起随机数,Sobol 序列着重于在概率空间中产生均匀的分布。但这并不是单纯的使用网格来填满,而是使用一个本质上随机,但是巧妙的方法去“填满”概率空间,即之后产生的随机数会分布到之前没有采样到的区域。
Radical Inversion 运算和 Van der Corput 序列
本节来自低差异序列(一)- 常见序列的定义及性质,欢迎读者移步查看原文。
Radical Inversion 运算的定义:
上式中,如果将任意一个整数 i 表示成
如果将上述过程简化一下,另 C 为单位矩阵,可以得到 Van der Corput 序列:
举个例子,正整数8以2为底数的radical inverse的计算过程如下。首先算出8的2进制表示,1000。此处假设 C 单位矩阵,所以直接将1000镜像到小数点右边,0.0001。这个二进制数的值就是最终结果,把它转换回10进制就得到1/16, 即
Van der Corput序列有几个属性:
- 每一个样本点都会落在当前已经有的点里“最没有被覆盖”的区域。例如 Φ2(3)=3/4 是刚好落在了 [0,1) 区间中被覆盖最少的区域( Φ2(0)=0,Φ2(1)=1/2,Φ2(2)=1/4 )。
- 样本个数达到 bm 个点时对 [0,1) 会形成 uniform 划分。
- 很多时候并不能够代替伪随机数,因为点的位置和索引有很强的关系,例如在以2为底的Van der Corput序列中,索引为基数时候序列的值大于等于0.5,偶数时则小于0.5。
产生 Sobol 序列
本节来自低差异序列(二)- 高效实现以及应用,由于原文写很棒,我也就不献丑了,直接贴过来。
Sobol序列的每一个维度都是由底数为 2 的 radical inversion 组成,但每一个维度的radical inversion都有各自不同的矩阵
C
。
因为完全以2为底数,所以Sobol序列的生成可以直接使用bit位操作实现radical inversion,非常高效。Sobol序列的分布具有不仅均匀,而且当样本的个数为2的整数次幂时,在
Radical Inversion 实现
目前介绍过的所有序列都基于Radical Inversion这个操作,而这个操作尽管公式略微复杂,但它的实现非常直观和简单,下面贴上代码以及注。
double IntegerRadicalInverse(int Base, int i)
{
int numPoints, inverse;
numPoints = 1;
// 此循环将i在"Base"进制下的数字左右Flip
for(inverse = 0; i > 0; i /= Base)
{
inverse = inverse * Base + (i % Base);
numPoints = numPoints * Base;
}
// 除以Digit将这个数镜像到小数点右边
return inverse / (double) numPoints;
}
Sobol序列的所有维度都是基于2为底数的Radical Inversion。只不过每个维度有自己不同的生成矩阵(Generator Matrix)
C
。因为是以2为底数,将数字从二进制中提取每一位的数字,以及和矩阵
double Sobol(uint64 i, uint Dimension)
{
double r;
// 将i依次右移,提取2进制里的每一位
for (uint k = 0; i; i >>= 1, k++)
if (i & 1) // 若当前位为1,则用异或和矩阵相乘
r ^= C[Dimension][k];
return r / (double) (1 << M); // 除以2^M,移到小数点右边
}
那么如何计算出能生成如此高质量分布的矩阵呢?Quasi Monte Carlo的学者们已经花了数10年的时间搜索这种矩阵,现在我们可以在这个网页(Sobol sequence generator )找到可以生成21201维度的Sobol序列的矩阵。