【RAY TRACING THE REST OF YOUR LIFE 超详解】 光线追踪 3-2 蒙特卡罗(二) 重要性采样

 

书本内容：见相册

注：鉴于不少网站随意爬取数据，可能致使内容残缺以及引用失效等问题，影响阅读，请认准原创网址：

https://www.cnblogs.com/lv-anchoret/category/1368696.html

 preface

还记的咱们上一篇说的Monte Carlo 维度诅咒吗

上一篇算是二维的例子吧，你们看了以后是否想着写一个一维的Monte Carlo模拟积分？（我想了，没写出来）

为何要整这个嘞

光照渲染中多处涉及重积分，最终结果是要求取一个近似值，于是须要对其值进行数值估计，Monte Carlo方法就是一个较为理想的方案。

其实咱们的光线追踪器不只用了不少向量运算，还用了不少数值分析的知识，好比以前的背景用的是最简单的插值，柏林噪声纹理用的是三线性插值等

 

 ready

你须要对如下知识有了解：

二重定积分

几率密度函数

反函数

 

 Chapter 2：One Dimensional MC Integration

积分是用来求取面积或者体积的东东，常见形式为：I=∫x² dx

假设积分区间为0->2,咱们用计算机符号表示为I = area(x²,0,2)

那么咱们来模拟如下这个积分

咱们先来看一个普通的：

 

那么，咱们就能够写代码了（程序2-1）

void MC_integration(double(*f)(double x), double low, double high)
{
    size_t all = 1000000;
    double sum{ 0 };
    for (int i = 0; i < all; ++i)
    {
        double x = (high - low) * rand01() + low;
        sum += x*x;
    }
    stds cout << "I = " << (high - low) * sum / all << stds endl;
}

int main()
{
    MC_integration([](double x) {return x*x; }, 0, 2);
}

正如咱们所指望的，但上述方法也能够计算那些咱们用解析法没法解决的积分，如被积函数为log（sin（x）。

在图形学中，咱们常常有一些咱们能够评估但不能明确写下来的函数，或者咱们只能几率评估的函数。 事实上，前面两本书的光线追踪代码中的lerp函数，其实咱们不知道在每一个方向上看到的是什么颜色，但咱们能够在任何给定的维度上统计估算它。

 

我 又想起来上本书刚刚写过的区域光照渲染图片，里面的黑色噪点真的把好精美的Cornell box 盒子给糟蹋了，那时由于咱们对全部的东西作了统一的采样，而没有对光源进行足够的光源采样，因此咱们能够对光源进行更多的随机采样，可是咱们又须要控制这个比重，以防止过采样的发生，因此，如何控制这个比例，已达到更好的效果，这就须要引入一个很是重要的概念——重要性采样

可是在讲重要性采样以前，咱们须要了解一些关于几率密度函数的东东

　　

 

先看书上这张图 

若是咱们添加更多的tree，那么这张图的长条将会增高（变长），由于他们的数量更多了，若是Height轴细分为更多的小块，那么，每一个长条的高度会下降（变短）。

离散密度函数与直方图的不一样之处在于它将频率y轴归一化为分数或百分比。

而连续直方图中，咱们将Height轴细分为无数个小块，那么，纵轴将不会是一个分数，由于全部小格对应的长条基本没有高度（每一个小格表明的树高区间基本只有一棵树或者没有树，出现频率几乎为0）

因此对于密度函数，咱们调整小格对应的影响因素，使得咱们添加更多的小格时，长条不会过低

对于上述条形图，咱们采起以下：

小格的高 = 树高介于H~H'的比例 / (H - H')

它将是颇有用的！ 咱们能够将其解释为树高的统计预测器：

即：一个介于H~H'的随机树高的几率 = 小格的高 * (H - H')

若是咱们须要知道多个小格对应的几率，那么加和便可

其实上述公式就构成了几率密度函数，即，面积表明着对应横坐标区间的几率，也就是上面的第二个公式

而几率密度函数就是一种连续的分数直方图，简称为pdf(probability density function)

 

让咱们来整一个pdf，或许会有更深刻的理解。若是我想要一个0~2的随机数r它的出现几率和r本身成比例，咱们将指望pdf为以下图像：

 

 

r位于（x0,x1）的几率为C*area(p(r),x0,x1),其中C为常量,为了简便，C通常取1

而area(p(r),0,2) = 1（几率密度函数的总面积就是总几率）

由于p(r)和r成比例，因此，设另外一个常量C',则p(r) = C'r

故咱们解积分：

1 = ∫_0->2 C'r dr = C'r²/2|_0->2 = 2C' =>C' = 1/2  公式2-1

因此 p(r) = r/2

咱们如何用pdf p(r)生成一个随机数？为此咱们须要更多的机器。不要担忧这不会永远持续下去！（不会像上一篇最后那个程序，永无止境。。）

 

给一个随机数d = rand01(),咱们知道咱们写的rand01()产生的随机数对于0~1的每一个浮点数的几率都是均匀一致的，咱们须要找个函数f(d)获得咱们想要的。假定e = f(d) = d²,此次将再也不是一个均匀一致的pdf了，之前是的，以后咱们讲（把此处记为pos1遗留问题）

咱们为了观察这个函数还须要累积几率分布函数Cpdf(Cumulative probability distribution function)

记为P(x),则P(x) = area(p,-inf,x),几何意义就是从负无穷到x的累积几率和，也就是x左边部分的积分造成的面积

并且咱们规定在没有定义的定义域内p(x) = 0,若采用pdf函数p(r) = r/2，那么，P(x)以下：

P(x) = 0,　　x < 0

P(x) = x²/4, 0 < x < 2

P(x) = 1,　　x > 2

当x∈（0,2）,P(x) = ∫_0->x p(r) dr = ∫_0->x r/2 dr = x²/4

那么，x和r是什么关系？其实就是两个虚拟变量，相似于程序中的参数，若是咱们把x调整为有效区间的中点，那么

P(1.0) = 1/4

这就是说基于咱们的pdf几率密度函数p(r) = r/2的设定下产生一个小于1的随机数的几率为25%

这产生了一种巧妙的观察结果，这种观察结果是产生非均匀随机数的许多方法的基础。

咱们想要一个函数f（），它知足：

当咱们调用f（drand48（））时，咱们获得一个基于设定的Cpdf（x²/4）的一个返回值。 咱们不知道那是什么，但咱们知道当x为25％时应该小于1，而参数为75％应该高于1。 若是f（）为增函数，那么咱们指望f（0.25）= 1.0。 

上面说了一堆，意思就是x和r无非就是函数参数，即P(1.0) = 0.25的时候，咱们一样想获得一个函数f()，它知足f(0.25)=1.0

即：根据Cpdf得知随机数小于1.0的几率为0.25，而存在一个f()，使得以0.25为参数的时候可以获得那个1.0（即累积几率分布函数右端的边界线）

这能够推广到每一个可能的输入

即：f(P(x)) = x

也就是说f(x) = P^-1(x)

也就是P(x)的反函数。对于咱们的目的而言，上述的意思就是：若是咱们有了pdf函数p()和对应的Cpdf函数P(),咱们为f()输入一个随机数，那么它就会给出咱们想要的：

e = P^-1(rand01())

对于咱们的pdf函数p(x) = x/2，以及对应的P(x),咱们须要计算P^-1

即：若是咱们有 y = x²/4

那么，x = sqrt(4y)

所以基于pdf,咱们输入随机数，获得的就是

e = sqrt(4*rand01())

那么e∈（0,2），正如咱们所指望的，咱们输入一个0.25，它就返回一个1.0

 

到了这里，其实上面重复了不少废话，上面一大堆，表面上干了一个什么事情呢，用一个0~1的随机数映射0~2这个积分区间

咱们想一下咱们要作的是什么，模拟x的二次方曲线的积分，咱们要用随机数模拟随机选取积分区间的值，而咱们的积分区间就是0~2，

可是不必这么麻烦啊，0~1随机数随机模拟0~2取值，直接2 * rand01()不就完啦，也就是程序2-1

 

不不不，它们有着天差地别，咱们上述一顿操做是为了引入一个概念——Important Sample

咱们经过设置pdf使得，对于积分区间的自变量取值有了权重比例，也就是说有些地方的采样要多一点，有些地方要少一点。

也就是：I = 1/n * Σ(F(x_i)/pdf(x_i))

由此可知，每一个积分采样 I_i = 采样高度（函数值）/该采样点占的权重，全部的积分采样取均值，实现了不一样采样点为整个积分作的贡献是不一样的

constexpr double pdf(const double x) {return 0.5 * x; }

void pdf1()
{
    size_t all{ 100000 };
    double sum{ 0. };
    for (int i = 0; i < all; ++i)
    {
        double x = sqrt(4 * rand01());
        sum += x*x / pdf(x);
    }
    stds cout << "I = " << sum / all << stds endl;
}

 

咱们来解释pos1处的遗留问题

以前的程序2-1所指的积分方案中，是均匀一致采样，即对于任意采样点同等对待

即,设p(x) = C

则公式2-1 写为 

1 = ∫_0->2 C dr = Cr|_0->2 = 2C =>C = 1/2

即，开篇的方法中的pdf函数为p(x) = 1/2

因此程序2-1能够改写为：

constexpr double pdf(const double x) {return 0.5; }

void pdf1()
{
    size_t all{ 100000 };
    double sum{ 0. };
    for (int i = 0; i < all; ++i)
    {
        double x = 2 * rand01();
        sum += x*x / pdf(x);
    }
    stds cout << "I = " << sum / all << stds endl;
}

 

在important sample中咱们以为函数峰高的部分对于整个积分的贡献更多，低峰对于积分贡献不大，因此咱们对于x²积分函数采起的pdf函数为x/2

其实，pdf应该采起的就是积分函数自己，在知道答案的状况下，咱们采起以下pdf函数最佳：

p(x) = (3/8)x²

则P(x) = x³/8

P^-1(x) = pow(8x,1/3)

则函数为：

constexpr double pdf(const double x) { return 3 * x * x / 8; }

void pdf1()
{
    size_t all{ 1 };
    double sum{ 0. };
    for (int i = 0; i < all; ++i)
    {
        double x = pow(8 * rand01(), 1. / 3);
        sum += x*x / pdf(x);
    }
    stds cout << "I = " << sum / all << stds endl;
}

 

至此，全部的内容都已经阐述完毕

咱们回顾一下，由于这个是光线追踪蒙特卡罗中的大部分概念

1. 你有一个被积函数f(x)的积分域【a,b】

2. 在域【a,b】中选择一个非零的pdf函数p()

3. 对全部的积分采样I_i = f(r)/p(r)取平均，pdf函数p()，随机数r

这老是会收敛到正确答案，p越接近f，收敛越快

 

感谢您的阅读，生活愉快~