PBRT阅读：第十四章 蒙特卡罗积分I：基本概念

转自：http://www.opengpu.org/forum.php?mod=viewthread&tid=6768

第14章
蒙特卡罗积分I：基本概念

在介绍计算沿光线到达相机的辐射亮度的SurfaceIntergrator类和VolumeIntegrator类之前，我们先做一些关于求解散射积分方程的基本工作。这些积分方程通常没有解析解，所以我们必须求助于数值方法。虽然诸如梯形积分法和高斯积分法这类标准的数值积分技术在求低维光滑积分时很有效，但它们对于在渲染过程中所遇到的高维不连续积分而言，其收敛速度还是太慢了。

蒙特卡罗积分技术为这个问题提供了一个解决方案。这类方法使用了随机性做积分求值，且收敛速度跟被积函数的维的个数无关。在本章中，我们将回顾一些概率论中的重要概念，并为使用蒙特卡罗解决渲染过程中的积分问题打下一个基础。下一章将描述提高收敛速度的计算，并描述pbrt所用到的技术。

对随机性的明智运用使算法设计领域发生了革命性的变化。随机性算法大致分两大类：拉斯维加斯算法和蒙特卡罗算法。拉斯维加斯算法使用了随机性但最终会产生相同的结果（如在快速排序算法中随机地选择一个数组项作为主元素）。而蒙特卡罗算法依赖于在算法过程中所使用的不同的随机数列而得到不同的结果，但结果大致上（在平均意义上）是正确的。所以，对蒙特卡罗算法的多次运行结果做平均，就有可能找到跟正确答案很接近的统计意义上的结果。

蒙特卡罗积分是一项用随机采样来估算积分值的技术。蒙特卡罗积分的一个非常有用的特性是：为了估算积分值，我们只需保证能够在定义域中的任意点对被积函数求值。这个特性不仅使得蒙特卡罗技术易于实现，而且可以应用于非常广泛的被积函数类型，包括那些不连续函数。

在渲染过程中所遇到的积分是很难或不可能直接求积分的。例如，为了计算表面上某一个点上的反射光的量，我们要用方程(5.6)来对入射辐射亮度和BSDF的乘积做单位球面上的积分。因为物体的可见性使得在复杂场景中的入射辐射亮度函数会有难以预测的变化，对于这个乘积中的所有项就很难找到封闭形式的表达式，即使表达式存在，对之进行解析方式的积分也是通常不可能的。因为有了蒙特卡罗积分，计算反射辐射亮度就成为可能，因为我们只需选择球面上的一组方向，计算相应的入射辐射亮度，分别乘以相应方向上的BSDF值，再用上一个权值项。这样一来就可以处理任意的BSDF、光源和几何体；所要做的只是在任意点上对这些函数求值。

蒙特卡罗积分的主要缺点是，如果用了n个采样点来估算积分，算法收敛到正确解的速度是O(n-1/2)。换句话来说，如果将误差减半，就要使用4倍的采样点。在渲染过程中，为了计算被积函数的值，每个采样点需要追踪一条或多条光线，当用蒙特卡罗积分做图像合成时，这真是很痛苦的代价。在图像中，蒙特卡罗采样所产生的人为缺陷会以噪声的形式出现，即像素随机性地出现太黑或太亮的情况。当前许多关于蒙特卡罗积分的研究着重于如何尽可能地减少误差又要减少采样点的数量。

14.4 背景和概率论知识的回顾

我们先定义一些基本的术语并回顾一下概率论的基本知识。我们先假定读者已经熟悉了概率论的基本概念；需要对这一领域做近一步的全面了解的读者可以参考相应的教科书，例如Sheldon Ross的Introduction to Probability Models(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-1 Pj <
ξ≤ Σj=1,i Pj

注意所有的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)对此有正规的表述：它描述了随机变量在某个值上取值的相对概率。PDF p(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,b p(x) dx

上式可以根据微积分第一基本定理和PDF的定义得出。

14.1.2 期望值和方差

函数f的期望值Ep[f(x)]被定义成该函数在其定义域的上某种p(x)分布的平均值。在下一节中我们将会看到如何用蒙特卡罗积分计算任意积分的期望值的。在域D中的期望值定义如下：

Ep[f(x)] = ∫D f(x)p(x)
dx

举个例子，考虑一下如何求余弦函数在[0,π]上的期望值，其中p是均匀的。因为PDF p(x)在定义域上的积分必须为1，所以 p (x) = 1 / π。故有：

E[cosx] = ∫0, π p(x) dx = (-sin π + sin 0) / π = 0

这正是所期望得到的值，考虑一下cosx在[0,π]的曲线就明白了。

一个函数的方差是函数值跟其期望值的误差函数。方差是用来量化蒙特卡罗算法估算值的误差的基本概念。它提供了误差量化的精确方法，并可以测量蒙特卡罗算法为减少误差所做的改进程度。第15章的大部分内容将用来介绍减少方差从而改进pbrt的计算结果的技术。一个函数f的方差定义为：

V[f(x)] = E [ (f(x) - E [f(x)]2 ]

从期望值和方差定义中，我们马上可以看到它们有三个重要的性质。

E [ a f(x)] = a E [f(x)]

E [Σi f(Xi)]
= Σi E [f(Xi)]

V [ a f(x)] = a2 V [f(x)]

利用这些性质做些简单的代数推导，可以得到方差的一个更简单的表达式：

V[f(x)] = E [(f(x)) 2] - E [f(x)] 2

这样，方差只是函数平方的期望值减去期望值的平方。如果随机变量是相互独立的，那么方差还有这样的性质，即方差的和等于和的方差：

Σi V[f(Xi)]
= V [Σi f(Xi)]

14.2 Monte Carlo估计量

现在我们可以定义基本的Monte Carlo估计量了，它可以近似地估算任意积分的值。这是第16，17章中光传输算法的基础。

假定我们想要对一个一维积分∫a, b f(x) dx求值。对于一个均匀分布的随机变量Xi∈[a,b]，Monte Carlo估计量用下式给出期望值：

FN = (b - a) / N Σ1,n f(Xi)

E(FN)就是这个积分值。这可以用几步推导来得出。注意随机变量Xi所对应的PDF p(x)必须等于1/(b-a)，因为p必须是常量，并且在[a,b]上的积分为1。推导如下：

我们做一些一般化推广，就可以去掉随机变量必须是均匀的限制。这是极其重要的一步，因为对生成采样的PDF的选取是减少蒙特卡罗方差（第15.4节）的重要技术。如果随机变量是某个任意PDF p(x)来抽取的，那么下面的估计量可用来估算积分：

FN = 1/ N Σ1,Nf(Xi) / p(Xi)

p(x)的唯一的限制是它必须在所有 满足F(x) > 0的x为非零。类似地，我们可以看出这个估计量的期望值计算所要求的积分：

我们可以直接将这个估计量推广到高维或复杂的积分域中。从多维PDF生成的N个采样也可以照样使用这个估计量。例如，考虑一个三维积分：

∫xo,x1 ∫yo,y1 ∫yo,y1 f(x,y,z)
dxdydz

如果采样Xi = (xi, yi, zi)从三维盒(x0,y0,z0)-(x1,y1,z1)中均匀地取出，则PDF p(X)就是常量：

1/(x1-x0)(y1-y0)(z1-z0)

估计量为：

(x1-x0)(y1-y0)(z1-z0) / N Σ f(Xi)

注意采样个数N是任意选取的，跟被积函数的维数无关。这是另一个Monte Carlo对传统积分的优势。蒙特卡罗所使用的采样个数完全独立于积分维数，而标准的数值积分计算所需要的采样个数跟维数呈指数级增长。

只指出蒙特卡罗估计量可以收敛于正确结果还不够，因为良好的收敛速度也是很重要的。我们这里不做推导，只指出其误差按照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如下：

为了在这个分布中取样，我们先求CDF P(x)。在连续的情况下，P是p的不定积分。在离散的情况下，我们可以直接从左边开始，将相应的直方条叠在一起，如图：（注意最右的一叠的高度为1，因为所有的概率之和为1）。

为了在分布中抽样，我们使用一个均匀随机数ξ，用它通过CDF来选择其中一个可能的结果，并且要求该结果的概率等于它自己的概率。如图，事件的概率被投影到垂直轴上，随机变量ξ从中选取概率值。应该清楚这个抽样的分布是正确的，即均匀采样的概率恰好等于它所碰到的特定直方条的高度。为了将这个技术推广到一般的连续分布上，我们可以考虑当离散的概率个数趋于无穷大的情况。则上面的PDF变成一个光滑的曲线，而CDF就是它的积分。上面所表示的投影过程仍然一样，这个投影有一个很方便的数学解释，即它表示求CDF的反函数，并在ξ处求反函数的值。这个技术被称为逆转法。

更精确地讲，我们可以从一个任意的PDF p(x)中按照下列步骤抽取采样Xi:

1. 计算CDF P(x) = ∫0, x f(x') dx'。

2. 求反函数P-1(x)。

3. 取一个均匀分布的随机数ξ。

4. 计算Xi = P-1(ξ)。

14.3.2举例：幂分布

作为一个该过程的一个例子，我们考虑从一个幂分布抽样的过程。我们在对Blinn微平面模型采样时就是这种情况。幂分布的PDF为：

p(x) = c x n

第一个任务是求PDF。在大多数情况下，这值涉及到求比例常数c。我们可以利用∫p(x)dx = 1就可以得到：

所以，p(x) = (n+1)xn。我们对之积分，就可以得到CDF:

P(x) = ∫0, x f(x') dx' = x n+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):

该函数的反函数为：

P-1(x) = - ln(1-x) / a

故我们可以得到下列的抽样：

X = - ln(1- ξ) / a

我们可以对该式做近一步地简化，因为ξ是均匀分布的随机数，故1- ξ也是均匀分布的，故而我们可以将1- ξ替换成ξ而不影响分布，故有：

X = - ln(ξ) / a

14.3.4 举例：一维分段函数

另一个很有趣味的练习是在一个一维分段函数上取样。不失一般性，我们只考虑定义在[0,1]上的分段函数

假定该函数的定义域被分为N等分，每等分为Δ=1/N。这些区间的端点为xi = i Δ，其中i的范围是0到N。在每个区间中，f(x)为一个常量。如图：

f(x) = vi , 其中 xi <=x < xi+1 (i = 0, N - 1)

积分∫f(x)dx为：

c = ∫0,1f(x)dx = Σi=1,N-1 vi
/ N

我们很容易构造出f(x)的PDF p(x)， 即f(x)/c。CDF P(x)是分段的线性函数（上面图b）。

在点xi和xi+1之间，CDF按照斜率vi/c线性增加

回忆一下为了对f(x)采样，我们需要求CDF的翻函数，并找到值x，使得：

ξ = ∫0,xf(x')dx' = P(x)

由于CDF是单调递增的，故x必在xi到xi+1之间，并满足P(xi) <= ξ<= P(xi+1)

为了有效地进行采样，我们提供一个函数，它从f(x)的vi中取值，计算在xi的CDF值，并返回f(x)的积分值。

<MC Function Definitions> =

void ComputeStep1dCDF(float *f, int nSteps, float *c, float *cdf) {

<Compute integral of step function at xi>

<Transform step function integral into cdf>

}

该函数先计算f(x)的积分。它将结果存放到cdf数组，故无需再申请额外的内存。调用者必须为cdf数组申请nStep+1个浮点数的内存空间，因为如果f(x)有N个阶梯值，则需要存放CDF在每个N+1个xi处的CDF值。在数组后面存放CDF值1是冗余的，但有利于简化代码。

<Compute integral of step function at Xi>

int i;

cdf[0] = 0.;

for (i = 1; i < nSteps+1; ++i)

cdf[ i ] = cdf [ i - 1 ] + f [ i - 1 ] / nSteps;

由于在[0,1]上的积分值存放到在cdf[nSteps]中，CDF可以用这个值进行归一化：

<Transform step function integral into cdf> =

*c = cdf[nSteps];

for (i = 1; i < nSteps+1; ++i)

cdf[ i ] /= *c;

下列函数用来对函数采样：

<MC Function Definitions> +=

float SampleStepld(float *f, float *cdf, float c,

int nSteps, float u, float *pdf) {

<Find surrounding cdf segments>

<Return offset along current cdf segment>

}

首先我们需要找到覆盖ξ的一对CDF值。因为cdf数组是单调递增的，我们可以用二分查找来得到它。其中用到c++标准函数lower_bound()。

<Find surrounding cdf segments> =

float *ptr = std::lower_bound(cdf, cdf+nSteps+1, u);

int offset = (int) (ptr-cdf-1);

有了这一对CDF值，我们就可以计算x值了。我们确定ξ在cdf[offset]和cdf[offset+1]中的位置。因为CDF是线性的，x位于xi和xi+1之间。这个采样的PDF p(x)很容易计算，因为我们有归一化值c，它满足p(x) = f(x) / c。

<Return offset along current cdf segment> =

u = (u - cdf[offset]) / (cdf[offset+1] - cdf[offset]);

*pdf = f[offset] / c;

return (offset + u) / nSteps;

14.3.5 拒绝法

对于某些函数，我们无法使用积分的方法求它们的PDF，或者无法解析地得到CDF的反函数。拒绝法(rejection method)可以不做这两步也可以按照函数的分布来生成采样，这实际上是掷飞镖方法。假定我们要在函数f(x)上采样，但我们有一个PDF p(x)满足f(x) < cp(x)，c为某个常数，并且假定我们知道如何从p上采样。这时拒绝法就很简单：

loop forever:

sample X from p's direction

if ξ< f(X) / (cp(X)) then

return X

这个例程重复不断地选择一对随机变量(X, ξ)。如果点(X, ξcp(X))位于f(X)之下，则接受采样X，否则就拒绝它，并选择下一对采样。无需深入讨论，我们就可以明白该方法的效率依赖于cp(x)界定f(x)的紧密程度。该技术适用于任何维数。

在实际应用中， pbrt中的任何蒙特卡罗算法都没有使用拒绝法。我们通常需要找跟f(x)类似的可以直接进行采样的分布，下一章我们再介绍其中的原因。然而，拒绝法是一项值得关注的重要技术，特别是在调试蒙特卡罗的程序时。例如，如果我们怀疑使用逆转法抽样程序有错误，我们可以用更直接的基于拒绝法的实现代替之，再看蒙特卡罗估计量是否计算出相同的结果。当然，这要进行大量的采样。

14.3.6 举例：用拒绝法采样单位圆

假定我们想要在一个单位圆内选一个均匀分布的点。如果使用拒绝法，只需随机地在外切正方形内选择一个随机位置(x,y)，再看它是否落在圆内。如图：

下面的RejectionSampleDisk()函数实现了这个算法。我们可以用类似的方法来生成任意复杂形状之内的均匀分布采样，只要该形状有一个测试点在形状之内或之外的方法。

<Monte Carlo Function Definitions> +=

void RejectionSampleDisk(float *x, float *y)

{

float sx, sy;

do {

sx = 1.f - 2.f * RandomFloat();

sy = 1.f - 2.f * RandomFloat();

} while (sx*sx + sy*sy > 1.f);

*x = sx;

*y = sy;

}

一般地说，拒绝法的效率依赖于被拒绝的采样所占的比例。对于2D的均匀采样的情况，还是很容易计算的，即是圆面积除以外切正方形面积：π/4，约等于78.5%。对于一般的n维情况，在超球内采样的效率会随着n值的增大而降低（因为n维的超球体积趋于0）。

14.4 分布之间的变换

为了描述逆转法，我们介绍了一种按某种分布生成采样的技术，它将典型均匀分布的随机变量以特定的方式进行变换。这里我们将研究将采样从任意一种分布变换到另一种分布的更一般化的问题。

假定我们有一个从某个PDF Px(x)抽样的随机变量Xi。现在，如果我们要计算Yi = y(Xi)，我们想要找到新随机变量Yi的分布。这看上去是个很奇怪的问题，但我们将体会到理解这种变换对从多维分布函数中抽样至关重要。

函数y(x)必须是一一映射。如果多个x值被映射到同一个y值，那么就无法清楚地描述一个y值的概率密度。由于y必须是一一映射，那么它的导数必须严格大于0或严格小于0，这就意味着下式成立：

Pr {y <= y(x) } = Pr { X <= x}

所以：

Py(y) = Py(y(x))
= Px(x)

利用这个CDF之间的关系可以直接导出它们PDF之间的关系。如果我们假定y的导数大于0，通过求微分，可知：

py(y) dy/dx = px(x)

故有：

py(y) = (dy/dx) -1 px(x)

一般而言，y的导数要么严格为正，要么严格为负，它们的密度关系式如下：

py(y) = |dy/dx| -1 px(x)

怎么用这个公式呢？假定px(x) = 2x，定义域为[0,1]，令Y = sinX。随机变量Y的PDF是什么呢？因为我们知道dy/dx = cosx，

py(x) = px(x)
/ |cosx| = 2x / cosx = 2sin-1y/((1-y2)1/2)

这个例程似乎有些倒退---通常我们已知某个采样的PDF，而不是某个变换。例如，我们想从某个px(x)抽样的X，想要从某个分布py(y)计算出Y。那么我们用什么变换呢？我们只需要求它们的CDF相等，即Py(y)
= Px(x)， 那么就可以得到变换：

y(x) = Py-1(Px(x))

这是逆转法的一般化情况，因为如果X是在[0,1]均匀分布的，就有Px(x)=x，就会得到上一节的相同的例程。

14.4.1 多维情况下的变换

对于一般的n维情况，我们可以有类似的推导来得到不同密度之间的类似关系。这里不再列出推导过程了，因为它跟一维情况有相同的形式。假定我们有一个n维随机变量X，其密度函数为px(x)。另Y = T(X)， 其中T为一个双射函数。这时，密度关系式为：

py(y) = py(T(x))
= px(x) / |JT(x)|

其中|JT|是T的雅克比矩阵的行列式值的绝对值。

14.4.2 举例： 极坐标

极坐标变换公式为：

x = r cosθ

y = r sinθ

假定我们要从每个密度函数p(r, θ)抽样，相应的密度函数p(x,y)应该是怎样的？这个变换的雅克比矩阵为：

它的行列式值为r(cos2θ + sin2θ)
= r。所以p(x,y) = p(r, θ) / r。当然并非我们想要得到的结果--因为我们先从笛卡尔坐标系开始采样再将之变换到极坐标，这时我们有p(r, θ) = r p(x,y)。

14.4.3 举例：球面坐标

球面坐标变换公式为：

x = r cosθ cosΦ

y = r sinθsinΦ

z = r cosθ

这个变换的雅克比行列式值为| JT| = r2sinθ，所以相应的密度函数为：

p(r, θ, Φ) = r2sinθ p(x,y,z)

这个变换很重要，因为它可以将方向用单位球上的点(x,y,z)来表示。我们记得立体角是用单位球面上的点集面积来定义的。在球面坐标系中，我们在前面推导过：

dω = sinθ dθ dΦ

所以，如果我们有一个在立体角Ω上的密度函数，这意味着下式成立：

Pr {ω ϵ Ω} = ∫Ωp(ω)dw

我们就可以推导出关于θ 和Φ的密度函数:

p(θ, Φ) dθ d Φ = p(ω) d ω

p(θ, Φ) = sinθ p(ω)

14.5 使用多维变换的2D采样

假定我们有一个2D联合密度函数p(x,y)做(X,Y)的采样。有时多维函数是可分离的，可以被表示为多个一维密度函数的乘积，例如：

p(x,y) = px(x) . py(y)

在这种情况下，我们可以独立地从px采样X，从py中采样Y来得到(X,Y)。但是，许多有用的密度函数并不是可分离的，所以我们要介绍如何在一般情况下进行多维分布的采样理论。

给定一个2D密度函数，我们可以“积分掉”两维中的一个而得到边缘密度函数:

p(x) = ∫p(x,y)dy

这个函数可以视为X单独的密度函数。更精确的说法是，它是特定的x对所有可能的y值的平均密度。

条件密度函数p(y|x)是选定某个x的关于y的密度函数：

p(y|x) = p(x,y) / p(x)

从联合分布中进行2D采样的基本思想是，先计算出边缘密度函数，将一个变量分离出来，再用标准的一维技术利用该密度函数采样。一旦得到这个采样，再计算这个采样值所对应的条件密度函数，然后再利用标准的一维技术继续采样。

14.5.1 举例：在半球上的均匀采样

做为一个例子，我们在半球上根据立体角均匀地对方向采样。我们还记得均匀分布的意义是密度函数是一个常量，所以我们有p(ω) = c。由于密度函数在其定义域上积分必须为1，我们就得到：

故而p(ω) = 1/(2π)，或者p(θ,Φ) = sinθ/（2π）。注意这个密度函数是可分离的。然而，我们仍使用边缘密度和条件密度来演示多维采样技术。

我们先来采样θ。为此，我们需要知道θ的密度函数p(θ):

p(θ) = ∫0,2π p(θ,Φ)d Φ = ∫0,2πsin
θdΦ = sinθ

现在我们计算Φ的条件密度：

p(Φ|θ) = p(θ,Φ)/p(Φ) = 1/(2π)

注意到Φ的密度函数本身是均匀的，这可以从半球的对称性上直观地看出来。现在我们用一维逆转技术分别对这些PDF进行采样：

P(θ) = ∫0,θ p(θ')dθ' = 1 - cosθ

P(Φ|θ) =∫0, 1/2 π d Φ ' = Φ/2 π

我们可以很方便地得到它们的反函数，注意到ξ和1-ξ在[0,1]都是均匀分布的随机变量，我们可以用ξ代替1-ξ，就得到：

θ = cos -1 ξ1

Φ = 2π ξ2

将这些值变换会笛卡尔坐标系中，就得到最后的采样公式：

我们用下面的代码来实现这个采样策略。我们传入两个均匀分布的随机数u1，u2，得到一个半球上的向量：

<MC Function Definitions>　＋＝

COREDLL Vector UniformSampleHemisphere(float u１, float u2) {

float z = u１;

float r = sqrtf(max(0.f, １.f - z*z));

float phi = 2 * M_PI * u2;

float x = r * cosf(phi);

float y = r * sinf(phi);

return Vector(x, y, z);

}

在pbrt中，对于每个类似这样的采样例程，都有一个对应的返回特定采样的PDF值的函数。对于这样的函数，我们需要搞清楚到底是用哪一个PDF求值---例如，在对半球上方向上采样时，我们已经看到这些密度函数即可以用立体角表示，也可以用（θ，Φ）表示。对半球（包括所有其它关于方向的采样）而言，这些函数返回关于立体角的值，立体角的PDF是一个常数 p(ω) = 1/(2π)。

<MC Function Definitions> +=

COREDLL float UniformHemispherePdf(float theta, float phi) {

return INV TWOPI;

}

在整个球面上的均匀采样公式可以用相同的方式推导出来，这里略去了。公式如下：

<MC Function Definitions> +=

COREDLL Vector UniformSampleSphere(float u1 , float u2)

float z = 1.f - 2.f * u1;

float r = sqrtf(max(0.f, 1.f - z*z));

float phi = 2.f * M_PI * u2;

float x = r * cosf(phi);

float y = r * sinf(phi);

return Vector(x, y, z);

}

<MC Function Definitions>+=

COREDLL float UniformSpherePdf() {

return 1.f / (4.f * M_PI);

}

14.5.2 举例：在单位圆盘上的采样

虽然圆盘看上去比半球更简单，但要进行均匀采样却更麻烦一些。一个直观但不正确的解决方案是： r = ξ1, θ = 2π ξ2。虽然所得到的结果是随机的且落在圆盘之内，但它并不是均匀分布的，而是采样点挤在靠近圆心的位置。

因为我们是根据面积来进行均匀采样，PDF p(x,y)必须是个常数。通过归一化约束我们得到 p(x,y) = 1 / π。如果我们变换到极坐标系中，则有p(r, θ) = r / π。现在我们使用前面介绍的方法计算边缘密度和条件密度：

跟半球上采样一样，由于圆的对称性，我们知道p(θ | r)是一个常数。通过积分和求反函数，就得到了P(r), P-1(r), P(θ) 和p-1(θ)。我们可以求得在圆盘上均匀采样的正确解：

r = ξ11/2

θ = 2π ξ2

<Monte Carlo Function Definitions> +=

void UniformSampleDisk(float u1, float u2, float *x, float *y) {

float r = sqrtf(u1);

float theta = 2.0f * M_PI * u2;

*x = r * cosf(theta);

*y = r * sinf(theta);

}

虽然解决了目前的问题，但圆盘上的面积却产生了变形（如图），第15.2.3节再介绍这种变形所带来的坏处。

Peter Shirley(1997)开发出一种从单位正方形到单位圆的“同心”映射来避免这个问题（如图）。

同心映射将正方形[-1,1]x[-1,1]上的点均匀地按照“同心正方形--同心圆”的方式进行映射，如图：正方形上的三角楔形被映射为饼条形状，映射公式为：

r = x

θ = y / x

<Monte Carlo Function Definitions> +=

void ConcentricSampleDisk(float u1, float u2, float *dx, float *dy) {

float r, theta;

<Map uniform random numbers to [−1, 1]x[-1,1]>

<Map square to (r , θ)>

*dx = r * cosf(theta);

*dy = r * sinf(theta);

}

<Map uniform random numbers to [−1, 1]x[-1,1]>

float sx = 2 * u1 - 1;

float sy = 2 * u2 - 1;

<Map square to (r , θ)> =

<Handle degeneracy at the origin>

if (sx >= -sy) {

if (sx > sy) {

<Handle first region of disk>

}

else {

<Handle second region of disk>

}

}

else {

if (sx <= sy) {

<Handle third region of disk>

}

else {

<Handle fourth region of disk>

}

}

theta *= M_PI / 4.f;

<Handle first region of disk> =

r = sx;

if (sy > 0.0) theta = sy/r;

else theta = 8.0f + sy/r;

14.5.3 举例：带余弦权值的半球采样

我们将在下一章看到，我们常常要从跟被积函数形状相似的分布函数进行采样。例如，散射方程用一个余弦项作为BSDF和入射辐射亮度的乘积的权值，我们希望找到这样的方法：它倾向于在靠近半球顶部的地方生成方向采样，在这个地方的余弦项要比靠近底部的余弦项大。

在数学意义上，这意味着我们在一个PDF采样方向ω，而p(ω)正比于cosθ。

下面是归一化推导，从而得到密度函数p:

我们可以象前面所做的那样计算出边缘密度和条件密度，但这里我们使用一种称为Malley法的技术来生成这些带余弦权值的点。其背后的思想是，如果我们从单位圆盘上均匀地采样，再将点投影到半球上去，就可以得到所要求的带有余弦分布的方向。

为什么可以这样做呢？令(r, Φ)是在圆盘上的极坐标。根据前面的计算，我们知道点在圆盘上的联合密度函数是p(r, Φ) = r / π。

现在我们将这个点映射到球面上，由于是垂直投影，故而sinθ = r。然后计算从(r, Φ)到(sinθ, Φ)的变换，我们可以得到雅克比行列式值：

|JT| = cosθ

故而p(θ, Φ) = |JT| p(r, Φ) = cosθ r / π，这正是我们想要得到的结果。

<Monte Carlo Utility Declarations> +=

inline Vector CosineSampleHemisphere(float u1, float u2) {

Vector ret;

ConcentricSampleDisk(u1, u2, &ret.x, &ret.y);

ret.z = sqrtf(max(0.f, 1.f - ret.x*ret.x - ret.y*ret.y));

return ret;

}

注意pbrt中所有的PDF求值例程都是用立体角定义的（而不是用球面坐标），故而PDF函数返回cosθ r/ π。

<Monte Carlo Utility Declarations> +=

inline float CosineHemispherePdf(float costheta, float phi) {

return costheta * INV_PI;

}

14.5.4 举例： 三角形上的采样

虽然在三角形上均匀采样看上去很简单，但实际上比前面所见到的采样都复杂。为了简化问题，我们值考虑在一个面积为1/2的等腰直角三角形上采样。采样例程的输出用的是重心坐标，这里的方法应适用于所有形状的三角形。

我们记点的重心坐标为(u,v)，因为我们根据面积采样，那么 PDF p(u,v)必然是个常数，且等于面积的倒数，即p(u,v)=2。

我们先球边缘密度p(u)：

条件密度为：

通过积分我们可以得到CDF：

对这些函数求反函数，在使用均匀的随机变量，得到：

u = 1 – (ξ1)1/2

v = ξ2 (ξ1)1/2

注意两个变量是并不是相互独立的。

<Monte Carlo Function Definitions> +=

void UniformSampleTriangle(float u1, float u2, float *u, float *v) {

float su1 = sqrtf(u1);

*u = 1.f - su1;

*v = u2 * su1;

}