【RAY TRACING THE REST OF YOUR LIFE 超详解】 光线追踪 3-4 基于重要性采样的材质初探

 

 Preface

咱们今天来把第三本书从开局到如今讲的一大堆理论运用到咱们的框架中，那么今天咱们首先将原始的材质改成基于重要性采样原理的材质

这一篇是代码工程中进行MC理论应用的初步尝试篇

 

 Ready

咱们须要这几个重要的文件，我担忧你们手上的文件可能不太对，因此再贴一下

/// rectangle.hpp

// -----------------------------------------------------
// [author]        lv
// [begin ]        2019.1
// [brief ]        the rectangle-class for the ray-tracing project
//                from the 《ray tracing the next week》
// -----------------------------------------------------

#pragma once


namespace rt
{

//the statement of xy_rect class

class xy_rect : public intersect
    {
public:
    xy_rect() {  }

    /*
    @brief: the rectangle in the x-y plane
    @param: the boundary of x axis
            the boundary of y axis
            the value of z axis
            the material of rectangle
    */
    xy_rect(rtvar x1, rtvar x2, rtvar y1, rtvar y2, rtvar z0, material* mat);

    virtual bool hit(const ray& sight, rtvar t_min, rtvar t_max, hitInfo& info)const override;

    virtual aabb getbox()const override;

private:
    material * _materialp;

    rtvar _x1, _x2;

    rtvar _y1, _y2;

    rtvar _other;
    };

//the statement of xz_rect class

class xz_rect : public intersect
    {
public:
    xz_rect() {  }

    /*
    @brief: the rectangle in the x-z plane
    @param: the boundary of x axis
            the boundary of z axis
            the value of y axis
            the material of rectangle
    */
    xz_rect(rtvar x1, rtvar x2, rtvar z1, rtvar z2, rtvar y0, material* mat);

    virtual bool hit(const ray& sight, rtvar t_min, rtvar t_max, hitInfo& info)const override;

    virtual aabb getbox()const override;

private:
    material * _materialp;

    rtvar _x1, _x2;

    rtvar _z1, _z2;

    rtvar _other;
    };

//the statement of yz_rect class

class yz_rect : public intersect
    {
public:
    yz_rect() {  }

    /*
    @brief: the rectangle in the y-z plane
    @param: the boundary of y axis
            the boundary of z axis
            the value of x axis
            the material of rectangle
    */
    yz_rect(rtvar y1, rtvar y2, rtvar z1, rtvar z2, rtvar x0, material* mat);

    virtual bool hit(const ray& sight, rtvar t_min, rtvar t_max, hitInfo& info)const override;

    virtual aabb getbox()const override;

private:
    material * _materialp;

    rtvar _z1, _z2;

    rtvar _y1, _y2;

    rtvar _other;
    };



// the implementation of xy_rect class

inline xy_rect::xy_rect(rtvar x1, rtvar x2, rtvar y1, rtvar y2, rtvar z0, material* mat)
    :_x1(x1)
    , _x2(x2)
    , _y1(y1)
    , _y2(y2)
    , _other(z0)
    , _materialp(mat)
    {
    }

bool xy_rect::hit(const ray& sight, rtvar t_min, rtvar t_max, hitInfo& info)const
    {
    rtvar t = (_other - sight.origin().z()) / sight.direction().z();
    if (t < t_min || t > t_max)return false;

    rtvar x = sight.origin().x() + t*sight.direction().x();
    rtvar y = sight.origin().y() + t*sight.direction().y();
    if (x < _x1 || x > _x2 || y < _y1 || y > _y2)
        return false;

    info._u = (x - _x1) / (_x2 - _x1);
    info._v = (y - _y1) / (_y2 - _y1);
    info._t = t;
    info._materialp = _materialp;
    info._p = sight.go(t);
    info._n = rtvec(0, 0, 1);

    return true;
    }

aabb xy_rect::getbox()const
    {
    return aabb(rtvec(_x1, _y1, _other - 0.0001), rtvec(_x2, _y2, _other + 0.0001));
    }


// the implementation of xz_rect class

inline xz_rect::xz_rect(rtvar x1, rtvar x2, rtvar z1, rtvar z2, rtvar y0, material* mat)
    :_x1(x1)
    , _x2(x2)
    , _z1(z1)
    , _z2(z2)
    , _other(y0)
    , _materialp(mat)
    {
    }

bool xz_rect::hit(const ray& sight, rtvar t_min, rtvar t_max, hitInfo& info)const
    {
    rtvar t = (_other - sight.origin().y()) / sight.direction().y();
    if (t < t_min || t > t_max)return false;

    rtvar x = sight.origin().x() + t*sight.direction().x();
    rtvar z = sight.origin().z() + t*sight.direction().z();
    if (x < _x1 || x > _x2 || z < _z1 || z > _z2)
        return false;

    info._u = (x - _x1) / (_x2 - _x1);
    info._v = (z - _z1) / (_z2 - _z1);
    info._t = t;
    info._materialp = _materialp;
    info._p = sight.go(t);
    info._n = rtvec(0, 1, 0);

    return true;
    }

aabb xz_rect::getbox()const
    {
    return aabb(rtvec(_x1, _other - 0.0001, _z1), rtvec(_x2, _other + 0.0001, _z2));
    }

// the implementation of yz_rect class

inline yz_rect::yz_rect(rtvar y1, rtvar y2, rtvar z1, rtvar z2, rtvar x0, material* mat)
    :_y1(y1)
    , _y2(y2)
    , _z1(z1)
    , _z2(z2)
    , _other(x0)
    , _materialp(mat)
    {
    }

bool yz_rect::hit(const ray& sight, rtvar t_min, rtvar t_max, hitInfo& info)const
    {
    rtvar t = (_other - sight.origin().x()) / sight.direction().x();
    if (t < t_min || t > t_max)return false;

    rtvar y = sight.origin().y() + t*sight.direction().y();
    rtvar z = sight.origin().z() + t*sight.direction().z();
    if (y < _y1 || y > _y2 || z < _z1 || z > _z2)
        return false;

    info._u = (y - _y1) / (_y2 - _y1);
    info._v = (z - _z1) / (_z2 - _z1);
    info._t = t;
    info._materialp = _materialp;
    info._p = sight.go(t);
    info._n = rtvec(1, 0, 0);

    return true;
    }

aabb yz_rect::getbox()const
    {
    return aabb(rtvec(_other - 0.0001, _y1, _z1), rtvec(_other + 0.0001, _y2, _z2));
    }

}// rt namespace 

rectangle.hpp

/// box.hpp
// https://www.cnblogs.com/lv-anchoret/p/10307569.html
// -----------------------------------------------------
// [author]        lv
// [begin ]        2019.1
// [brief ]        the box-class for the ray-tracing project
//                from the 《ray tracing the next week》
// -----------------------------------------------------

#pragma once


namespace rt
{

// the statement of box class

class box: public intersect
    {
public:
    box() {  }

    box(const rtvec& pointmin, const rtvec& pointmax, material * mat);

    virtual bool hit(const ray& sight, rtvar t_min, rtvar t_max, hitInfo& info)const override;

    virtual aabb getbox()const override;
    /*
public:
    // the normal of xy_rect, the normal's direction is +
    inline const rtvec& xy_n()        const    { return _list[0].Getnormal(); }

    // the normal of xy_rect, the normal's direction is -
    inline const rtvec& xy_nflip()    const    { return _list[1].Getnormal(); }
    
    inline const rtvec& xz_n()        const    { return _list[2].Getnormal(); }

    inline const rtvec& xz_nflip()    const    { return _list[3].Getnormal(); }
        
    inline const rtvec& yz_n()        const    { return _list[4].Getnormal(); }

    inline const rtvec& yz_nflip() const    { return _list[5].Getnormal(); }


    inline void set_xy_material(material* m) { _list[0].setmaterial(m); }

    inline void set_xy_flipmaterial(material* m) { _list[1].setmaterial(m); }

    inline void set_xz_material(material* m) { _list[2].setmaterial(m); }

    inline void set_xz_flipmaterial(material* m) { _list[3].setmaterial(m); }

    inline void set_yz_material(material* m) { _list[4].setmaterial(m); }

    inline void set_yz_flipmaterial(material* m) { _list[5].setmaterial(m); }
    */

private:
    rtvec _min;
    
    rtvec _max;
    intersections* _list;
    //rectangles* _list;
    };



// the implementation of box class

inline     box::box(const rtvec& pointmin, const rtvec& pointmax, material * mat)
    :_min(pointmin)
    ,_max(pointmax)
    {
        intersect ** list = new intersect*[6];
        list[0] = new xy_rect(_min.x(), _max.x(), _min.y(), _max.y(), _max.z(), mat);
        list[1] = new flip_normal(new xy_rect(_min.x(), _max.x(), _min.y(), _max.y(), _min.z(), mat));
        list[2] = new xz_rect(_min.x(), _max.x(), _min.z(), _max.z(), _max.y(), mat);
        list[3] = new flip_normal(new xz_rect(_min.x(), _max.x(), _min.z(), _max.z(), _min.y(), mat));
        list[4] = new yz_rect(_min.y(), _max.y(), _min.z(), _max.z(), _max.x(), mat);
        list[5] = new flip_normal(new yz_rect(_min.y(), _max.y(), _min.z(), _max.z(), _min.x(), mat));
        _list = new intersections(list, 6);
        /*
    rectangle ** list = new rectangle*[6];
    list[0] = new xy_rect(_min.x(), _max.x(), _min.y(), _max.y(), _max.z(), mat);
    list[1] = new xy_rect(_min.x(), _max.x(), _min.y(), _max.y(), _min.z(), mat, true);
    list[2] = new xz_rect(_min.x(), _max.x(), _min.z(), _max.z(), _max.y(), mat);
    list[3] = new xz_rect(_min.x(), _max.x(), _min.z(), _max.z(), _min.y(), mat, true);
    list[4] = new yz_rect(_min.y(), _max.y(), _min.z(), _max.z(), _max.x(), mat);
    list[5] = new yz_rect(_min.y(), _max.y(), _min.z(), _max.z(), _min.x(), mat, true);
    _list = new rectangles(list, 6);*/
    }

bool box::hit(const ray& sight, rtvar t_min, rtvar t_max, hitInfo& info)const
    {
    return _list->hit(sight, t_min, t_max, info);
    }

aabb box::getbox()const
    {
    return aabb(_min, _max);
    }

} // rt namespace 

box.hpp

为了统一测试效果，因此把main文件代码也贴在下面了

#define LOWPRECISION
//#define uvwtest
//#define listtest
//#define accumulatetest
//#define attenuationtest
//#define colortest
#define STB_IMAGE_IMPLEMENTATION

#include <fstream>
#include <stb\stb_image.h>
#include "include\texture\RTtexture.hpp"
#include "include\material\RTmaterial.hpp"
#include "include\hit\RThit.hpp"
#include "camera.hpp"
using namespace rt;


rtvec lerp(const ray& sight, intersect* world, int depth)
{
    hitInfo info;
    if (world->hit(sight, (rtvar)0.001, rtInf(), info))
    {
        ray scattered;
        rtvec attenuation;
        rtvec light = info._materialp->emitted(info._u, info._v, info._p);
        if (depth < 50 && info._materialp->scatter(sight, info, attenuation, scattered))
            return light + attenuation * lerp(scattered, world, depth + 1);
        else
            return light;
    }
    else
        return rtvec();
}

void Cornell(intersections** scene, camera** cam, rtvar aspect)
{
    intersect ** list = new intersect*[9];
    size_t cnt = 0;
    material * red = new lambertian(new constant_texture(rtvec(0.65, 0.05, 0.05)));
    material * blue = new lambertian(new constant_texture(rtvec(0.05, 0.05, 0.73)));
    material * white = new lambertian(new constant_texture(rtvec(0.88, 0.88, 0.88)));
    material * green = new lambertian(new constant_texture(rtvec(0.12, 0.45, 0.15)));
    material * light = new areaLight(new constant_texture(rtvec(20, 20, 20)));

    list[cnt++] = new flip_normal(new yz_rect(0, 555, 0, 555, 555, green));
    list[cnt++] = new yz_rect(0, 555, 0, 555, 0, red);
    list[cnt++] = new xz_rect(200, 350, 220, 340, 550, light);
    list[cnt++] = new flip_normal(new xz_rect(200, 350, 220, 340, 550, light));
    list[cnt++] = new flip_normal(new xz_rect(0, 555, 0, 555, 555, white));
    list[cnt++] = new xz_rect(0, 555, 0, 555, 0, white);
    list[cnt++] = new flip_normal(new xy_rect(0, 555, 0, 555, 555, blue));

    list[cnt++] = new translate(new rotate_x(new rotate_y(new box(rtvec(), rtvec(165, 165, 165), white), -45), -30), rtvec(130, 200, 65));
    list[cnt++] = new translate(new rotate_z(new rotate_y(new box(rtvec(), rtvec(165, 330, 165), white), 30), 20), rtvec(265, 60, 295));

    *scene =  new intersections(list, cnt);

    rtvec lookfrom(278, 278, -800);
    rtvec lookat(278, 278, 0);
    rtvar dist_to_focus = 10.0;
    rtvar aperture = 0.;
    rtvar vfov = 40.0;
    *cam = new camera(lookfrom, lookat, rtvec(0, 1, 0),
        vfov, 1, aperture, dist_to_focus, 0., 1.);
}

void build_1_1()
{
    stds ofstream file("graph-1-1.ppm");
    size_t W = 200, H = 200, sample = 200;

    if (file.is_open())
    {
        file << "P3\n" << W << " " << H << "\n255\n" << stds endl;

        intersections * world;
        camera * cma;

        Cornell(&world, &cma, (rtvar)W / H);

        for (int y = H - 1; y >= 0; --y)
            for (int x = 0; x < W; ++x)
            {
                rtvec color;
                for (int cnt = 0; cnt < sample; ++cnt)
                {
                    lvgm::vec2<rtvar> para{
                        (lvgm::rand01() + x) / W,
                        (lvgm::rand01() + y) / H };
                    color += lerp(cma->get_ray(para), world, 0);
                }
                color /= sample;

#ifdef colortest
                if (color.x() < 0)stds cout << "color.x < 0 " << stds endl;
                if (color.y() < 0)stds cout << "color.y < 0 " << stds endl;
                if (color.z() < 0)stds cout << "color.z < 0 " << stds endl;
#endif

                color = rtvec(sqrt(color.r()), sqrt(color.g()), sqrt(color.b()));    //gamma 校订
                int r = int(255.99 * color.r());
                int g = int(255.99 * color.g());
                int b = int(255.99 * color.b());
                file << r << " " << g << " " << b << stds endl;
            }
        file.close();

        if (world)delete world;
        if (cma)delete cma;

        stds cout << "complished" << stds endl;
    }
    else
        stds cerr << "open file error" << stds endl;
}

int main()
{
    build_1_1();
}

 

 Chapter4 Important Sampling Materials

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

让咱们假定咱们发送的一束朝向光源的光线的pdf为pLight(direction)

让咱们假定一个和s有关的pdf函数，咱们叫它为pSurface(direction)

咱们能够把二者进行线性组合，获得一个综合的pdf，最简单的形式以下

p(direction) = 0.5*pLight(direction) + 0.5*pSurface(direction)

其实，只要咱们的权重是正的，而且加起来等于1,也就是咱们常常提到的——pdf积分必定要能积到1，任意pdf函数组合造成的pdf只要知足这两点都是合理的，因此咱们能够将几个份量因素的pdf组合控制，使咱们的pdf函数更加贴近真实效果。虽然咱们可使用任何的pdf，可是它必定不能改变咱们收敛的答案，这是最重要的一点

所以，这就是一个关于寻找pdf的小游戏，要弄清楚如何使s（direction）*color（direction）大的pdf更大，也就是找到合适的pdf，创造更真实的颜色比例效果

对于漫反射（diffuse or Lambertian）表面，咱们猜想它更偏重color（direction）因子。 对于镜面（metal）材质表面来讲，s（）因子只在一个方向附近很大，因此它更重要。

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

 

咱们今天先来测试一下漫反射表面，咱们用以前的Cornell盒子作测试，为了方便，咱们把相机的设置也放在里面，由于不少时候放在外面，常常换场景，没法记录每一个场景采用的合适的摄像机参数，因此就把它和场景放在一块儿了。

already中的main函数中有Cornell函数，这里就不贴了

 

咱们先用以前的方法跑一个200*200的图，sample为200，光线渲染路径计算递归上限50次

可能效果不是很好，可是噪声多一点，方便看效果

 

用时如图 

 

效果如图

咱们的目的是减小噪声，咱们将构造一个pdf，使得更多的光线反射到光源

 

趁着等待渲染的时间，再叨叨两句咱们的核心思想，咱们的渲染路径是这样计算的

首先咱们从eye（或者 camera）发射采样光线，它指向投影屏幕的一个像素位置，也就是当前采样光线正在计算的那个像素位置。

而后光线以射线的形式计算eye与投影屏幕之间全部的几何体交点，获得与eye最近的那个，而后咱们以那个交点为落脚点，做为新的eye，按照光学表面散射原理计算新的射线的方向（结合具体的材质），而后继续测交。

若是在递归深度上限范围内没有找到光源，那么说明该像素位置不会有光传入眼睛，也就是当前像素位置是黑色的；若是它找到光源了，也就是通过屡次散射最后指向了光源，那么就说明光源发出的光线能够沿着咱们的路径逆向进入眼睛，那么咱们就看到了这个像素，那么如何计算像素值呢？

这就涉及到了路径渲染方程，里面有一个rgb份量衰减比例控制参数，它依据材质和纹理计算得出，这是以前的，咱们如今要加上pdf函数，帮助咱们更逼真地计算像素值

好了，唠完了，上图片

貌似没区别，可是咱们下面还会继续优化

刚刚那个图跑了53min，接下来你就会知道为何时间还多了

 

咱们要把重要性采样嵌入进去，首先MC模拟f(x)的积分形式为

　　　　f(r)/p(r)

对于Lambertian材质，咱们用上一篇的pdf函数：

　　　　p(direction) = cosθ/π

因此咱们修改material基类以下： 

/// material.hpp

// -----------------------------------------------------
// [author]        lv
// [begin ]        2018.12
// [brief ]        the material-class for the ray-tracing project
//                from the 《ray tracing in one week》
// -----------------------------------------------------

#pragma once

namespace rt
{

// the statement of material class

class material
    {
public:

    /*
    @brief: produce a scattered ray
    @param: InRay -> Incident light
            info -> the information of intersect-point(hit-point)
            attenuation -> when scattered, how much the ray should be attenuated by tis reflectance R
            scattered -> as we talk, it is a new sight; or
                         it is the scattered ray with the intersect-point
            pdf -> the Important Sample's pdf function
    @retur: the function calculate a scattered ray or not
    */
    virtual bool scatter(const ray& InRay, const hitInfo& info, rtvec& attenuation, ray& scattered, rtvar& pdf)const
        {
        return false;
        }

    /*
    @brief: produce a scattered ray
    @param: 
    InRay -> Incident light
    info -> the information of intersect-point(hit-point)
    scattered -> as we talk, it is a new sight; or
                it is the scattered ray with the intersect-point
    @retur: the function calculate a scattered ray or not
    */
    virtual rtvar scatter_pdf(const ray& InRay, const hitInfo& info, const ray& scattered)const
        {
        return false;
        }

    virtual bool scatter(const ray& rIn, const hitInfo& info, rtvec& attenuation, ray& scattered)const
        {
        return false;
        }

    /*
    @brief: 自发光
    @param: 纹理所需信息
    @retur: 纹理像素值
    */
    virtual rtvec emitted(rtvar u, rtvar v, const rtvec& p)const 
        {
        return rtvec(); 
        }

    };

}

 

/// diffuse.hpp
// https://www.cnblogs.com/lv-anchoret/p/10198423.html
// -----------------------------------------------------
// [author]        lv
// [begin ]        2018.12
// [brief ]        one of the materials
// -----------------------------------------------------

#pragma once

namespace rt
{

class texture;

// the statement of lambertian class

class lambertian : public material
    {
public:
    lambertian(texture* _tex);

    rtvar scatter_pdf(const ray& rIn, const hitInfo& info, const ray& scattered)const;

    bool scatter(const ray& rIn, const hitInfo& info, rtvec& attenuation, ray& scattered, rtvar& pdf)const;

public:
    const texture* get_texture()const { return _albedo; }

    void set_texture(texture* tex) { _albedo = tex; }

protected:

    texture* _albedo;
    };



// the implementation of lambertian class

inline lambertian::lambertian(texture* _tex)
    :_albedo(_tex) 
    {
    }

rtvar lambertian::scatter_pdf(const ray& rIn, const hitInfo& info, const ray& scattered)const
    {
    rtvar cosine = dot(info._n, scattered.direction().ret_unitization()); 
    if (cosine < 0)cosine = 0;
    return cosine / π;
    }

bool lambertian::scatter(const ray& rIn, const hitInfo& info, rtvec& attenuation, ray& scattered, rtvar& pdf)const
    {
    rtvec target = info._p + info._n + lvgm::random_unit_sphere();
    scattered = ray{ info._p, (target - info._p).ret_unitization(), rIn.time() };
    attenuation = _albedo->value(info._u, info._v, info._p);
    pdf = dot(info._n, scattered.direction()) / π;
    return true;
    }

} // rt namespace

 

以及主函数中的lerp函数

rtvec lerp(const ray& sight, intersect* world, int depth)
{
    hitInfo info;
    if (world->hit(sight, (rtvar)0.001, rtInf(), info))
    {
        ray scattered;
        rtvec emitted = info._materialp->emitted(info._u, info._v, info._p);
        rtvar pdf;
        rtvec albedo;
        if (depth < 50 && info._materialp->scatter(sight, info, albedo, scattered, pdf))
            return emitted + albedo *info._materialp->scatter_pdf(sight, info, scattered)*lerp(scattered, world, depth + 1) / pdf;
        else
            return emitted;
    }
    else
        return rtvec();
}

 

 因此，它运行时间多了3min，新代码的运算只增无减

解释一些东西：

scatter_pdf函数里面的cosine值是散射光线和表面法线的夹角，因此大于90°无效，表示采样失败，上一篇说过

scatter函数仍是原来的差很少

 

上述是基于上一本书的第一次修改

下面来进行第二次修改

咱们尝试一下不一样的采样策略，好比，咱们选择从表面的上半球部分随机采样（能够把物体放入球坐标系中，咱们只取上半球部分），那么咱们的p（direction） = 1/（2π）

由于球体的球面度是4πsr，因此半球就是2πsr，采起表面均匀采样的话，就是一立体角对应的球面度量1sr，因此就是1/（2π）

bool lambertian::scatter(const ray& rIn, const hitInfo& info, rtvec& attenuation, ray& scattered, rtvar& pdf)const
    {
    rtvec direction;
    do{
        direction = lvgm::random_unit_sphere();
    } while (dot(direction, info._n) < 0);
    scattered = ray{ info._p,direction.ret_unitization(),rIn.time() };
    attenuation = _albedo->value(info._u, info._v, info._p);
    pdf = 0.5 / π;
    return true;
    }

 

按理说应该获得同上述两张图片相同的图片，可是可能事与愿违

 

书做者采用的是我们2-7教程里面的以下图左上角那张，而咱们今天采用的是右下角那张，由于它让问题看的更清楚

 

下面是书中的结尾

 

其实你看一下这张图顶部的白光（这一变化更容易让咱们捕捉到），就能够验证做者所提出的问题

 

 

原本是周六就写了这篇的草稿的，可是由于抄错代码了，出不来效果，白天调了很久。。。可能不适合抄东西。。

 

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