Please enable Javascript to view the contents

采样与低差异序列

 ·  ☕ 8 分钟

最近在开始将PBRT中一些算法搬到Ilum渲染器上，边读边写边学边练

1. 采样理论

采样理论是大二信号与系统的内容了，考虑一个一维函数$f(x)$，假设我们可以采样其定义域内任意一点$x'$的值$f(x')$，希望通过一系列的采样点和采样值，重建出$f(x)$的估计$\tilde f(x)$，如下图所示

image-20220324213403574

讲到采样，就不得不提频域和傅里叶变换，傅里叶变换通过表达式：
$$
F(\omega)=\int_{-\infty}^{+\infty}f(x)e^{-i2\pi\omega x}\mathrm dx
$$
将时域信号$f(x)$变换到频域信号$F(\omega)$，同时，频域信号也可以通过逆变换：
$$
f(x)=\int_{-\infty}^{+\infty}F(\omega)e^{i2\pi\omega x}\mathrm d\omega
$$
变换到时域。

考虑冲激信号函数$\delta(x)$，满足$\int\delta (x)\mathrm dx=0$，且$\forall x\neq 0$，有$\delta (x)=0$，则有
$$
\int f(x)\delta (x)\mathrm dx=f(0)
$$
对于以采样率$T$的均匀连续采样，shah采样函数$III_T(x)$定义为：
$$
III_T(x)=T\sum_{i=-\infty}^{\infty}\delta(x-iT)
$$
将shah函数与待采样函数$f(x)$相乘可得：
$$
III_T(x)f(x)=T\sum_i \delta(x-iT)f(iT)
$$
shah函数相当于对$f(x)$以一定的采样率进行采样，如下图所示

image-20220324215323827

通过与一个滤波函数$r(x)$进行卷积来重建函数$\tilde f$：
$$
(III_T(x)f(x))\otimes r(x)
$$
卷积$f(x)\otimes g(x)$定义为：
$$
f(x)\otimes g(x)=\int_{-\infty}^{+\infty}f(x')g(x-x')\mathrm dx'
$$
因此，
$$
\tilde f(x)=T\sum_{i=-\infty}^{\infty}f(iT)r(x-iT)
$$
从上式找到一个合适的滤波函数$r(x)$具有一定的难度，考虑傅里叶变换的卷积性质：
$$
\begin{aligned}
\mathcal F\{f(x)\otimes g(x)\}&=F(\omega)G(\omega)\\\\
\mathcal F\{f(x)g(x)\}&=F(\omega)\otimes G(\omega)
\end{aligned}
$$
因此
$$
\begin{aligned}
\mathcal F\{III_T(x)f(x)\}&=\int_{-\infty}^{+\infty}F(\omega')III_{1/T}(\omega-\omega')\mathrm dx'\\\\
&=\int_{-\infty}^{+\infty}F(\omega')\frac{1}{T}\sum_{i=-\infty}^{\infty}\delta (\omega-\frac{i}{T}-\omega')\mathrm dx'\\\\
&=\frac{1}{T}\sum_{i=-\infty}^\infty F(\omega-\frac{i}{T})
\end{aligned}
$$
可以看到，$F(\omega)$与shah频域函数的卷积为$F(\omega)$的一个周期形式，如下图所示：

image-20220324220658526

因此我们的滤波函数只需能够截选其中一段即可，这里选用盒函数$\Pi_T(x)$：
$$
\Pi_T(x)=\begin{cases}
1/(2T)&|x|<T\\\\
0&\mathrm{otherwise}
\end{cases}
$$
因此重建函数的频域形式为：
$$
\tilde F(\omega)=(F(\omega)\otimes III_{1/T}(\omega))\Pi_T(\omega)
$$
逆变换得到：
$$
\tilde f(x)=(f(x)III_T(x))\otimes \mathrm{sinc}(x)=\sum_{-\infty}^{\infty}\mathrm{sinc}(x-i)f(i)
$$
当采样率$T$太低时，由奈奎斯特采样定理，容易发生走样现象：

image-20220324221305235

常见的反走样方法有：

  • 非均匀采样，例如采用采样函数：
    $$
    \sum_{i=-\infty}^{\infty}\delta\left(x-\left(i+\frac{1}{2}-\zeta\right)T\right)
    $$
    其中$\zeta$为$[0,1]$的随机数

  • 自适应采样

  • 预滤波

2. 随机采样序列

2.1. 随机采样模式的差异性

在基于光线追踪的渲染方法中，本质也是对材质、对光源进行采样，来重建出场景的光照结果，因此采样时使用的随机序列对采样结果至关重要。

一个随机采样模式的好坏可以通过差异性进行度量，设有$n$维空间$[0,1)^n$，考虑其子集合$B$满足：
$$
B=\{[0, v_1]\times[0, v_2]\times\cdots\times[0, v_n]\}
$$
其中$0\leq v_i<1$，对给定采样点序列$P=x_1,\cdots,x_N$，$P$对于$B$的差异性定义为：
$$
D_N(B,P)=\sup_{b\in B}\left|\frac{\#\{x_i\in b\}}{N}-V(b) \right|
$$
其中$\#\{x_i\in b\}$为$b$中点的数量，$V(b)$为$b$的体积

上式能够用于度量随机采样模式好坏的一大原因便是$\#\{x_i\in b\}/N$的值可以作为给定点集$P$下$b$的体积的近似，因此差异性表示为体积近似的最大误差。若$B$的集合覆盖了空间的边界时，差异性表示为$D^\ast_N(P)$

通常来说，分布越均匀的点集，任意区域内的点集数量占点总数量的比例也会越接近于这个区域的体积，因此低差异序列往往是我们在渲染中需要用到的随机采样序列。

生成随机序列的方法很多，这里参考了https://github.com/nvpro-samples/vk_raytrace/blob/master/shaders/random.glsl中的pcg随机数实现，可以实现$[0,1)$范围内的随机数：

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72

// https://dl.acm.org/doi/10.5555/1921479.1921500 
/*GPU Random Numbers via the Tiny Encryption Algorithm*/
uint tea(in uint val0, in uint val1)
{
    uint v0 = val0;
    uint v1 = val1;
    uint s0 = 0;

    for (uint n = 0; n < 16; n++)
    {
        s0 += 0x9e3779b9;
        v0 += ((v1 << 4) + 0xa341316c) ^ (v1 + s0) ^ ((v1 >> 5) + 0xc8013ea4);
        v1 += ((v0 << 4) + 0xad90777d) ^ (v0 + s0) ^ ((v0 >> 5) + 0x7e95761e);
    }

    return v0;
}

uint init_random(in uvec2 resolution, in uvec2 screen_coord, in uint frame)
{
    return Tea(screen_coord.y * resolution.x + screen_coord.x, frame);
}

// https://www.pcg-random.org/
uint pcg(inout uint state)
{
    uint prev = state * 747796405u + 2891336453u;
    uint word = ((prev >> ((prev >> 28u) + 4u)) ^ prev) * 277803737u;
    state     = prev;
    return (word >> 22u) ^ word;
}

uvec2 pcg2d(in uvec2 v)
{
    v = v * 1664525u + 1013904223u;
    v.x += v.y * 1664525u;
    v.y += v.x * 1664525u;
    v = v ^ (v >> 16u);
    v.x += v.y * 1664525u;
    v.y += v.x * 1664525u;
    v = v ^ (v >> 16u);
    return v;
}

uvec3 pcg3d(in uvec3 v)
{
    v = v * 1664525u + uvec3(1013904223u);
    v.x += v.y * v.z;
    v.y += v.z * v.x;
    v.z += v.x * v.y;
    v ^= v >> uvec3(16u);
    v.x += v.y * v.z;
    v.y += v.z * v.x;
    v.z += v.x * v.y;
    return v;
}

float rand(inout uint seed)
{
    uint r = pcg(seed);
    return uintBitsToFloat(0x3f800000 | (r >> 9)) - 1.0f;
}

vec2 rand2(inout uint prev)
{
    return vec2(rand(prev), rand(prev));
}

vec3 rand3(inout uint prev)
{
    return vec3(rand(prev), rand(prev), rand(prev));
}

生成结果如下图所示：

uniform

2.2. Stratified采样

Stratified采样的核心思想是将采样域分为不重叠区域，然后对每一个区域进行采样。这样能够使得生成的随机点覆盖更广、更均匀，如下图所示：

image-20220325102508236

  • 随机模式(a)是一种无效的的采样模式，因为大量区域没有被采样到，图形采样不足
  • 统一分层模式(b)分布较好，但会加剧混叠效应
  • 分层抖动模式(c)在保持分层优势的同时，将均匀模式的混叠转化为高频噪声

生成结果如下图所示：

统一分层模式

stratified

分层抖动模式

stratified_jitter

2.3. Latin Hypercube采样

Latin Hypercube采样（LHS）是分层采样的一种优化，用于避免直接分层采样可能带来的样本在某一维上发生聚集到局部的情况，如下图所示：

image-20220325153523228

对于$n$个采样点，LHS方法过程如下：

  1. 先将采样区域$[0,1)^2$分成$n\times n$部分
  2. 采样对角线部分区域
  3. 沿y$方向随机排序
  4. 沿$x$方向随机排序

相当于对每一维上进行分层采样，如下图所示：

image-20220325153930056

2.4. Halton采样

Halton采样序列是一种常用的低差异序列，在介绍Halton采样以及其他低差异序列之前需要了解一个基本运算Radical Inversion，Radical Inversion先将一个正整数$a$表达为以$b$的幂为基底的加权组合形式：
$$
a=\sum_{i=1}^m d_i(a)b^{i-1}
$$
则Radical Inversion函数$\Phi_b$表示为：
$$
\Phi_{b,C}(a)=(b^{-1}\cdots b^{-m})\left[C(d_0(a)\cdots d_{m-1}(a))^T\right]
$$
直观上来说，$a$为一个正整数，先将$a$表示为$b$进制的数，再把每一位上的数$d_i(a)$排成一个向量，和一个生成矩阵$C$相乘得到一个新的向量，最后将这个新向量镜像到小数点右边就能得到$a$以$b$为底数，以$C$为生成矩阵的Radical Inversion函数$\Phi_b$

当生成矩阵$C$为单位阵时，将退化为Van der Corput序列，即：
$$
\Phi_{b}(a)=(b^{-1}\cdots b^{-m})(d_0(a)\cdots d_{m-1}(a))^T
$$
此时则是直接把$b$进制的数镜像到小数点右边即可。

以$b=2$的Van der Corput序列$\Phi_2$为例，不同的$a$取值结果如下：

image-20220325132055685

可以看到底为2的Van der Corput序列递归地将一维空间进行等距分割，其差异度为：
$$
D^\ast_N(P)=O\left(\frac{\log N}{N}\right)
$$
代码实现：

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15

float radical_inverse(uint base, uint i)
{
    float digit, radical, inverse;
    digit = radical = 1.0 / float(base);
    inverse         = 0.0;

    while (i != 0)
    {
        inverse += digit * float(i % base);
        digit *= radical;
        i /= base;
    }

    return inverse;
}

此时我们可以定义Halton序列：
$$
x_a=(\Phi_2(a), \Phi_3(a), \Phi_5(a),\cdots,\Phi_{p_n}(a))
$$
其中，$p_1,p_2,\cdots$互为质数，Halton序列的差异度为：
$$
D^\ast_N(x_a)=O\left(\frac{(\log N)^n}{N}\right)
$$
达到渐近最优。

当采样点数量$N$固定时，可以使用Hammersley序列来进一步降低差异度，Hammersley序列定义为：
$$
x_a=(\frac{a}{N}, \Phi_{p_1}(a), \Phi_{p_2}(a), \Phi_{p_3}(a),\cdots,\Phi_{p_n}(a))
$$
代码实现中，我们可以预存一个质数数组Primes，在计算Halton序列时只需简单地查询和计算Van der Corput序列即可：

1
2
3
4

float halton(uint dim, uint index)
{
    return radical_inverse(Primes[dim % 1000], index);
}

一些采样结果如下图所示：

$(\Phi_2(i),\Phi_3(i))_{i=0}^{3000}$ $(\Phi_{15}(i),\Phi_{16}(i))_{i=0}^{3000}$ $(\Phi_{29}(i),\Phi_{30}(i))_{i=0}^{3000}$
halton2_3 halton15_16 halton29_30

可以看到当Halton序列的底比较大时，当序列的点的数量不够多时，点会趋于落在一条直线上分布，以29和31为底的序列为例，一开始的点为$1/29,2/29,3/29\cdots$，落在一条直线上。通常的解决方法是进行Scrambling，一个常用的方法是Faure Scrambling：
$$
\Phi_b(a)=\sum_{i=0}^{m-1}\sigma_b(d_i(a))b^{-i-1}
$$
在Radical Inverse的时候不直接将数字镜像到小数点右边，而是先通过一个permutation操作$\sigma_b$先转换为另一个数字，不同底数$b$下有不同的permutation，例如$\sigma_4=[0,2,1,3]$、Scrambling本质上在交换均等划分的部分，因此处理后的序列依然具有Radical Inverse的性质。

在实际实现尤其是在着色器实现中，按公式计算需要有大量的二元运算和循环，效率低下，Leonhard的网页https://gruenschloss.org/给出了高效的求解方法https://gruenschloss.org/halton/halton.zip，通过查表可以得到Scrambling结果，例如$\Phi_5(a)$的实现如下：

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14

const uint FaurePermutation5[125] = {
    0, 75, 50, 25, 100, 15, 90, 65, 40, 115, 10, 85, 60, 35, 110, 5, 80, 55,
        30, 105, 20, 95, 70, 45, 120, 3, 78, 53, 28, 103, 18, 93, 68, 43, 118, 13, 88, 63, 38, 113, 8, 83, 58, 33, 108,
        23, 98, 73, 48, 123, 2, 77, 52, 27, 102, 17, 92, 67, 42, 117, 12, 87, 62, 37, 112, 7, 82, 57, 32, 107, 22, 97,
        72, 47, 122, 1, 76, 51, 26, 101, 16, 91, 66, 41, 116, 11, 86, 61, 36, 111, 6, 81, 56, 31, 106, 21, 96, 71, 46,
        121, 4, 79, 54, 29, 104, 19, 94, 69, 44, 119, 14, 89, 64, 39, 114, 9, 84, 59, 34, 109, 24, 99, 74, 49, 124};

float halton5(uint index)
{
    return (FaurePermutation5[index % 125u] * 1953125u + FaurePermutation5[(index / 125u) % 125u] * 15625u +
            FaurePermutation5[(index / 15625u) % 125u] * 125u +
            FaurePermutation5[(index / 1953125u) % 125u]) *
        4.096e-09;
}

而$\Phi_2(a)$可以直接通过位运算求得：

 1
 2
 3
 4
 5
 6
 7
 8
 9
10

float halton2(uint index)
{
    index = (index << 16) | (index >> 16);
    index = ((index & 0x00ff00ff) << 8) | ((index & 0xff00ff00) >> 8);
    index = ((index & 0x0f0f0f0f) << 4) | ((index & 0xf0f0f0f0) >> 4);
    index = ((index & 0x33333333) << 2) | ((index & 0xcccccccc) >> 2);
    index = ((index & 0x55555555) << 1) | ((index & 0xaaaaaaaa) >> 1);
    index = 0x3f800000u | (index >> 9);
    return uintBitsToFloat(index) - 1.f;
}

$(\Phi_2(a),\Phi_5(a))$序列如下图所示：

halton

2.5. Sobel采样

Sobel序列的每一维均由底数为2的Radical Inversion组成，但每一维的Radical Inversion均有不同的生成矩阵$C$，即：
$$
x_a=(\Phi_{2,C_1}(a), \Phi_{2,C_2}(a), \Phi_{2,C_3}(a),\cdots,\Phi_{2,C_n}(a))
$$
Sobel序列不需要像Halton序列一样进行Scrambling，高质量分布的矩阵$C$可以在网站https://web.maths.unsw.edu.au/~fkuo/sobol/找到，因为矩阵$C$实在太大了，因此也没有在Shader里实现了。

其他的采样序列还有0-2序列、最大最小距离序列以后有空再补

分享
Wenbo Chen
作者
Wenbo Chen
CG Student

相关内容

目录