PBRT阅读：第七章 采样和重构 第7.1节

http://www.opengpu.org/forum.php?mod=viewthread&tid=4805&fromuid=10107

第七章 采样和重构

虽然象pbrt这样的渲染器的最终输出结果是二维的彩色像素阵列，但实际上入射辐射亮度是定义在胶片平面上的连续函数。从这个连续函数计算出离散像素值的方式对最终的图像质量有着显著的影响；如果这个计算过程有纰漏，图像上就会出现人为的误差。幸运地是，只需做一些额外的计算，所渲染的图像质量会有实质性的提高。

本章介绍采样理论，即对于在连续域上定义的函数，讨论如何获取函数上离散的采样值，并用这些采样值重构出和原函数很接近的新函数。以采样理论为基础，本章所要介绍的采样器(Samplers)在胶片平面上的入射辐射亮度函数上选取采样点。在本章里，我们要介绍三个Sampler类的实现，它们覆盖了许多关于采样问题的解决方法。在这一章的最后还介绍了Filter类。Filter类用于确定如何把每个像素相邻的采样进行混合，从而得到最终的像素值。

7.1 采样理论

一幅数字图像通常是由一个长方形的像素值阵列来表示的。当一个物理设备显示一个数字图像时，这些值用来设置显示器上的像素的强度和颜色。当我们研究数字图像时，把图像像素和显示像素区分开来十分重要：图像像素代表在一个特定采样位置上的函数值，而显示像素是可以放射出某种光谱分布的光的物理实体。例如，在CRT显示器上，每个萤光点放光，形成沿其中心向外衰减的区域分布。像素强度也有一个角向分布(angular distribution),对于CRT显示器而言它是均匀分布的，而对LCD显示器而言它却是有方向性的。显示器用图像像素值来构造一个在显示器表面上的新图像函数。这个新函数对显示器上所有的点都有定义，而并非局限于像素上的点。这个获得一组采样值并把它们转换成一个连续函数的过程就是是重构(reconstruction)。

为了计算数字图像上的离散像素值，有必要对原连续图像函数进行采样。跟其它基于光线追踪的渲染器一样，pbrt获取图像函数信息的唯一方式是用光线追踪算法来采样。例如，并没有通用的方法来计算胶片平面上两个点的函数值差异的界限值。虽然我们可以只在像素位置上采样来生成图像，但是更好的方法是在不同的位置上再采集额外的多个采样值，然后再合成为最终的像素值。实际上，为了取得最佳效果，像素值的计算结果要使得重构出的在显示器上的图像和虚拟相机胶片上的原场景图像尽可能地一致。要注意的是，这个目标并不等同于仅仅令显示器上的像素取图像函数在像素位置上的实际值。 本章所介绍的算法的主要目标就是处理这个不同之处。

因为采样和重构过程使用了近似技术，这就会出现所谓“走样”(aliasing)的失真现象，这种失真会有多种情况，其中包括锯齿边和动画中的闪烁现象。失真现象的出现是因为采样过程并不能够获取连续图像函数的全部信息。

作为一个例子，考虑一个一维函数f(x)，我们在函数的定义域中求位置x'的函数值f(x')。x'即被称为采样位置，f(x')被称为采样值。下图显示了一个平滑的一维函数，还有一个被重构出的函数f~用来拟合原函数f。在这个例子中，f~是一个由折线构成的函数，是通过对相邻的采样值进行线性插值的方法来拟合原函数f的。因为我们对原函数f的信息只局限于这些采样点，f~就不太可能对f有完美的匹配。

我们可以用傅立叶分析的理论来评估重构函数和原函数之间的匹配质量。本节将介绍采样和重构过程所需要的关于傅立叶分析的主要思想，并略去了许多证明和无关乎采样算法的细节。
 
7.1.1 频域和傅立叶变换

傅立叶分析的基础之一是傅立叶变换，它是一个定义在频域(frequency domain)上的函数(我们说函数一般地是在空间域spatial domain上定义的）。下图的(a)显示了一个关于x的变化缓慢的函数，而(b)显示了一个变化快得多的函数。
 
而下图是上面两个函数在频率空间中的表示。我们可以看到低频率的那个函数比高频率的函数更快地趋向零。
 
绝大多数函数可以被分解为具有相移的正弦函数的加权和。这个重大发现是首先由Joseph Fourier做出的，傅立叶变换就是把一个函数转换为这种表示。一个函数在频域中的表示可以显现出一些重要特征：正弦函数的频率分布对应于原函数的频率分布。利用这种变换，就可以更深入地了解采样和重构过程所产生的失真现象，并研究如何减少这种失真现象的方法。
一维函数f(x)的傅立叶变换为：
           
(其中e ^ix = cosx + i sinx)为了简单起见，我们只考虑偶函数，即f(-x)=f(x)，这时f的傅立叶变换没有虚数项。新函数F是关于频率ω的函数。我们用符号 F来表示傅立叶变换算子，即 F{f(x)} = F(ω)。显然 F是线性的，即 F{a f(x)} = a  F{f(x)}， F{f(x) + g(x) } =  F{f(x)} +  F{g(x)}。

上述等式称为傅立叶分析式，又称傅立叶变换。我们也可以利用傅立叶合成式（或称傅立叶逆变换）从频域变换会空间域：
               
下表显示了一些重要函数在频率空间中的表示：

（表7.1）

有许多函数基于Dirac  δ分布（又称单位脉冲函数），它是满足ʃ δ(x)dx = 1 （x不等于零，且δ(x)=0）的特殊函数。一个重要的性质为：
           
delta函数无法用标准的数学函数来表示，但是它通常可以被想象成一个中心在原点的单位方盒函数当宽度趋向于零时的面积的极限。

7.1.2 理想的采样和重构

利用频率空间分析，我们可以比较正式地研究采样的性质。采样过程需要我们选择一组等间距的采样位置并计算出这些位置上的函数值。正式的说法是用一个“shah”或“脉冲序列(impulse train)”函数（即一个等距的delta函数的无限项的求和函数）来乘以这个函数。shah函数III _T(x)定义为：
               
其中T定义了周期，或采样速率。这个采样过程的正式定义如图所示。
 
相乘的结果就是一个在等距采样点上的函数值所构成的无限序列：
 
我们选择一个滤波函数r(x)，就可以利用这些采样值求如下的卷积还定义重构函数f ^~:
 
其中卷积算子的定义如下：
 
对于重构而言，卷积给出了以采样点为中心的重构滤波器做过比例变换之后的实例加权和：
 
例如，对于本章开始的图，使用了三角形滤波器f(x) = max (0, 1 - |x|)。下图表示了比例变换后的三角滤波函数：
 

我们经历了看似没有缘由的复杂过程却得到一个很直观的结果：被重构成的函数f~也可以由某种插值方法得到。然而，建立好这样的背景知识后，我们可以更容易地使用傅立叶分析进行重构。
我们可以在频域中分析被采样的函数来得到对采样过程的更深入的理解。特别地，我们可以得到能从采样值精确地复原出原函数的条件----这是一个很有力的结果。在以后的讨论中，我们假定f(x)是有带宽局限(band-limited)的，即存在某个频率ω0，使得f(x)不超过大于ω0的频率。根据定义，有局限带宽的函数在频率空间中有紧支撑(compact support)的表示，使得对于所有的|ω| > ω0, 有F(ω) = 0。

在傅立叶分析中，一个重要的思想是两个函数乘积的傅立叶变换可以用两个函数的各自的傅立叶变换的卷积来表示：
 
类似地，在空间域中的卷积等价于频域中的乘积：
 
在关于傅立叶分析的标准参考文献中都有对这些性质的推导。利用这些性质我们可以知道，在空间域中采样这一步骤里，我们可以得到shah函数和原函数的乘积，这一步骤也可以描述为在频域中的F(ω)和另一个shah函数的卷积。

从表7.1中我们可以得知shah函数III _T(x)的频谱；周期为T的shah函数的傅立叶变换是另一个以1/T为周期的shah函数。我们要记住这一重要的倒数关系：这意味着如果在空间域中采样点相距很远，那么它们在频域中就相距很近。

所以，采样信号的频域中的表示就是F(ω)和这个新shah函数的卷积。求一个函数跟一个delta函数的卷积就得到这个函数的一个拷贝，所以，跟一个shah函数求卷积就得到原函数的拷贝的无限序列，其间距就是shah的周期。下图就是一个例子：
 
我们得到了这个函数的拷贝的无限序列的频谱以后，怎么来重构原函数呢？答案很明显：只需要留下位于中心的拷贝，丢弃其它所有的拷贝就可以了。为了丢掉除了中心拷贝以外的其它所有拷贝，我们要乘上一个具有适当宽度的方盒函数。宽度为T的方盒函数Π _T(x)定义如下：

        Π _T(x) = 1 /(2T),     当|x| < T
        Π _T(x) = 0,               当|x| >= T

（如图）这个乘法相当于在空间域中跟重构滤波器做卷积。这是理想状态下的采样和重构过程，总结如下：
      
这是一个很有意义的结果：我们仅仅靠一组等距采样值，就可以确定f(x)的精确的频率空间中的表达。

同样地，在空间域中使用相同的过程也以精确地恢复f(x)。因为方盒函数的傅立叶逆变换是sinc函数，在空间域中的理想重构函数为：
      
不幸地是，因为sinc函数的取值范围是无限的，为了计算f ^~(x)在空间域中的某个值就必须用到所有的采样值f(i)。在具体实现中，我们要采用有限空间范围的滤波器，虽然这样做就无法构造出完美的原函数。

在计算机图形学中，常用的方法是利用方盒函数进行重构，在x附近的某个区域取其中所有采样值的平均值。这是一个很糟糕的选择，考虑一下方盒滤波器在频域中的表现就可以看出这一点：这个技术试图利用“乘以一个sinc函数”来孤立出函数频谱的中心拷贝，这样做不仅没有很好地选择出中心拷贝，还把其它的无限个拷贝的高频贡献也包括进来了。

7.1.3 走样

除了sinc函数的无限取值范围的问题，另一个严重的实际问题是理想的采样和重构过程假定信号是有带宽局限的。对于没有带宽局限的信号、或者使用了采样速率不够高的那些信号而言，上述过程所重构的函数就跟原信号不同。

重构的成功关键在于：能够对被采样的频谱乘以适当宽度的方盒函数而精确地得到原频谱F(ω)。注意在图7.6中，信号的频谱拷贝是被空的间隔分离开来的，这样使得完美的重构成为可能。考虑一下用低一些的采样速率对原函数采样的情况。回忆一下：周期为T的shah函数的傅立叶变换是一个周期为1/T的新的shah函数。这意味着当在空间域中的采样间隔增大时，频域中的采样间隔就减小，使得F(ω)的拷贝更加彼此靠近。当它们靠得太近了，就出现重叠。
 
由于拷贝出现叠加现象，所得到的频谱就不再像是原频谱的多个拷贝了。当这个新频谱乘上一个方盒函数，结果是一个跟原F(ω)相似但不相等的频谱：在原信号中的高频率的细节被泄漏到重构信号中的低频区域中。这个新的低频的误差就被称为“走样”（因为高频信息被视为低频的了），下图显示了对一维函数f(x) = 1 + cos(4x ²)进行低速率采样而出现的重构的走样现象：
 
    （图7.9）

一个对重叠频谱问题的可能解决方案是增加采样速率，直到频谱拷贝相距足够远，没有了重叠，也就摆脱了走样。实际上，采样定理告诉我们需要多大的采样速率。这个定理说只要均匀采样点ωs的频率大于信号中最高频率ω0的两倍，就可以利用采样值完美地重构出原信号来。这个最小的采样频率被称为Nyquist（奈奎斯特）频率。

对于那些没有带宽局限的信号(ω0为无穷大)而言，就无法通过高速率采样完成完美重构。无带宽限制的信号有无限支撑的频谱，所以不管它们的频谱拷贝相距多远，总会产生重叠。不幸地是，计算机图形学中所用的函数很少是有带宽局限的。特别是任何不连续的函数都不是带宽局限的，因此也无法进行完美的采样和重构。我们可以理解这一点，因为采样点总会跳过不连续点，因此采样无法提供不连续性的位置。因此，除了增加采样速率以外，还必须使用其它不同的方法来减少因走样而带来的渲染图像的失真。

7.1.4 反走样技术

如果我们对采样和重构不够仔细，最终的图像就会出现奇怪的人为的误差。有时把因采样而产生的误差跟因重构产生的误差区分开来是很有用处的，我们称前者为前走样(prealiasing)，后者为后走样(postaliasing)。消除这些走样的技术统称“反走样”。这一节将简介几个反走样技术。

非均匀采样

虽然我们知道图像函数有无限大的频率分量，故不能从采样点得到完美的重构，但是我们有可能通过非均匀地采样来减少走样所产生的视觉影响。令ξ为在0到1之间的随机数，一组基于脉冲序列(impulse train)的非均匀采样为：
      
对于固定的采样速率而言，这并不足以捕捉到整个函数信息，无论是均匀采样还是非均匀采样都会产生不正确的重构信号。但是非均匀采样可以把走样产生的误差转换成立了噪声，对人的视觉系统的不良效果也就降低了。在频率空间中，被采样信号的拷贝被随机地移位了，所以重构后的结果是随机的失真而不是一致性的走样。

自适应的采样

另一种要推荐的反走样技术是自适应超采样(adaptive supersampling)。如果我们识别出高于Nyquist频率的信号区域，就可以只在这些区域增加采样，而无需对所有区域增加采样速率而招致不必要的计算开销。不幸地是，在实际应用中这个方法很难得到有效的利用，因为找到需要超采样的区域是很困难的。大多数的技术要检查相邻采样值，找到那些变化值很大的地方，并假设这个区域有高频的信号。

一般情况下，相邻的采样值无法确定地告诉我们真实的情况：即使这些值是相同的，函数在它们中的变化可能仍然很大。还有，在没有任何走样的情况下，相邻采样值仍然可以有很大的变化。例如，在第11章中的纹理滤波算法要很艰难地对付表面上因图像映射和过程纹理而产生的走样；我们并不想在纹理值变化很大但并无极高频率信号的地方进行自适应采样。

自适应方法总会错过一些需要超采样的区域，使得增加采样速率成为唯一的选择。自适应反走样技术很擅长于把走样非常严重的图像转换成走样不太严重的图像，但是对于复杂的场景而言，它跟只增加采样速率的方法比较，并不是更有效率地产生无视觉瑕疵的图像。

前置滤波

采样理论提供的另一项去除走样的方法是对原函数滤波，使得在所使用的采样速率下不再存在无法精确捕捉到的高频信号。这种方法将在第11章中的纹理函数中用到。虽然这种方法去掉了被采样函数中的信息，从而改变了函数的特性，但是比起走样而言，还是更令人接受些。

回忆一下，我们用一个方盒滤波器乘以原函数的频谱，我们选择方盒的宽度，使得高于Nyquist极限的频率刚好被去掉。在空间域中，这相当于用sinc滤波器对原函数做卷积。

在实际应用中，我们可以用一个效果良好的具有有限空间范围的滤波器。这个滤波器在频率空间的表示有助于理解它对理想的sinc滤波器的模拟情况。

图7.10表示了函数1+cos(4x ²)跟一个具有有限空间范围的sinc函数的变种（7.6节有介绍）进行卷积的情况，这个函数可以用图7.9中的采样速率进行采样和重构而不会出现走样。
 

7.1.5 图像合成

我们可以把前面所介绍的理论直截了当地应用到所渲染的二维图像的采样和重构中：我们有一个图像，可以把它想象为一个把二维图像位置(x,y)映射为辐射亮度L的函数：

            f(x,y) ->L.

好消息是我们可以利用我们的光线跟踪器在任何(x,y)点上对函数求值。坏消息是在一般情况下不可能对f进行前置滤波，在采样之前去掉高频信号。因此，本章所介绍的采样器会使用两种策略：第一个是增加采样速率，第二个是进行非均匀分布的采样，把走样转换成噪声。

为了方便，我们把场景函数的定义一般化为更高维的函数，使它也依赖于时间t和所采样的镜头位置(u,v)。因为相机光线是基于这5个参数的，改变其中任何一个变量就会得到一条不同的光线，从而f函数值也会有所不同。对于一个特定的图像位置，辐射亮度会随时间的变化（场景中有移动的物体）和镜头上的位置的变化（相机镜头的光圈有限）而变化。

对于更一般的情况，例如第16章和第17章定义的积分器使用了统计技术来估算沿光线的辐射亮度，它们可能会对用一条光线返回不同的辐射亮度值。如果我们把场景辐射亮度函数做进一步地扩充，把积分器所使用的采样值也包括进去（例如用于在面光源上选择点的那些值），我们就得到一个更高维的图像函数：

            f(x,y,t,u,v,i1, i2, ...) ->L

对所有这些维度进行采样是有效地生成高质量图像的重要一步。例如，如果我们确保(x,y)附近上的图像位置上有不同的镜头(u,v)位置，所得到的图像结果就会有少一些的失真，因为每个采样比其邻近的采样更能够代表场景中的信息。下面几节中的Sampler类将尽可能地解决所有这些维度上的采样问题。

7.1.6 渲染过程中的走样的原因

在渲染图像过程中，几何形状是最容易产生走样的原因之一。当把物体投影到图像平面中，物体的边界就产生一个阶梯函数：图像的函数值从一个值马上跳到另一个值。阶梯函数不仅有前面所述的无限频率内容，更糟糕的是，对走样的采样使用完美的重构滤波器时会出现人为的误差：即在被重构出的函数中出现振铃现象(ringing artifact)，又称Gibbs现象。下面显示了一个1维函数的这种现象。选择对付走样的重构滤波器需要科学、艺术和个人风格的一种组合，我们会在本章的后面见到这一点。

 

场景中非常小的物体也会产生几何走样。如果几何体过小，就会落在图像平面上的采样点中间，在动画的多个帧中，就会无法预测地时隐时现。

另一个走样的根源是物体上的纹理和材料。没有正确地滤波的纹理映射，或者光亮表面上细小的高光也会产生明暗处理走样(shading aliasing)。如果采样速率对这些特征而言不够高，就会出现走样现象。还有，一个物体的很深的阴影也会在最终的图像中产生阶梯函数。虽然有可能在图像平面上的几何边上确定阶梯函数的位置，但是在阴影边界上探测出阶梯函数要困难得多。
一个关于走样的关键的洞见是，我们绝不可能去掉所有这些走样的根源，所以我们必须开发降低走样对图像结果的影响的技术。

7.1.1 理解像素

在读本章时要注意关于像素的两点知识：首先，必须牢记图像上的像素是在图像平面上的离散点上的点采样(point samples)，跟“区域”没有任何干系。正如Alvy Ray Smith(1995)所强调的，把像素看成一个带有限面积的小方块是导致一系列错误的思维模型。本章采用了信号处理的方法来介绍这个论题，就为更精确的思维模型打下了基础。

第二个问题是，虽然在最终的图像上像素在像素阵列中很自然地用离散的整数(x,y)来定义，但是本章中的Sampler使用了连续的浮点数位置(x,y)。这两种表达方式的一个很自然的转换方法就是把浮点数坐标取整为最近的离散坐标。这个方法颇有吸引力，因为它把跟离散坐标相同的浮点坐标转换为这个离散坐标，但是，对于在[x0, x1]上的离散坐标集合而言，覆盖这个范围的连续坐标集合的范围是[x0-.5, x1+.5)。因此，对于一个给定的离散像素范围而言，所生成的连续采样位置都有0.5的偏移量。在编写程序当中我们就很容易地忘记这些，从而导致不易发现的错误。

如果我们用下列公式把连续坐标c变换为离散坐标d：
            d = floor(c)

然后用下列公式把离散值变换为连续值：

            c = d + .5

那么离散范围[x0,x1]所对应的连续范围很自然地就是[x0, x1+1)， 程序代码就简单很多(Heckbert 1990a)。我们在pbrt中采用这个约定，如图所示：