《用余生学习光线追踪》2.光照散射与重要性采样

本项目参考自《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.生成球坐标随机方向
《用余生学习光线追踪》4.直接对光照采样

球面各方向上的蒙特卡洛积分

Ray Tracer需要选取随机方向，方向可以定义为单位圆上的一个点。我们需要定义一个二维的PDF，假设我们要对任意方向求解以下积分：
$$\int co{s}^2(\theta )$$

如果要用蒙特卡洛算法来求解积分，则我们要计算$co{s}^2(\theta )/p(direction)$，可以用基于极坐标的参数$(\theta ,\varphi )$来定义方向。

创建一个单位球上的点的函数如下：

vec3 random_on_unit_sphere() {
    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\pi )$，如果$\theta$是球心到球面上一点连线与z轴正方向的夹角，则$cos(\theta )$为该单位向量的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_on_unit_sphere();
            float cosine_squared = d.z()*d.z();
            sum += cosine_squared / pdf(d);
    }
    std::cout << "I =" << sum/N << "\n";
}

上述积分的解析解为$\frac{4}{3}\pi$

备注：solid angle 和 area 和 direction 是一回事。

光照散射

接下来，我们会大幅度修改光照的算法。之前的光照都是直接光线和物体表面散射得到的，并没有用到概率论，接下来把概率论用上。首先，我们去假设光线是否被吸收了：

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

• 光被吸收的概率为$1-A$

A代表albedo，也就是我们所说的反射率。
我们将会使用波长来代替RGB，RGB对应长、中、短波长。如果光线确实发生散射了，那就会有一个方向性的概率分布，我们将其定义成关于立体角(solid angle)的散射pdf（scattering pdf），数学表示为$s(direction)$，散射pdf会随着方向而变化。根据散射概率、散射pdf和颜色函数，颜色函数定义如下：
$$Color=\int A\cdot s(direction)\cdot color(direction)$$

注意，A、s()可能会随着观察的方向和物体表面的位置不同而变化。
套用蒙特卡洛积分公式，得到估算公式如下：
$$Color=(A\cdot s(direction)\cdot color(direction))/p(direction)\text{\hspace{1em}(0)}$$
其中，$p(direction)$为随机生成的方向的pdf。
对于兰伯特材质来说，我们已经隐性地套用过这个公式了：我们假设其$p()$是一个cosine函数。s()跟$cos(\theta )$成正比，其中$\theta$ 为相对于法线表面的夹角。注意，所有的pdf在积分区间内的积分结果为1。当 $\cos (\theta )<0$的时候，$s(direction)=0$，半球的cos的积分值为 $\pi$。
在球坐标系中，有公式：
$$\delta A=sin(\theta )\delta \theta \delta \varphi$$
于是有：
$$Area={\int }_0^{2\pi }{\int }_0^{\pi /2}cos(\theta )sin(\theta )\delta \theta \delta \varphi =2\pi \cdot \frac{1}{2}=\pi$$
因此，对于一个兰伯特表面，散射pdf为：
$$s(direction)=cos(\theta )/\pi \text{\hspace{1em}(1)}$$
如果采样pdf也是用这个函数，即
$$p(direction)=cos(\theta )/\pi \text{\hspace{1em}(2)}$$
则s(direction)和p(direction)相互抵消后有：
$$Color=A\ast color(direction)$$
这就是原来的color函数。
现在有了概率论，我们可以朝着更加重要的方向（例如光源的方向）额外发射光线。
有些教材里面会通过双向反射分布函数（bidirectional reflectance distribution function，简称BRDF）描述反射：
$$BRDF=A\ast s(direction)/cos(\theta )$$
以兰伯特材质为例，BRDF公式为：
$$BRDF=A/\pi$$

光照散射与重要性采样

上一节的蒙特卡洛算法求解积分的步骤如下：

1. 找到在区间$[a,b]$内被积函数$f(x)$
2. 在区间$[a,b]$内挑选一个合适的非零的概率密度函数$p$
3. 生成一堆满足上述概率密度函数$p$分布的随机数$r$，从而得到一堆$f(r)/p(r)$，求这堆$f(r)/p(r)$的平均值

接下来两节，目标都是向着光源的位置，发送一堆额外的光线，从而降噪。

假设我们可以用一个概率密度函数$pLight(direction)$生成一堆额外向着光源的光线。

假设我们有一个跟s有关的概率密度函数$pSurface(direction)$
两个pdf的线性求和，仍然是pdf，例如：
$$p(direction)=0.5\cdot pLight(direction)+0.5\ast \cdot pSurface(direction)$$
只要权重是整数且权重之和为1，任意的pdf求和结果仍然是pdf。
我们希望能够找到一个pdf，使得$s(direction)\ast color(direction)$增大的同时，这个pdf也同样增大。对于漫反射表面来说，这也是在猜测哪里会有更大的$color(direction)$值。

为了降噪，我们需要会构造一个能发射更多光线到光源的pdf。
回忆一下蒙特卡洛积分公式：
$$\int f(x)\approx f(r)/p(r)$$
对于Lambertian材质，我们定义：
$$p(direction)=cos(\theta )/\pi$$

接下来，修改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);
}

其中，color的计算来源于公式一和公式二。
渲染结果和之前一样，没有区别。

接下来，根据经验，换一种采样策略。选取命中表面上的一个半球的某个方向作为散射方向，则有：
$$p(direction)=1/2\pi$$

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;
}

结果如下：

左box前面的颜色似乎更加均匀了。这究竟是bug还是feature？

参考资料：
1.Ray Tracing: The Rest of Your Life