采样与重建

渲染一张2D图片，实际上就是使用离散的像素点去描述真实世界的场景。在光线追踪算法中，需要通过光线去“采样”场景中的信息，然后”重建“在计算机的屏幕上。一般来说，“采样”指的是使用一些样本去获得某个连续函数的特征，而重建指的是通过采样获得的特征去还原原来的连续函数。

下面是一个采样一维函数并重建的示意图（PBRT）：

在图(a)中采样了9个样本并计算了它对应的样本值，然后在图(b)中使用这些样本值来尝试重建原来的函数，这个重建并不精确，因为样本信息不够完善。

采样与重建更深入的理论需要了解傅里叶变换（Fourier Transform），简单来说傅里叶变换可以把一个函数变换为多个三角函数来表示，这样就可以得到函数的频率信息，然后根据时域的卷积等于频域的乘积：( ，就可以根据频域去制定样本。
这里暂时不做深入的研究，简单了解一下傅里叶变换的公式：

$$F(\omega )={\int }_{-\infty }^{\infty }f(x)e^{-i2\pi \omega x}dx$$

它的逆变换为：

$$f(x)={\int }_{-\infty }^{\infty }F(\omega )e^{-2\pi \omega x}dx$$

一般来说，采样就是将脉冲函数与被采样函数$f(x)$相乘，脉冲函数的定义为：

$$H_T(x)=T\sum _{i=-\infty }^{\infty }\delta (x-iT)$$

对函数$f(x)$进行采样，可以表示为：

$$H_T(x)f(x)=T\sum _{i=-\infty }^{\infty }\delta (x-iT)f(iT)$$

使用图像表示为：

(a)为被采样函数$f(x)$，(b)为脉冲函数，(c)为采样结果。

而重建的过程可以表示为将采样结果与滤波函数$r(x)$求卷积，卷积的定义为：

$$f(x)\otimes g(x)={\int }_{-\infty }^{\infty }f(\tau )g(x-\tau )d\tau$$

其中$\tau$被称为约束变量。卷积运算的定义有点复杂，这里暂时记住有这么一个运算就行。
对函数$f(x)$重建的过程可以表示为：

$$f^{{}^{\prime }}(x)=H_T(x)f(x)\otimes r(x)$$

其中滤波函数$r(x)$有多种，常用的有盒型滤波和高斯滤波等。

采样与积分

在渲染中，通常使用蒙特卡洛积分来进行渲染方程的求解，在求解的过程中，一个好的采样分布直接影响到积分的收敛速度，所以在设计渲染算法时，一个可以生成优秀分布的样本点的采样算法是非常重要的。

Discrepancy

什么是一个优秀的样本分布，通常使用Discrepancy来衡量，它的定义是：

这个公式不需要做太深入的理解，简单来说就是在单位区域P点集中随机选择一个区域B，B中的采样点数量与总采样点数量的比值与B的体积相差得越小，那么这些样本的Discrepancy就越小，这些样本的分布就越优秀。通俗的理解就是优秀的样本分布应该是样本均匀地分布在整个样本空间。如果一个序列的Discrepancy很小，那么就称它是一个低差序列（Low Discrepancy Sequence）。
在程序中，采样时一般都是从一个伪随机数序列中选取样本，所以，一个好的采样算法应该是可以生成低差序列的。

采样器的设计（PBRT）

采样器（Sampler）用来生成样本点，在渲染中可以使用这些样本进行蒙特卡洛积分的计算，不同的生成算法将会影响到渲染的最终效果。
采样器需要根据不同的情况生成不同维度的采样点，因为在渲染方程中，不同阶段的光线可能需要不同维度的样本。

采样器的基类

定义一个采样器的基类：

class Smapler
{
...
};

首先，它需要提供每个像素的采样数，这个变量可以作为成员变量保存起来，一般在渲染系统初始化的时候确定。

int SamplesPerPixel() const;

然后，采样器需要有一个StartPixelSample方法，在最开始生成camera ray时由Integrator调用，Integrator需要提供采样像素的坐标和采样的索引，采样的索引标识当前像素正在进行第几个样本的采样，它应该大于0但是小于等于每个像素的采样数（SamplesPerPixel的返回值），Integrator也可以指定从第几个样本进行采样，

void StartPixelSample(Point2i p, int sampleIndex, int dimension = 0);

这样做的目的有两个，第一，让采样器知道此次采样的索引，可以利用算法对样本进行优化，防止临近的两次采样取得相近的样本值；第二，可以在最开始进行初始化时就确定采样器的状态，采样器将为固定的采样索引提供相同的样本值，这样方便程序的debug并且不影响渲染结果。

接下来，采样器需要提供两个接口来返回采样获得的一维和二维样本，更高维的样本在渲染算法中一般不需要，需要时也可以将这两个低维的样本组合，

Float Get1D();
Point2f Get2D();

最后是一个可选的接口GetPixel2D，它专门用来生成flim平面上的坐标点，有的子类会实现它，有的子类只在其中调用Get2D方法，

Point2f GetPixel2D();

Independent Sampler

这是一个最简单的采样器，它生成一些独立均匀的样本，不考虑任何优化采样的措施，在渲染时一般只用于测试。

class IndependentSampler public Smapler
{
...
};

IndependentSampler中有一个伪随机的随机数生成器，设置不同的seed它就会生成不同的随机数，

class IndependentSampler public Smapler
{
public:
    IndependentSampler(int samplesPerPixel, int seed = 0) : samplesPerPixel(samplesPerPixel), seed(seed) {}  // 构造函数
    ...
private:
    int samplesPerPixel, seed;
    RNG rng; // 随机数生成器
};

其中RNG的状态要在每次使用完毕后重置，这样就可以保证每一次生成都是独立的。

接下来实现StartPixelSample方法，制定一个像素点的坐标和该像素的采样索引，“激活”这个采样器，

void StartPixelSample(Point2i p, int sampleIndex, int dimension) {
    rng.SetSequence(Hash(p, seed));
    rng.Advance(sampleIndex * 65536ull + dimension);
}

RNG使用SetSequence生成一个随机数序列，使用Advance指定这个序列的偏移，这样就可以保证同一个像素点的同一个采样索引生成的样本值是相同的，同时sampleIndex可以保证不同的样本值在序列中不会位于相近的地方。

接下来获取样本值的方法就比较简单了，直接使用rng生成随机数即可，

Float Get1D() { return rng.Uniform<Float>(); }
Point2f Get2D() { return {rng.Uniform<Float>(), rng.Uniform<Float>()}; }
Point2f GetPixel2D() { return Get2D(); }

这个采样器采样的效果非常不好，独立均匀的样本往往会在图片中产生很多噪点，但它已经是一个可以正常工作的采样器了，它的实现也非常清晰可观。

Halton Sampler

Halton采样器可以通过生成低差序列来得到一个更好的样本分布，它的采样效率会比Independent更高。

Radical Inversion

这是构造Halton采样器的基础思想，它的定义为：

$$\begin{array}{r}i=\sum _{l=0}^{M-1}{a}_l(i)b^l\\ {\Phi }_{b,C}(i)=(b^{-1}...b^{-M})[C(a_0(i)...a_{M-1}(i){)}^T]\end{array}$$

它的含义是，对于一个正整数b，将任何一个整数i表示为b进制的数，然后把每一位上的数字组成一个向量，和一个矩阵C相乘得到一个新的向量，再把这个新向量中的数字镜像到小数点右边。这样得到的小数称为以b为底数，以C为生成矩阵的radical Inversion${\Phi }_{b,C}(i)$。

Van der Corput序列

为了方便，先假设C是一个单位矩阵，这样乘以矩阵的过程就可以先忽略，整个过程就变成了把b进制的各个位上的数镜像到小数点右边。这就是Van der Corput序列。
举个例子，假设i = 8，b = 2，将8以二进制表示，得到1000，然后镜像小数点到右边，得到0.0001，转换为10进制就是0.0625，所以${\Phi }_2(8)=0.0625$。下表为以而为底的Van der Corput序列例子。

i	Base2	${\Phi }_2(i)$
0	0	0
1	1	0.1 = 1/2
2	10	0.01 = 1/4
3	11	0.11 = 3/4
4	100	0.001 = 1/ 8
…	…	…
Van der Corput有以下属性：

1. 每一个新的样本点都会落在当前样本点“最稀疏的区域”。
2. 样本点个数达到$b^m$个时会对$[0,1)$形成均匀的划分，比如上面样本点个数为4时，刚好4等分单位区间。
3. 很多时候并不能够替代伪随机数，因为样本值与索引有很大关系。

Halton序列和Hammersley点集

Halton序列定义如下：

$$X_i:=({\Phi }_{b1}(i),{\Phi }_{b2}(i),...,{\Phi }_{bn}(i))$$

确定样本点索引i，n维的Halton序列的每一个维度都是基于不同的底数$bn$的Van der Corput序列，通常来说，所使用的底数互质，比如$b1=2,b2=3,...,b5=11,...$。

Hammersley点集定义如下：

$$X_i:=(\frac{i}{N},{\Phi }_{b1}(i),{\Phi }_{b2}(i),...,{\Phi }_{bn-1}(i))$$

与Halton序列唯一不同的是将第一位的序列改成了$\frac{i}{N}$，N为样本的个数。

下图为Halton序列和Hammersley点集分别生成216个2维样本点的图（选择的radical Inversion低为2和3）：

图(a)为Halton序列，$({\Phi }_2(i),{\Phi }_3(i)){)}_{i=0}^{215}$
图(b)为Hammersley点集，$(\frac{i}{216},{\Phi }_2(i){)}_{i=0}^{215}$
从上图中可以看出Hammersley点集比Halton序列的样本分布得更加均匀，但是Hammersley点集生成必须预先知道N的大小（即总样本个数），而Halton序列没有这个限制。

Radical Inversion算法的实现（PBRT）

Float RadicalInverse(int baseIndex, uint64_t a) 
{ 
    // base通过索引从质数列表中取
    int base = Primes[baseIndex]; 
    Float invBase = (Float)1 / (Float)base, invBaseM = 1; 
    uint64_t reversedDigits = 0; 
    while (a) {
        // 将a看作base进制的数的话，digit是最后一位的数字，next是去掉最后一位剩下的数字
        // 即digit = a % base;
        uint64_t next = a / base; 
        uint64_t digit = a - next * base; 
        // 将数字累加到reversedDigits
        reversedDigits = reversedDigits * base + digit; 
        // 每取一位，invBaseM就会缩小Base倍，用来将数字镜像到小数点后边
        // 每乘以一次invBase，最终结果的小数点就会向左移动一位
        invBaseM *= invBase; 
        a = next;
    } 
    // 保证结果不会大于1
    return std::min(reversedDigits * invBaseM, OneMinusEpsilon); 
}

Randomization via Scrambling（通过扰乱进行随机化）

Halton序列和Hammersley点集固然是非常优秀的低差序列，但是它们有一个缺点，当使用的底数b增大时，样本的分布并不是特别均匀，事实上，只有点的数量接近于b的幂时样本分布才是最均匀的。（上面以2和3为底数时生成了216个样本，因为$2^3\ast 3^3=216$。)
如果底数b太大，而生成的样本点数不多，那么它们随机均匀分布的性质就会大打折扣。为了解决这个问题，需要使用扰乱（Scrambling）的方法来保证样本的随机性。扰乱的方法有许多中，对于Halton序列和Hammersley点集来说，一般是在将每一位数字镜像到小数点右边时，通过一组排列（permutations）先将其转换为另一个数字，这种转换并不会影响到样本的随机性，因为Halton序列是将一个域进行均匀的划分，而这种扰乱只是影响划分的顺序。
permutations的生成方法也很多，构造好的permutations也可以减少积分的错误，为了简单起见，先只考虑随机生成的permutations。
下面是PBRT中的算法实现：
对于一个底数base，DigitPermutation类用于生成和保存它的permutations，

class DigitPermutation
{
public:
    DigitPermutation(int base, uint32_t seed, Allocator alloc)
    : base(base)
    {
        // 先计算出此base的Permutation长度
        nDigits = 0;
        Float invBase = (Float)1 / (Float)base, invBaseM = 1;
        while (1 - (base - 1) * invBaseM < 1) {
             ++nDigits;
            invBaseM *= invBase;
        } 
    }
private:
    int base, nDigits; 
    uint16_t *permutations;
};

Sobol Samplers

sobol

Reference

https://zhuanlan.zhihu.com/p/20197323