PBRT读书笔记10-Monte Carlo Integration

这个系列是Physically Based Rendering: From Theory To Implementation的读书笔记，本节主要介绍Monte Carlo积分的相关知识。

在渲染过程中我们经常需要去计算各种各样的高维数值积分。这些被积函数往往十分复杂，甚至没有解析形式。计算这类函数的积分一般需要在被积范围内进行采样，因此这种方法称为Monte Carlo积分(Monte Carlo integration)。

Background and Probability Review

在正式介绍Monte Carlo积分前首先简要复习一下概率论和随机变量的相关知识。对于随机变量我们使用斜体$X$来表示，把某个函数$f$作用在$X$上时我们会得到一个新的随机变量$Y = f(X)$。离散型随机变量的累计概率分布函数(cumulative distribution function, CDF)定义为：

Continuous Random Variables

在渲染中更为常用的是连续型随机变量，其中[0, 1)区间上均匀分布的随机变量记为$\xi$。$\xi$是一个非常重要的随机变量，利用$\xi$可以生成更加复杂分布的随机变量。我们定义连续型随机变量的概率密度函数(probability density function, PDF)为CDF的微分：

对于均匀分布的随机变量$\xi$，它的PDF为

需要注意的是任意PDF必须满足非负而且积分为1。在给定区间[a, b]上PDF的积分即为随机变量落在该范围上的概率：

Expected Values and Variance

我们定义函数$f$在给定概率分布$p(x)$上的期望(expected value)为$f(x)$与$p(x)$的积分：

定义$f$的方差(variance)为$f(x)$偏离自身期望的程度：

期望和方差的常用性质如下：

对于相互独立的随机变量，它们的方差具有可加性：

The Monte Carlo Estimator

对于一元函数$f(x)$，它在[a, b]上的积分$\int_a^b f(x) dx$可以通过[a, b]上均匀分布随机变量$X_i$的函数来估计：

实际上$F_N$的期望恰好等于积分：

实际上我们并不需要将随机变量限制为均匀分布，对于任意分布$p(x)$只需要对估计进行一定的变形：

同样可以证明上式的期望等于积分值：

Monte Carlo积分同样可以拓展到多元函数上。假设我们想要计算一个三重积分：

为了计算积分值需要从$[x_0, x_1] \times [y_0, y_1] \times [z_0, z_1]$进行均匀采样，每个样本的PDF均为：

因此使用Monte Carlo积分得到的估计值为：

对于Monte Carlo积分而言样本数$N$可以是任意的。这样的好处在于无论被积函数的维度如何我们都可以很好地控制计算复杂度，这也是Monte Carlo积分优于很多其它数值积分方法的原因之一。

除此之外还需要考虑的是Monte Carlo积分的精度，实际上可以证明Monte Carlo积分误差的收敛速度为$O(\sqrt{N})$且与维数无关。对于一维的情况很多数值积分方法都有着更好的收敛速度，但随着维数的增加这些数值积分的精度会显著下降。因此在多维的情况下Monte Carlo积分几乎是唯一可行的积分方法。

Sampling Random Variables

Monte Carlo积分的难点在于如何对各种不同类型的分布进行采样，在本节中我们暂时只介绍最简单的两种采样方法，更多复杂的采样方法会留到后面的章节再进行介绍。

The Inversion Method

最基本的采样方法是逆变换采样(inverse transform sampling)。对于任意分布的随机变量$X$可以证明把CDF作用到$X$上后会得到均匀分布，即$P(X) \sim U(0, 1)$。因此我们只需要从均匀分布进行采样，然后利用$P(x)$的反函数即可获得来自$X$分布的样本。

逆变换采样的步骤如下：

• 计算随机变量$X$的CDF：\(P(x) = \int_0^x p(x') \ dx'\)
• 计算CDF的反函数$P^{-1}(x)$
• 从均匀分布$U(0, 1)$进行采样，得到样本$\xi$
• 利用CDF的反函数得到来自$X$分布的样本$X_i = P^{-1}(\xi)$

Power Distribution

对于幂分布(power distribution)，它的PDF可以表示为：

利用PDF的归一化条件可以求得常数$c$为：

因此幂分布的CDF为：

对应的反函数为：

因此要获得幂分布的样本只需要对均匀分布的样本$\xi$进行逆变换即可：

Exponential Distribution

对于参数为$a$的指数分布，其CDF为：

对应的反函数为：

因此获取指数分布的样本只需要利用变换：

Piecewise-Constant 1D Functions

很多时候我们需要考虑分段的概率分布。假设每一段的大小是均匀的$\Delta = \frac{1}{N}$，此时未归一化的PDF可以表示为：

利用PDF的归一化条件可以得到归一化常数$c$：

因此CDF在不同间隔上的值为：

进行逆变换采样时需要计算$\xi$位于的区间使得$P(x_i) \leq \xi$且$\xi \leq P(x_{i+1})$。在已知CDF的情况下可以使用二分查找的方法来加速计算过程。

The Rejection Method

对于某些分布的CDF直接进行计算可能非常复杂甚至没有解析形式，此时就不能使用逆变换采样。在这种情况下可以使用接受-拒绝采样(acceptance-rejection sampling)来生成样本。接受-拒绝采样不需要计算概率分布的CDF，只需要目标分布$f(x)$和另一个已知分布$p(x)$的PDF即可，同时两个分布的PDF需要满足$f(x) \leq c p(x)$，其中$c$为常数。

接受-拒绝采样的步骤如下：

• 从已知分布$p(x)$上生成样本$X$
• 从均匀分布$U(0, 1)$上采样一个随机数$\xi$
• 如果$\xi \leq \frac{f(X)}{c p(X)}$则接受$X$作为样本，否则回到第一步重新进行采样

从几何意义上讲接受-拒绝采样会生成随机变量对$(X, \xi)$，如果点$(X, \xi \cdot c p(X))$位于曲线$f(X)$的下方则接受样本。接受-拒绝采样的一个优势在于它与随机变量的维数无关，只要能够计算PDF即可。

Rejection Sampling a Unit Circle

接受-拒绝采样的经典应用是生成圆内均匀分布的样本。我们只需要在正方形区域内进行均匀采样，然后拒绝掉落在圆外的样本即可。

在PBRT中使用RejectionSampleDisk()函数来描述单位圆内进行均匀采样的过程：

Point2f RejectionSampleDisk(RNG &rng) {
    Point2f p;
    do {
        p.x = 1 - 2 * rng.UniformFloat();
        p.y = 1 - 2 * rng.UniformFloat();
    } while (p.x * p.x + p.y * p.y > 1);
    return p;
}

Metropolis Sampling

Metropolis采样(Metropolis sampling)是一种针对非负函数$f(x)$进行采样的方法，其中可以认为$f(x)$是归一化前的概率密度。和前面介绍过的采样方法相比，Metropolis采样不需要积分也不需要计算函数的反函数，只要能对$f(x)$求值即可。而且Metropolis采样不会拒绝生成的样本，在每一轮的采样过程中都可以得到新的数据。因此Metropolis采样比前面介绍过的采样方法都要高效。

Metropolis采样的缺陷在于它是一个迭代的算法，换句话说连续的样本之间是存在一定关联的。它的另一个缺陷在于当采样数比较少时生成的样本可能无法覆盖整个定义域，因此很多降低方差的技术不能直接用在Metropolis采样中。

Basic Algorithm

Metropolis采样的目标是从给定的函数$f: \Omega \rightarrow \mathbb{R}$上采样出一系列样本$X_i$。设初始的样本为$X_0$，在每一轮采样过程中样本$X_i$都是由前一个样本$X_{i-1}$经过一个随机突变(mutation)来获得，突变的结果记为$X’$。我们可以接受或拒绝突变的样本，因此$X_i$的取值为$X_{i-1}$或$X’$。在这种采样模式下$X_i$会达到一种平衡态分布，这个概率分布称为平稳分布(stationary distribution)。在极限情况下$X_i$的概率分布即为$f(x)$定义的概率分布：\(p(x) = \frac{f(x)}{\int_\Omega f(x) \ d\Omega}\)。

显然Metropolis采样的核心是设计突变函数以及接受突变的策略。突变函数可以有不同的形式，比如说对当前样本$X$进行扰动或是直接生成一个新的样本都可以。假设我们已经有了某种突变策略，接下来需要考虑从$X$突变为$X’$的概率$T(X \rightarrow X’) = P(X’ \vert X)$，称为状态转移函数(transition function)。在此基础上我们可以定义一个接受概率$a(X \rightarrow X’)$表示我们是否要接受突变的样本，从而保证样本$X$的分布正比于$f(x)$。当系统达到平衡态时任意两个状态相互转换的概率必须相等，即：

上式称为细致平衡条件(detailed balance)。当$f$和$T$确定时我们可以得到接受概率$a$的表达式：

容易验证上式定义的接受概率$a$满足细致平衡条件。如果进一步假设$X$和$X’$相互突变的概率是相等的，则可以得到更简单的接受概率：

这样我们就得到了最基本的Metropolis采样算法，它的伪代码形式如下：

X = X0
for i = 1 to n
    X' = mutate(X)
    a = accept(X, X')
    if (random() < a)
        X = X'
    record(X)

当$f(X’)$的值比较小时Metropolis采样会有很大的概率拒绝突变，这会使得只有很少的样本会分布在函数值比较小的区域。为了克服这个问题我们可以使用求期望的方法来改进标准的Metropolis采样算法。此时在每一轮采样中需要同时记录$X$和$X’$，并且为它们赋予相应的概率(权重)。改进后的Metropolis采样算法伪代码如下：

X = X0
for i = 1 to n
    X' = mutate(X)
    a = accept(X, X')
    record(X, 1 - a)
    record(X', a)
    if (random() < a)
        X = X'

Choosing Mutation Strategies

只要能够计算转移概率$T(X \rightarrow X’)$，理论上突变函数可以是任意的。实际应用中为了化简Metropolis采样的过程一般会要求$X$与$X’$之间相互突变的概率是相同，在这个前提下人们总结出了许多常用的突变策略。

一般来说我们希望突变后的样本能够迅速转移到函数值较大的区域而不是停留在函数值较小的区域，但是当$f(X)$很大时Metropolis采样会以很小的概率转移到突变后的样本。我们希望能够避免这种情况，并且使样本尽可能分散在整个定义域上。因此一种突变方法是对$X$进行随机扰动：假设$X$是一个向量且第$i$维是$x_i$，我们通过对$x_i$加减一个随机值实现扰动：

其中$s$表示扰动的尺度；对扰动结果取余数用来进行归一化。显然这种扰动方法满足对称性$T(X \rightarrow X’) = T(X’ \rightarrow X)$。

一种类似的突变方法是直接使用随机数来代替$x_i’$：

这种突变方式同样是对称的，而且它还可以保证样本在突变过程中能够覆盖到整个定义域而不会局限在某些区域上。

除此之外我们还可以利用已知的PDF来进行突变。假设已知一个概率分布$p(x)$，我们可以直接从遮盖概率分布中进行采样作为突变后的样本。此时转移概率可以表示为：

此时突变的结果与当前状态$X$无关。

Start-up Bias

在Metropolis采样中还需要注意初值$X_0$的选取，如果$X_0$选择的不好会导致一些初值误差(start-up bias)。要解决这个问题一种常用的做法是从一个随机初值开始进行采样，抛弃刚开始生成的样本然后把当前状态作为初值重新进行采样。这种做法的缺陷在于抛弃刚开始产生的样本会浪费计算资源，同时我们无法得知要进行多少次迭代才能开始真正的采样。

因此实践中更常用的方式是首先从另一个分布$p(x)$中采样得到初值$X_0$，然后从$X_0$开始进行Metropolis采样。同时这些样本还具有权重：

如果$f(X_0)=0$，此时所有的样本权重均为0。为了解决这个问题我们可以从$p(x)$中采样出一系列初值的候选值$Y_1, Y_2, …, Y_N$，每个候选值的权重为：

然后根据$Y_i$的权重抽样出初始值$X_0$进行Metropolis采样即可，此时采样出的样本权重均为$w_i$的平均值。

1D Setting

接下来我们看一个实际的例子。假设$f(x)$的表达式为：

$f(x)$的函数图像可见下图。为了展示Metropolis采样的能力，这里假装$f(x)$是一个黑盒：我们可以对$f(x)$进行求值，但不知道具体的解析式。

下面考虑两种不同的突变策略。第一个策略是直接从[0, 1]均匀分布上采样一个新的样本，在这种情况下我们有：

第二个策略是对$X$进行一个随机扰动，扰动的大小是0.05：

对于初值$X_0$，我们直接从[0, 1]上的均匀分布进行采样：

这样每个样本的权重均为$w = f(X_0)$。

然后我们就可以通过Metropolis采样来重建$f(x)$。由于在每次采样中我们都记录了两个样本$X$和$X’$以及它们对应的权重，在重建时可以通过直方图的方式来估计样本的概率密度，此时只需要对每个小区间内的样本权重进行求和即可。重建结果如下图所示，其中左图为单纯使用第一种突变策略的重建结果，而右图为混合使用两种突变策略的结果，其中有10%的概率选择第一个策略90%的概率选择第二个策略。

直观来看直接在均匀分布上进行采样并不是一个很好的策略。它的缺陷在于当函数值$f(x)$比较大时不能很好地利用该点的概率密度，仍会随机跳跃到其它区域上。尽管如此，在这种策略下Metropolis采样仍然能够收敛到正确的分布上。而混合策略则可以更好地近似$f(x)$对应的分布，这是由于此时的采样策略可以利用$f(x)$的函数值来移动到高概率的区间上。

如果单纯使用随机扰动作为突变则会得到下图所示的重建结果。其中左图是10,000个样本的重建结果，而右图是300,000个样本的结果。产生这种现象的原因在于当突变的步长比较小时样本很难从一侧跨越到曲线的另一侧，因此需要更多的采样次数才能得到合理的结果。

Estimating Integrals with Metropolis Sampling

利用Metropolis采样我们可以估计形如$\int f(x) g(x) \ d\Omega$的积分。利用Monte Carlo积分公式可以得到估计值：

根据Metropolis采样，我们可以把$p(x)$取为$f(x)$。这样积分就等于：

Transforming between Distributions

在逆变换采样一节中我们介绍过将均匀分布转换为其它分布的方法，本节中我们将讨论把函数$f$作用在某个已知分布后会得到怎样的分布。假设已知随机变量$X$服从某个概率分布，其PDF为$p_x(x)$，我们的目标是计算随机变量$Y = y(X)$的分布。当函数$y(x)$是一个一对一映射时，利用概率的定义可以得到：

因此CDF满足：

我们进一步假设$y(x)$的导数严格大于0，可以得到PDF的关系：

即

更宽泛地来说，如果$y(x)$导数严格为正或是严格为负，$Y$和$X$之间的PDF满足：

利用上式就可以计算出$Y$的分布。举个简单的例子，设$X$是定义在[0, 1]区间上的随机变量且PDF为$p_x(x)=2x$。令$Y=\sin(X)$，则$Y$的PDF为：

在实际中更常见的情况是从样本分布来估计函数$y(x)$。比如说我们有一些样本$X$，目标是计算来自$p_y(y)$分布的样本$Y$。那么如何取建立$X$和$Y$之间的联系呢？实际上非常简单，只需要令$X$和$Y$的CDF相等即可构造出映射：

不难发现当$X$为[0, 1]区间上的均匀分布时上式即为逆变换采样，因此上式可以看做是对逆变换采样的推广。

Transformation in Multiple Dimensions

在很多情况下我们需要计算多元随机变量之间PDF的关系，其中$Y=T(X)$且$T$是一个双射。此时$X$和$Y$的PDF满足：

式中$\vert J_T(x) \vert$表示映射$T$ Jacobian矩阵行列式的绝对值，即

Polar Coordinates

对于极坐标的情况，首先考虑坐标转换公式：

因此坐标变换的Jacobian矩阵为：

它的行列式为$r (\sin^2(\theta) + \cos^2(\theta)) = r$，因此极坐标和直角坐标的PDF转换公式为：

在实际应用中一般是从直角坐标进行采样然后变换到极坐标上，因此上式更常见的形式为：

Spherical Coordinates

在三维空间中球坐标有非常多的应用，它和直角坐标的转换公式为：

容易验证球坐标转换到直角坐标的Jacobian矩阵满足：

因此两个坐标系下PDF的转换公式为：

回忆在单位球上的微分立体角满足公式：

假设我们在立体角上定义了某个概率分布$p(\omega)$，利用上式可以得到它和单位球参数平面上PDF的转换关系：

2D Sampling with Multidimensional Transformations

假设我们想要从二维联合概率$p(x, y)$进行采样。当$X$和$Y$相互独立时，我们可以对联合概率进行分解：

此时只要分别对$X$和$Y$进行采样即可。然而在大部分情况下$X$和$Y$都不是相互独立的，因此我们需要介绍对多维随机变量进行采样的相关理论。

我们定义随机变量的边缘概率密度函数(marginal density function)为对其它变量积分后的结果：

显然边缘概率密度仅仅是单个随机变量的函数。类似地，我们可以定义条件概率密度函数(conditional density function)为联合概率与边缘概率的商：

因此对二维随机变量进行采样的基本思路是首先从边缘概率中采样一个样本，然后通过条件概率采样另一个维度上的随机变量。

Uniformly Sampling a Hemisphere

接下来看一个简单的例子。假设我们想要对半球面上的均匀分布进行采样，利用归一化条件可以得到立体角的PDF：

转换到参数平面上有：

尽管$\theta$和$\phi$是相互独立的，我们仍然可以利用边缘概率和条件概率来进行采样。首先考虑$\theta$，通过对$\phi$积分可以得到$\theta$的边缘概率密度：

然后计算$\phi$所需的条件概率密度：

注意到条件概率密度为[0, $2\pi$]区间上的均匀分布。利用逆变换方法即可完成对边缘概率和条件概率的采样：

因此只需要采样两个均匀分布的随机数即可：

注意这里在采样$\theta$时使用了$1-\xi_1$来代替$\xi_1$。如果想要得到直角坐标系下的表示则可以利用球坐标和直角坐标的转换公式：

在PBRT中使用UniformSampleHemisphere()函数来对半球进行采样，其输入u为两个均匀分布随机数组成的Point2f实例：

Vector3f UniformSampleHemisphere(const Point2f &u) {
    Float z = u[0];
    Float r = std::sqrt(std::max((Float)0, (Float)1. - z * z));
    Float phi = 2 * Pi * u[1];
    return Vector3f(r * std::cos(phi), r * std::sin(phi), z);
}

每个样本的PDF均为$\frac{1}{2 \pi}$，在PBRT中可以使用UniformHemispherePdf()函数来计算：

Float UniformHemispherePdf() {
    return Inv2Pi;
}

如果需要在整个球面上进行采样，只需对公式进行简单变形：

对应的采样函数为：

Vector3f UniformSampleSphere(const Point2f &u) {
    Float z = 1 - 2 * u[0];
    Float r = std::sqrt(std::max((Float)0, (Float)1 - z * z));
    Float phi = 2 * Pi * u[1];
    return Vector3f(r * std::cos(phi), r * std::sin(phi), z);
}

Float UniformSpherePdf() {
    return Inv4Pi;
}

Sampling a Unit Disk

接下来考虑如何在单位圆内进行均匀采样。一个直观的想法是直接生成两个随机数然后映射到单位圆内：

然而这种做法是完全错误的，实际上它会导致样本集中在圆心的位置，如下图中左边的结果。正确的采样结果应该是样本均匀分散在圆内，如下图中右边的结果。

要得到正确的结果我们还是需要从边缘概率以及条件概率进行采样。首先考虑单位圆内的PDF为：

利用直角坐标和极坐标PDF的转换公式得到：

因此得到边缘概率和条件概率为：

再根据逆变换方法可以得到采样公式：

在PBRT中使用UniformSampleDisk()函数来表示对单位圆盘进行均匀采样的过程：

Point2f UniformSampleDisk(const Point2f &u) {
    Float r = std::sqrt(u[0]);
    Float theta = 2 * Pi * u[1];
    return Point2f(r * std::cos(theta), r * std::sin(theta));
}

尽管这种采样方法是正确的，但有时会造成采样结果的扭曲。下图中每块圆环和扇形的面积都是相等的，但随着半径的增加圆环的区域会被挤压。

更好的采样方法是保证每个小区域的形状在对应圆环上都是相似的，如下图所示。

要改进原来的采样方法也非常简单，只需要把单位正方形区域中的点压缩映射到单位圆上即可。变换公式如下：

在PBRT中使用ConcentricSampleDisk()函数来描述上述采样过程：

Point2f ConcentricSampleDisk(const Point2f &u) {
    // Map uniform random numbers to [-1, 1]^2
    Point2f uOffset = 2.f * u - Vector2f(1, 1);

    // Handle degeneracy at the origin
    if (uOffset.x == 0 && uOffset.y == 0)
        return Point2f(0, 0);

    // Apply concentric mapping to point
    Float theta, r;
    if (std::abs(uOffset.x) > std::abs(uOffset.y)) {
        r = uOffset.x;
        theta = PiOver4 * (uOffset.y / uOffset.x);
    } else {
        r = uOffset.y;
        theta = PiOver2 - PiOver4 * (uOffset.x / uOffset.y);
    }
    
    return r * Point2f(std::cos(theta), std::sin(theta));
}

Cosine-Weighted Hemisphere Sampling

渲染中的很多函数都包含一个余弦项$\cos(\theta)$，因此我们希望进行采样时能够更多地采样到接近法线的方向上。更严格的说法是我们希望能够构造一个立体角上的PDF满足：

利用归一化条件我们可以得到：

因此得到参数平面上的PDF：

我们可以利用前面用过的边缘概率以及条件概率进行采样，但这里再介绍一种新的采样方法，Malley方法(Malley’s method)。它的思想是首先在单位圆上进行均匀采样，然后把这些采样得到的点映射到单位球的上半球上，这个过程如下图所示。

接下来证明Malley方法的正确性。假设单位圆盘内的点坐标为$(r, \phi)$，对应的PDF为：

然后我们把这个点映射到单位球的上半球上。根据球坐标可以得到：

因此单位圆到半球面的映射可以表示为：

对应Jacobian矩阵的行列式为：

因此半球和单位圆盘上PDF的转换关系为：

上式即为根据法线夹角余弦在上半球进行采样时的PDF。

在PBRT中使用CosineSampleHemisphere()函数来描述使用Malley方法对上半球进行采样的过程：

inline Vector3f CosineSampleHemisphere(const Point2f &u) {
    Point2f d = ConcentricSampleDisk(u);
    Float z = std::sqrt(std::max((Float)0, 1 - d.x * d.x - d.y * d.y));
    return Vector3f(d.x, d.y, z);
}

inline Float CosineHemispherePdf(Float cosTheta) {
    return cosTheta * InvPi;
}

这里需要说明的是在PBRT中我们默认使用立体角而不是球面坐标来描述方向，因此这里PDF的返回值为$\frac{\cos{\theta}}{\pi}$。

Sampling a Cone

在某些情况下我们需要对锥面进行均匀采样。此时的PDF为：

因此，

上式说明$\theta$和$\phi$是相互独立的，我们只需要分别进行采样即可。同时在直角坐标系下我们更常用的是$\cos(\theta)$而不是$\theta$，因此可以构造出$\cos(\theta$的采样公式：

Vector3f UniformSampleCone(const Point2f &u, Float cosThetaMax) {
    Float cosTheta = ((Float)1 - u[0]) + u[0] * cosThetaMax;
    Float sinTheta = std::sqrt((Float)1 - cosTheta * cosTheta);
    Float phi = u[1] * 2 * Pi;
    return Vector3f(std::cos(phi) * sinTheta, std::sin(phi) * sinTheta,
                    cosTheta);
}

Float UniformConePdf(Float cosThetaMax) {
    return 1 / (2 * Pi * (1 - cosThetaMax));
}

Sampling a Triangle

对三角形进行采样时需要使用重心坐标(barycentric coordinate)。为了便于推导，我们假设三角形内的中心坐标为$(u, v)$且三角形的面积为$\frac{1}{2}$，那么在三角形内进行均匀采样时的PDF为：

因此边缘概率密度为：

对应的条件概率密度为：

分别进行积分得到CDF：

利用逆变换采样得到样本的重心坐标：

在PBRT中对应的采样函数为UniformSampleTriangle()：

Point2f UniformSampleTriangle(const Point2f &u) {
    Float su0 = std::sqrt(u[0]);
    return Point2f(1 - su0, u[1] * su0);
}

Sampling Cameras

这一小节中介绍了对真实相机进行采样的方法。由于我们不考虑真实相机模型，此处暂时略过。

Piecewise-Constant 2D Distributions

对于二元离散型随机变量，其PDF会通过一个矩阵来进行存储。在这种情况下进行采样需要一些额外的处理，首先我们需要把函数$f(u_i,v_j)$规范到$[0, 1]^2$区间上并对积分进行归一化：

这样可以得到PDF：

其中$u$、$v$为规范化平面上的坐标，$\tilde{u}$、$\tilde{v}$分别为矩阵元素的坐标：

因此边缘概率密度为：

对应的条件概率密度为：

在PBRT中使用Distribution2D类来存储二元离散随机变量的PDF，并实现了采样功能。

class Distribution2D {
public:
    Distribution2D(const Float *data, int nu, int nv);
    Point2f SampleContinuous(const Point2f &u, Float *pdf) const {
        Float pdfs[2];
        int v;
        Float d1 = pMarginal->SampleContinuous(u[1], &pdfs[1], &v);
        Float d0 = pConditionalV[v]->SampleContinuous(u[0], &pdfs[0]);
        *pdf = pdfs[0] * pdfs[1];
        return Point2f(d0, d1);
    }

    Float Pdf(const Point2f &p) const {
        int iu = Clamp(int(p[0] * pConditionalV[0]->Count()),
                       0, pConditionalV[0]->Count() - 1);
        int iv = Clamp(int(p[1] * pMarginal->Count()),
                       0, pMarginal->Count() - 1);
        return pConditionalV[iv]->func[iu] / pMarginal->funcInt;
    }

private:
    std::vector<std::unique_ptr<Distribution1D>> pConditionalV;
    std::unique_ptr<Distribution1D> pMarginal;
};

Distribution2D::Distribution2D(const Float *func, int nu, int nv) {
    for (int v = 0; v < nv; ++v) {
        // Compute conditional sampling distribution for v
        pConditionalV.emplace_back(new Distribution1D(&func[v * nu], nu));
    }

    // Compute marginal sampling distribution p[v]
    std::vector<Float> marginalFunc;
    for (int v = 0; v < nv; ++v)
        marginalFunc.push_back(pConditionalV[v]->funcInt);
    
    pMarginal.reset(new Distribution1D(&marginalFunc[0], nv));
}

Russian Roulette and Splitting

我们定义估计值$F$的效率(efficiency)为：

其中$V[F]$是$F$的方差，而$T[F]$是计算$F$所需的时间。为了提高Monte Carlo积分的效率，人们开发出了俄罗斯轮盘赌(Russian roulette)以及splitting两种技术。其中俄罗斯轮盘赌可以解决某些样本对积分贡献小的问题，而splitting则可以在采样时让样本更关注相对重要的维度。

Russian Roulette

下面以反射方程为例来介绍俄罗斯轮盘赌。设光源为$L_d$，则$p$点接收到来自光源的直接光照所产生的反射radiance为：

如果我们只有2个样本，则反射radiance的估计值为：

在进行计算时大部分的计算资源会用在测试光源与$p$点之间是否存在遮挡上，如果存在遮挡则来自光源的radiance为$L_d(p, \omega_i) = 0$。显然对于那些存在遮挡的方向$\omega_i$，我们希望能够直接跳过它们因为它们对积分值的贡献为0。而俄罗斯轮盘赌则可以实现这样的功能，不仅如此它还可以跳过那些$f_r(p, \omega_o, \omega_i)$比较小或是$\vert \cos(\theta_i) \vert$接近于0的光线样本。

使用俄罗斯轮盘赌时首先需要设置一个终止概率$q$，它表示对于某些样本会以$q$的概率进行跳过。同时我们需要设置一个常数$c$，它表示跳过这些样本时使用$c$来代替，一般可以将$c$设置为0。而对于那些没有跳过的样本，我们为它们赋予权重$\frac{1}{1-q}$。这样我们可以利用$F$构造一个新的估计$F’$：

实际上可以证明$F’$是无偏的：

需要说明的是俄罗斯轮盘赌不能降低Monte Carlo积分的方差。实际上除非$c=F$，俄罗斯轮盘赌反而会增加方差。而且当$q$设置的不好时，俄罗斯轮盘赌可能会使方差急剧上升进而降低估计值的效率。

Splitting

splitting的思想是通过更好地分配样本来提高性能。以计算直接光照为例，图像上$(x, y)$像素区域内的直接光照可以表示为：

要计算这个积分我们可以直接采样$N$个样本然后使用Monte Carlo积分来进行估计，其中每个样本包含两道光线：从$(x, y)$到场景的光线以及从场景到光源的光线。

假设$N=100$，此时我们需要100道从相机发出的光线和100道从场景中发出的光线。显然100道从相机发出的光线有一点多余了，就算是考虑反走样也不需要这么多的光线。更合理的策略是只从相机发射出少量的光线，更多的计算资源用在场景向光源发射光线上。

splitting的思想就在于此，我们可以把积分估计值分解为：

然后考虑5条从相机发出的光线以及20条从场景发出的光线。这样我们就只要考虑105条光线而不是200条光线，从而可以极大地提高计算效率。

Careful Sample Placement

接下来我们介绍降低Monte Carlo积分方差的相关技术。在PBRT中大量使用了经过设计的样本而不是完全随机的样本，而这些设计过的样本可以大大减少方差。

Stratified Sampling

在Stratified Sampling一节中介绍过分层抽样的概念，接下来我们把分层抽样应用在Monte Carlo积分中。假设整个积分域为$[0, 1]^s$，在每一个区域$\Lambda_i$内$f(x)$的估计值为：

其中$X_{i,j}$表示在区域$\Lambda_i$上通过分布$p_i$进行抽样得到的第$j$个样本。因此整个定义域上Monte Carlo积分可以表示为：

其中$v_i$为每个区域占积分域总体的比例。

对每个区域单独分析可以得到区域$\Lambda_i$上的均值和方差为：

因此当区域内的样本数为$n_i$时，$F_i$的方差为$\frac{\sigma_i^2}{n_i}$。把所有区域加起来可以得到分层抽样情况下总体的方差为：

如果我们进一步假设每个区域内的样本数正比于区域的大小，即$n_i = v_i N$，则可以得到化简后分层抽样的方差：

作为对比，我们考虑不使用分层抽样时的方差。

其中$Q$为$F$在整个积分域$\Lambda$上的均值。具体的推导过程可以参考Veach的博士论文。

把两种采样方法的方差进行对比不难发现使用分层抽样时的方差是要小于等于正常抽样的结果的，换句话说分层抽样永远不会增大方差。如果想要最大化分层抽样的性能则需要让$\sum v_i (\mu_i - Q)^2$取最大值，即让每一个积分区域内的均值尽可能远离$Q$。也因此增大分层的数量减少每个区域的大小有利于提升分层抽样的性能。

当然分层抽样也有自身的缺点，比如说它会受到维数灾难(curse of dimensionality)的影响很难直接应用在高维积分中。要缓解这个问题可以使用前面介绍过的随机关联方法，同时在进行分层时要尽可能选择和积分比较相关的维度。除此之外还可以考虑使用Latin超立方抽样来降低所需的样本数，尽管Latin超立方抽样可能不如进行分层抽样有效但也要好于直接抽样的效果。

Quasi Monte Carlo

前面介绍过的低差异序列并不是真正的随机数，使用它们进行采样的方法称为伪Monte Carlo(quasi Monte Carlo)。和标准Monte Carlo方法相比，伪Monte Carlo往往有着更好的收敛速度，而且很多Monte Carlo方法中使用的技术可以直接应用在伪Monte Carlo中。在渲染中使用伪Monte Carlo方法一般会比Monte Carlo方法有更好的表现。

Warping Samples and Distortion

在PBRT中进行采样一般是在标准区间如$[0, 1)^2$上进行的，而在实际计算积分时往往需要会对区间进行缩放，这就导致了原始区间上均匀分布的样本在缩放后产生一定的扭曲。以下图为例，(a)图表示从$4\times4$区间上进行分层抽样然后缩放到长方形区域的样本，而(b)图则是直接在长方形区域上使用低差异序列进行采样的样本。显然原本在正方形区域上均匀分布的样本经过缩放后的质量下降了，这个例子也说明了低差异序列和分层抽样的样本相比有着更好的适应性。

Bias

减少方差的另一种方法是引入偏差(bias)。我们定义估计值的偏差为：

当$\beta=0$时我们称估计值时无偏(unbiased)的。

接下来我们用一个例子来说明为什么引入偏差可以降低方差。假设目标是估计[0, 1]区间上均匀分布的均值，那么一个无偏的估计是：

同时我们考虑一个有偏的估计：

显然$F_1$是一个无偏的估计，它的方差为$O(N^{-1})$。而$F_2$的期望为：

因此$F_2$是一个有偏的估计，但可以证明它的方差为$O(N^{-2})$。这个简单的例子说明通过引入偏差可以降低估计的方差。

我们再看一个图像滤波的例子。当渲染完成图像后一般需要通过滤波来处理噪声，此时重建的图像可以表示为：

其中$I(x, y)$为重建后的图像，$f(x, y)$为滤波函数，$L(x, y)$为渲染得到的图像。为了便于讨论我们假设了$f(x, y)$的积分为1。

如果使用均匀采样和Monte Carlo积分来求解，我们可以得到估计：

其中$p_c$为均匀采样的PDF。

而一个有偏的估计为：

实际上第二种估计有更小的方差，也是PBRT中使用的滤波计算方法。为了证明这一点可以假设$L(x, y)$恒为1，此时有偏的估计会同样给出1而无偏估计则会由于采样无法给出稳定的值。对于更加复杂的图像，有偏估计引入的偏差会远小于无偏估计自身的方差，因此推荐使用有偏估计的公式来对图像进行滤波。

Importance Sampling

重要性采样(importance sampling)是最常用的方差缩减技术。它的思想是当分布$p(x)$的形状接近于积分函数$f(x)$时，我们只需要少量的样本就可以近似积分值。想要严格证明这个性质比较复杂，这里我们用一个简单的例子来进行说明。

假设待计算的积分为$\int f(x) \ dx$，我们可以构造一个分布$p(x)$使得$p(x) \propto f(x)$，即：

当我们从$p(x)$进行采样并且计算Monte Carlo积分时，每一个样本都会直接得到正确的积分：

此时我们的估计值没有任何方差。当然在大多数情况下我们没有办法构造出这样的概率分布$p(x)$，不过只要$p(x)$的形状接近于$f(x)$就可以减少方差。

重要性采样的缺陷在于当$p(x)$和$f(x)$的形状差异很大时它反而会增大估计值的方差。不过对于渲染中的积分我们总是可以找到形状类似于被积函数的概率分布，关于如何计算这些积分的方法我们留到后面的章节再进行介绍。

Multiple Importance Sampling

很多情况下被积函数可以分解成几个函数乘积的形式，即$\int f(x) g(x) \ dx$。此时我们既可以对$f(x)$进行重要性采样，也可以对$g(x)$进行重要性采样。如何设计采样的策略会对最终的渲染效果有重要的影响。

接下来看一个具体的例子，直接光照的反射方程为：

我们既可以对BSDF项$f_r$进行采样，也可以对光源项$L_d$进行采样，不过这种直接进行采样的方法往往不会得到非常好的渲染结果。比方说当$p$点位于镜面上时$f_r$只在很少的方向上有非零值，如果此时对$L_d$进行采样会导致大部分的样本光线都是0，从而产生非常大的方差；类似地，如果$p$点位于粗糙平面上对$f_r$进行采样也会导致比较大的方差。

为了处理这样的问题人们提出了多重重要性采样(multiple importance sampling, MIS)的方法。对于形如$\int f(x)g(x) \ dx$的积分使用MIS时需要同时对$f$和$g$进行采样，此时的估计值为：

其中$X_i$和$Y_j$为采样的样本，$p_f$和$p_g$为样本的概率，$n_f$和$n_g$为采样的数目，而$w_f$和$w_g$则是样本的权重。MIS的核心在于如何设计权重使得计算结果是无偏的。一种常用的权重函数是根据概率密度和样本数量进行平衡加权：

为了证明MIS的效果我们考虑这样一种情况：假设从$p_f(x)$中进行采样得到了样本$X$，此时$p_f(X)$和$f(X)$的值都比较小但$g(X)$的值比较大。如果直接使用Monte Carlo积分会得到：

由于分母$p_f(X)$很小会得到一个很大的样本值，这说明此时积分的方差会很大。而如果使用MIS则会得到：

此时只要$p_g(x)$和$g(x)$比较接近就可以避免分母很小导致估计值偏大的问题。从这个例子可以发现MIS可以通过平衡不同分布的权重来地减少Monte Carlo积分的方差。

在PBRT中使用BalanceHeuristic()函数来计算平衡加权的权重：

inline Float BalanceHeuristic(int nf, Float fPdf, int ng, Float gPdf) {
    return (nf * fPdf) / (nf * fPdf + ng * gPdf);
}

在实际应用中还可以根据幂次进行加权：

参数$\beta$根据经验可以取$\beta = 2$。在PBRT中使用PowerHeuristic()函数来计算按幂次加权的权重：

inline Float PowerHeuristic(int nf, Float fPdf, int ng, Float gPdf) {
    Float f = nf * fPdf, g = ng * gPdf;
    return (f * f) / (f * f + g * g);
}

Reference

• 13 Monte Carlo Integration
• Chapter 2: Monte Carlo Integration
• Monte Carlo Integration