Sampling and Reconstruction(采样与重建)
Sampling Theory
数字图像使用一组图像像素(image pixels)进行表示,通常排列在矩形网格上,关于数字图像的函数可以是连续的(渲染器渲染的场景),也可以是离散的(jpg,png 图像);而显示器使用一系列显示像素(display pixels)值,调整发射的光谱功率,从而在显示表面上构造出显示图像。当显示器呈现数字图像时,需要使用图像像素在显示表面上构造一个新的图像函数时,该函数在显示器上的所有点上定义,而不是数字图像像素上的点上定义,因此我们需要从数字图像的连续函数(仅考虑连续函数,因为离散数字图像也可以重构成连续函数)上取一组样本值并将其转换回连续函数,这个过程称为重构(reconstruction)。
为了获得最佳质量的结果,重构的图像要尽可能接近场景的原始图像,但整个采样和重建的过程涉及近似,会引入误差,这种误差被称为 aliasing(走样),例如锯齿状边缘或动画中的闪烁。这些错误的发生是因为采样过程不能从连续定义的图像函数中捕获所有信息。
傅里叶分析(Fourier analysis)可以用来评价重建函数与原始函数之间的匹配质量。
Fourier Transform(傅里叶变换)
傅里叶分析中一个重要的概念是傅里叶变换(Fourier transform),它表示一个频域(frequency domain)中的函数(其他函数一般被定义在时域 spatial domain 上)。
如下图所示,值变化较慢的函数(a)相比于值变化较快的函数(b)被认为具有较低频率(lower-frequency)的内容。
下图是这两个函数的频域空间表示,低频函数趋向于 0 的速度比高频函数快。
约瑟夫·傅里叶(Joseph Fourier)提出,大多数函数可以分解成移位正弦曲线的加权和。使用傅里叶变换可以将函数转换成这种形态。函数的这种性质可以让我们深入了解它的一些特性——函数对应的正弦函数中的频率分布对应于原始函数中的频率分布。使用这种形式,可以使用傅立叶分析来深入了解由采样和重建过程引入的误差,以及如何减少这种误差的感知影响。
一维函数 f ( x ) f(x) f(x)的傅里叶变换为:
F ( ω ) = ∫ − ∞ ∞ f ( x ) e − i 2 π ω x d x = ∫ − ∞ ∞ f ( x ) cos ( 2 π ω x ) d x − i ∫ − ∞ ∞ f ( x ) sin ( 2 π ω x ) d x \begin{equation} F(\omega)=\int_{-\infin}^\infin{}f(x)\mathrm{e}^{-\mathrm{i}2\pi\omega{}x}\,\mathrm{d}x=\int_{-\infin}^\infin{}f(x)\cos(2\pi\omega{x})\,\mathrm{d}x-\mathrm{i}\int_{-\infin}^\infin{}f(x)\sin(2\pi\omega{x})\,\mathrm{d}x \end{equation} F(ω)=∫−∞∞f(x)e−i2πωxdx=∫−∞∞f(x)cos(2πωx)dx−i∫−∞∞f(x)sin(2πωx)dx
其中,后一个形式是由 e i x = cos x + i sin x \mathrm{e}^{\mathrm{i}x}=\cos{}x+\mathrm{i}\sin{}x eix=cosx+isinx得到的,余弦部分被称为 F ( ω ) F(\omega) F(ω)的实部,记为 R ( ω ) R(\omega) R(ω),正弦部分为虚部,记为 I ( ω ) I(\omega) I(ω)。
频谱函数的振幅和相位分别为: ∣ F ( ω ) ∣ = R ( ω ) 2 + I ( ω ) 2 |F(\omega)|=\sqrt{R(\omega)^2+I(\omega)^2} ∣F(ω)∣=R(ω)2+I(ω)2, φ ( ω ) = arctan ( I ( ω ) R ( ω ) ) \varphi(\omega)=\arctan(\frac{I(\omega)}{R(\omega)}) φ(ω)=arctan(R(ω)I(ω))。
ω \omega ω是频率值,将 ω \omega ω推广到多维我们就能得到多维的傅里叶变换。同时,函数 F F F是线性的,有 F { a f ( x ) } = a F { f ( x ) } F\{af(x)\}=aF\{f(x)\} F{af(x)}=aF{f(x)}, F { f ( x ) + g ( x ) } = F { f ( x ) } + F { g ( x ) } F\{f(x)+g(x)\}=F\{f(x)\}+F\{g(x)\} F{f(x)+g(x)}=F{f(x)}+F{g(x)}。
式(1)被称为傅里叶分析方程(Fourier analysis equation)或傅里叶变换(Fourier transform),我们也可以用傅里叶合成方程(Fourier synthesis equation),或者说傅里叶反变换(inverse Fourier transform)从频域变换回空间域:
f ( x ) = ∫ − ∞ ∞ F ( ω ) e i 2 π ω x d ω \begin{equation} f(x)=\int_{-\infin}^\infin{}F(\omega)\mathrm{e}^{\mathrm{i}2\pi\omega{}x}\,\mathrm{d}\omega \end{equation} f(x)=∫−∞∞F(ω)ei2πωxdω
若一个函数满足 Dirac delta distribution,即 ∫ δ ( x ) d x = 1 \int\delta(x)\mathrm{d}x=1 ∫δ(x)dx=1且 ∀ x ≠ 0 , δ ( x ) = 0 \forall{}x\neq0, \delta(x)=0 ∀x=0,δ(x)=0,则有性质:
∫ f ( x ) δ ( x ) d x = f ( 0 ) \begin{equation} \int{}f(x)\delta(x)\,\mathrm{d}x=f(0) \end{equation} ∫f(x)δ(x)dx=f(0)
delta distribution 不能表示为标准的数学函数,它们通常被认为是以原点为中心,宽度接近 0 的单位面积框函数的极限。
下表展示了一些重要的函数及其频域表示,一些函数用到了 delta distribution。
Spatial Domain | Frequency Space Representation |
---|---|
Box: f ( x ) f(x) f(x) if $ | x |
Gaussian: f ( x ) = e − π x 2 f(x)=\mathrm{e}^{-\pi{}x^2} f(x)=e−πx2 | Gaussian: F ( ω ) = e − π ω 2 F(\omega)=\mathrm{e}^{-\pi{}\omega^2} F(ω)=e−πω2 |
Constant: f ( x ) = 1 f(x)=1 f(x)=1 | Delta: F ( ω ) = δ ( ω ) F(\omega)=\delta(\omega) F(ω)=δ(ω) |
Sinusoid: f ( x ) = cos x f(x)=\cos{}x f(x)=cosx | Translated delta: F ( ω ) = π ( δ ( 1 − 2 π ω ) + δ ( 1 + 2 π ω ) ) F(\omega)=\pi(\delta(1-2\pi\omega)+\delta(1+2\pi\omega)) F(ω)=π(δ(1−2πω)+δ(1+2πω)) |
Shah: f ( x ) = I I I T = T ∑ i δ ( x − T i ) f(x)={III}_T=T\sum_i\delta(x-Ti) f(x)=IIIT=T∑iδ(x−Ti) | Shah: F ( ω ) = I I I 1 / T ( ω ) = ( 1 / T ) ∑ i δ ( ω − T i ) F(\omega)={III}_{1/T}(\omega)=(1/T)\sum_i\delta(\omega-Ti) F(ω)=III1/T(ω)=(1/T)∑iδ(ω−Ti) |
Ideal Sampling and Reconstruction(理想采样和重构)
利用频率空间分析,我们现在可以正式研究采样的性质。抽样过程要求我们选择一组等间隔的样本位置,并计算这些位置上的函数值。形式上,这对应于将函数乘以一个 shah 或脉冲序列(impulse train)函数,即一个等间隔的 delta 函数的无限和:
I I I T = T ∑ i = − ∞ ∞ δ ( x − T i ) \begin{equation} {III}_T=T\sum_{i=-\infin}^{\infin}\delta(x-Ti) \end{equation} IIIT=Ti=−∞∑∞δ(x−Ti)
其中, T T T表示周期(period)或采样率(sampling rate)
shah 与原始函数的乘积可以得到采样点在原始函数上的值:
I I I T ( x ) f ( x ) = T ∑ i = − ∞ ∞ δ ( x − T i ) f ( i T ) \begin{equation} {III}_T(x)f(x)=T\sum_{i=-\infin}^{\infin}\delta(x-Ti)f(iT) \end{equation} IIIT(x)f(x)=Ti=−∞∑∞δ(x−Ti)f(iT)
如下图所示,(a)表示原始函数,(b)是 shah 函数,( c )是对应的乘积
有了这些采样值,我们再选择一个重构滤波函数 r ( x ) r(x) r(x),就可以通过卷积(convolution)来得到重构函数 f ~ \tilde{f} f~:
( I I I T ( x ) f ( x ) ) ⊗ r ( x ) ({III}_T(x)f(x))\otimes{}r(x) (IIIT(x)f(x))⊗r(x)
卷积运算 ⊗ \otimes ⊗的定义为:
f ( x ) ⊗ g ( x ) = ∫ − ∞ ∞ f ( x ′ ) g ( x − x ′ ) d x ′ f(x)\otimes{}g(x)=\int_{-\infin}^\infin{f(x')g(x-x')\,\mathrm{d}x'} f(x)⊗g(x)=∫−∞∞f(x′)g(x−x′)dx′
对于重建,卷积给出每个采样点的权值,然后通过加权和的形式得到重建函数:
f ~ ( x ) = T ∑ i = − ∞ ∞ f ( i T ) r ( x − T i ) \begin{equation} \tilde{f}(x)=T\sum_{i=-\infin}^{\infin}f(iT)r(x-Ti) \end{equation} f~(x)=Ti=−∞∑∞f(iT)r(x−Ti)
例如,使用三角重建滤波器 r ( x ) = max ( 0 , 1 − ∣ x ∣ ) r(x)=\max(0,1-|x|) r(x)=max(0,1−∣x∣)可以得到上文原始函数的重建函数:
为方便为函数的分析,我们假设 f ( x ) f(x) f(x)是 band limited,即存在 ω 0 \omega_0 ω0使得 ∀ ∣ ω ∣ > ω 0 , F ( ω ) = 0 \forall{|\omega|>\omega_0}, F(\omega)=0 ∀∣ω∣>ω0,F(ω)=0。
傅里叶分析推导出了原始函数和频域上函数间的关系,时域上的乘积是频域上的卷积:
F { f ( x ) g ( x ) } = F ( ω ) ⊗ G ( ω ) \begin{equation} F\{f(x)g(x)\}=F(\omega)\otimes{}G(\omega) \end{equation} F{f(x)g(x)}=F(ω)⊗G(ω)
相似的,时域上的卷积也是频域上的乘积:
F { f ( x ) ⊗ g ( x ) } = F ( ω ) G ( ω ) \begin{equation} F\{f(x)\otimes{}g(x)\}=F(\omega){}G(\omega) \end{equation} F{f(x)⊗g(x)}=F(ω)G(ω)
现在考虑 shah 函数 I I I T ( x ) {III}_T(x) IIIT(x),从上表可知,shah 函数的傅里叶变换也是一个 shah 函数,但是周期为 1 / T 1/T 1/T。
由上面的性质可知,使用脉冲函数 shah 进行采样时,
I
I
I
T
(
x
)
f
(
x
)
{III}_T(x)f(x)
IIIT(x)f(x)相当于
F
(
ω
)
F(\omega)
F(ω)和一个 shah 函数的卷积,而
F
(
ω
)
F(\omega)
F(ω)和 shsh 函数卷积的结果,是
F
(
ω
)
F(\omega)
F(ω)的位移和拷贝,即产生的函数图像是一个周期函数,如下图所示:
要使用这个周期频谱来重构原始函数,我们可以使用一个盒函数(Box function)与其相乘,丢弃除 0 点处的频谱副本,如下图所示,将(a)乘以(b)得到我们重构的周期频谱©:
将宽度为 T T T的盒函数定义为 Π T ( x ) \Pi_T(x) ΠT(x):
Π T ( x ) = { 1 / T ∣ x ∣ < T / 2 0 o t h e r w i s e \begin{equation} \Pi_T(x)=\begin{cases} 1/T&|x|<T/2\\ 0&\mathrm{otherwise} \end{cases} \end{equation} ΠT(x)={1/T0∣x∣<T/2otherwise
我们可以写出重构后的频谱:
F ~ = ( F ( ω ) ⊗ I I I 1 / T ( ω ) ) Π 1 / T ( ω ) \begin{equation} \tilde{F}=(F(\omega)\otimes{}{III}_{1/T}(\omega))\Pi_{1/T}(\omega) \end{equation} F~=(F(ω)⊗III1/T(ω))Π1/T(ω)
等式中的卷积象征着采样的过程,与盒函数的乘积是重构的过程,这个重构过程是非常理想的,只要 f ( x ) f(x) f(x)有 band limited,我们即可以找到需要表示这个 f ( x ) f(x) f(x)需要的最低采样频率,由此即可完美还原原始函数。
这个过程也可以在时域上进行,由于盒函数的逆傅里叶变换是 sinc 函数,因此有:
f ~ = ( f ( x ) I I I T ( x ) ) ⊗ s i n c T ( x ) \begin{equation} \tilde{f}=(f(x){III}_{T}(x))\otimes{}\mathrm{sinc}_T(x) \end{equation} f~=(f(x)IIIT(x))⊗sincT(x)
同时,由于 s i n c T ( x ) = s i n c ( T x ) \mathrm{sinc}_T(x)=\mathrm{sinc}(Tx) sincT(x)=sinc(Tx),可以将其写成离散形式:
f ~ ( x ) = ∑ i = − ∞ ∞ s i n c ( x − T i ) f ( T i ) \begin{equation} \tilde{f}(x)=\sum_{i=-\infin}^\infin \mathrm{sinc}(x-Ti)f(Ti) \end{equation} f~(x)=i=−∞∑∞sinc(x−Ti)f(Ti)
但由于 sinc 函数的范围是无限的,我们需要计算所有 f ( T i ) f(Ti) f(Ti)点,在实际应用中,一般会框定一定的范围,虽然这样会丢失部分高频信息,但效果还是较为理想的。在图形学应用中,还可能会使用盒函数代替 sinc 函数在时域上进行卷积。
Aliasing(混叠)
在理想的采样和重建中,我们做出了 band limited 的假设,但对于没有 band limited 的信号,或者对其频率内容没有以足够高的采样率采样的信号,前面描述的过程将重建一个与原始信号不同的函数。我们可以通过傅里叶分析来深刻理解这个问题。
考虑在采样频率较低的情况下,
F
(
ω
)
F(\omega)
F(ω)和
I
I
I
1
/
T
(
ω
)
{III}_{1/T}(\omega)
III1/T(ω)卷积的结果的周期太短,如下图(a)所示,
F
(
ω
)
F(\omega)
F(ω)的低频信息会泄露到高频率上,仅有(b)上显示的频率信息是真实有效的。
这种 artifacts 被称为 aliases(混叠/走样),当混叠出现在采样过程中,则被称为 prealiasing,若出现在重构过程中,则被称为 postaliasing,任何用来修复这些错误的方法都被称为 antialiasing(抗混叠/反走样)。
采样定理告诉我们,对于一个最高频率为 ω 0 \omega_0 ω0的频谱,完美重建需要 F ( ω ) ⊗ I I I 1 / T ( ω ) F(\omega)\otimes{}{III}_{1/T}(\omega) F(ω)⊗III1/T(ω)的最低频率 ω s > 2 ω 0 \omega_s>2\omega_0 ωs>2ω0,这个最低的采样频率被称为奈奎斯特频率(Nyquist frequency)。
在渲染中,采样数量和时间成本是成正比的,同时需要渲染的图像也不是 band limited 的,因此无法简单通过提高采样率来解决混叠问题。
渲染中的混叠
本书在介绍的是渲染相关的知识,因此我们必须了解在渲染中,又有哪些混叠会产生,一般又用什么方法解决
Source of aliasing
- 几何体:在渲染中,几何体的存在是最常见的混叠原因之一,对象的边界引入了阶跃函数(step function),图像函数的值立即从一个值跳转到另一个值。步进函数不仅具有的无限频率内容,更糟糕的是,当应用于混叠的样本时,重构滤波器会产生伪影,即在重构函数中出现环形伪影,这种效应被称为吉布斯现象(Gibbs phenomenon)。此外,太小的几何体也可能导致混叠,它可能会在某几帧出现,又在其他帧消失。
- 纹理:混叠的另一个来源可能来自物体的纹理和材质。没有正确的采样材质可能导致 Shading aliasing,例如如果采样率不够高,无法捕捉到闪亮表面的小高光,就会产生混叠。此外,硬阴影使用的 shadow map 引入了一个新的阶跃函数,因此也会造成混叠。
Strategy for aliasing
- Adaptive Sampling:如果我们能识别频率高于奈奎斯特极限的信号区域,我们就能使用自适应超采样,即在这些区域内进行额外的采样
- Prefiltering:我们也可以在采样前对原始函数进行模糊,抹去高频信息
Sampling Patterns
在给定固定采样率时,提高图像质量的方法仅能通过改变样本位置的分布来实现。在实践中,我们发现考虑随机采样(stochastic sampling)是值得的,即使用一个或多个随机数来确定样本位置。因此,我们需要区分算法可能生成的所有样本集的统计属性(采样模式,sample pattern)和它生成的点集。
功率谱密度(power spectral density, PSD)的概念有助于完成这项任务,对于 f ( x ) f(x) f(x)在傅里叶基中表示的函数 F ( x ) F(x) F(x),PSD 定义为:
P f ( ω ) = F ( ω ) F ( ω ) ‾ \begin{equation} P_f(\omega)=F(\omega)\overline{F(\omega)} \end{equation} Pf(ω)=F(ω)F(ω)
其中 F ( ω ) ‾ \overline{F(\omega)} F(ω)是 F ( ω ) F(\omega) F(ω)的共轭复数函数(假设 F ( ω ) F(\omega) F(ω)只有实部,则 P f ( ω ) = F ( ω ) 2 P_f(\omega)=F(\omega)^2 Pf(ω)=F(ω)2),PSD 丢弃了频谱函数的相位信息,因此无法恢复原始傅立叶系数。
PSD 的一个有用的性质是两个函数 f , g f,g f,g在空间域中的乘积的 PSD 由它们的 PSD 在傅里叶域中的卷积给出:
P f g ( ω ) = P f ( ω ) ⊗ P g ( ω ) \begin{equation} P_{fg}(\omega)=P_{f}(\omega)\otimes{}P_{g}(\omega) \end{equation} Pfg(ω)=Pf(ω)⊗Pg(ω)
由这个性质可知,对一个函数 f f f进行采样,即用 f f f与一个脉冲函数(就像 shah 函数一样)相乘可以直接表示成他们频谱函数的卷积。由于采样模式函数上的点都是 delta distribution 的,采样模式的 PSD 可以通过对样本点的求和得出。
一种简单的随机采样的方法是通过抖动(jittering),将随机偏移添加到均匀间隔的样本点上。使用 ξ ∈ [ 0 , 1 ] \xi\in[0,1] ξ∈[0,1],可以将随机样本集表示为:
s T ( x ) = ∑ i = − ∞ ∞ δ ( x − ( i + 1 2 − ξ ) T ) \begin{equation} {s}_T(x)=\sum_{i=-\infin}^\infin\delta(x-(i+\frac{1}{2}-\xi)T) \end{equation} sT(x)=i=−∞∑∞δ(x−(i+21−ξ)T)
这个采样策略的 PSD 的期望为(注意这里 PSD 的期望值直接用 PSD 函数来表示了):
P s ( ω ) = 1 − s i n c 2 ( T ω 2 ) + δ ( ω ) \begin{equation} P_s(\omega)=1-\mathrm{sinc}^2(\frac{T\omega}{2})+\delta(\omega) \end{equation} Ps(ω)=1−sinc2(2Tω)+δ(ω)
T
=
1
T=1
T=1时抖动采样的 PSD 如下图所示:
使用 T = 1 T=1 T=1时的 shah 函数和抖动采样,采样原始函数为(a)的函数得到的重构函数如下图(b), ©所示:
(b)红线表示的是 shah 函数采样重构的结果,可以看到效果非常差
©红线表示的是抖动采样函数采样重构的结果,可以看到函数图像在
∣
ω
∣
<
1
/
2
|\omega|<1/2
∣ω∣<1/2时基本完美重合,效果非常好
一般来说,采样函数的 PSD 在低频能量最小时,可以最有效的减少混叠(前提是函数
f
f
f的能量主要集中在低频上)。
shah 函数采样会导致结构化误差(structured error),而抖动采样在其 PSD 的所有较高频率中具有大致相同的能量,因此它将被采样函数 f f f中的高频能量传播到许多频率上,将混叠转换为高频噪声,这比低频噪声更能在视觉上更符合人类观察者。
PSD 有时用颜色来描述。例如,白噪声(white noise)分布在所有频率上都具有相等的功率,就像白光在所有可见频率上(或多或少)具有相等的功率一样,如下图1所示。蓝色噪声对应于一个分布,它的功率集中在较高的频率,而在较低的频率功率较少,再次对应于蓝光所表现出的功率与频率的关系,如下图2所示。
Sampling and Integration
前文讨论的重点在于采样和重构的数学解释,在实际的渲染应用中,我们需要通过蒙特卡洛积分来进行采样,前文提到过的蒙特卡洛积分的误差和混叠的概念并不相同。然而,蒙特卡罗与傅里叶分析和其他分析点采样算法的方法之间存在许多联系。例如,抖动抽样是分层采样的一种形式,即一种方差减小技术,因此,抖动采样从两个角度都是有利的。
我们希望找到一种蒙特卡罗积分的最佳抽样方法,但这个问题没有简单的答案,这是由于:
- 分析采样技术的各种数学方法经常不一致
- 渲染的图像通常是供人类使用的,我们需要考虑到人类视觉的系统
在实践中,通过采样模式生成具有蓝噪声特征的图像在视觉上是可取的,但此图像可能不会比那些没有蓝噪声特征的图像具有更低的数值误差。
Fourier Analysis of Variance
傅里叶分析也可以应用于评估蒙特卡洛积分背景下的采样模式,从而深入了解各种采样算法的方差和收敛率。为方便理解,我们在分析时做出如下假设:
- 样本点均匀分布,权值相等(即不使用重要性抽样)
- 所使用的蒙特卡罗估计是无偏的
- 采样点的性质相对于采样域上的环面平移(toroidal translation)是齐次的(如果不是,则分析有效地覆盖所有可能的样本点随机平移)
假设第三点的齐次性出现是因为我们的采样模式有很多都会进行分层,并在每个分层中放置单个样本,分层会出现不连续性和边界问题,需要齐次性来解决。
我们首先引入傅里叶级数(Fourier series)的概念,傅里叶级数使用一系列系数 f j f_j fj来表示有限域上的函数 f f f, j ≥ 0 , j ∈ Z j\geq0, j\in\mathbb{Z} j≥0,j∈Z,对于域 [ 0 , 1 ) [0,1) [0,1),傅里叶级数的参数可以表示为:
f j = ∫ [ 0 , 1 ) f ( x ) e − i 2 π j x d x \begin{equation} f_j=\int_{[0,1)}f(x)\mathrm{e}^{-\mathrm{i}2\pi{}jx}\,\mathrm{d}x \end{equation} fj=∫[0,1)f(x)e−i2πjxdx
用傅里叶级数系数 f j f_j fj表示原始函数 f f f,有:
f ( x ) = ∑ j ∈ Z f j e − i 2 π j x \begin{equation} f(x)=\sum_{j\in\mathbb{Z}}f_j\mathrm{e}^{-\mathrm{i}2\pi{}jx} \end{equation} f(x)=j∈Z∑fje−i2πjx
可以证明,连续傅里叶变换对应于无限范围取傅里叶级数的极限。
函数在傅里叶级数基中的 PSD 由每个系数和它的复共轭的乘积给出的:
P f ( j ) = f j f j ‾ \begin{equation} P_f(j)=f_j\overline{f_j} \end{equation} Pf(j)=fjfj
为了在频域上分析蒙特卡洛积分,我们将采样函数 s ( x ) s(x) s(x)设定为 n n n个采样点 x i x_i xi的平均值,采样点用 delta distribution 表示:
s ( x ) = 1 n ∑ i = 1 n δ ( x − x i ) \begin{equation} s(x)=\frac1n{}\sum_{i=1}^n\delta(x-x_i ) \end{equation} s(x)=n1i=1∑nδ(x−xi)
有了采样函数后,我们可以将蒙特卡洛估计器写成积分形式:
∫ [ 0 , 1 ) f ( x ) d x ≈ 1 n ∑ i = 1 n f ( x i ) = ∫ [ 0 , 1 ) f ( x ) s ( x ) d x \begin{equation} \int_{[0,1)}f(x)\,\mathrm{d}x\approx\frac1n{}\sum_{i=1}^nf(x_i)=\int_{[0,1)}f(x) s(x)\,\mathrm{d}x \end{equation} ∫[0,1)f(x)dx≈n1i=1∑nf(xi)=∫[0,1)f(x)s(x)dx
将(18)式带入上式,此积分形式便可以用傅里叶级数表示,有:
∫ [ 0 , 1 ) f ( x ) s ( x ) d x = ∑ j ∈ Z f j ‾ s j \begin{equation} \int_{[0,1)}f(x) s(x)\,\mathrm{d}x=\sum_{j\in\mathbb{Z}}\overline{f_j}s_j \end{equation} ∫[0,1)f(x)s(x)dx=j∈Z∑fjsj
我们知道 f 0 = f 0 ‾ = ∫ f ( x ) d x f_0=\overline{f_0}=\int{}f(x)\,\mathrm{d}x f0=f0=∫f(x)dx,同时由于我们假设了 s ( x ) s(x) s(x)是均匀无权值采样,因此有 s 0 = 1 s_0=1 s0=1,因此,蒙特卡洛积分器的误差为:
∣ ∫ [ 0 , 1 ) f ( x ) d x − ∫ [ 0 , 1 ) f ( x ) s ( x ) d x ∣ = ∣ f 0 − ∑ j ∈ Z f j ‾ s j ∣ = ∑ j ∈ Z ⋆ f j ‾ s j \begin{equation} |\int_{[0,1)}f(x) \,\mathrm{d}x-\int_{[0,1)}f(x) s(x)\,\mathrm{d}x|=|f_0-\sum_{j\in\mathbb{Z}}\overline{f_j}s_j|=\sum_{j\in\mathbb{Z^\star}}\overline{f_j}s_j \end{equation} ∣∫[0,1)f(x)dx−∫[0,1)f(x)s(x)dx∣=∣f0−j∈Z∑fjsj∣=j∈Z⋆∑fjsj
其中 Z ⋆ \mathbb{Z}^\star Z⋆是所有非 0 的整数
深入理解蒙特卡洛积分器的误差,假如 f f f是 band limited 的,那么 ∃ j m a x , ∀ j > j m a x , f j = 0 \exists j_{max}, \forall j>j_{max}, f_j=0 ∃jmax,∀j>jmax,fj=0,在这种情况下,如果采样频率够高,使得 ∀ 0 < j < j m a x , s j = 0 \forall 0<j<j_{max},s_j=0 ∀0<j<jmax,sj=0,那么估计器便是零方差的。
由此式便可以使用 f ( x ) f(x) f(x)和 s ( x ) s(x) s(x)的 PSD 来表示估计器的方差:
V [ 1 n ∑ i = 1 n f ( x i ) ] = ∑ j ∈ Z ⋆ P f ( j ) P s ( j ) \begin{equation} V[\frac1n\sum_{i=1}^nf(x_i)]=\sum_{j\in\mathbb{Z^\star}}P_f(j)P_s(j) \end{equation} V[n1i=1∑nf(xi)]=j∈Z⋆∑Pf(j)Ps(j)
这个公式很好理解,我们希望在函数 f ( x ) f(x) f(x)的 PSD 高的时候采样模式的 PSD 低一些,但函数 f ( x ) f(x) f(x)通常是不可解析的,我们一般认为函数 f ( x ) f(x) f(x)在低频率下功率较高,在这个假设下,我们的采样模式需要尽可能在低频率下有较低的功率,这个假设和我们之前看到的蓝噪声相同。
特别的,在均匀随机采样(白噪声)下, P s = 1 / n P_s=1/n Ps=1/n常数,方差便会等于:
1 n ∑ j ∈ Z ⋆ P f ( j ) = O ( 1 n ) \begin{equation} \frac1n\sum_{j\in\mathbb{Z^\star}}P_f(j)=O(\frac1n) \end{equation} n1j∈Z⋆∑Pf(j)=O(n1)
这与蒙特卡洛章节导出的方差相同。
更进一步,如果采样技术的 PSD 可以渐近有界,则可以证明该技术在给定合适的函数进行积分时表现出更高的方差缩减率,比如在 2D 中,给定一个平滑的被积函数,抖动采样模式可以实现 O ( n − 2 ) O(n^{-2}) O(n−2)方差。
傅里叶分析还揭示了泊松盘采样模式的渐近收敛性。泊松盘点集的构造使得两点之间的距离不能小于某一最小距离,多年来,它们被认为优于抖动采样的。泊松盘采样模式禁止多个样本聚集在一起,如下图(b)所示:
比较抖动采样和泊松盘采样的的 PSD 图像和曲线可以发现,泊松盘采样在 PSD 较低的原点周围具有更大的频率范围,象征着具有更优越的蓝噪声特征。
Low Discrepancy and Quasi Monte Carlo
在傅里叶分析之外,另一种评估样本点质量的有用方法是基于一个叫做差异(discrepancy)的概念。分布良好的采样模式具有较低的差异,因此样本模式生成问题可以看作是寻找一个合适的低差异点模式的问题。
我们首先需要区分样本集(sample sets)和样本序列(sample sequences):
- 样本集:特定数量的点集,位置确定
- 样本序列:生成任意数量样本点的算法
随后引入差异的概念,为了计算一组点的差异,我们首先选择一系列形状 B B B,它们是 [ 0 , 1 ) d [0,1)^d [0,1)d的子集。例如,使用一个角位于原点的 box 盒子。这对应于 B = { [ 0 , v 1 ] × [ 0 , v 2 ] × . . . × [ 0 , v d ] } B=\{[0,v_1]\times[0,v_2]\times...\times[0,v_d]\} B={[0,v1]×[0,v2]×...×[0,vd]},其中 0 ≤ v i < 1 0\leq{}v_i<1 0≤vi<1。给定 n n n个采样点 P = { x 1 , . . . , x n } P=\{x_1,...,x_n\} P={x1,...,xn},相对于 B B B的差异 P P P 表示为:
D n ( B , P ) = sup b ∈ B ∣ ♯ { x i ∈ b } n − V ( b ) ∣ \begin{equation} D_n(B,P)=\sup_{b\in{}B}|\frac{\sharp\{x_i\in{b}\}}{n}-V(b)| \end{equation} Dn(B,P)=b∈Bsup∣n♯{xi∈b}−V(b)∣
其中 ♯ { x i ∈ b } \sharp\{x_i\in{b}\} ♯{xi∈b}是 b b b中点的数量, V ( b ) V(b) V(b)是 b b b的体积
差异的定义符合直觉的理解是,我们在使用一个盒子内点的数量 ♯ { x i ∈ b } / n \sharp\{x_i\in{b}\}/n ♯{xi∈b}/n来近似盒子的体积 V ( b ) V(b) V(b),差异是在所有可能的方框中,得到误差最大的近似方案。
当 B B B是一个角位于原点的盒子是,这个值被称为星型差异(star discrepancy) D n ⋆ ( P ) D^\star_n(P) Dn⋆(P)。此外还有 axis-aligned boxes 等等。
对于某些点集,差异有解析解,比如在一维情况下 x i = i n x_i=\frac{i}{n} xi=ni,则 x i x_i xi的星型差异为 D n ⋆ = 1 n D_n^\star=\frac1n Dn⋆=n1( b = [ 0 , 1 / n ) b=[0,1/n) b=[0,1/n)时)。当 x i = i − 1 / 2 n x_i=\frac{i-1/2}{n} xi=ni−1/2时, D n ⋆ = 1 2 n D_n^\star=\frac1{2n} Dn⋆=2n1。
我们可以看到,当样本均匀分布时,差异最小,但傅里叶分析表明,抖动优于均匀采样。幸运的是,高维中的低差异模式远不如一维中的均匀,因此在实践中通常作为示例模式工作得相当好。
我们说 d 维上的点序列具有低差异(low discrepancy)需要满足差异的数量级为:
O ( ( log n ) d n ) \begin{equation} O(\frac{(\log{n})^d}{n}) \end{equation} O(n(logn)d)
低差异点集和序列通常使用特定的算法生成,使用这样的点来采样函数进行积分,我们就得到了准蒙特卡罗(quasi–Monte Carlo, QMC)方法。在常规蒙特卡罗算法中使用的许多技术都可以很好地处理这种准随机样本点。
Koksma–Hlawka inequality 用于将积分的一组点的差异与函数积分的估计误差联系起来,写作:
∣ ∫ f ( x ) d x − 1 n ∑ i f ( x i ) ∣ ≤ D n ( B , P ) V f \begin{equation} |\int{f(x)\,\mathrm{d}x}-\frac1n\sum_if(x_i)|\leq{}D_n(B,P)V_f \end{equation} ∣∫f(x)dx−n1i∑f(xi)∣≤Dn(B,P)Vf
其中, V f V_f Vf是被积分函数 f f f的总变分(total variation),有 V f = sup 0 = y 1 < y 2 < . . . < y m = 1 ∑ i = 1 m ∣ f ( y i ) − f ( y i + 1 ) ∣ V_f=\sup_{0=y_1<y_2<...<y_m=1}\sum_{i=1}^m|f(y_i)-f(y_{i+1})| Vf=sup0=y1<y2<...<ym=1∑i=1m∣f(yi)−f(yi+1)∣
从本质上讲,总变化表示函数值在点之间变化的速度有多快,差异表示用于积分的点在捕捉函数变化方面的有效性。TODO
从 Koksma-Hlawka inequality 可以看出,随着被积量维数的增加,低差异的积分误差越来越接近 O ( n − 1 ) O(n^{-1}) O(n−1),它的渐进性比蒙特卡洛积分的 O ( n − 1 / 2 ) O(n^{-1/2}) O(n−1/2)好。实践中 QMC 的收敛速度会比较好。
由于 QMC 积分器是确定性的,因此不可能使用方差作为估计器质量的度量,尽管当然仍然可以计算均方误差。另外,可以使用不损害其差异的方法对样本点进行随机化。随机化甚至可以提高收敛速度,这种方法被称为随机拟蒙特卡罗(randomized quasi–Monte Carlo, RQMC)。
Reference
https://www.pbr-book.org/4ed/Sampling_and_Reconstruction
PS
本文主要是对PBRT内容的翻译总结,欢迎交流学习~