TC130：游戏渲染进阶​zhuanlan.zhihu.com

蒙特卡洛积分是图形学中常用的数学工具, 这里就来总结下蒙特卡洛积分的原理和使用方式. 很多教程中把概率分布和积分是混在一起讲的, 个人觉得分开讲比较合适. 这篇文章就先来讲下概率分布变换和随机采样的部分.

概率论基础

这里快速回顾下概率论的基础, 这里不会特别深入精确地描述. 需要的朋友可以参考概率论相关的教材.

设

是随机变量,

是任意实数, 称函数

为随机变量

的分布函数/CDF, 或称

服从

, 记为

.

分布函数满足:单调不减;

右连续, 即

;

.

从分布函数求概率:

(1) 离散型变量

随机变量

只能取有限个可能的值, 称

为离散型随机变量, 称

为

的概率分布, 记为

.

可以用矩阵形式表示为:

.

离散型随机变量的概率分布满足

.

比如记掷骰子的点数为

, 得

的分布函数为

(2) 连续型变量

如果随机变量

的分布函数可以表示为:

称

为连续型随机变量, 称

为

的概率密度函数/概率密度/PDF, 记为

.

概率密度函数满足

.

对任意实数

, 有

.

取在某个区间的概率为:

比如现在假设

是一个均匀随机分布在

上的连续随机变量, 得

(3) 如果

是定义在样本空间上n个随机变量, 则称

为n维随机变量.

n维随机变量的分布函数定义为

n维随机变量的性质和上面类似, 这里不再一一描述.

随机值采样

在计算机中, 得到一个均匀随机分布在

上的随机数是很简单的, 我们这里用

表示服从

均匀分布的随机变量. 现在我们就来用

来得到我们想要的服从特定概率分布的随机变量.

(1) 离散型随机变量

对于离散型随机变量, 计算过程比较简单, 已知

, 假设从

推导出

的函数为

.

考虑到

在

上均匀分布, 只需要将

依概率映射到

样本空间中每个值即可. 得到

比如现在要得到

易得

(2) 连续型随机变量

比如现在要得到概率密度为

, 概率分布为

的随机变量

. 设变换函数为

, 即

. 为了下面计算方便, 我们先假设

是一个单调递增函数.

由概率分布定义可知:

已知

是单调递增函数, 可得:

已知

在

上均匀分布, 可得

可以得到

互为反函数, 即

我们平时遇到的概率分布函数都是不满足单调递增的, 只需要去掉概率密度为0的部分即可.

现在举两个例子:

A. 次方分布

设

在

上服从n次方分布, 即概率密度满足

, 设

, 由概率密度性质可知

可以解得:

由此算出

的概率分布函数:

将

限制在

上, 可得

B. 指数分布

设

的概率密度满足

.

推导的部分和上面相同, 可得:

这里的

的概率分布和

是相同的, 因此也可以写成

(3) 拒绝式随机

对于一些无法求出解析解的概率分布函数, 或者无法得到

反函数的概率分布函数, 可以用拒绝式随机方法.

假设我们现在想得到一个概率密度为

的随机变量

.

现在我们已有一个概率密度为

的随机变量 , 我们可以任意次得到一个服从

分布的随机变量

, 且其概率密度满足

.

这样, 我们就可以通过下面的方法来随机得到

:取一个服从

分布的随机变量值

;

在

上随机得到一个变量值

;

如果

, 则该次随机结果被接受, 返回

. 否者该次随机被拒绝, 重新执行第一步.

拒绝式随机方法的效率取决于

和

之间的贴合程度, 如果二者之间空隙很大, 就可能需要多次随机, 效率会比较低.

一个常见的拒绝式随机法的应用场景就是随机在一个圆中取一个点, 大致过程为:随机在单位圆中取一点

point p;

do {

p.x = rand() * 2 - 1;

p.y = rand() * 2 - 1;

} while(p.x * p.x + p.y + p.y > 1);

return p;

概率分布变换

现在已知一个随机变量

的概率密度为

, 现在我们令

, 现在我们要尝试求出

的概率密度函数

. 为了计算方便, 我们只考虑函数

是严格单调递增的情况, 平时我们需要求解的函数大部分都是满足严格单调递增的.

由概率分布函数定义可知:

对两边一起求导得:

这样, 我们成功计算出了

的概率密度函数.

比如现在有

, 令

, 可算出

的概率密度为:

现在, 让我们来考虑多维随机变量, 设

和

都是n维的随机变量,

和

之间的转换关系为

.

为矩阵函数, 即

,

.

可以推导得出:

表示

的雅可比矩阵的行列式的绝对值,

的雅可比矩阵为:

现在来看下实际应用的例子:

A. 极坐标系

极坐标系的变换为

假设我们现在已知关于极坐标的概率密度函数

, 现在来计算直角坐标系的概率密度.

对应的雅可比矩阵为:

求得行列式值为

. 这样, 我们得到两种坐标系之间的变换公式为:

B. 球坐标系

球坐标系到直角坐标系变换为:

可解得雅可比矩阵行列式值为

, 相应的概率密度为:

现在来考虑在单位球面上的情况. 在球坐标系中, 我们从立体角的定义可以得到:

立体角在某个

范围内的概率为:

得到概率密度的转换为:

二维随机变量采样

现在可以来尝试从二维随机变量中采样.

(1) 联合概率密度

在开始之前, 我们还需要来简单回顾下联合概率密度的概念.

设现在有二维连续型随机变量

, 二维随机变量的联合概率密度为

,

的联合分布函数为

二维随机变量的概率密度满足

的边缘概率密度为:

在

的条件下,

的条件概率密度为:

(2) 单位半球面采样

在单位半球面上均匀采样时, 每个立体角上都是等可能的. 由此得关于立体角的概率密度

是常数, 令其为

, 得

解得

, 由前面得到的结论可知

.

先来计算

, 得到

的边缘概率密度为:

再得到

的条件概率密度为:

的概率密度在

确定时是固定的, 这和我们的直觉是相同的. 接下来来计算相应的概率分布函数:

求相应的反函数, 并将

替换为

, 得到:

将结果用直角坐标系来表示:

(3) 随机单位球面采样

推导过程和上面的几乎一模一样, 这里不再赘述. 最终结果为:

(4) 随机单位圆采样

一个常见的错误是随机取半径, 随机取角度, 使用

来采样. 这样得到的结果会使得在圆的中心区域概率比边缘部分要高.左边是错误的采样方式得到的分布, 右边是正确的方式得到的分布

在单位圆上均匀采样时, 关于面积的概率密

是个常数, 可解得

. 转换为极坐标系下得表示为

. 使用和前面一样得推导过程得:

在

确定时, 因为圆的对称性,

是个固定常数. 进一步计算分布函数并取反函数可求得:

另外一种方式是使用正方形随机采样, 然后同心映射到圆上.

其中一个1/8部分的映射公式为:

另外七个部分的映射公式可用相似的方式得到.

(5) 单位半球面余弦权重采样

求解图形学中的渲染方程时, 许多BRDF方程都是和夹角余弦相关的, 因此按照余弦采样是很有必要的. 即

, 求解概率密度为:

这样, 我们就可以继续使用上面的方式来推导出结果.

不过这里要介绍下Malley方法的实现, Malley方法就是先在单位圆上随机采样, 然后将单位圆上的点作为半球面上点的投影, 来得到半球面上的点.

下面我们来验证一下这种方式的正确性:Malley方法得到半球面上按余弦权重采样点

已知单位圆上随机采样的点极坐标为

, 概率密度为

. 单位半球面上对应的点极坐标系为

, 两个坐标的关联为

. 这样得到雅可比矩阵为:

行列式值为

, 变换概率密度分布得:

刚好符合上面我们想要得概率密度函数, 这样就可以从单位圆采样得到单位半球面上得采样.

其余的在锥形区域, 三角形, 长方形区域随机采样的过程和结果都是类似的, 这里不再给出. 这样我们可以随意按照自己想要的概率密度进行随机数采样, 下一篇会讲述如何使用随机数来实现蒙特卡洛积分.