【光线追踪系列十五】MC光照散射与基于重要性采样的材质

本文主要参照 Ray Tracing: The Rest of Your Life，其中只是主要精炼光追相关理论，具体实现可参照原文。

本章内容之前推荐最好先预备以前文章【PBR系列二】辐射度量学理论相关知识。

我们之前的光线追踪器中有一个过程是获取一个三维空间的随机方向，我们利用单位球体内的随机点来代表该方向，那么我们今天用Monte Carlo 来实现相同的效果

一、重要性采样在光追中作用

主要作用就是加速函数收敛，可以使用较小的样本数获取最优的趋近值。

结合光追来看就是：小光源会产生太多的噪声，因为我们的均匀采样不够频繁。如果我们发送更多的随机样本到光照下，我们可以减少这个问题，但是我们需要降低这些样本的权重，以适应过采样。我们如何调整？为此，我们需要概率密度函数的概念。

通俗的解释来说就是：光源很小，我们发射的采样光线大多数根本打不到光源上，比如挨着的两个像素，而且这两个像素显示的内容并不是某物体的边界处，两个像素各自发射100个采样光线，在随机散射下，其中一个像素有80个采样光线最终打中了光源，比较亮，但是另一个像素很不巧，只打中了20个，就很暗，但按理来说它们挨着，又不是边界，明暗应该很接近才对，所以产生很强烈的噪声（当然如果你每个像素发射10000个光线，随机噪声就会大大降低，但是也更浪费时间，是只发射100个光线的渲染速度的100倍！）。所以我们应该让像素发射的光线尽可能在能打到光的方向进行散射，换成这里的积分的例子，即我们应该要让更多的样本产生在x比较大的地方，因为在积分中，大的x对积分的准确性影响更大，但是我们又不能将所有x都取大的，不然得到的积分值肯定大于正常值，所以需要有一个权重来衡量：即我们要产生更多对最后结果有用的样本，同时还要让这些样本通过某种权重防止其过采样。

再举个例子详细说明一下，比如，我要计算x^2从0到2的积分，我们真实地计算一下发现，从0到1积分起来的值只有1/3，而从1到2积分的值是7/3，所以说，从1-2的积分的值不就占比重更大吗？当我们要通过估计的方法来获得值的时候，假如我们随机产生200个数，其中100个在（0-1）之间，100个在（1-2）之间，假设通过100个数的估计偏差在2%左右，则从0-1积分的估计值与真实值的偏差是0.0067（即(1/3)*0.02），而从1-2积分值的估计值与真实值的偏差为 0.047（即(7/3)*0.02），两个偏差加起来的值是0.053左右。而如果我们进行一下分配，让随机数在1-2之间产生的更多，这样就会降低1-2之间积分的偏差。根据刚才的计算，我们知道1-2之间占积分值比重很大，所以这样就能很好的降低总的偏差。

二、球面蒙特卡洛积分

首先我们需要定义一个2D的pdf函数，我们假定下面这个积分适用于任意方向：

如果要用蒙特卡洛算法来求解积分，则我们要计算：

这里的direction在我们的环境中的定义是：基于极坐标的参数 (θ,ϕ) 来定义方向。

不管你怎么整，一定要记住：**设计一个pdf函数，它的积分必须是1，并且pdf代表着该方向被采样的相对概率。

首先我们来看一下之前利用单位球体内的随机点来代表该方向的代码：

vec3 random_in_unit_disk() {
    vec3 p;
    do {
        p = 2.0*vec3(random_double(),random_double(),random_double()) - vec3(1,1,1);
    } while (dot(p,p) >= 1.0);
    return p;
}

如果我们要得到球体表面的三维空间点，那么我们只需要单位化一下即可，也就是return p的单位化向量。

vec3 random_in_unit_disk() {
    vec3 p;
    do {
        p = 2.0*vec3(random_double(),random_double(),random_double()) - vec3(1,1,1);
    } while (dot(p,p) >= 1.0);
    return unit_vector(p);
}

这些单位球上的随机点的PDF，就是1/area，其中area为单位球的面积，也就是1/(4π)，如果θ是球心到球面上一点连线与z轴正方向的夹角，则cos(θ) 为该单位向量的z轴上的分量。

写成代码如下：

inline float pdf(const vec3& p) {
    return 1 / (4*M_PI);
}

int main() {
    int N = 1000000;
    float sum;
    for (int i = 0; i < N; i++) {
            vec3 d = random_in_unit_disk();
            float cosine_squared = d.z()*d.z();
            sum += cosine_squared / pdf(d);
    }
    std::cout << "I =" << sum/N << "\n";
}

上述积分的解析解为 (4/3)π。

三、散射光照

在之前光追中，我们基于表面或者次表面实现散射光线，这是一个普适模型。但是，更自然的方式是概率论。

首先，我们去假设光线是否被吸收了：

• 光被散射（没被吸收）的概率为A
• 光被吸收的概率为1−A

这里的A其实代表之前材质中的albedo(latin for whiteness)，Albedo代表着某些反射形式的反射比例。

当我们实现玻璃材质的时候，albedo伴随着入射光线方向的变化而变化，进而计算相应变化的像素值。在大多数的基于物理的渲染器，我们将用一组波长表示浅色而不是RGB，我们用长中短波长来代表RGB。如果实现光散射，那么我们可以基于立体角设计一个pdf来描述光线散射方向分布。我将其称为散射pdf：s(direction)，这个散射pdf也会随着入射方向的变化而变化。

根据散射概率、散射pdf和颜色函数，颜色函数定义如下：

需要注意的是：A 和 s()都是依赖于观察方向的，当然，color也随着观察方向的变化而变化,A 和 s()也可能随着表面或体内部的位置而变化。套用蒙特卡洛积分公式，得到估算公式如下：

其中，p(direction) 是方向参数随机模拟生成的一个pdf函数。

对于Lambertian材质表面，我们假定其密度是呈余弦规律的。所以, 一个Lambertian表面的p()是正比于cosθ，而θ是光线方向和表面法线的夹角，还需要注意一点：所有设计的pdf函数在定义域内的积分必须是1。我们还知道，对于磨砂表面来说，入射光线和表面法线的夹角不会超过90°，所以，我们规定：

• 对于 cosθ < 0 的情况，我们设定s(direction) = 0

所以有效区间为半球面，而基于半球的cos积分为π。

具体推导如下：

1. 在球坐标系中，有公式：

1. 对其求半球积分，于是有：
   
   所以，Lambertian表面散射pdf为：
   
   如果我们采样也使用同样的pdf函数：
   
   那么分子分母可以消去，就得到了：
   
   也有相关资料中会用BRDF来表示反射，你可能会看到双向反射分布函数（BRDF）描述的反射。 它的表示形式非常简单：
   
   以Lambertian材质为例，BRDF公式为：
   

四、基于重要性采样的材质

我们的目标就是向光源发送额外的光线使得我们的图片噪声更少。

1. 假设我们可以用一个概率密度函数pLight(direction)生成一堆额外向着光源的光线。
2. 假设我们有一个跟s有关的概率密度函数pSurface(direction)。

我们可以把两者进行线性组合，得到一个综合的pdf，最简单的形式如下：

其实，只要我们的权重是正的，并且加起来等于1,也就是我们经常提到的：pdf积分一定要能积到1，任意pdf函数组合形成的pdf只要满足这两点都是合理的，所以我们可以将几个分量因素的pdf组合控制，使我们的pdf函数更加贴近真实效果。

因此，寻找最佳pdf便成了光追渲染中的主要问题，要弄清楚如何使s(direction)*color(direction)增大的同时，这个pdf也同样增大，也就是找到合适的pdf，创造更真实的颜色比例效果。

对于漫反射（diffuse or Lambertian）表面，我们猜测它更偏重color（direction）因子。 对于镜面（metal）材质表面来说，s()因子只在一个方向附近很大，所以它更重要。

事实上，大多数渲染器将镜面设计成一种特殊情况，只是隐含地使用s / p 进行计算，所以，我们的代码目前也是这样。

为了降噪，我们需要会构造一个能发射更多光线到光源的pdf。
回忆一下蒙特卡洛积分公式：

对于Lambertian材质，我们定义：

接下来，修改material类的scatter函数，新增函数scattering_pdf():

class material  {
    public:
        virtual bool scatter(const ray& r_in,
            const hit_record& rec, vec3& albedo, ray& scattered, float& pdf) const {
            return false;
        }
        virtual float scattering_pdf(const ray& r_in, const hit_record& rec,
                                     const ray& scattered) const {
            return 0;
        }
        virtual vec3 emitted(float u, float v, const vec3& p) const {
            return vec3(0,0,0);
        }
};

修改lambertian类，scatter函数新增了计算pdf的功能。

class lambertian : public material {
    public:
        lambertian(texture *a) : albedo(a) {}
        float scattering_pdf(const ray& r_in,
            const hit_record& rec, const ray& scattered) const {
            float cosine = dot(rec.normal, unit_vector(scattered.direction()));
            if (cosine < 0)
                return 0;
            return cosine / M_PI;
        }
        bool scatter(const ray& r_in,
            const hit_record& rec, vec3& alb, ray& scattered, float& pdf) const {
            vec3 target = rec.p + rec.normal + random_in_unit_sphere();
            scattered = ray(rec.p, unit_vector(target-rec.p), r_in.time());
            alb = albedo->value(rec.u, rec.v, rec.p);
            pdf = dot(rec.normal, scattered.direction()) / M_PI;
            return true;
        }
        texture *albedo;
};

color函数修改如下：

vec3 color(const ray& r, hittable *world, int depth) {
    hit_record rec;
    if (world->hit(r, 0.001, MAXFLOAT, rec)) {
        ray scattered;
        vec3 attenuation;
        vec3 emitted = rec.mat_ptr->emitted(rec.u, rec.v, rec.p);
        float pdf;
        vec3 albedo;
        if (depth < 50 && rec.mat_ptr->scatter(r, rec, albedo, scattered, pdf)) {
            return emitted + albedo*rec.mat_ptr->scattering_pdf(r, rec, scattered)
                *color(scattered, world, depth+1) / pdf;
        }
        else
            return emitted;
    }
    else
        return vec3(0,0,0);
}

渲染结果和之前一样，没有区别（每个像素采样200）。

接下来，根据经验，换一种采样策略。选取命中表面上的一个半球的某个方向作为散射方向，则有：

scatter函数换成：

bool scatter(const ray& r_in,
    const hit_record& rec, vec3& alb, ray& scattered, float& pdf) const {
    vec3 direction;
    do {
        direction = random_in_unit_sphere();
    } while (dot(direction, rec.normal) < 0);
    scattered = ray(rec.p, unit_vector(direction), r_in.time());
    alb = albedo->value(rec.u, rec.v, rec.p);
    pdf = 0.5 / M_PI;
    return true;
}

渲染结果如下：

五、总结

我们的核心思想、渲染路径是这样计算的：

• 首先我们从eye（或者camera）发射采样光线，它指向投影屏幕的一个像素位置，也就是当前采样光线正在计算的那个像素位置。
• 然后光线以射线的形式计算eye与投影屏幕之间所有的几何体交点，得到与eye最近的那个，然后我们以那个交点为落脚点，作为新的eye，按照光学表面散射原理计算新的射线的方向（结合具体的材质），然后继续测交。
• 如果在递归深度上限范围内没有找到光源，那么说明该像素位置不会有光传入眼睛，也就是当前像素位置是黑色的；如果它找到光源了，也就是经过多次散射最后指向了光源，那么就说明光源发出的光线可以沿着我们的路径逆向进入眼睛，那么我们就看到了这个像素，那么如何计算像素值呢？
• 这就涉及到了路径渲染方程，里面有一个rgb分量衰减比例控制参数，它依据材质和纹理计算得出，这是之前的，我们现在要加上pdf函数，帮助我们更逼真地计算像素值。