采样和重建(Sample And Reconstruction)

Motivation: 渲染的输出是一张2D的像素网格, 是离散的. 入射辐射亮度实际上是定义在胶片平面的连续函数. 从连续到离散就必定会带来质量下降. 下降的多少决定了最终渲染输出的效果的好坏. 而从连续信号到离散信号的过程一般称为采样. 良好的采样方法应该产生低差异(low-discrepancy)采样点. 采样结束之后我们还希望从采样结果中能恢复原始信号, 这个过程叫做重建.

##采样理论
数字图像是一系列像素值的网格. 像素值用来决定显示的物理设备发射的光谱能量. 首先需要区分图像像素和显示像素. 图像像素是特定采样位置的函数值, 而显示像素是以一定分布发射光的物理器件. 物理器件使用图像像素值来构建显示设备表面的新的图像函数, 这个函数在显示设备的所有点都有定义, 而不仅仅是图像像素的点. 因而, 我们生成数字图像的过程实际上是将物理器件的连续图像函数用数字图像的像素离散化采样得到的.
采样中不得不提到的一个重要技术就是傅里叶分析. 傅里叶分析从频率与考虑信号的采样和重建的相关问题, 并可以用来评估重建信号和原始信号的匹配度.

Fig.1: 采样和重建

频率域和傅里叶变换

傅里叶分析的基础之一就是傅里叶变换, 他是频率域中的一个函数. 事实上大多数函数都可以分解成偏移的正弦曲线的加权和. 这个惊人的事实第一次被Joseph Fourier发现, 而傅里叶变换就是进行这个分解的过程. 频率域可以给我们一个原始信号某些特点的更深入的理解.1D函数\(f(x)\)的傅里叶变换是:
$$ F(\omega) = \int_{-\infty}^\infty f(x) e^{- \boldsymbol{i} 2\pi \omega x} dx$$
其中i是复数项, i =\(\sqrt{ -1}\). 并且
$$ e^ {\boldsymbol{i} x} = cosx + \boldsymbol{i} sinx$$
从频率域当然也可以变换回空间域:
$$ f(x) = \int_{-\infty}^\infty F(\omega) e^{ \boldsymbol{i} 2\pi \omega x} d \omega$$

如果想要更深刻的理解, 建议去查阅更详细的傅里叶分析介绍.
我们考虑空间域常数为1的函数
$$ f(x) = 1$$

这个函数可以认为频率为0, 也就是没有任何变化. 但是整个基本区域的积分值为1, 所以频率为0的函数值可以认为无限大(并不是无穷大, 但是可以简单理解为无穷大). 这个函数是傅里叶分析中一个特殊的函数, 叫做Dirac delta函数
$$
\begin{cases}
\int \delta(x)dx =1 \\
\delta(x) =0, x \ne 0
\end{cases}
$$

Dirac delta函数的一个重要性质是:
$$\int f(x) \delta(x) dx = f(0)$$

理想采样和重建

我们将采样过程用方程描述.采样实际上是将均匀分布的delta函数的无穷累加和一个叫做脉冲序列(“shah” or “impulse train”)函数相乘得到的结果.脉冲序列函数定义如下:
$$\bigvee_T(x) = T \sum_{i= {-\infty}}^{\infty} \delta(x- i T)$$

其中T叫做采样频率,或者周期.将脉冲序列函数和原始信号相乘,就得到了均匀的无穷序列采样值.
$$S(x) = \bigvee_T(x)f(x) = T \sum_i \delta(x - i T)f(iT)$$

上面的采样信号想要重建的话,就是运行一个重建滤波函数并进行卷积操作:

$$\tilde{f}(x) = S(x) \otimes r(x)$$, 而\(\otimes \)的运算是:
$$f(x) \otimes g(x) = \int_{- \infty}^{\infty} f(x’) g(x - x’) dx’$$
所以可以得到:
$$
\tilde{f}(x) = T \sum_i \delta(x - i T)f(iT) \otimes r(x) \\
= \int_{- \infty}^{\infty} T \sum_i \delta(x’ - i T)f(iT) r(x - x’) dx’ \\
= T \sum_{i= {-\infty}}^{\infty} f(iT) r(x - iT)
$$
傅里叶变化的一个重要性质就是:两个函数的乘积的傅里叶变换等于每个函数单独傅里叶变换的卷积, 即
$$\mathcal{F} \{f(x)g(x)\} = F(\omega) \otimes G(\omega)$$
反之亦然.即空间域和频率域的卷积和相乘可以互相转化.因而上面的重建操作的卷积可以用频率域的相乘代替. shah函数及其对应的频率域表达分别是:
$$ \bigvee_T(x)= T \sum_i \delta(x - i T)\Longleftrightarrow \bigvee_{1/T}(\omega) = (1/T) \sum_i \delta (\omega - i / T)$$

这里值得注意的一点是,shah函数的空间域和频域函数的周期是倒数关系,这意味着如果空间域的样本距离很大,那么频率域的样本距离很小.
因而采样信号的频率域表达\(\mathcal{F} (S(x))\)可以由\(F(\omega)\)和上面的shah的频率域函数的卷积得到.而将函数与delta函数进行卷积得到的就是原函数的若干复制,所以shah函数卷积会得到原始函数的无穷多的序列复制,等距分布, 周期和shah的周期一样.这就是样本序列的频率空间的表达!!

Fig.2: 样本频率空间和Shah函数卷积结果

很显然,想要得到原始信号,应该将除了原点周围的信号全部去除掉.我们可以采用一个箱函数box function来进行这个操作. 假定箱函数宽度为\(T\),则箱函数定义是:
$$
\Pi_T(x)=\begin{cases}
1 / 2T, & \text{如果$ \vert{x}\vert <T $} \\
0, & \text{其他情况}
\end{cases}
$$
频率域的相乘意味着空间域的卷积,所以这里我们相乘得到重建信号的频率表示:
$$\tilde{F} = (F(\omega) \otimes \bigvee_{1/T}(\omega)) \Pi_T(\omega)$$

我们重新梳理一遍.

空间域: 原始信号\(f(x)\), 周期为T的采样函数shah函数\(\bigvee_T(x)\), 采样信号\(S(x)\), 重建信号\(\tilde{f}(x)\).
频率域:原始信号\(F(\omega)\), shan函数\(\bigvee_{1/T}(\omega)\), 采样信号\(S(\omega)\),重建信号\(\tilde{F}(\omega)\)

则空间域
$$S(x) = \bigvee_T(x)f(x)$$
$$\tilde{f}(x) = S(x) \otimes r(x)$$
频率域:
$$S(\omega) = \bigvee_{1/T}(\omega) \otimes F(\omega)$$
$$\tilde{F}(\omega) = S(\omega) \Pi_T(\omega)$$

所以要想完美重建信号,需要找到完美箱函数\(\Pi_T(x)\)的空间域表达.箱函数的逆傅里叶变换是\(sinc\)函数,所以理想的空间域重建是:
$$\tilde{f}(x) = S(x) \otimes r(x) = \sum_{i= - \infty}^ \infty sinc(x - i)f(i)$$

悲剧的是,因为sinc函数具有无限的边界,所以必须使用\(f(i)\)的所有采样值来计算空间域的某个特定值\(\tilde{f}(x)\). 实际使用只会采取一定范围之内.

Fig.3: 理想重建过程

混叠(aliasing)

sinc函数的无边界特性只是一个问题.另外一个理想采样和重建过程中的严重问题是其对于信号都是假设有边界的.对于无限边界的信号,如果采样频率不够的话,将会重建出不同于原信号的结果.
成功重建信号的关键在于精确覆盖原信号频率域\(F(\omega)\)以一定宽度的箱函数,如图三所示.但是多种情况会导致无法完美重建.第一种情况,原始信号的采样频率过低.我们之前说过空间域和频率域的shah函数的周期是倒数关系,所以空间域样本间距过大会导致频率域样本间距过小,从而使得\(F(\omega)\)的复制样本之间重叠.这种情况下无论如何也不可能完美恢复原信号.这会造成原始信号的高频部分信息泄漏到低频区域,从而遗失高频信号,并且扰乱低频信号.这些新的低频部分的干扰叫做混叠.

这种情况导致的混叠很容易去除:增大采样频率.采样理论告诉我们,如果原始信号的最大频率是\(\omega_0\),那么均匀采样\(\omega_s\)的频率必须大于\(2\omega_0\).\(2\omega_0\)的频率是一个分界线,叫做奈奎斯特频率(Nyquist frequency)

对于无边界的信号,他们在频率域无穷范围都有值,所以不可能不重叠.而不幸的是,计算及图形学中只有很少的感兴趣函数是有边界的.因而,除了单纯增大采样频率,我们必须有更好的方法来抵消混叠效果.

Fig.3: 采样频率过低

混叠消除技术(消除走样技术, anti-aliasing)

经过前面的讨论我们知道,混叠会来源于两方面:采样频率过低和重建过程.我们将前者称为prealiasing, 后者称为postaliasing. 消除二者的技术统称为antialiasing.下面介绍一些除了提高采样频率之外的方法

非均匀采样(NonUniform Sampling)

尽管自然图像函数具有无穷的频率部分, 无法从点集中完美重建,但是我们可以通过以非均匀的方式变化样本之间的间隔来减少混叠带来了的视觉影响. 用\(\xi \)表示0~1之间的随机数,则基于脉冲序列的非均匀采样点集是:
$$ \sum_{i=-\infty}^{\infty}(x-(i+1/2 - \xi)T)$$

当采样频率不足时,均匀采样和非均匀采样都不能完美重建信号,但是费均匀采样将混叠转变成了噪声,后者比前者对于人的视觉系统更容易接受.从频率域考虑,非均匀采样信号的复制的偏移更加随机,因而重建的时候会导致随机的噪声而不是相关的混叠

自适应采样(Adaptive Sampling)

这种方法的思想是:找到信号频率高于奈奎斯特频率的地方,在这些区域进行额外的采样.这样只会在需要的地方进行高频率采样,其他地方保持低频率采样.这种方法实际很难实现,因为找到高频率区域很困难.大多数实现方法是检查邻接采样值,找到变化很大的区域,并认为这样的区域是信号中的高频区域.
但是实际上邻接采样值并不能确定性的告诉我们真正发生了什么.即使两个采样值一样,中间也可能发生很大变化;而不同的值可能也不是混叠造成的.

前置过滤(Prefiltering)

该方法过滤掉原始信号的高频部分.例如一些材质函数会用到.当然这会带来原始信号的信息丢失.但是通常这种缺失比混叠更不容易发现.
另外一种方法,会像空间域的理想重建过程
$$S(x) \otimes r(x) = S(x) \otimes sinc(2\omega_sx))$$
理想的sinc对于低频率采样会带来混叠,我们可以认为将其截断,采用一个近似sinc的低频信号来卷积.这种方法会导致原始信号的高频信号丢失,但是低频信号的混叠也会消失

图像生成的应用

2D图像中的采样和重建也可以应用上面方法. 图像可以认为是位置到辐射亮度L的函数:
$$f(x,y) \rightarrow L$$
在光线追踪中,图像中的任意位置\((x,y)\)都可以得到\(L\), 但是一般不太可能对\(f\)运用前置过滤的方法来过滤掉高频部分.因此, 将重点关注两种策略:增大采样频率, 非均匀采样.
更一般的,从相机发射的光线基于五个量:时间\(t\), 采样的透镜位置\((u,v)\), 图像位置\((x,y)\).对于特定的图像位置,该点的辐射亮度随着时间(场景中的运动物体)和镜头位置(相机是有限孔径的透镜)而改变.
更一般的,许多积分器都会采用统计方法来估计特定光线的辐射亮度.如果我们将图像函数拓展到统计中使用的样本值的话, 会更加高维:
$$ f(x,y,t,u,v,i_1,i_2,…) \rightarrow L$$

渲染中混叠的来源

几何体是渲染图像中混叠来源的最通用来源之一.当投影到图像平面,物体的边界是一个阶梯函数(step function),其具有无穷高的频率.当完美重建滤波器进行重建时,会产生振铃效应(ringing artifacts),这种现象叫做吉布斯现象(Gibbs phenomenon).
场景中非常小的物体也会造成几何混叠.如果几何体足够小以至于在图像平面上位于采样样本之间,那么在动画里面他就会不可预期地出现和消失.
另外一种混叠来源于物体的贴图和材质.着色混叠(Shading aliasing)可能来源于没有正确滤波的贴图, 或者光滑表面的小的高光区域.如果采样频率没有高到可以采样这些特征,就会产生混叠.
物体产生的尖锐的阴影边界也会产生阶梯函数.这更加难以检测.
渲染图像中的混叠永远不可能消失.我们能做的只是降低混叠的影响.

采样接口

前面我们提到过,图像函数是多维输入到辐射亮度的函数.采样就是产生多维输入的过程.我们将输入归一化到[0,1] 范围, 则采样就是产生\({[0,1]}^n\)的n维样本序列.最简单的采样实现就是[0,1]的均匀采样.这种采样需要更多的样本,更长的时间才能产生和其他复杂采样器相同效果的图片.简单来说,采样器越简单,需要的采样数就越多.而采样器本身的运行消耗差不多是一样的,所以更好的采样器会有效降低渲染时间.

延伸阅读:评估采样模式的指标: 差异(Discrepancy)

傅里叶分析给了我们评估2D采样模式的一种方法.但是我们只能从增加更多采样样本中提高效果.在面对图像中的无穷频率和高维采样向量,傅里叶分析往往独木难支.
除了傅里叶分析,数学家发明了一个叫做差异的概念来评估n维采样模式的质量.分布良好的样本往往都是低差异的.幸运的是,高维数据的低差异采样我们有许多确定性的技术可以生成.
差异性可以这样理解:在\([0,1]^n\)范围内,取某个特定的区域, 在这个区域内计算样本数量,并计算样本占总样本的比例.这个比例和该区域面积占总面积的比例认为是一样的.如果不一样,就产生了差异. 所以这样的区域中差异最大的, 就是整个采样的差异.一般来说,样本分布越均匀,差异性越小
如果你还记得,前面说过不均匀的扰动模式在视觉感知上要优于均匀模式,因为他们将混叠转变成了噪声.而现在的低差异则要求采样尽量均匀.确实如此,地差异性采样视觉上更容易受到混叠的影响.幸运的是,高维里面的低差异采样通常会”不是特别均匀”.
一种称为泊松盘(Possion disk)的算法可以用来产生比较好的低差异采样.泊松盘采样中任意两个点的距离都大于d.研究表明人眼的锥细胞和杆细胞具有类似的分布.实际中采样2D图像泊松盘非常有效,但是对于更高维的采样, 其效果不是特别好.

分层采样（Stratified Sampling）

分层采样是将原始图像区域分为多个矩形区域, 在每个区域内部产生一个采样样本.这些区域通常称为”地层”(strata).为什么这么做?这样的话可以确保不会有采样距离过近,因为我们就不太可能丢失图像的重要部分信息.另一方面,距离过近的采样有的时候并不是必要的,因为可能信息接近造成冗余.从信号处理角度,我们隐式定义了一个采样频率,地层面积越小,个数越多,则采样频率越高.
分层采样在每个地层中随机点生成采样,随机点的位置是用地层的宽和高的一半的随机数来扰动地层中心位置.这里的扰动造成的非均匀性有利于将混叠转化成噪声.当然也可以产生均匀采样.这种时候往往用来进行不同渲染算法的比较,而不是生成高质量的图片.
直接进行高维采样的分层会导致样本数量的爆炸.例如图像函数认为是5D函数的话, 如果我们每个维度都生成四个地层,那么每个像素的采样总数就是\(4^5 = 1024\)!!我们可以在某些维度降低采样个数,但是这样会丢失那些维度的分层采样的好处.这个问题叫做维度灾难(curse of dimensionality). 实际使用并不是每个分层相乘,而是相加. 也就是每个维度单独进行分层采样,然后将不同维度的分层进行随机组合,这样的总个数就是\(4^5 = 20\)

Fig.4: 多维度分层采样

延伸阅读:哈尔顿采样(Halton Sampler)

分层抽样还是有可能出现不同地层之间的采样距离过近的情况.哈尔顿采样是一种低差异采样,其会确保产生的采样点不会距离过近,同时又可以在所有维度均匀分布.

哈默斯利序列和哈尔顿序列(Hammersley and Halton Sequence)

这两个序列都都是低差异的,都基于”倒根函数”(radical inverse).

延伸阅读:(0,2)序列采样

这是另一种产生高质量采样的低差异序列采样方法.

延伸阅读:最大最小距离采样

延伸阅读:苏泊尔采样(Sobol Sample)

图像重建

无论采用什么采样方法,我们现在得到了图像样本,现在需要将样本值转变成最终像素的辐射亮度值.根据信号处理理论,我们需要做三个步骤:

• 从采样数据中重建连续信号 \(\tilde{F}\)地方
• 对\(\tilde{F}\)进行前置滤波,移除超过奈奎斯特频率的而部分
• 对\(\tilde{F}\)进行像素位置的采样,得到最终像素值

因为我们知道我们队\(\tilde{F}\)在像素位置进行重新采样,所以没必要构建\(\tilde{F}\)的显示表达.可以采用一个滤波函数实现这个过程.
再次回顾理想重建的过程:均匀以高于奈奎斯特频率对原信号进行采样,并以sinc滤波器进行重建,那么重建信号将与原信号完美匹配.但是,由于图像函数总是有更高的频率,所以这个过程不太可能实现.因为我们采样非均匀采样,将混叠转变成噪声.
尽管人们试了很多方法来试图用非均匀采样完美重建,但是还没有一个可接受的方法实现.现实中完美重建基本不可能实现.因而,显示中的重建在于重建函数和原始函数的误差最小,而不局限于原始信号是否是有界的.
重建像素值可以用下面公式来表达:
$$I(x,y) = \frac{\sum_i f(x- x_i, y -y_i)\omega(x_i, y_i) L(x_i, y_i)}{\sum_if(x-x_i, y -y_i)}$$
其中\(f\)是滤波函数,\(\omega\)是采样分布权重, \(L\)是采样辐射亮度值. sinc滤波器在这里并不是很好的选择:一方面,它会产生吉布斯现象;另一方面,他是无边界的.事实上,滤波器的选择并没有最好一说.每个场景可能都有最适合的滤波器.

盒子滤波器 和三角滤波器(Box Filter And Triangle Filter)

最常用的滤波器之一.最简单,也可能是效果最差的.前面我们说过盒子滤波器会允许采样数据的高频信息泄漏到重建值中.这导致了postaliasing:即使原始采样值频率足够高,但是滤波器引入了新的误差.
三角滤波器要好一点点,但是本质还是存在上面问题.

高斯滤波器(Gaussian Filter)

这可能是实际应用中效果比较好的一种了.高斯分布将造成轻微的模糊.

米切尔滤波器(Mitchell Filter)

滤波器的设计是很困难的事情,混合了数学分析和感知试验.Mitchell在1988年进行了一系列试验,并发明了一种在振铃效应和模糊之间平衡的滤波器:米切尔滤波器.米切尔滤波器很像高斯滤波器,但是其具有负数的末端分布.这些负数部分造成了边缘部分的尖锐效果,会减少模糊.其缺点是如果负数部分过大,振铃效应就会出现,并且可能需要阶段来保证没有负数输出.

Fig.4: 高斯和米切尔滤波器

窗口sinc滤波

将原始无边界的sinc滤波乘以某个窗口函数变成有边界的,这回降低振铃效应