本项目参考自《Ray Tracing in One Weekend》系列。
项目链接:https://github.com/maijiaquan/ray-tracing-with-imgui
目录:
《用两天学习光线追踪》1.项目介绍和ppm图片输出
《用两天学习光线追踪》2.射线、简单相机和背景输出
《用两天学习光线追踪》3.球体和表面法向量
《用两天学习光线追踪》4.封装成类
《用两天学习光线追踪》5.抗锯齿
《用两天学习光线追踪》6.漫反射材质
《用两天学习光线追踪》7.反射向量和金属材质
《用两天学习光线追踪》8.折射向量和电介质
《用两天学习光线追踪》9.可放置相机
《用两天学习光线追踪》10.散焦模糊
《用一周学习光线追踪》1.动态模糊
《用一周学习光线追踪》2.BVH树、AABB相交检测
《用一周学习光线追踪》3.纯色纹理和棋盘纹理
《用一周学习光线追踪》4.柏林噪声
《用一周学习光线追踪》5.球面纹理贴图
《用一周学习光线追踪》6.光照和轴对齐矩形
《用一周学习光线追踪》7.长方体和平移旋转
《用余生学习光线追踪》1.蒙特卡洛积分与概率密度函数PDF
《用余生学习光线追踪》2.光照散射与重要性采样
《用余生学习光线追踪》3.生成球坐标随机方向
离散型随机变量
如果随机变量X只可能取有限个或至多可列个值,则称X为离散型随机变量。例如:
X
=
{
1
,
硬
币
正
面
0
,
硬
币
反
面
X =\left\{ \begin{aligned} 1, 硬币正面 \\ 0, 硬币反面 \end{aligned} \right.
X={1,硬币正面0,硬币反面
X = { 1 , 骰 子 点 数 为 1 2 , 骰 子 点 数 为 2 3 , 骰 子 点 数 为 3 4 , 骰 子 点 数 为 4 5 , 骰 子 点 数 为 5 6 , 骰 子 点 数 为 6 X =\left\{ \begin{aligned} 1, 骰子点数为1 \\ 2, 骰子点数为2 \\ 3, 骰子点数为3 \\ 4, 骰子点数为4 \\ 5, 骰子点数为5 \\ 6, 骰子点数为6 \end{aligned} \right. X=⎩⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎧1,骰子点数为12,骰子点数为23,骰子点数为34,骰子点数为45,骰子点数为56,骰子点数为6
PMF, PDF, CDF 的区别
概率质量函数(probability mass function,简写为PMF),是离散随机变量在各特定取值上的概率。
概率密度函数(probability distribution function,简写为PDF) 定义了一个连续的概率分布,例如正态分布函数。
维基百科中PMF和PDF的区别:
PMF是对离散随机变量定义的,本身代表该值的概率。
PDF是对连续随机变量定义的,本身不是概率,只有对连续随机变量的PDF在某区间内进行积分后才是概率。
PDF可以求一个区间内的概率:
P
r
(
a
≤
X
≤
b
)
=
∫
a
b
p
d
f
(
x
)
d
x
.
Pr( a \le X \le b) = \int_a^b pdf(x)\:dx.
Pr(a≤X≤b)=∫abpdf(x)dx.
PDF最重要的特点就是,整个概率区间的概率之和必定为1,即:
∫ − ∞ ∞ p d f ( x ) d x = 1. \int_{-\infty}^{\infty} pdf(x) dx = 1. ∫−∞∞pdf(x)dx=1.
下图展示了标准正态分布函数,这是高斯分布函数的一个特例。

CDF (Cumulative Distribution Function) 累积分布函数

CDF将PDF的当前样本及其之前的所有样本都累加了起来。
因为PDF是一个连续函数,所以CDF也是一个连续的函数。
不难发现:
1.CDF是PDF的积分,积分区间为 − ∞ -\infty −∞到 x x x
c d f ( x ) = ∫ − ∞ x p d f ( x ) d x cdf(x) = \int_{-\infty}^{x} pdf(x) dx cdf(x)=∫−∞xpdf(x)dx
2.PDF是CDF的导数
p d f ( x ) = d d x c d f ( x ) pdf(x) = { d \over {dx} } cdf(x) pdf(x)=dxdcdf(x)
用均匀函数生成任意PDF的随机值
PDF可能是一个正态分布函数,也可能是其他形态。
我们的目标就是要使用C++提供的均匀分布随机数的drand48,根据给定的任意一个PDF,产生符合这个PDF形态的随机值。例如正态分布,应该会有更多的样本落在中部。
先上结论:
如果要根据给定的PDF,生成满足该PDF分布的随机值,则可以对PDF求原函数得到CDF,再求CDF的反函数,该函数即为满足PDF概率分布的随机数函数。
这个结论,稍后在代码例2和例3中会用到。
蒙特卡洛积分简介
假设被积函数 f ( x ) f(x) f(x)图像如下,积分区间为 [ a , b ] [a,b] [a,b],我们要求的积分为 F = ∫ a b f ( x ) d x F = \int_a^b f(x)\;dx F=∫abf(x)dx:

( b − a ) f ( x ) (b-a)f(x) (b−a)f(x)几何上可表示为一个宽为 ( b − a ) (b-a) (b−a),高为 f ( x ) f(x) f(x)的矩形面积:

下图中,假设我们采样4次,则相当于计算了4个矩形的面积,并求平均值:

定义:
⟨
F
N
⟩
=
1
N
∑
i
=
0
N
−
1
(
b
−
a
)
f
(
x
)
\langle F^N\rangle = \dfrac{1}{N } \sum_{i=0}^{N-1} (b-a)f(x)
⟨FN⟩=N1i=0∑N−1(b−a)f(x)
变形得:
⟨
F
N
⟩
=
(
b
−
a
)
1
N
∑
i
=
0
N
−
1
f
(
X
i
)
(1)
\langle F^N\rangle = (b-a) \dfrac{1}{N } \sum_{i=0}^{N-1} f(X_i) \tag{1}
⟨FN⟩=(b−a)N1i=0∑N−1f(Xi)(1)
其中, ⟨ F N ⟩ \langle F^N \rangle ⟨FN⟩为样本数量为 N N N时的 F F F的平均值,也就是积分结果的近似值,以上公式叫做Basic Monte Carlo Estimator(基本蒙特卡洛估算器),这个公示在下面的代码例1中会用到。
由于随机数是等概率均匀产生的,上述方程的pdf为
p
d
f
(
X
i
)
=
1
(
b
−
a
)
pdf(X_i) = \frac{1}{(b - a)}
pdf(Xi)=(b−a)1
上述公式(1)仅在随机数是等概率均匀生成的情况下才生效。
如果将蒙特卡洛估算器扩展到任意的PDF,一个更加通用的
⟨
F
N
⟩
\langle F^N \rangle
⟨FN⟩可以写成:
⟨
F
N
⟩
=
1
N
∑
i
=
0
N
−
1
f
(
X
i
)
p
d
f
(
X
i
)
(2)
\langle F^N \rangle = \dfrac{1}{N} \sum_{i=0}^{N-1} \dfrac{f(X_i)}{pdf(X_i)} \tag{2}
⟨FN⟩=N1i=0∑N−1pdf(Xi)f(Xi)(2)
除以PDF,是为了抵消掉不均匀带来的影响。例如某个区域会有很多样本,那就除以该位置的PDF,以削弱其权重。
用PDF求解蒙特卡洛积分
在没有学习公式一之前我们用积分可以用来算面积的思想来用蒙特卡洛算法求解积分。下设我们要求以下积分:
I
=
∫
0
2
x
2
d
x
I = \int_{0}^{2} x^2 dx
I=∫02x2dx
本质上就是求0到2区间内的函数曲线下方和x轴所围成的面积,即:
I
=
a
r
e
a
(
x
2
,
0
,
2
)
I = area( x^2, 0, 2 )
I=area(x2,0,2)
本质上,也是求区间的长度和区间内函数平均值的乘积:
I
=
2
⋅
a
v
e
r
a
g
e
(
x
2
,
0
,
2
)
I = 2 \cdot average(x^2, 0, 2)
I=2⋅average(x2,0,2)
其实这就是公式一:
⟨
F
N
⟩
=
(
b
−
a
)
1
N
∑
i
=
0
N
−
1
f
(
X
i
)
\langle F^N\rangle = (b-a) \dfrac{1}{N } \sum_{i=0}^{N-1} f(X_i)
⟨FN⟩=(b−a)N1i=0∑N−1f(Xi)
即
⟨
F
N
⟩
=
2
1
N
∑
i
=
0
N
−
1
x
2
\langle F^N\rangle = 2 \dfrac{1}{N } \sum_{i=0}^{N-1} x^2
⟨FN⟩=2N1i=0∑N−1x2
写成代码如下:
//例1
int main() {
int inside_circle = 0;
int inside_circle_stratified = 0;
int N = 1000000;
float sum = 0;
for (int i = 0; i < N; i++) {
float x = 2*random_double();
sum += x*x;
}
std::cout << "I =" << 2*sum/N << "\n"; //输出得到 I =2.66534
}
上述代码中,我们通过暴力增加样本数量的方式,使得平均值尽量接近真实答案。但是如果要结果误差尽量小,就需要很多样本,才能够模拟出正确的结果。接下来,问题变成了能否用最少的样本,模拟出最接近真实值的结果。
第二种方式:借助PDF。
先人为给
f
(
x
)
=
x
2
f(x)= x^2
f(x)=x2 选择一个PDF,例如取
p
(
r
)
=
r
/
2
p(r) = r/2
p(r)=r/2,则可以用CDF反函数的方式来生成满足PDF概率分布的随机数。
根据PDF为
p
(
r
)
=
r
/
2
p(r) = r/2
p(r)=r/2,CDF为
P
(
x
)
=
x
2
4
P(x) = \frac{x^2}{4}
P(x)=4x2
如果x是均匀生成的[0,1]区间内的随机数,求CDF的反函数,就可得到符合PDF的概率分布的随机数生成函数:
e
=
4
∗
randomdouble
(
)
e = \sqrt{4*\text{randomdouble}()}
e=4∗randomdouble()
接下来套用公式二:
⟨
F
N
⟩
=
1
N
∑
i
=
0
N
−
1
f
(
X
i
)
p
d
f
(
X
i
)
(2)
\langle F^N \rangle = \dfrac{1}{N} \sum_{i=0}^{N-1} \dfrac{f(X_i)}{pdf(X_i)} \tag{2}
⟨FN⟩=N1i=0∑N−1pdf(Xi)f(Xi)(2)
写成代码如下:
//例2
inline float pdf(float x) {
return 0.5*x;
}
int main() {
int inside_circle = 0;
int inside_circle_stratified = 0;
int N = 1000000;
float sum = 0;
for (int i = 0; i < N; i++) {
float x = sqrt(4*random_double());
sum += x*x / pdf(x);
}
std::cout << "I =" << sum/N << "\n";
}
上述pdf似乎并不能帮我们快速求解积分。接下来我们人为定义一个更加合适的PDF:
p ( x ) = 3 8 x 2 p(x) = \frac{3}{8}x^2 p(x)=83x2
对应有CDF为
P
(
x
)
=
x
3
8
P(x) = \frac{x^3}{8}
P(x)=8x3
CDF的反函数为:
P
−
1
(
x
)
=
8
x
1
3
P^{-1}(x) = 8x^\frac{1}{3}
P−1(x)=8x31
这是一个堪称完美的pdf,因为我们已经预先算出了 I = ∫ 0 2 x 2 d x I = \int_{0}^{2} x^2 dx I=∫02x2dx 的值是8/3,从而使得一次蒙特卡洛积分采样就能给出结果:
inline float pdf(float x) {
return 3*x*x/8;
}
int main() {
int inside_circle = 0;
int inside_circle_stratified = 0;
int N = 1;
float sum;
for (int i = 0; i < N; i++) {
float x = pow(8*random_double(), 1./3.);
sum += x*x / pdf(x);
}
std::cout << "I =" << sum/N << "\n";
}
选择一个非均匀的pdf,使得能够在被积函数较大的地方采样更多的样本,这样能够减少噪声、快速收敛,这叫做重要性采样。
综上,使用蒙特卡洛算法求解积分的步骤如下:
1.找到在区间 [ a , b ] [a,b] [a,b]内被积函数 f ( x ) f(x) f(x)
2.在区间 [ a , b ] [a,b] [a,b]内挑选一个合适的非零的概率密度函数 p p p
3.生成一堆满足上述概率密度函数 p p p分布的随机数 r r r,从而得到一堆 f ( r ) / p ( r ) f(r)/p(r) f(r)/p(r),求这堆 f ( r ) / p ( r ) f(r)/p(r) f(r)/p(r)的平均值
参考资料: