通俗讲解清楚光线追踪中的蒙特卡洛积分

其实应该说是路径追踪中的蒙特卡洛积分，一般蒙特卡洛积分是在路径追踪这个计算方法里使用的。这篇文章将从最基本的地方开始，逐步追踪内容去讲解清楚路径追踪中的蒙特卡洛积分，给在图形渲染中遇到蒙特卡洛积分而且看不懂的伙伴梳理清楚——什么是蒙特卡洛积分，在路径追踪中用来干嘛？

主要网上关于蒙特卡洛积分的介绍五花八门的，这回我必须要把它说清楚，让小白也能清楚看明白。

---

1.从最基本的，我们平时是怎么求解积分的

一个积分公式，最基础是长什么样子的：

$$I={\int }_a^{b}f(x)dx$$

很基础是吧，那么我给它套成一个很简单的积分公式：
$$I={\int }_0^1{x}^{2}dx$$

看着就很舒服很亲切是吧，我们甚至不需要画函数图去看这个函数长什么样子，感觉我们可以很轻易把这个积分求出来，那么我们是怎么求解这个积分的呢？我们平常数学课学习的都是直接用微积分的方式去计算出来。如果你忘记了或者没反应出来，没关系，我带你回顾一下。

如果这个函数 $f(x)=x^2$在区间0到1内连续(显然连续），且存在原函数$F(x)$， 那么：
$$I={\int }_a^{b}f(x)dx=F(b)-F(a)$$

我不讲解为什么这样计算，这个是定理，是学校知识，你不需要问到底，你现在应该能够天然地回忆起就是这么计算的。

那么相应的，我上面给出简单的积分公式，原函数$F(x)=\frac{1}{3}{x}^3$，因此：
$$I=F(1)-F(0)=\frac{1}{3}\cdot 1^3-\frac{1}{3}\cdot 0^3=\frac{1}{3}=0.3333$$

我们成功地算出来了这个函数在区间[0,1]的定积分是0.3333。

---

2.如果用蒙特卡洛积分方法去计算这个积分的话呢？

我现在同样要计算这个积分，但是这次用蒙特卡洛方法去计算。
$$I={\int }_0^1{x}^{2}dx$$

我先不解释蒙特卡洛方法是怎么样的，我先直接讲一下我现在怎么去算出来。

我其实能够算出来0到1区间中任意一个点的函数值，我现在通过采样的想法，在0到1区间中均匀地生成N个点，$x_1,x_2,...,x_N$，我计算他们的函数值的平均值：
$$\overline{f(x)}=\frac{1}{N}\sum _{i=1}^{N}f(x_i)dx$$

再乘以区间长度，将最终积分近似为：
$$I\approx （b-a)\times \overline{f(x)}$$

看起来有点像计算长方形面积了，这样算出来的结果是可以去近似$f(x)$ 的积分最终值的，当然采样数量N越大，结果约接近0.3333。

先说这种计算方式的缺点：

• 我要计算那么多个样本值，计算量好大，没有刚刚一步到位那么快
• 我算的结果还没刚刚算的那么准确，只能算出来一个近似值

那么假如我给的函数没有原函数的话呢，你怎么求它的定积分

例如，我给你一个没有原函数的$f(x)$:

$$I={\int }_0^1{e}^{-x^2}dx$$

你还怎么算它的定积分，但是用我的方法还是可以如法炮制地把近似的最终积分值算出来，对吧。

再例如，我给你一个二维积分和三维积分

$$I={\int }_0^1{\int }_0^1(x^2+y^2)dxdy$$

和
$$I={\int }_0^1{\int }_0^1{\int }_0^1(x^2+y^2+z^2)dxdydz$$
我让你算定积分，是不是看着头都大，让你用代码去算你都想直接放弃了。但是用我刚刚的方法，感觉还是能够把近似的结果给算出来。

显然这样算有什么优点：

• 代码实现挺方便的，用离散计算的方式就能算出一个近似该积分值
• 就算函数没有原函数，我也能算出来结果，相当于限制变小了
• 对高维积分的计算非常友好

这个计算方法其实就是用的蒙特卡洛积分，你可以稍微摸清楚蒙特卡洛用的什么想法去计算的。但是这里介绍的只是一个最简单的特例，我们接下来要对这个方法规范化，介绍蒙特卡洛积分的全貌。你或许好奇上面哪一部分是属于不规范的地方，那就是在做采样的那部分，我们很直观地做了均匀采样，但是我们要数学地表示我们怎么做均匀采样的，才能说是规范。

---

3.规范地表达我们均匀采样的方法

我们要求的积分是
$$I={\int }_a^{b}f(x)dx$$
这句话的意思是：
“对每一个$x$，把$f(x)$乘上一个很小的$dx$，加起来”。

我们上面采用的是均匀采样，很自然的采样方式，现在要数学地表示这种采样方式，就用到概率密度函数，我记为p(x)，那么我再[a,b]区间内均匀撒点，那么概率密度是:$$p(x)=\frac{1}{b-a}$$

我用两个例子让你理解这个$p(x)$的意思：

假如，积分区间是[0, 1]，那么$p(x)=1$，意味着在0至1之间采样的概率是1，或者说100%。

再假如，积分区间是[0, 2]， 那么算出来$p(x)=\frac{1}{2}$，意味着在0至1之间采样概率是0.5，或者说50%，而在1至2之间采样的概率是0.5，也是50%。

所以它不是「直接的概率」，而是「单位区间的概率密度」，反映的是在某一个点附近取到样本的可能性有多大。

因为$p(x)$的结果是一个不随$x$变化的量，所以不管你$x$是多少（即在采样哪个点），结果概率都是一样的，所以就是均匀采样。

回到我们上一节演示蒙特卡洛积分用到的例子：$$I={\int }_0^1{x}^{2}dx$$我强调这个例子很特殊，是因为区间是[0, 1]，我们屏蔽了概率密度函数。甚至我用了理所当然的方式解释了该积分的计算：
$$I\approx （b-a)\times \overline{f(x)}$$

我已经用概率密度函数去表达了均匀采样。我们现在开始规范一下，怎么把概率密度函数加入到积分中，或者怎么从蒙特卡洛积分的角度去理解上面这个积分的计算。

---

4.规范地表达蒙特卡洛积分

我们的概率密度函数，规定满足：
$${\int }_a^{b}p(x)dx=1$$
整个区间[a, b]是你可能取到的所有结果的集合，你不管怎么取点，都落在区间[a, b]内，那么所有可能发生的事情的概率总和，必须是100%。

这个规定是很自然的，不然就会出现「概率大于100%」或者「概率加起来小于100%」的奇怪现象。

那么我们可以把任何函数$f(x)$写成：
$$I={\int }_a^{b}f(x)dx=I={\int }_a^b\frac{f(x)}{p(x)}p(x)dx$$

我们从这个顺序开始看这条修改后的公式：

• $p(x)dx$是新的采样方式
• 这里$\frac{f(x)}{p(x)}$是调整后的函数，这个是要关键理解的地方
• 我是先按$p(x)$的方式采样，再用$\frac{f(x)}{p(x)}$调整回来！

先说新的采样方式，我们一直介绍的是均匀采样，那么我们是通过$p(x)$的概率去得到了$x_1,x_2,...,x_N$生成了这么多个均匀采样点。

那么如果我们不是在做均匀采样呢，而是在做非均匀采样呢，想象你有一块蛋糕，蛋糕有不同的甜度密度$p(x)$，你要评价这个蛋糕是什么味道的，但是你又喜欢吃甜的，你只能吃十口蛋糕，那么你会更多地选甜的地方吃，而味道淡的地方你会吃少一点。这就是非均匀采样，我们先不介绍非均匀采样对计算这个积分$f(x)$有什么助益，你现在只需要先知道我们的采样点可以不均匀地，而且采样概率密度还是能够用$p(x)$表示出来。我们不需要去想象这个$p(x)$长什么样子，只要知道我们是通过$p(x)$的概率去得到了$x_1,x_2,...,x_N$生成了这么多个非均匀采样点。

我们吃蛋糕时更多地挑在甜的地方吃，那么我们最后对蛋糕的评价也会认为这个蛋糕整体很甜，那么这个评价是有偏颇的，同样非均匀采样去对$f(x)$积分后结果也是有偏差的，所以我们要用补偿这个变化，必须要除以$p(x)$。

直白点说：本来一个$x$应该平均出现一次，结果因为你的$p(x)$偏好这个点，它出现了10次，那么这个点的贡献就被多算了10倍，所以要除回来，除以10，也就是除以$p(x)$。所以我们变成了对$\frac{f(x)}{p(x)}$积分。

$${\int }_a^b\frac{f(x)}{p(x)}p(x)dx$$
分成$p(x)dx$和$\frac{f(x)}{p(x)}$两部分去看，这就是蒙特卡洛积分的全貌！

关键思想：用「随机采样」代替「穷举积分」！

---

5.再次用蒙特卡洛积分方法去计算这个积分

我们用回第2节用到的例子：
$$I={\int }_0^1{x}^{2}dx$$

我们重新用规范的蒙特卡洛方法去算这个积分

1. 区间是[0, 1]，我们做均匀采样，那么$p(x)=1$
2. 我们用$p(x)dx$去采样，得到$x_1,x_2,...,x_N$这些采样点。
3. 而这个积分计算变成，因为$p(x)=1$：
   $$I=\frac{1}{N}\sum _{i=1}^N\frac{f(x_i)}{p(x_i)}=\frac{1}{N}\sum _{i=1}^{N}f(x_i)$$

你会思考，第2节中不是用到这个计算吗？
$$I\approx （b-a)\times \overline{f(x)}$$

其实$p(x)=\frac{1}{b-a}$，那么代入上面蒙特卡洛积分中的:
$$I=\frac{1}{N}\sum _{i=1}^N\frac{f(x_i)}{p(x_i)}=（b-a)\times \frac{1}{N}\sum _{i=1}^{N}f(x_i)$$

其实是一样的。

如果积分区间改为[0, 2]那么$p(x)=\frac{1}{2}$，积分值变成：
$$I=\frac{1}{N}\sum _{i=1}^N\frac{f(x_i)}{\frac{1}{2}}=2\times \frac{1}{N}\sum _{i=1}^{N}f(x_i)=(2-0)\times \frac{1}{N}\sum _{i=1}^{N}f(x_i)$$

---

6.重要性采样(非均匀采样)+蒙特卡洛积分

要彻底弄明白蒙特卡洛积分是怎么计算的，那就不能只看均匀采样，最重要还要看明白非均匀采样时是怎么计算的。

我们再来看重要性采样，从刚刚吃蛋糕的例子你也大概知道是怎么样的采样方式了。我们要从计算效率看这件事，用重要性采样是为了提高蒙特卡洛积分的效率，降低误差！

我们刚刚均匀采样，有个天然缺点:

如果函数 f(x) 在某些地方变化剧烈、在其他地方很平稳，均匀采样会很浪费！平稳地方，取很多样本也没啥用；剧烈地方，取样本太少，导致误差大！结果就是需要非常非常多的样本，才能得到准确的估计。

而重要性采样的核心是：「让重要的地方取更多点，不重要的地方取少点」

也就是说：根据$f(x)$的特点来定制采样方式！如果某个地方特别大或变化特别快，就多撒点。如果某个地方很小或很平滑，就少撒点。这样，有限的样本资源就能用在最有影响力的地方，整体估计误差就可以大大减少！（而且使用更少的样本数会减少计算资源耗费）

想象你要调查一个地方的收入水平：

• 如果你随机找人问，可能80%的人收入很普通，10%的人非常有钱，但你随机问的话，可能问不到有钱人！
• 而如果你「刻意」多去金融圈、商业区调查，就能更准确了解高收入群体对整体收入的贡献。这就相当于做了「重要性采样」！

现在我们要更聪明地去采样，我们先知道如果$p(x)$越接近$f(x)$，效果越好，如果是非均匀采样，那么要用$\frac{1}{p(x)}$来修正采样偏差，即上文蒙特卡洛积分中$\frac{f(x)}{p(x)}$，就是$f(x)$乘了修正项$\frac{1}{p(x)}$。

同样用回刚才的积分，看看用重要性采样会怎么计算，这回我用$p(x)=2x$来作为概率密度函数，这下我们可以看看图：

$f(x)$在靠近1的区域变化越剧烈，而黄色曲线是$p(x)$，能够让采样更集中在1的区域。

把$f(x)$和$p(x)$代入刚刚蒙特卡洛积分的公式中，实际结果都是相等的。（代码实现，结果是可以用更少的样本数量更快地得到结果），我直接说结果，分别用均匀采样和重要性采样去得到10000个样本去计算刚刚的积分，重要性采样的标准误差约小60%，意味同样的样本数，重要性采样给出的结果更稳定、更可靠。那更少的样本数就能够获得等用于均匀采样的计算效果，而更少的样本数在光线追踪中意味着更少的计算量，所以重要性采样是路径追踪中最核心的加速技术之一。

---

7.高维蒙特卡洛积分（二维和三维）

蒙特卡洛积分在处理高维积分是特别好用，刚刚大篇幅的介绍都是围绕一维展开的，现在要延伸至二维和三维的蒙特卡洛积分是怎么样的。

7.1 回顾一维蒙特卡洛积分

目标积分：
$$I={\int }_b^{a}f(x)dx\approx \frac{1}{N}\sum _{i=1}^N\frac{f(x_i)}{p(x_i)}$$
其中采样点用$p(x)$采样，且满足:
$${\int }_a^{b}p(x)dx=1$$

别忘了$N$是采样点个数，是根据$p(x)dx$得来的。

7.2 二维蒙特卡洛积分

二维即多了y轴的维度，顶点坐标变成了(x, y)，相应的重要性采样用的概率密度函数也是二维的。

选一个二维概率密度函数p(x, y)，满足
$$\int {\int }_{D}p(x,y)dxdy=1$$

这里的$D$当然是指(x, y)的区间，只是缩小写成$D$。

那么目标积分是二维的，写成：

$$I=\int {\int }_{D}f(x,y)dxdy$$

那么用蒙特卡洛积分：

$$I\approx \frac{1}{N}\sum _{i=1}^N\frac{f(x_i,y_i)}{p(x_i,y_i)}$$

举例（二维）：
假设$f(x)=x^2+y^2$ 在单位圆内$x^2+y^2<1$（即区间）

我们用极坐标$r,\theta$来表达$f(x,y)$和$p(r,\theta )$。而改为$f(r,\theta )=r^2$，接着可以用$p(r,\theta )=\frac{r}{\pi }$这个概率密度函数，在单位圆中采样$r,\theta$时常用的一个分布，我们不需要考虑为什么和曲线怎么样的，我们将这个现成的概率密度函数拿来即用。

然后积分为：
$$I\approx \frac{1}{N}\sum _{i=1}^N\frac{f(r_i,{\theta }_i)}{p(r_i,{\theta }_i)}$$

生成了采样点后，我们可以写成程序计算出结果。

7.3 三维蒙特卡洛积分

三维正是我们图形世界里面常用到的，是路径追踪中要用到的。三维顶点(x, y, z)对应我们的空间坐标点。

目标函数是：
$$I=\int \int {\int }_{V}f(x,y,z)dxdydz$$
这里的$V$是指(x, y, z)的区间。

而在路径追踪中，常用的是方向半球上的积分：

$$I={\int }_{\Omega }f(\omega )d\omega$$

这里$\Omega$指半球空间，$\omega$是一个方向(单位向量)，在$\Omega$这个半球表面上取得，也就是说积分区域是半球表面。

我们继续看路径追踪中该积分的计算，我们用到的概率密度函数是$p(\omega )$，同样先不管具体公式是怎么样的，直接看积分值的计算：

$$I\approx \frac{1}{N}\sum _{i=1}^N\frac{f({\omega }_i)}{p({\omega }_i)}$$

是不是发现如出一辙，实际上你不断往高维拓展也是一样的。你现在已经清楚高维蒙特卡洛积分的计算方式了。

---

8. 路径追踪中的蒙特卡洛积分

在实际应用场景中，$f(x)$和$p(x)$都不会简单到如上文般能够画出曲线图，所以我们的理解要摆脱曲线图。我们理解$p(x)$只需要理解它要尽量接近$f(x)$且满足积分等于1。当然不会他俩长什么样子，或者是几维的，计算起来也是如上般

说说路径追踪中的蒙特卡洛积分，关于蒙特卡洛积分是什么和用来做什么的，看完上文应该很清晰了。

我以问答的形式帮你梳理路径追踪中的蒙特卡洛积分：

问：为什么路径追踪中会用到蒙特卡洛积分？
答：因为路径追踪中出现了高维的积分计算，用蒙特卡洛积分是高效简单的计算方法，而且将积分计算转变为离散点的计算，用代码实现也很方便。

问：那么路径追踪在哪里用到了蒙特卡洛积分？
答：在渲染方程中，有一部分计算是积分计算。在光线碰撞点处，要在其半球区域上计算来自其他地方的光的贡献，这就是在半球区域上做三维的积分，那么就用蒙特卡洛积分转变成在半球区域上做采样。

我们看看路径追踪中使用的渲染方程：

我们先忽略每个符号具体所指，我们直接看该公式是结果 = A + B的形式。

一句话总结就是：出射光(量) = 自发光(量) + 来自其他方向的光贡献量。其中来自其他方向的光贡献量便是渲染方程中的积分部分，这部分我们便用蒙特卡洛积分去计算。

那么相当于我们现在的积分是
$$I={\int }_{\Omega }{f}_r(x,{\omega }_i,{\omega }_o)\cdot L_i(x,{\omega }_i)\cdot ({\omega }_i,n)d{\omega }_i$$

我们选择$p(\omega )$为：

• 均匀半球采样： $p(\omega )=\frac{1}{2\pi }$
• 非均匀采样（常用的余弦加权采样）：
  $$p(\omega )=\frac{\cos \theta }{\pi }$$

当面$p(\omega )$是二选一，那么积分结果为：
$$I\approx \frac{1}{N}\sum _{i=1}^N\frac{f_r(x,{\omega }_i,{\omega }_o)\cdot L_i(x,{\omega }_i)\cdot ({\omega }_i,n)}{p({\omega }_i)}$$

看到这里，你应该了解清楚了蒙特卡洛积分在路径追踪中是怎么用的了。你也可以止步于此，看路径追踪的内容也可以很好理解积分部分的计算了。

---

9.More：在深入地了解蒙特卡洛积分在路径追踪中怎么计算

我不打算展开介绍渲染方程，但是继续往下看需要了解渲染方程，这一节是介绍路径追踪渲染方程中最关键的部分——递归部分。

先回顾一下渲染方程：

其中：

• $L_o$：出射光部分（要求的结果）
• $L_i$：入射光部分（未知量）
• $f_r$：材质反射函数（已知量——根据采样量得到）
• $w_i$：来自哪个方向的光向量（采样量）
• $w_o$：出射光向量（已知量）
• $x$：光线碰撞点（已知量）

关键点：
问题关键出现在$L_i(x,{\omega }_i)$这里了，它表示从方向$w_i$过来并命中了$x$处的光，所以它本身又是另一个碰撞点$x‘$处的出射光，即$L_o(x‘,-{\omega }_i)$。

所以，从数学上看，渲染方程是递归的，我想知道$L_i$，那就的回溯到它上一个碰撞顶点计算它的出射光量，这个点在计算时，我还得递归去往回找，直至找到光源。那从数学上就在不断嵌套积分公式。

这个渲染方程的递归展开，用数学公式表达就是：
$$L_o=L_e+{\int }_{\Omega }{f}_r\cdot (L_e+{\int }_{\Omega }{f}_r\cdot (L_e+\int ...))$$
实现起来基本无法解了。

而用蒙特卡洛积分，则不需要这样显示去嵌套积分。还记得蒙特卡洛是通过计算多个采样值后平均，从而近似积分吧。蒙特卡洛没法帮助我们逃脱不了递归，但是可以摆脱这个积分。

路径追踪中还是需要递归，但是蒙特卡洛积分让我们可以在不求解积分的前提下，间接估计出积分值，从而帮助我们从无法解变成可以解。

是的，计算量依然很大，你把上面积分公式换成用蒙特卡洛积分层层嵌套，就是实际计算方式了，手写会很复杂，换成代码就是写成递归函数去计算。这正是路径追踪计算量大的核心原因之一。

现在举例介绍一下两次递归用蒙特卡洛积分怎么表示出来：

1. 第一次碰撞： 对渲染方程的积分部分使用蒙特卡洛估计，只使用1个采样方向，即N = 1：
   
   采样光线方向是“根据材质的 概率密度函数来随机生成”的。另外我们知道：
   $$L_i(x,{\omega }_1)=L_o(x_1,-{\omega }_1)$$

2. 那么单独看第二次碰撞的展开：
   

3. 把第二次碰撞的结果嵌套进第一次碰撞的公式中：
   

这就是路径追踪中「两次采样」的完整公式结构！这计算量别手工推了，还是交给程序去计算吧！

• 注意N = 1这一点，在实际计算时都是采用N = 1，即一次采样，避免递归时指数爆炸。

---

总结：

从一维出发解释蒙特卡洛积分，再延申至三维的理解。讲解了路径追踪中的蒙特卡洛积分用途，并用两次递归的公式结构展示了一下计算量。希望讲解能够面面俱到，认真看完能对蒙特卡洛积分彻底清晰。如果内容有疏漏或者错误的地方，可以给我留言私信改进一下。