前言

开年第一篇博客~ 整理了三四个小时才整理完orz。
这一部分是光线追踪三部曲的最后一部，主要介绍了蒙特卡洛积分、重要性采样等内容。场景上没有什么大的改变，基本上就是在Cornell Box中渲染的，本篇主要在加速收敛，提升渲染效率上下功夫。
参考链接：https://raytracing.github.io/books/RayTracingTheRestOfYourLife.html

什么是蒙特卡洛积分

蒙特卡洛方法是一类通过随机采样来求解问题的算法的统称，要求解的问题是某随机事件的概率或某随机变量的期望。通过随机抽样的方法，以随机事件出现的频率估计其概率，并将其作为问题的解。

蒙特卡洛的基本做法是通过大量重复试验，通过统计频率，来估计频率，从而得到问题的求解。举个例子，如下图所示，矩形内不规则图案的面积$A_{shape}$可由$A\ast p$得出，$p$是矩形内的点在不规则图案中的概率。

现在重复往矩形范围内随机投射点，样本点有一定概率会落在不规则图形内，重复n次试验，落在不规则图形内的次数为k，则频率为k/n，若样本数量较大，根据伯努利大数定律，则有：
$$p=\frac{A_{shape}}{A}\approx \frac{k}{n}$$
由此，$A_{shape}$的值可由$\frac{kA}{n}$估计得出。

用蒙特卡洛方法来计算函数积分，这里给出一般的定义，设$X_1,X_2,...,X_n$是独立同分布的一组样本，它们的概率密度函数为$pdf(x)$，则函数的积分可以表示为：
$$F(X)=\frac{1}{n}\sum _{k=1}^n\frac{f(X_k)}{pdf(X_k)}$$
这就是蒙特卡洛积分的一般等式，其中除以概率密度函数的做法，可以理解为是对样本的统计处理。我们可以证明蒙特卡洛法的积分估计量的正确性：
$$\begin{gathered} E(F(X))=E[\frac{1}{n}\sum _{k=1}^n\frac{f(X_k)}{pdf(X_k)} \\ =\frac{1}{n}\sum _{k=1}^n\int \frac{f(x)}{pdf(x)}\cdot pdf(x)dx \\ =\frac{1}{n}\sum _{k=1}^n\int f(x)dx \\ =\int f(x)dx \end{gathered}$$
蒙特卡洛法的积分估计值得数学期望等于被积函数的积分真值，证明$F(X)$是无偏估计量。

什么是重要性采样

采样函数$pdf(x)$的变化趋势与$f(x)$越相近，蒙特卡洛积分的收敛速度就会越快。因此，我们可以把采样点集中在被积函数数值较高的部分，以及来得到更准确更有效率的采样估算。这种非均匀的采样方式被称为重要性采样。

正文记录

接下来，我将顺序罗列正文中我认为重要的部分。

一个简单的蒙特卡洛程序

在做采样时，我们一般可以采用两种随机算法：Monte Carlo以及Las Vegas。

Las Vegas算法可以保证得到正确的结果，但是时间可能会很长，不能保证什么时候可以获得结果。一段采用Las Vegas随机的代码片段如下所示，在一个单位球内随机选取一个点：

inline vec3 random_in_unit_sphere() {
    while (true) {
        auto p = vec3::random(-1,1);
        if (p.length_squared() < 1)
            return p;
    }
}

用Monte Carlo算法可能会得到正确的结果，但是也可能出错，出错时我们将该样本点剔除就好，时间上不需要做过多假设。一个Monte Carlo算法的例子，估计$\pi$的值：

int main() {
    int N = 100000;
    int inside_circle = 0;
    for (int i = 0; i < N; i++) {
        auto x = random_double(-1,1);
        auto y = random_double(-1,1);
        if (x*x + y*y < 1)
            inside_circle++;
    }
    std::cout << std::fixed << std::setprecision(12);
    std::cout << "Estimate of Pi = " << (4.0 * inside_circle) / N << '\n';
}

Monte Carlo算法一个糟糕的点是，随着采样进行，每一个样本对结果的帮助都小于上一个样本。我们可以通过一种称为分层样本（Jittering）的方式来改进采样。分层法在采样时会在细分的一个个小网格里采样，进行了一部分约束，示意图如下：

用分层的方式估计$\pi$值与不用分层方式进行对比，代码和结果如下：

int main() {
    int inside_circle = 0;
    int inside_circle_stratified = 0;
    int sqrt_N = 1000;
    for (int i = 0; i < sqrt_N; i++) {
        for (int j = 0; j < sqrt_N; j++) {
            auto x = random_double(-1,1);
            auto y = random_double(-1,1);
            if (x*x + y*y < 1)
                inside_circle++;
            x = 2*((i + random_double()) / sqrt_N) - 1;
            y = 2*((j + random_double()) / sqrt_N) - 1;
            if (x*x + y*y < 1)
                inside_circle_stratified++;
        }
    }

    std::cout << std::fixed << std::setprecision(12);
    std::cout
        << "普通采样估计的Pi值 = "
        << (4.0 * inside_circle) / (sqrt_N*sqrt_N) << '\n'
        << "分层采样估计的Pi值 = "
        << (4.0 * inside_circle_stratified) / (sqrt_N*sqrt_N) << '\n';
}

结果：

普通采样估计的Pi值 = 3.143392000000

分层采样估计的Pi值 = 3.141532000000

分层采样可以以更好的渐进率收敛。不幸的是，这种优势随着问题维度的增大而减弱（例如，用3D球体体积与正方体体积之比来估计$\pi$的值，两种方法的差距会变小）。这就是所谓的维度诅咒。

以下是采用分层采样与不进行分层采样渲染Cornell Box场景的对比：

上图为不进行分层采样，下图为进行分层采样，这里的分层采样只在光线从相机第一次向场景发射时使用，后续的散射与分层采样无关。教程中说分层采样在高频信息（物体边缘）上的精度比不进行分层采样更高。（虽然我觉得两张图结果差不多…）

如何生成随机变量

这里补充一下概率论中有关累积分布函数(Cumulative Distribution Function, CDF)和概率密度函数(Probability Distribution Function, PDF)的知识。

设$X$是一个随机变量，$x$是任意实数，函数 c d f ( x ) = P { X ≤ x } , − ∞ < x < + ∞ cdf(x) = P\left\{ X \leq x\right\}, -\infty < x < +\infty cdf(x)=P{X≤x},−∞<x<+∞称为$X$的累计分布函数。

如果对于随机变量$X$的累积分布函数$pdf(x)$，存在非负函数$pdf(x)$，使对于任意实数$x$，有：
$$cdf(x)={\int }_{-\infty }^{x}pdf(t)dt$$
则称$X$为连续型随机变量，其中函数$pdf(x)$称为$X$的概率密度函数，简称概率密度。

那么，如何生成一个随机变量呢？

设$X$是一个随机变量，它的概率密度函数为$pdf(x)$，它的累积密度函数可以表示为：
$$cdf(x)={\int }_{-\infty }^{x}pdf(t)dt$$
计算符合该概率分布的随机数方法如下所示：

• 对于概率密度函数$pdf(x)$，计算它的累积分布函数$cdf(x)$，如上面的等式所示；
• 计算$cdf(x)$的反函数$cd{f}^{-1}(x)$；
• 对于一个范围在$[0,1]$之间均匀分布的随机数$\xi$，则$X=cd{f}^{-1}(\xi )$，就是符合该概率分布的随机数。

举个例子，区间$[0,1]$之间的概率密度函数$pdf(x)=(n+1)x^n,x\in [0,1]$，计算其累积分布函数：
$$cdf(x)={\int }_0^x(n+1)t^{n}dt=x^{n+1}$$
其反函数为
$$cd{f}^{-1}(x)=\sqrt[n+1]{x}$$
则符合该概率分布的随机数为
$$X=\sqrt[n+1]{\xi }$$

重要性采样

正如前言中提到的，利用蒙特卡洛法进行积分的计算公式如下：
$$F(X)=\frac{1}{n}\sum _{k=1}^n\frac{f(X_k)}{pdf(X_k)}$$
当进行采样的分布律$pdf$与被采样函数越接近时，这个采样便会越快收敛。我们可以把采样点集中在被积函数数值较高的部分，以及来得到更准确更有效率的采样估算。这种非均匀的采样方式被称为重要性采样。

举个例子，我们对积分${\int }_0^2{x}^{2}dx$进行估计，采用均匀采样和重要性采样两种方式来进行蒙特卡洛积分。这个积分的值为$\frac{8}{3}$，换成小数为2.666666666666…

均匀采样

此时，我们令$pdf(x)=\frac{1}{2},x\in [0,2]$，计算其累积分布函数：
$$cdf(x)={\int }_0^{x}pdf(t)dt={\int }_0^x\frac{1}{2}dt=\frac{x}{2}$$
其反函数为
$$cd{f}^{-1}(x)=2x$$
因此符合该概率分布的随机数为
$$X=2\xi ,\xi \in [0,1]$$
这里的计算代码如下：

double f(double d) {
    return 2.0 * d;
}

double pdf(double x) {
    return 0.5;
}

int main() {
    int N = 1000000;
    auto sum = 0.0;
    for (int i = 0; i < N; i++) {
        auto x = f(random_double());
        sum += x*x / pdf(x);
    }
    std::cout << std::fixed << std::setprecision(12);
    std::cout << "I = " << sum / N << '\n';
}

运算结果：

I = 2.666942705288

重要性采样

我们让进行采样的分布律$pdf$与被采样函数接近，设$pdf(x)=c{x}^2,x\in [0,2]$，于是有
$$\begin{gathered} {\int }_0^{2}c{x}^{2}dx=1 \\ c{\int }_0^2{x}^{2}dx=c\ast \frac{8}{3}=1 \\ c=\frac{3}{8} \end{gathered}$$
此时我们得到一个与被采样函数接近的分布函数$pdf(x)=\frac{3}{8}{x}^2,x\in [0,2]$，计算其累积分布函数
$$cdf(x)={\int }_0^x\frac{3}{8}{t}^{2}dt=\frac{x^3}{8}$$
其反函数为
$$cd{f}^{-1}(x)=2\sqrt[3]{x}$$
因此符合该概率分布的随机数为
$$X=2\sqrt[3]{\xi },\xi \in [0,1]$$
这里的计算代码如下：

double f(double d) {
    return 2.0 * pow(d, 1.0/3.0);
}

double pdf(double x) {
    return (3.0/8.0) * x*x;
}

int main() {
    int N = 1;
    auto sum = 0.0;
    for (int i = 0; i < N; i++) {
        auto x = f(random_double());
        sum += x*x / pdf(x);
    }
    std::cout << std::fixed << std::setprecision(12);
    std::cout << "I = " << sum / N << '\n';
}

运算结果：

I = 2.666666666667

在结果中，我们可以明显感受到重要性采样的威力。（何止更快，只需要迭代1次就能得到正确结果）

生成半球面上的随机方向

这一小节是生成半球面上的随机方向。

首先是生成相对于Z轴的随机方向，为了简单起见，假设z是平面的法线。

首先我们考虑半球上的均匀随机方向生成，设此时随机方向的概率密度函数为$pdf(w)=c$，其中$w$是立体角，极坐标与立体角的对应关系有
$${\int }_{{\Omega }^2}p(w)dw={\int }_{{\Omega }^2}sin(\theta )d\varphi d\theta$$
于是有
$${\int }_{{\Omega }^2}p(w)dw=c{\int }_0^{\frac{\pi }{2}}sin(\theta ){\int }_0^{2\pi }d\varphi d\theta =2\pi c=1$$
则$c=1/2\pi$，那么
$$pdf(\theta ,\varphi )=sin\theta /2\pi$$
可以计算出
$$pdf(\theta )={\int }_0^{2\pi }p(\theta ,\varphi )d\varphi =sin\theta$$
再根据条件概率的公式，可得
$$pdf(\varphi |\theta )=\frac{pdf(\theta ,\varphi )}{pdf(\theta )}=\frac{1}{2\pi }$$
分别计算累积分布函数
$$\begin{gathered} cdf(\theta )={\int }_0^{\theta }sintdt=1-cos\theta \\ cdf(\varphi |\theta )={\int }_0^{\varphi }\frac{1}{2\pi }dt=\frac{\varphi }{2\pi } \end{gathered}$$
根据之前介绍的随机变量生成，设${\xi }_1,{\xi }_2$是$[0,1]$之间均匀分布的随机数，可以用$1-{\xi }_1$替换${\xi }_1$，我们能够求出：
$$\theta =co{s}^{-1}{\xi }_1,\varphi =2\pi {\xi }_2$$
代入极坐标，可得：
$$\begin{gathered} x=sin\theta cos\varphi =cos(2\pi {\xi }_2)\sqrt{1-{\xi }_1^2} \\ y=sin\theta sin\varphi =sin(2\pi {\xi }_2)\sqrt{1-{\xi }_1^2} \\ z=cos\theta ={\xi }_1 \end{gathered}$$
除此之外，还有一种用cos加权的随机生成方法，此时采用的是Lambertian光照模型。这时$pdf(w)$与$cos\theta$正相关，计算可得$pdf(\theta ,\varphi )=\frac{1}{\pi }cos\theta sin\theta$，用相似的方法，可以计算出
$$\theta =co{s}^{-1}{\sqrt{\xi }}_1,\varphi =2\pi {\xi }_2$$
代入极坐标，可得：
$$\begin{gathered} x=sin\theta cos\varphi =cos(2\pi {\xi }_2)\sqrt{1-{\xi }_1} \\ y=sin\theta sin\varphi =sin(2\pi {\xi }_2)\sqrt{1-{\xi }_1} \\ z=cos\theta =\sqrt{{\xi }_1} \end{gathered}$$

正交基

这一部分主要是让上一章随机生成的方法支持生成任意表面的法向量。

相对坐标系由一个原点$O$和三个基方向$u、v、w$来定义的。

给定一个表面法线向量$n$，如何构建这个相对坐标系？

我们不需要原点，因此只需要求两外两个互相垂直的向量。

1. 选择一个辅助向量$a$，这个向量不能与$n$垂直，通常可以用y轴正方向，如果$n$接近与y轴平行，则用将$a$设置为x轴正方向，即：

if (fabs(n.x()) > 0.9)
    a = vec3(0, 1, 0);
else
    a = vec3(1, 0, 0);

2. 得出第二个轴向量$s$，$s$的方向可由$n$与$a$叉乘得出。最后一个轴向量$t$，可由$n$与$s$叉乘得出，即：

vec3 s = unit_vector(cross(n, a));
vec3 t = cross(n, s);

一旦我们得到了正交基$s$、$t$和$n$，我们在获取一个相对z轴随机的向量$(x,y,z)$之后，便可以获得一个相对于$n$轴的随机向量：
$$Randomvector=x\text{s}+y\text{t}+z\text{n}$$
可以单独构造一个表示正交基的类，这个类的核心代码如下：

#ifndef ONB_H
#define ONB_H

#include "rtweekend.h"

class onb {
  public:
    onb() {}

    vec3 operator[](int i) const { return axis[i]; }
    vec3& operator[](int i) { return axis[i]; }

    vec3 u() const { return axis[0]; }
    vec3 v() const { return axis[1]; }
    vec3 w() const { return axis[2]; }

    vec3 local(double a, double b, double c) const {
        return a*u() + b*v() + c*w();
    }

    vec3 local(const vec3& a) const {
        return a.x()*u() + a.y()*v() + a.z()*w();
    }

    void build_from_w(const vec3& w) {
        vec3 unit_w = unit_vector(w);
        vec3 a = (fabs(unit_w.x()) > 0.9) ? vec3(0,1,0) : vec3(1,0,0);
        vec3 v = unit_vector(cross(unit_w, a));
        vec3 u = cross(unit_w, v);
        axis[0] = u;
        axis[1] = v;
        axis[2] = unit_w;
    }

  public:
    vec3 axis[3];
};


#endif

生成指向光源的随机方向

我们可以让所有的散射光线都努力指向光源的方向，在指向光源方向途中，如果折射的光线射入物体表面之下，或者方向和光源平面平行，则丢弃。

设光源表面的面积为$A$，我们假设光线打到光源上的概率是均匀的，则打在平面上的每一点的概率密度为$1/A$。

对于光线打向的这一小块区域$dA$，采样这块区域的概率为$p_q(q)\cdot dA$，其中$p_q(q)=1/A$。而在光线的发射球面上，采样球表面小片区域$dw$的概率为$p(w)\cdot dw$，$dw$和$dA$有如下的几何关系：
$$dw=\frac{dA\cdot cos(\theta )}{distanc{e}^2(p,q)}$$
由于在采样中$dw$和$dA$的概率必须相等，因此有
$$\begin{gathered} p(w)\cdot dw=p_q(q)\cdot dA \\ p(w)\cdot \frac{dA\cdot cos(\theta )}{distanc{e}^2(p,q)}=p_q(q)\cdot dA \end{gathered}$$
其中$p_q(q)$的值为$\frac{1}{A}$，因此我们有
$$p(w)=\frac{distanc{e}^2(p,q)}{cos(\theta )\cdot A}$$
如何生成这样的射线

随机在光源平面上取一个点，与当前表面的点构成射线，如果射线在表面点和光源之间没有别的遮挡，那么就可以成功发出。这条射线的$pdf$已由上式推导计算出。

直接生成指向光源的射线，其渲染结果如下：

这样做的效果是，牺牲了折射的次数，好处是噪声比较少。

构建射向球面的随机方向

这一部分构建从一个点到球面的射线的均匀随机采样。

参考生成半球面上的随机方向，我们要进行均匀采样，设$pdf(w)=C$因此有
$${\int }_{{\Omega }^2}pdf(w)dw=C{\int }_0^{{\theta }_{max}}sin(\theta ){\int }_0^{2\pi }d\varphi d\theta =2\pi C(1-cos{\theta }_{max})=1$$
其中$w$是立体角，$\theta ,\varphi$为极坐标，${\theta }_{max}$为光线能够打到球面上时，$\theta$的最大值，如下图所示：

接上面的公式，因此有
$$\begin{gathered} pdf(w)=C=\frac{1}{2\pi \cdot (1-cos{\theta }_{max})} \\ pdf(\theta ,\varphi )=C\cdot sin\theta \end{gathered}$$
那么
$$\begin{gathered} pdf(\theta )={\int }_0^{2\pi }pdf(\theta ,\varphi )d\varphi =2\pi C\cdot sin\theta \\ pdf(\varphi |\theta )=\frac{pdf(\theta ,\varphi )}{pdf(\theta )}=\frac{1}{2\pi } \end{gathered}$$
分别计算累积分布函数
$$\begin{gathered} cdf(\theta )={\int }_0^{\theta }2\pi C\cdot sintdt=2\pi C(1-cos\theta ) \\ cdf(\varphi |\theta )={\int }_0^{\varphi }\frac{1}{2\pi }dt=\frac{\varphi }{2\pi } \end{gathered}$$
设${\xi }_1,{\xi }_2$是$[0,1]$之间均匀分布的随机数，我们有
$$\begin{gathered} \theta =co{s}^{-1}(1-\frac{{\xi }_1}{2\pi C})=co{s}^{-1}(1+{\xi }_1\cdot (cos{\theta }_{max}-1)) \\ \varphi =2\pi {\xi }_2 \end{gathered}$$
因此我们有下面的式子：
$$\begin{gathered} z=cos\theta =1+{\xi }_1\cdot (cos{\theta }_{max}-1) \\ x=cos\varphi \cdot sin\theta =cos(2\pi \cdot {\xi }_2)\cdot \sqrt{1-z^2} \\ y=sin\varphi \cdot sin\theta =sin(2\pi \cdot {\xi }_2)\cdot \sqrt{1-z^2} \end{gathered}$$
对玻璃球采用重要性采样，渲染结果如下：

密度混合

我们可以创建任何PDF的线性混合，以生成随机变量。PDF的任何加权平均值也是PDF。只要权重为正并且加起来为1，我们就得到了一个新的PDF分布。

例如，在生成随机散射光线时，我们提到了两种方式，一种是让散射光线随机射向光源，另一种是让散射光线根据表面法线随机散射（Lambertian光照模型）。我们也可以将两种方式结合，让光线既有可能随机射向光源，又可能按照Lambertian光照模型进行散射。这种思路的伪代码如下：

if (random_double() < 0.5)
    pick direction according to pSurface
else
    pick direction according to pLight

实现这种密度混合的代码如下：

class mixture_pdf : public pdf {
  public:
    mixture_pdf(shared_ptr<pdf> p0, shared_ptr<pdf> p1) {
        p[0] = p0;
        p[1] = p1;
    }

    double value(const vec3& direction) const override {
        return 0.5 * p[0]->value(direction) + 0.5 *p[1]->value(direction);
    }

    vec3 generate() const override {
        if (random_double() < 0.5)
            return p[0]->generate();
        else
            return p[1]->generate();
    }

  private:
    shared_ptr<pdf> p[2];
};

最后，展示一下最终的渲染结果：

ps：文中没有讲到的有利于采样的内容：低差异序列。

完整代码

链接：https://pan.baidu.com/s/1n3EcaPQoZcG6dyRyNWuhUA?pwd=itd6
提取码：itd6
–来自百度网盘超级会员V6的分享

参考

https://raytracing.github.io/books/RayTracingTheRestOfYourLife.html

https://zhuanlan.zhihu.com/p/146144853

https://dreamerchen.com/post/ImportanceSampling.html