第14章蒙特卡罗积分I

第14章蒙特卡罗积分I：基本概念在介绍计算沿光线到达相机的辐射亮度的SurfaceIntergrator类和VolumeIntegrator类之前，我们先做一些关于求解散射积分方程的基本工作。这些积分方程通常没有解析解，所以我们必须求助于数值方法。虽然诸如梯形积分法和高斯积分法这类标准的数值积分技术在求低维光滑积分时很有效，但它们对于在渲染过程中所遇到的高维不连续积分而言，其收敛速度还是太慢了。蒙特卡罗积分技术为这个问题提供了一个解决方案。这类方法使用了随机性做积分求值，且收敛速度跟被积函数的维的个数无关。在本章中，我们将回顾一些概率论中的重要概念，并为使用蒙特卡罗解决渲染过程中的积分问题打下一个基础。下一章将描述提高收敛速度的计算，并描述pbrt所用到的技术。对随机性的明智运用使算法设计领域发生了革命性的变化。随机性算法大致分两大类：拉斯维加斯算法和蒙特卡罗算法。拉斯维加斯算法使用了随机性但最终会产生相同的结果（如在快速排序算法中随机地选择一个数组项作为主元素）。而蒙特卡罗算法依赖于在算法过程中所使用的不同的随机数列而得到不同的结果，但结果大致上（在平均意义上）是正确的。所以，对蒙特卡罗算法的多次运行结果做平均，就有可能找到跟正确答案很接近的统计意义上的结果。蒙特卡罗积分是一项用随机采样来估算积分值的技术。蒙特卡罗积分的一个非常有用的特性是：为了估算积分值，我们只需保证能够在定义域中的任意点对被积函数求值。这个特性不仅使得蒙特卡罗技术易于实现，而且可以应用于非常广泛的被积函数类型，包括那些不连续函数。在渲染过程中所遇到的积分是很难或不可能直接求积分的。例如，为了计算表面上某一个点上的反射光的量，我们要用方程(5.6)来对入射辐射亮度和BSDF的乘积做单位球面上的积分。因为物体的可见性使得在复杂场景中的入射辐射亮度函数会有难以预测的变化，对于这个乘积中的所有项就很难找到封闭形式的表达式，即使表达式存在，对之进行解析方式的积分也是通常不可能的。因为有了蒙特卡罗积分，计算反射辐射亮度就成为可能，因为我们只需选择球面上的一组方向，计算相应的入射辐射亮度，分别乘以相应方向上的BSDF值，再用上一个权值项。这样一来就可以处理任意的BSDF、光源和几何体；所要做的只是在任意点上对这些函数求值。蒙特卡罗积分的主要缺点是，如果用了n个采样点来估算积分，算法收敛到正确解的速度是O(n-1/2)。换句话来说，如果将误差减半，就要使用4倍的采样点。在渲染过程中，为了计算被积函数的值，每个采样点需要追踪一条或多条光线，当用蒙特卡罗积分做图像合成时，这真是很痛苦的代价。在图像中，蒙特卡罗采样所产生的人为缺陷会以噪声的形式出现，即像素随机性地出现太黑或太亮的情况。当前许多关于蒙特卡罗积分的研究着重于如何尽可能地减少误差又要减少采样点的数量。14.4背景和概率论知识的回顾我们先定义一些基本的术语并回顾一下概率论的基本知识。我们先假定读者已经熟悉了概率论的基本概念；需要对这一领域做近一步的全面了解的读者可以参考相应的教科书，例如SheldonRoss的IntroductiontoProbabilityModels(2002)。一个随机变量X是在某个随机过程中所选定的一个值。我们通用用大写字母来代表随机变量，特殊情况下我们使用希腊字母代码一些特殊的随机变量。随机变量总是在某些定义域中抽取出来的，它既可以是离散的（例如一个固定的概率集合），也可以是连续的（如实数集R）。对随机变量X使用函数f就得到另一个新的随机变量Y=f(X)。例如，随机掷一个骰子的结果就是一个离散的随机变量，它被采样于一个事件集合Xi={1,2,3,4,5,6｝。每个事件有一个概率pi=1/6，这样概率总和Σpi一定是1。我们可以取一个连续的均匀分布的随机变量ξ∈[0,1]，并将之映射为一个随机变量，如果下式满足，则取值Xi：Σj=1,i-1Pjξ≤Σj=1,iPj注意所有的Pi之和为1。一个随机变量的累积分布函数（CDF）P(x)是变量分布中的一个值小于或等于某个值x的概率。P(x)=Pr{X≤x}以掷骰子为例，P(2)=1/3，因为小于或等于2的情况占总数6个中的2个。14.1.1连续随机变量在渲染过程中，跟离散随机变量相比，会更多地用到连续随机变量，这些变量在一个连续的定义域内取值（如实数域，或单位球上的方向）。其中一个非常重要的随机变量是典型均匀随机变量，记为ξ。该变量之所以重要有两个原因。首先，在软件中生成这个分布的变量很容易，因为大多数运行库中的伪随机数发生器就可以做到这一点。第二，我们可以先生成典型均匀随机变量再施加上某种变换来得到任意的随机分布。前面所提到的将ξ映射到骰子中6个面就是这个技术在离散情况中的应用。我们来看另一个定义在[0,2]上的连续随机变量的例子：它在值x的概率值跟2-x成正比：即在0附近取值的可能性是在1处附近取值时的两倍。概率密度函数(PDF)对此有正规的表述：它描述了随机变量在某个值上取值的相对概率。PDFp(x)是CDF的导数：p(x)=dP(x)/dx均匀随机变量的p(x)是个常数。对于ξ，我们有：p(x)=1（当x∈[0,1])或0(其它情况）PDF必须是非负的，并在其定义域上的积分是1。给定一个定义域中的一个区间[a,b]，PDF给出了随机变量处于该区间的概率：P(x∈[a,b])=∫a,bp(x)dx上式可以根据微积分第一基本定理和PDF的定义得出。14.1.2期望值和方差函数f的期望值Ep[f(x)]被定义成该函数在其定义域的上某种p(x)分布的平均值。在下一节中我们将会看到如何用蒙特卡罗积分计算任意积分的期望值的。在域D中的期望值定义如下：Ep[f(x)]=∫Df(x)p(x)dx举个例子，考虑一下如何求余弦函数在[0,π]上的期望值，其中p是均匀的。因为PDFp(x)在定义域上的积分必须为1，所以p(x)=1/π。故有：E[cosx]=∫0,πp(x)dx=(-sinπ+sin0)/π=0这正是所期望得到的值，考虑一下cosx在[0,π]的曲线就明白了。一个函数的方差是函数值跟其期望值的误差函数。方差是用来量化蒙特卡罗算法估算值的误差的基本概念。它提供了误差量化的精确方法，并可以测量蒙特卡罗算法为减少误差所做的改进程度。第15章的大部分内容将用来介绍减少方差从而改进pbrt的计算结果的技术。一个函数f的方差定义为：V[f(x)]=E[(f(x)-E[f(x)]2]从期望值和方差定义中，我们马上可以看到它们有三个重要的性质。E[af(x)]=aE[f(x)]E[Σif(Xi)]=ΣiE[f(Xi)]V[af(x)]=a2V[f(x)]利用这些性质做些简单的代数推导，可以得到方差的一个更简单的表达式：V[f(x)]=E[(f(x))2]-E[f(x)]2这样，方差只是函数平方的期望值减去期望值的平方。如果随机变量是相互独立的，那么方差还有这样的性质，即方差的和等于和的方差：ΣiV[f(Xi)]=V[Σif(Xi)]14.2MonteCarlo估计量现在我们可以定义基本的MonteCarlo估计量了，它可以近似地估算任意积分的值。这是第16，17章中光传输算法的基础。假定我们想要对一个一维积分∫a,bf(x)dx求值。对于一个均匀分布的随机变量Xi∈[a,b]，MonteCarlo估计量用下式给出期望值：FN=(b-a)/NΣ1,nf(Xi)E(FN)就是这个积分值。这可以用几步推导来得出。注意随机变量Xi所对应的PDFp(x)必须等于1/(b-a)，因为p必须是常量，并且在[a,b]上的积分为1。推导如下：我们做一些一般化推广，就可以去掉随机变量必须是均匀的限制。这是极其重要的一步，因为对生成采样的PDF的选取是减少蒙特卡罗方差（第15.4节）的重要技术。如果随机变量是某个任意PDFp(x)来抽取的，那么下面的估计量可用来估算积分：FN=1/NΣ1,Nf(Xi)/p(Xi)p(x)的唯一的限制是它必须在所有满足F(x)0的x为非零。类似地，我们可以看出这个估计量的期望值计算所要求的积分：我们可以直接将这个估计量推广到高维或复杂的积分域中。从多维PDF生成的N个采样也可以照样使用这个估计量。例如，考虑一个三维积分：∫xo,x1∫yo,y1∫yo,y1f(x,y,z)dxdydz如果采样Xi=(xi,yi,zi)从三维盒(x0,y0,z0)-(x1,y1,z1)中均匀地取出，则PDFp(X)就是常量：1/(x1-x0)(y1-y0)(z1-z0)估计量为：(x1-x0)(y1-y0)(z1-z0)/NΣf(Xi)注意采样个数N是任意选取的，跟被积函数的维数无关。这是另一个MonteCarlo对传统积分的优势。蒙特卡罗所使用的采样个数完全独立于积分维数，而标准的数值积分计算所需要的采样个数跟维数呈指数级增长。只指出蒙特卡罗估计量可以收敛于正确结果还不够，因为良好的收敛速度也是很重要的。我们这里不做推导，只指出其误差按照O(N1/2)的速度收敛(N为采样个数)。虽然标准积分技术在一维的情况下的收敛速度要高于这个速度，但随着被积函数的维数增加，其性能会有指数级的下降，而蒙特卡罗的收敛速度跟维数无关，这样蒙特卡罗就成为唯一可行的高维积分的数值积分技术。我们已经遇到一些高维积分的例子，在第16章的路径跟踪公式实际上是个无限维的积分！14.3随机变量的采样为了对蒙特卡罗估计量求值，必须能够从所选定的概率分布中抽取随机采样。本节介绍这个过程的基本知识，并举一些简单的例子。下一节将介绍一般性的多维情况，下一章的第15.5节和15.6节将利用这些技术在BSDF和光源的分布中生成采样。14.3.1逆转法逆转法使用一个或多个均匀随机变量，并将之映射到所需分布的随机变量中。为了解释这个过程，我们先看一个简单的离散的例子我们有一个4种结果的过程，每个结果的概率分别为p1,p2,p3,p4，并且p1+p2+p3+p4=1。相应的PDF如下：为了在这个分布中取样，我们先求CDFP(x)。在连续的情况下，P是p的不定积分。在离散的情况下，我们可以直接从左边开始，将相应的直方条叠在一起，如图：（注意最右的一叠的高度为1，因为所有的概率之和为1）。为了在分布中抽样，我们使用一个均匀随机数ξ，用它通过CDF来选择其中一个可能的结果，并且要求该结果的概率等于它自己的概率。如图，事件的概率被投影到垂直轴上，随机变量ξ从中选取概率值。应该清楚这个抽样的分布是正确的，即均匀采样的概率恰好等于它所碰到的特定直方条的高度。为了将这个技术推广到一般的连续分布上，我们可以考虑当离散的概率个数趋于无穷大的情况。则上面的PDF变成一个光滑的曲线，而CDF就是它的积分。上面所表示的投影过程仍然一样，这个投影有一个很方便的数学解释，即它表示求CDF的反函数，并在ξ处求反函数的值。这个技术被称为逆转法。更精确地讲，我们可以从一个任意的PDFp(x)中按照下列步骤抽取采样Xi:1.计算CDFP(x)=∫0,xf(x')dx'。2.求反函数P-1(x)。3.取一个均匀分布的随机数ξ。4.计算Xi=P-1(ξ)。14.3.2举例：幂分布作为一个该过程的一个例子，我们考虑从一个幂分布抽样的过程。我们在对Blinn微平面模型采样时就是这种情况。幂分布的PDF为：p(x)=cxn第一个任务是求PDF。在大多数情况下，这值涉及到求比例常数c。我们可以利用∫p(x)dx=1就可以得到：所以，p(x)=(n+1)xn。我们对之积分，就可以得到CDF:P(x)=∫0,xf(x')dx'=xn+1求反函数很简单：P-1(x)=x(1/n+1)。这样，给定一个均匀随机变量ξ，从幂分布抽样所得到的结果是：X=ξ(1/n+1)14.3.3举例：指数分布当我们渲染有参与介质的图像时，经常要从一个指数分布中抽样。第一步是对该分布进行归一化，使其积分值为1。在这种情况下，我们希望所生成的采样能够覆盖[0,infnite)，而不是[0,1]，故有：所以c=a，而PDF为p(x)=ae-ax。对之积分得到P(x):该函数的反函