【RAY TRACING THE REST OF YOUR LIFE 超详解】 光线追踪 3-5 random direction & ONB

 

 Preface

日后看了几章，对这本书有了新的理解

上一篇，咱们第一次尝试把MC积分运用到了Lambertian材质中，固然，第一次尝试是失败的，做者发现它的渲染效果和现实有些出入，因此结尾处声明要经过实践，改进当前的效果

因而乎，就有了后面的章节，几乎整本书都在讲，如何一步一步地改进上一篇的画质，使其更加符合现实，上一篇实际上是抛砖引玉

这本书的小标题名为the rest of your life

经过前面几章，咱们能够更好地理解这句话：咱们经过MC积分优化效果，采用的是pdf函数，以前说过，这就是一场游戏：寻找一个pdf函数，使得使用它进行重要性采样获得的渲染图形更加贴合实际，其实它是没有止境的，好比pdf是一次曲线、二次曲线、高次曲线、正态分布、高斯分布等等，对应的研究方法也是没有止境的，好比：你能够经过对光源进行pdf采样实现最终目的（好比在双向追踪中，光源也要发射光线），你也能够经过对不一样材质表面的反射状态进行pdf采样，进而使得表面颜色变化更光滑更柔和更贴合实际。

上述为我的理解，可能有些出入，吾姑妄言之，汝姑妄听之，便罢。

 

 Ready

上述说到抛砖引玉，可是好像咱们用的不是一张图，思量再三，仍是先把砖整一个，毕竟以后都是围绕那块砖评说效果的，另辟蹊径可能不是明智之举

因此，咱们先把砖搞到手

造砖的代码：

void Cornell(intersections** scene, camera** cam, rtvar aspect)

{

      intersect ** list = new intersect*[8];

      size_t cnt = 0;

      material * red = new lambertian(new constant_texture(rtvec(0.65, 0.05, 0.05)));

      material * white = new lambertian(new constant_texture(rtvec(0.73, 0.73, 0.73)));

      material * green = new lambertian(new constant_texture(rtvec(0.12, 0.45, 0.15)));

      material * light = new areaLight(new constant_texture(rtvec(15, 15, 15)));

 

      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(213, 343, 227, 332, 554, 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, white));

      list[cnt++] = new translate(new rotate_y(new box(rtvec(), rtvec(165, 165, 165), white), -18), rtvec(130, 0, 65));

      list[cnt++] = new translate(new rotate_y(new box(rtvec(), rtvec(165, 330, 165), white), 15), rtvec(265, 0, 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.);

}

 

为了清晰点，sample改成了250，图片为200*200，为了方便重复作实验，因此参数就这样吧，能看清便可，"高清大图"实在是熬不起，都是夜啊，各位见谅~

上一篇结束以后的代码均不变，只是改一下Cornell 函数

咱们获得

 

这就是上一篇最后获得的效果，没错

咱们评说一下，这张图不只没有减小噪点，并且高的长方体表面的颜色趋于均匀一致，与实际有误差

评说好坏固然要有个参考，咱们放上两张第二本书中获得的图形（未加入MC积分）

   

 

上面三张图，不管你看哪张图都能发现，正对咱们的那个长方体表面是从上到下又黑到白渐变的，并且在正方体上表面高度处对应的长方体表面有一抹正方体上表面反射的白色光，而MC图形基本上呈均匀色调，甚至可能上面部分还稍白一些

 

因此咱们进入今天的这一篇

 

 Content

今天咱们讲两章，第一章是讲三维随机方向向量的生成，这有什么用呢，它给咱们的三维空间内进行重要性采样用，MC须要它！！

Chapter5 Generating Random Direction

在这一章和后续的两章，咱们须要增强咱们的理解和咱们手中的工具，搞明白什么才是正确的Cornell Box

让咱们首先从如何生成随机方向开始提及。

为了方便，咱们假定z轴为表面法线，以后再转换坐标系，与此同时，咱们规定θ为从法线张开的角度

咱们将仅仅处理关于z轴旋转对称的分布，因此其余相关的量，均设为均匀分布

 

给定一个和方向相关的pdf，p(direction) = f(θ)，一维pdf中θ 和 φ以下：

g(φ) = 1/(2π)　　　(均匀分布)

h(θ) = 2πf(θ)sinθ

 

对于两个随机生成的均匀变量r1和r2，咱们在第三篇内容中推导出

r₁ = ∫_0->φ 1/(2π) dθ = φ/(2π)

得   φ = 2πr₁

r₂ =  ∫_0->θ 2πf(t)sin(t) dt

t只是一个虚拟的量，用以代替变化的θ，而后辅助实现从0~θ对咱们以前的被积函数f(θ)积分（变限积份量不能相同，因此引入t）

 

咱们从下面开始尝试不一样的f()

首先，咱们考虑球体上采起均匀密度采样，由于整个球面为4πsr，故单位球体表面均匀密度采样的pdf函数为：p(direction) = 1/(4π)，因此获得：

r₂ = ∫_0->θ sin(t)/2 dt

 　= （-cos(t)/2）|_0->θ

　 = （1-cosθ）/2

cosθ = 1 - 2r₂

通常cosθ更经常使用一些，就不进一步求θ了，用的时候再acos()

为了在笛卡尔坐标系下生成一个指向（θ，φ）的单位向量，这就涉及到咱们的球坐标代换方程了：

x = cosφ sinθ

y = sinφ sinθ

z = cosθ

它对应的球坐标以下

此处r为1

 

咱们开始推导

把关于θ和φ的式子带入x,y,z中

x = cos(2πr₁)*sqrt(1-cos²θ)

y = sin(2πr₁)*sqrt(1-cos²θ)

z = 1-2r₂

解得：

x = cos(2πr₁)*sqrt(r₂*(1-r₂))

y = sin(2πr₁)*sqrt(r₂*(1-r₂))

z = 1-2r₂

 

而后咱们来作个图，验证一下它是否如咱们所愿，生成的点聚拢为单位球面

做者给出的画图方法是plot.ly，能够在线画，可是它要求注册帐号，流程还挺麻烦的，GitHub受权帐号总是没响应（可能我网很差），若是有兴趣的能够直接到这里

直接能够载入数据使用的网址

其实matlab最好使了，对不对啊

因而乎，咱们先写入txt,千万不要写入xls，仍是比较麻烦的

 

    stds ofstream outfile("random_direction.txt");

    for (int i = 0; i < 200; ++i)
    {
        double r1 = lvgm::rand01();
        double r2 = lvgm::rand01();
        double x = cos(2 * π * r1) * 2 * sqrt(r2 * (1 - r2));
        double y = sin(2 * π * r1) * 2 * sqrt(r2 * (1 - r2));
        double z = 1 - 2 * r2;
        outfile << x << "\t" << y << "\t" << z << stds endl;
    }
    
    outfile.close();

 

建立一个random_direction.xls的文件，点菜单栏中的数据->自文本，选择random_direction.txt，一路回车就加载完毕了，可见xls加载文本容易多了

假设你的路径是

E:\OpenGL\光线追踪\code\ray tracing 1-3\ray tracing 1-3\random_direction.xls

数据所在表的表名为sheet1

则matlab代码和效果以下图

 

 效果仍是很好的

咱们能够看到，它是均匀随机的，达到了咱们的预期

 

呐，咱们接下来试一下咱们第二经常使用的pdf

即p(direction) = cosθ/π

r₂ = ∫_0->θ 2π（cos(t)/π）sint dt

   = (-cos²t)|_0->θ

   = 1 - cos²θ

cosθ = sqrt(1-r₂)

因而乎咱们获得下述的x,y,z

x = cos(2πr₁)*sqrt(r₂)

y = sin(2πr₁)*sqrt(r₂)

z = sqrt(1-r₂₎

 

下面咱们就生成随机向量

inline rtvec random_cosine_direction()
{
    double r1 = lvgm::rand01();
    double r2 = lvgm::rand01();
    double z = sqrt(1 - r2);
    double φ = 2 * π*r1;
    double x = cos(φ) * 2 * sqrt(r2);
    double y = sin(φ) * 2 * sqrt(r2);
    return rtvec(x, y, z);
}

 

咱们按照原来的方法把第二个pdf函数产生的随机状况模拟一下：

 

 

那么咱们用pdf作一个数值模拟

 

void estimate2()
{
    int n = 1000000;
    double sum = 0.;
    for (int i = 0; i < n; ++i)
    {
        rtvec v = ::random_cosine_direction();
        sum += pow(v.z(), 3) / (v.z() / π);
    }
    stds cout << "π/2 = " << π / 2 << stds endl;
    stds cout << "Estimate = " << sum / n << stds endl;
}

 

还有一个模拟半球获得的近似，pdf = 1/（2π）

有兴趣的能够本身推一下，这里直接给代码，为了和上面作对比

void estimate()
{
    int n = 1000000;
    double sum = 0.;
    for (int i = 0; i < n; ++i)
    {
        double r1 = lvgm::rand01();
        double r2 = lvgm::rand01();
        double x = cos(2 * π * r1) * 2 * sqrt(r2 * (1 - r2));
        double y = sin(2 * π * r1) * 2 * sqrt(r2 * (1 - r2));
        double z = 1 - r2;
        sum += z*z*z / (1. / (2.*π));
    }
    stds cout << "π/2 = " << π / 2 << stds endl;
    stds cout << "Estimate = " << sum / n << stds endl;
}

 

主函数

int main()
{
    //build_1_1();

    estimate();
    estimate2();
}

 

结果

能够看到模拟半球的pdf函数效果没有p() = cos/π 这个好

这一章中全部的都是基于z轴的，而下一章咱们真正基于物体表面法线进行

 

Chapter6 Ortho-normal Bases

ONB是三个相互正交的单位向量集合，笛卡尔坐标系中的xyz轴也是一种ONB

咱们这一章的目的就是把上一章基于z轴的随机方向转换到基于表面法线的

 

假设，咱们有一个原点o和笛卡尔坐标系向量 x/y/z，当咱们描述一个位置

好比：（3，-2,7）时，咱们这样计算：

location = o + 3x - 2y +7z

假若有另一个坐标系，它的原点为o'，基向量为 u/v/w

咱们将要用以下方式表示（u，v，w）这个位置：

location = o' + uu + vv + ww

 

若是学过计算机图形学，那么你就能够用矩阵变换去完成坐标系，可是咱们这里不须要这种操做。

咱们须要的是生成具备相对于表面法线的集合分布的随机方向。 咱们不须要原点，由于向量无起点。 

咱们引用光线追踪1-8中所述的相机坐标来计算ONB的三个基向量

 

咱们如今拥有的是法向量n，咱们能够把它看作是lookfrom->lookat向量

咱们须要一个vup,假设法向量n自己几乎平行于特定轴，那么vup就取和n垂直的基向量，反之，咱们就使用特定轴做为vup

代码是描述思惟最好的方式之一:(假设特定轴为x轴)

if(fabs(n.x())>0.9)　//单位法向量n的x份量大于0.9，认为n//x

　　vup = (0,1,0)

else

　　vup = (1,0,0)

 

咱们设ONB基向量由 s,t,n 组成

那么 t = cross(vup,n).单位化

　　 s = cross(t,n)

理解可参考上图,t为图中的u, s为图中的v

 

因此咱们的坐标系中的任一点（x,y,z）表示以下
随机向量 = xs + yt + zn

 

至此，咱们上代码：

/// onb.hpp
// 
// -----------------------------------------------------
// [author]        lv
// [begin ]        2019.3
// [brief ]        ONB
// -----------------------------------------------------

#pragma once

namespace rt
{

class onb
    {
public:

    onb() {  }

    inline const rtvec& operator[](int index)const { return axis[index]; }
        
    inline const rtvec& u()const { return axis[0]; }
        
    inline const rtvec& v()const { return axis[1]; }

    inline const rtvec& w()const { return axis[2]; }

public:

    inline rtvec local(double a, double b, double c)const;

    inline rtvec local(const rtvec& v)const;

    void build_from_w(const rtvec&);

private:
        
    rtvec axis[3];

    };



inline rtvec onb::local(double a, double b, double c)const
    {
    return a*u() + b*v() + c*w();
    }

inline rtvec onb::local(const rtvec& v)const
    {
    return local(v.x(), v.y(), v.z());
    }

inline void onb::build_from_w(const rtvec& V)
    {
    axis[2] = V.ret_unitization();
    rtvec a;
    if (fabs(w().x()) > 0.9)
        a = rtvec(0, 1, 0);
    else
        a = rtvec(1, 0, 0);
    axis[1] = lvgm::cross(w(), a).ret_unitization();
    axis[0] = lvgm::cross(w(), v());
    }

} // rt namespace

 

咱们的Lambertian材质的scatter函数须要作相应的改动

原书没有do-while循环，可是可能会出现pdf为零的除零错误，因此咱们生成的随机方向必须可以算出一个非零的pdf，这并不妨碍咱们的算法，由于这一切都是随机的，随机错误就再随机一次

bool lambertian::scatter(const ray& rIn, const hitInfo& info, rtvec& alb, ray& scattered, rtvar& pdf)const
    {
    onb uvw;
    uvw.build_from_w(info._n);
    do {
        rtvec direction = uvw.local(random_cosine_direction());
        scattered = ray(info._p, direction.ret_unitization(), rIn.time());
        pdf = dot(uvw.w(), scattered.direction()) / π;
    } while (pdf == rtvar(0));
    alb = _albedo->value(info._u, info._v, info._p);
    return true;
    }

 

而后咱们跑一遍场景

感受还不是很真实，好像和之前的差很少（正常，这不能算鸽~，由于第八章之后差很少才能够。。）

探索过程也是很值得借鉴的嘛，毕竟这是在一步一步引入新的技术，算是循循善诱~

下面是原书结尾，下一篇咱们讲牛逼的技术——直接光源采样

 

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