第15章蒙特卡罗积分II

第15章蒙特卡罗积分II：提高效率蒙特卡罗光线追踪算法里的方差以图像噪声的形式体现出来。如何竭力地减小方差竟成了大多数优化蒙特卡罗算法的基本任务。回忆一下蒙特卡罗算法的收敛速度，这意味着为了对方差减半，就需要四倍的采样。由于估算例程的运行时间跟采样数目成比例关系，所以减少方差的花销可以是很高的。本章将介绍可以改进蒙特卡罗算法的效率而并不增加采样的理论和技术。一个估计量公式F的效率可以定义为：ε[F]=1/V[F]T[F]其中V[F]是它的方差，T[F],T[F]是求值所用的运行时间。根据这个度量方法，如果估计量F1跟另一个估计量F2相比，花费更少的时间而产生相同的方差，或者在相同时间内生成更小的方差，则称它更有效率。其中一项技术被认为是改进渲染效率的最为有效的技术，即重要性采样(importancesampling)。回忆一下，在蒙特卡罗估计量的实现中，我们有一定的自由来选择采样分布函数p(x)。如果我们选择跟被积函数很相像的采样分布，那么就可以减小方差。这项技术之所以被称为重要性采样，是因为我们趋向于在函数定义域的“重要”部分进行采样，而这些部分有相对大的函数值。本章的前半部分讨论重要性采样和几项其它改进蒙特卡罗算法效率的技术。本章的后半部分讨论根据BSDF、光源、体积散射来生成采样的技术，而这些技术可用做重要性采样中的采样分布。第16，17章将使用这些采样方法。15.1俄罗斯轮盘法和分割法俄罗斯轮盘法和分割法是两项可以改进蒙特卡罗算法效率的相关方法，它们的做法是要增加对结果有较大贡献值的采样的可能性。有些求值很耗时的采样却对最终结果的贡献值很小，俄罗斯轮盘法要解决就是这个问题。我们希望避免对这些采样进行求值却仍可以计算出正确的估计量。最为一个例子，我们考虑对直接光照积分的估算问题，该积分给出了场景中从光源直接照射到某点而产生的反射辐射亮度：假定我们从某个分布p(ω)取N=2个采样来计算估计量：其中累加式中绝大部分的计算开销来自于追踪起始于点p(看光源是否被遮挡)的阴影光线。对于那些被积函数值明显为0的那些方向（例如，对于fr(p,ωo,ωi)为0的那些方向）而言，我们可以跳过追踪阴影光线的工作，因为是否去追踪并不影响最终的计算结果。使用俄罗斯轮盘法可以让我们也跳过那些被积函数值很低但不一定为0的那些方向的光线追踪，而仍然大致可以计算出正确的结果。例如，在那些fr(p,ωo,ωi)很小或者ωi靠近水平面使得|cosθi|很小的方向上，我们可以避免进行光线追踪。当然，我们不能完全忽略这些采样，否则估算公式就会低估结果值。为了使用俄罗斯轮盘法，我们先要选定某个终止概率q。我们几乎可以用任意的方式来选择这个概率值；例如，我们可以根据一个选定的采样来估算出被积函数的值，并以此为基础确定一个概率值，当被积函数值减小时再对之增加。在使用概率q时，就不对被积函数求值，而是用一个常数c代替（一般用c=0）。在使用概率1-q时，就要对被积函数求值，但要乘以一个权值1/(1-q)，用来弥补那些被跳过的采样：F'=(F-qc)/(1-q),ξqF'=c,otherwise要注意的是它的期望值跟原估计量的期望值相同：E[F']=(1-q)((E[F]-qc)/(1-q))+qc=E[F]俄罗斯轮盘法绝不能减少方差。实际上，除非c=F，否则它总是加大方差的。然而，如果概率选择得当的话，就可以跳过那些对最终结果贡献很小的采样，从而提高了效率。其中一个陷阱是，如果俄罗斯轮盘法选择了一个很差的权值，就会大大地增加方差。想象一下我们用0.99做为终止概率对所有的相机光线使用俄罗斯轮盘法，我们即只追踪其中的1%的相机光线，而用1/0.1=100的权值。在严格的数学意义下结果仍然“正确”，但最终的视觉效果无比糟糕：大多数像素是黑色的，只有很少明亮的像素。15.1.1分割法俄罗斯轮盘法可以减少对次要的采样的求值，而分割法却以增加采样数的方式来提供效率。我们再考虑一下计算直接光照下的反射辐射亮度的问题。在忽略像素滤波的情况下，这个问题可以写成一个像素的面积A下以及方向球面δ2下的双重积分，其中Ld(x,y,ω)是图像坐标(x,y)处由方向ω上入射辐射亮度产生的出射辐射亮度：估算该积分的一个自然方式就是生成N个采样，每个采样包含一个(x,y)图像位置和一个朝向光源的方向ω，然后使用蒙特卡罗估计值公式。如果场景中有许多光源，或者有一个投射软阴影的面光源，为了使结果有一个可以接受的方差水平，就需要数十个或数百个采样。不幸地是，每个采样需要两条要被追踪的光线：第一条用来计算从图像平面(x,y)处看到的第一个可见表面，第二条是沿方向ω到光源的阴影光线。这个方法的问题是，如果使用N=100个采样来估算这个积分，那么就要追踪100条相机光线和100条阴影光线，而且100条相机光线对于像素反走样运算是过多了，因此对最终结果的方差的减少并没有什么贡献。分割法为了解决这个问题，将这个方法做了正规化处理：即在某些维上做了采样以后，对每一个采样再在其它维上做多个采样。利用分割法，这个积分的估计值公式可以写成（使用N个图像采样，每个图像采样使用M个光源采样）：这样一来，我们只使用5个图像采样，但对每个图像采样使用20个光源采样，就可以只追踪105条光线（而不是200条），我们仍旧使用100个面光源采样来得到高质量的软阴影效果。14.2精细的采样布置为了更好地扑捉到“重要的”被积函数的特征，有一类经典且有效的减少方差的技术是基于对采样的精细布置上的。这类技术是对重要性采样等技术的有益补充，并广泛地应用于pbrt中。实际上，第7章中的Sampler的一项任务就是为此目的而生成具有良好分布的采样。14.2.1分层采样我们在第7.3节介绍了分层采样，现在我们有了使用这项技术的工具。分层采样将积分定义域Λ分成n个互不重叠的区域Λ1,Λ2,...,Λn。每个区域称为一个层，它们的并集完全覆盖了Λ。为了从Λ中进行采样，我们将按照每个层中的密度pi来在Λi中抽取ni个采样。像素周围的区域被分成一个kxk网格，我们利用分层采样的方法在每个网格单元抽取一个采样。这样做要好于随机地抽取k2个采样，因为采样位置不会拥挤到一起。下面我们将表明为什么这种技术会减小方差。在一个层Λi中，蒙特卡罗估计量为：其中Xi,j是从密度pi中抽取的第j个采样。总体的估算值为F=Σi=1,nviFi，而vi是层i的体积所占的比例(0vi=1)。被积函数在层i的真值为：μi=E[f(Xi,j)]=1/vi∫Λif(x)dx故得在这个层中的方差：σi2=1/vi∫Λi(f(x)-μi)dx因此，由于层中的采样有ni个，故层的估计量的方差为σi2/ni。这表明总估计量的方差为：如果我们做这样一个比较合理的假定：采样个数ni正比于体积vi，则有ni=viN，那么总体估计量的方差公式为：V[FN]=1/NΣviσi2为了将这个结果跟没有分层的方差相比较，我们注意到选择一个没有分层的采样等价于这样的做法：即按照由体积vi定义的离散的概率分布选择一个分层I，然后在层Λi中选择一个随机采样X。在这种意义下，X是在I的条件下选择的，所以我们用条件概率公式有：V[F]=ExViF+VxEiF=1/N[Σviσi2+1/NΣvi(μi-Q)]这里Q是f在整个定义域上的平均值。对这个结果的推导见Veach(1997)。这个表达式有两个值得注意的地方。首先，第二项累加和必须是非负的，因为方差总是非负的。第二，它显示出分层采样绝不会增加方差。实际上，除非第二项的和为0，分层法总是可以减小方差的。只有当f在每个层有相同的平均值时，该式才为0。事实上，为了使分层采样效果达到最佳，我们希望第二项的和达到最大化，这样就使得各层上的平均值尽可能地不均等。如果我们不知道f的任何情况，就需要更“紧”(compact)的分层。如果分层很宽，就会有更大的方差，就会有更靠近真正平均值Q的μi。分层采样的主要缺点是它也受制于“维数灾难”(curseofdimensionality)。在D维空间中每个维上做S个分层，就会需要SD个采样，随着维数的增加，就很快地变为不可行了。幸运地是，我们常常可以只对某些维做独立的分层，然后随机地将来自不同维的采样联系起来，正如我们在第7.3节中所做的那样。我们应当这样选择需要分层的维：这些维在被积函数值的效果上要尽可能地高度相关联。例如，对于第15.1.1节上的直接光照的例子，对像素位置(x,y)和光线方向(θ,Φ)分层要比对(x,θ)和(y,Φ)要高效得多。另一个解决“维数灾难”的方法是前面介绍过的拉丁超立方体采样。不幸地是，拉丁超立方体采样在减少方差方面并不如分层采样那样高效。然而，拉丁超立方体采样的效果不比均匀随机采样差，而且经常要好上许多。15.2.2准蒙特卡罗算法在第7.4节介绍过的低差异采样计算是一类被称为“准蒙特卡罗”的算法的基础。准蒙特卡罗算法的核心部分就是将标准蒙特卡罗算法中所用的随机数用低差异点集来代替，而这些低差异点是由精细设计的确定性算法生成的。该方法的优点是，对于许多积分问题，准蒙特卡罗算法比基于标准蒙特卡罗的方法有更快的渐近收敛速度。许多用于常规蒙特卡罗算法的技术也可以很好地应用到准随机采样点上，但有些技术并非如此（例如俄罗斯轮盘法和采样拒绝法）。虽然渐近收敛速度通常无法适用于不连续的被积函数（因为收敛速率依赖于被积函数的光滑性），但准蒙特卡罗算法对这些被积函数而言仍然比常规蒙特卡罗算法要好一些，详见Keller2001。在pbrt中，我们忽略了这两种方法的差别，并在第7章中对Sampler做了适应性上的改动。这会引起细微错误的可能性，但只要Integrator使用这两类采样都适用的技术（如重要性采样或其它技术），就不会产生问题。15.2.3失真的采样和变形当用分层采样或低差异采样来处理诸如在光源采样等问题时，pbrt生成[0,1]x[0,1]上的一组采样(ξ1,ξ2)，然后用第14.4和14.5节中的变换算法将这些采样映射为光源上的点。我们期望在这个过程中[0,1]x[0,1]上的采样被变换后仍然保持分层性质---换句话说，邻近的采样在光源的表面上应该还是保持邻近关系，而彼此远离的采样在光源上应该还是彼此远离的。如果映射无法保留住这个性质，那么分层的好处就被丧失了。例如，这可以解释为什么Shirley的“正方形到圆”的映射要好于直接映射，因为直接映射在远离圆心的地方更有宽松些的分层。这也可以解释为什么(0,2)序列法在实际应用中要比分层法好上许多：因为通过定义域之间的变换后，它们比分层法相比更容易保持住良好的分布特性。15.3偏差另一个减少方差的方法是将偏差引入计算中：有时我们知道计算一个估计值并不能跟期望值相等，利用这一点却可以达到减小方差的目的。如果一个估计值的期望值等于正确结果，则称它是无偏差的。否则，它的偏差量等于：β=E[F]-FKalos和Whitlock（1986）给出了下面的一个利用偏差的例子。考虑一下如何估算一组在区间[0,1]之间的随机变量的平均值。我们可以用下面的估计值：1/NΣi=1,nXi或者用有偏差的估计值：1/2max(X1,X2,...,XN)第一个估计值实际上是无偏差的，但方差是O(1/N-1)。第二个估计值的期望值是：0.5N/(N+1)≠0.5所以它是有偏差的，虽然它的方差是O(N-2)，所以要好得多。对于很大的N值而言，我们可选用第二个估计值。在第7.6节中的像素重构方法也可以视为一个有偏差的估计值。作为一个蒙特卡罗估算问题，我们可以估算下面的公式：I(x,y)=∫∫f(x-x',y-y')L(x',y')dx'dy'其中I(x,y)是最终的像素值，f(x,y)是像素滤波函数（假定它已经被归一化，即积分为1），L(x,y)是图像辐射亮度函数。假定我们在图像平面上进行均匀采样，所有的采样有相同的概率密度（记为pc)，那么该方程的蒙特卡罗无偏差估计值为：I(x,y)≈1/NpcΣi=1,Nf(x-xi,y-yi)L(xi,yi)这个结果跟前面我们所用的像素滤波函数