3.1 Ray-Triangle and Ray-bounding box Intersection

1 蒙特卡洛基础

假设要求解积分 $\int f(x) \space dx$，往往没有解析解，直接求解困难。蒙特卡洛引入随机化，进而对积分进行估计。选择一个自变量 $X \sim P(x)$，概率密度为 $p(x)$，只要随机变量定义域能够覆盖被积函数定义域，就可以有如下恒等变换：
$$
\int f(x) \space dx = \int \frac{f(x)}{p(x)}\cdot p(x) \space dx=E\left[\frac{f(x)}{p(x)}\right]
$$
因此可以将求积分转为了求期望 $E\left[\frac{f(x)}{p(x)}\right]$。对于期望，可以使用其无偏估计(样本均值)，进而得到积分估计。对分布 $X \sim p(x)$ 进行多次采样，得到 $X_1, X_2, …, X_N$。使用这些采样点可以得到积分的无偏估计：
$$
E\left[\frac{f(x)}{p(x)}\right] \approx \frac{1}{N}\sum\limits_{i=1}^{N}\frac{f(X_i)}{p(X_i)}
$$

2 提高采样效率

蒙特卡洛积分理论建立的无偏估计能够通过提高采样数量而降低方差，即图像中的噪声。但实际中的样本数量是受限的，因此降低方差的唯一途径是最大化现有样本的利用率。

2.1 Stratified Sampling

分层采样是将积分区域划分为多个子区域，并在每个子区域分别放置样本，以这种方式对样本位置精心布局，更好地捕捉被积函数特征，从而降低方差。

2.2 低差异序列

Hammersley 和 Halton 并非随机采样生成的序列，而是通过数值计算得到的 low-discrepancy(低偏差) 点集。这两种序列生成原理相同，都是通过 radical inverse 操作，将一个输入值转为 low-discrepancy 点集中的一个点。这个过程不涉及随机性，相同输入得到相同的输出，但最终的点集整体表现出非均匀且 low-discrepancy 的分布特征，也就是说，具有无序性，同时任意相邻点不会相距太近或太远。

2.2.1 Radical Inverse

对于一个正整数 $a$ 可以使用一系列数字表示为以 $b$ 为基底的形式 $d_m(a)\cdots d_2(a)d_1(a)$，如下式
$$
a = \sum\limits^m_{i=1}d_i(a)\cdot b^{i-1} \space,\quad \mathrm{where} \space d_i(a) \in [0,b-1]
$$
Radical Inverse 函数 $\Phi_b$ 则通过将 $a$ 的每一位数字转换到小数点后，将 $a$ 转换到 $[0,1)$ 区间。$a$ 转换后的形式如下
$$
\Phi_b(a)=0.d_1(a)d_2(a)\cdots d_m(a)
$$
其数学计算形式为
$$
\Phi_b(a)=\sum\limits^m_{i=1} d_i(a) \cdot b^{-i}
$$
基底以 $2$ 为例，按照上述方式计算 $\Phi_2(a)$，即将 $a$ 进行 radical inverse，得到下表结果

$a$	Base 2	$\Phi_2(a)$
0	0	0
1	1	0.1=1/2
2	10	0.01=1/4
3	11	0.11=3/4
4	100	0.001=1/8
5	101	0.101=5/8
$\vdots$	$\vdots$	$\vdots$

可以看出 $\Phi_2(n)$ 是对 $\Phi_2(n-1)$ 划分后的区域再次进行均匀划分，这是 Radical Inverse 在基底 $2$ 的表现。

2.2.2 Radical Inverse 生成多维样本序列

假设样本序列中一个样本点的索引 $a$，对 $a$ 分别使用前 $n$ 个质数为基底，进行 Radical Inverse，则得到一个 $n$ 维样本点。假设前 $n$ 个质数为 $(p_1,\cdots,p_n)$，

2.2.2.1 Halton Sequence

对于样本总数未知的情况下，使用 Halton Sequence，则 $n$ 维样本为
$$
x_a = (\Phi_2(a),\Phi_3(a),\Phi_5(a),\cdots,\Phi_{p_n}(a))
$$
当样本总数为 $\prod p_i$ 的次幂时，Halton 序列达到最低差异性。

2.2.2.2 Hammersley Sequence

当样本总数已知时，Hammersley 序列差异性更小，假设样本总数为 $N$，那么 $n+1$ 维样本为
$$
x_a=(\frac{a}{N},\Phi_2(a),\Phi_3(a),\Phi_5(a),\cdots,\Phi_{p_n}(a))
$$

2.2.2.3 基底为 2 的高效实现

Radical Inverse 过程会将数字位序反转，即原位序 $d_m(a)\cdots d_2(a)d_1(a)$ 变为 $d_1(a)d_2(a)\cdots d_m(a)$ 。因此对于基底为 2 的 Radical Inverse，可以先使用位操作将位序反转，再转为小数。对于 32 为无符号整数 $n$ ，其位序反转过程为：

uint32_t ReverseBits32(uint32_t n){
    n = (n << 16) | (n >> 16);	// 将高 16 位与低 16 位对换
    // 分为两个16位，分别在每个16位中进行高8位与低8位对换
    n = ((n & 0x00ff00ff) << 8) | ((n & 0xff00ff00) >> 8);
    // 分为四个8位，分别在每个8位中进行高4位与低4位对换
	n = ((n & 0x0f0f0f0f) << 4) | ((n & 0xf0f0f0f0) >> 4);
    // 分为八个4位，分别在每个4位中进行高2位与低2位对换
	n = ((n & 0x33333333) << 2) | ((n & 0xcccccccc) >> 2);
    // 分为十六个2位，分别在每个2位中进行高1位与低1位对换
	n = ((n & 0x55555555) << 1) | ((n & 0xaaaaaaaa) >> 1);
	return n;
}

得到位序反转结果，乘以 0x1p-32 即可变为 $[0,1)$ 的点。