【Ray Tracing The Next Week 超詳解】 光線追蹤2-2

 

Chapter 2：Bounding Volume Hierarchies

今天我們來講層次包圍盒，乍一看比較難，篇幅也多，但是咱們一步一步來，相信大家應該都能聽懂

 

BVH 和 Perlin textures是這本書中最難的兩章，為什么把BVH放在第二章講呢，據說，層次包圍盒的渲染效率是比較高的，考慮到大家的渲染時間開銷，所以先講這個

光線 - 物體相交是光線跟蹤器中的主要時間瓶頸，運行時間與物體數量成線性關系。但它是對同一模型的重復搜索，所以我們可以采用二分搜索的速度進行對數級復雜度搜索。因為我們在同一模型上發送數百萬到數十億的光線，我們可以對模型進行分類，然后每個光線交叉點可以是次線性搜索。兩個最常見的兩類是

1）划分空間

2）划分對象

后者通常更容易編碼，並且運行速度與大多數模型一樣快。   關鍵思想是找到一個完全包圍（邊界）所有對象的體積。例如，你計算得出了10個對象包圍盒。任何未與包圍盒相交的射線肯定不會和十個對象相交。如果射線擊中了包圍盒，那么它可能會擊中十個物體中的一個。所以偽代碼：

 

if (ray hits bounding object)
    return whether ray hits bounded objects
else
    return false

 

最關鍵的一件事情是我們如何將各個物體划分到每個子集中，單個子集為一個包圍盒

引用書上一張圖

藍色和紅色邊界包含在紫色邊界中，但它們可能重疊，並且它們不是有序的，它們只是單純地被包含在內部。

右邊是左圖的樹結構

對應的檢測偽代碼是：

if (hits purple)
    hit0 = hits blue enclosed objects
    hit1 = hits red enclosed objects
if (hit0 or hit1)
  return true and info of closer hit
return false  

好了，我們下面就來實現上面的偽代碼框架

 

我們需要將場景中的物體進行划分，且包圍盒需要非常緊湊，以及考慮光線與包圍盒相交的方法，計算量盡量少。 在大多數模型的實踐中，軸對齊的盒子（axis - aligned bounding box，即AABB）比較好一些,我們只需要知道光線是否擊中盒體，而無需操心撞擊點的任何信息

有一種常用的“slab”的方法，它是基於n個維度的AABB，就是取n個坐標軸上的區間表示，稱為“slabs”

一維空間的一個區間，比如：x∈【3,5】，它是一條線段

二維空間的一個區間，比如：x∈【3,5】，y∈【3,5】,它是一塊矩形區域

　　

我們來確定光線與區間的相交信息

　　

>假設，上述圖中的這種情況，我們的光線和 x = x0,  x = x1 相交於 t0 和 t1

回顧視線方程：p(t) = a + t *b，若為x區間相交，則方程寫為 x(t) = a + t * b

則 x0 = a.x + t0 * b.x  =>  t0 = (x0 - a.x) / b.x

同理可得 t1 = (x1 - a.x) / b.x

 

>如果是二維的情況，那么，就要加上y(t) = a + t * b

我們的工作就是：

計算Q1(tx0,tx1)

計算Q2(ty0,ty1)

Q1 和 Q2 是否有交集

大致幾種分類情況如下：（下面紅色代表ty,綠色代表tx）

 

 

 

 

 

 

綜上，我們發現：

如果　兩個區間的左端點最大值 小於 右端點的最小值

　　光線一定和區域有交點

反之

　　光線和區域相離

 

>如果是三維的,那么同理，步驟如下：

計算Q1（tx0,tx1）

計算Q2（ty0,ty1）

計算Q3（tz0,tz1）

Q1、Q2和Q3是否有交集

 

代碼簡單描寫為：

_min 和 _max 分別指的是包圍盒中三個維度區間的左端點集合和右端點集合

它們是aabb的兩個數據成員

（這個代碼只用來理解用）

 

 后來我們觀察發現，一般情況下滿足 t0 <= t1 ，但是有時候 t0 > t1

當且僅當視線的方向向量在當前計算維度的分量中為負，此時 t0 > t1

所以，我們改寫成如下形式：

inline bool aabb::hit(const ray& sight, rtvar tmin, rtvar tmax)const
    {
    for (int i = 0; i < 3; ++i)
        {
        rtvar div = 1.0 / sight.direction()[i];
        rtvar t1 = (_min[i] - sight.origin()[i]) / sight.direction()[i];
        rtvar t2 = (_max[i] - sight.origin()[i]) / sight.direction()[i];
        if (div < 0.)stds swap(t1, t2);
        if (stds min(t2, tmax) <= stds max(t1, tmin))
            return false;
        }
    return true;
    }

同時，我們的相交類也要改一下

 

我們從現在開始增加各個子類實現

 關於sphere，球體的三個維度的左端點集合和右端點集合分別為

aabb sphere::getbox()const
    {
    return aabb(_heart - rtvec(_radius, _radius, _radius), _heart + rtvec(_radius, _radius, _radius));
    }

對於moving_sphere我們需要綜合開始時刻和結束時刻兩個球的盒體的邊界

aabb _surrounding_box(aabb box1, aabb box2);

但是，出於某種考慮，我覺得把它放在aabb盒體類中作為靜態成員函數比較好

 

/// aabb_box.hpp

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

namespace rt
{

//the statement of aabb class

class aabb
    {
public:
    aabb() {  }

    aabb(const rtvec& a, const rtvec& b);

    inline bool hit(const ray& sight, rtvar tmin, rtvar tmax)const;

    static aabb _surrounding_box(aabb box1, aabb box2);
    
public:

    inline rtvec min()const { return _min; }

    inline rtvec max()const { return _max; }


private:
    rtvec _min;

    rtvec _max;
    };


//the implementation of aabb class

inline aabb::aabb(const rtvec& a, const rtvec& b)
    :_min(a)
    , _max(b)
    {
    }

inline bool aabb::hit(const ray& sight, rtvar tmin, rtvar tmax)const
    {
    for (int i = 0; i < 3; ++i)
        {
        rtvar div = 1.0 / sight.direction()[i];
        rtvar t1 = (_min[i] - sight.origin()[i]) / sight.direction()[i];
        rtvar t2 = (_max[i] - sight.origin()[i]) / sight.direction()[i];
        if (div < 0.)stds swap(t1, t2);
        if (stds min(t2, tmax) <= stds max(t1, tmin))
            return false;
        }
    return true;
    }

aabb aabb::_surrounding_box(aabb box1, aabb box2)
    {
    auto fmin = [](const rtvar a, const rtvar b) {return a < b ? a : b; };
    auto fmax = [](const rtvar a, const rtvar b) {return a > b ? a : b; };
    rtvec min{    fmin(box1.min().x(),box2.min().x()),
                fmin(box1.min().y(),box2.min().y()),
                fmin(box1.min().z(),box2.min().z()) };
    rtvec max{    fmax(box1.max().x(),box2.max().x()),
                fmax(box1.max().y(),box2.max().y()),
                fmax(box1.max().z(),box2.max().z()) };
    return aabb(min, max);
    }

}

 

aabb moving_sphere::getbox()const
    {
    rtvec delt{ _radius, _radius, _radius };
    return aabb::_surrounding_box(aabb(_heart1 - delt, _heart1 + delt), aabb(_heart2 - delt, _heart2 + delt));
    }

 

現在我們開始着手，划分物體，並解決“光線是否擊中了當前盒體”這個開篇的問題

首先，我們需要創建像開篇那張圖中的一顆盒體范圍樹

樹節點定義：

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

    bvh_node(intersect** world, const int n, const rtvar time1, const rtvar time2);

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

    virtual aabb getbox()const override;

private:
    intersect* _left;

    intersect* _right;

    aabb _box;

    };

 

aabb bvh_node::getbox()const
    {
    return _box;
    }

 

 構造函數中那兩個時間實在不知道有什么用（=.=）

 

之后我們就需要寫hit函數了，其實很好寫

樹結構，遍歷左子樹遍歷右子樹，返回離eye最近的撞擊點信息即可

bool bvh_node::hit(const ray& sight, rtvar t_min, rtvar t_max, hitInfo& info)const
    {
    if (_box.hit(sight, t_min, t_max))
        {
        hitInfo linfo, rinfo;
        bool lhit = _left->hit(sight, t_min, t_max, linfo);
        bool rhit = _right->hit(sight, t_min, t_max, rinfo);
        if (lhit && rhit)
            {
            if (linfo._t < rinfo._t)
                info = linfo;
            else
                info = rinfo;
            return true;
            }
        else if (lhit)
            {
            info = linfo;
            return true;
            }
        else if (rhit)
            {
            info = rinfo;
            return true;
            }
        else 
            return false;
        }
    else 
        return false;
    }

}

 

構造函數設計：

1）隨機選擇一個軸
2）使用庫qsort對物體進行排序
3）在每個子樹中放一半物體
 

並且特判了兩種情況（物體個數為1 或者 2）
如果我只有一個元素，我在每個子樹中復制它。兩個物體的話，一邊一個。

明確檢查三個元素並且只跟隨一個遞歸可能會有所幫助，但我認為整個方法將在以后進行優化。即：

inline bvh_node::bvh_node(intersect** world, const int n, const rtvar time1, const rtvar time2)
    {
    int axis = static_cast<int>(3 * lvgm::rand01());
    if (axis == 0)
        qsort(world, n, sizeof(intersect*), _x_cmp);
    else if (axis == 1)
        qsort(world, n, sizeof(intersect*), _y_cmp);
    else
        qsort(world, n, sizeof(intersect*), _z_cmp);

    if (n == 1)
        _left = _right = world[0];
    else if (n == 2)
        _left = world[0],
        _right = world[1];
    else
        _left = new bvh_node(world, n / 2, time1, time2),
        _right = new bvh_node(world + n / 2, n - n / 2, time1, time2);

    aabb lbox = _left->getbox();
    aabb rbox = _right->getbox();

    _box = aabb::_surrounding_box(lbox, rbox);
    }

 

比較函數以_x_cmp為例：

inline    int _x_cmp(const void * lhs, const void * rhs)
    {
    intersect * lc = *(intersect**)lhs;
    intersect * rc = *(intersect**)rhs;
    aabb lbox = lc->getbox();
    aabb rbox = rc->getbox();

    if (lbox.min().x() - rbox.min().x() < 0.)
        return -1;
    else
        return 1;
    }

 

整個方法在之后可能會優化，但是目前確實不咋好：

 

測試代碼

#define LOWPRECISION

#include <fstream>
#include "RTmaterial.hpp"
#include "RThit.hpp"
#include "camera.hpp"
using namespace rt;

int Cnt;

rtvec lerp(const ray& sight, intersect* world, int depth)
{
    hitInfo info;
    if (world->hit(sight, (rtvar)0.001, rtInf(), info))
    {
        ray scattered;
        rtvec attenuation;
        if (depth < 50 && info.materialp->scatter(sight, info, attenuation, scattered))
            return attenuation * lerp(scattered, world, depth + 1);
        else
            return rtvec(0, 0, 0);
    }
    else
    {
        rtvec unit_dir = sight.direction().ret_unitization();
        rtvar t = 0.5*(unit_dir.y() + 1.);
        return (1. - t)*rtvec(1., 1., 1.) + t*rtvec(0.5, 0.7, 1.0);
    }
}

intersect* random_sphere()
{
    int cnt = 50000;
    intersect **list = new intersect*[cnt + 1];
    list[0] = new sphere(rtvec(0, -1000, 0), 1000, new lambertian(rtvec(0.5, 0.5, 0.5)));
    int size = 1;
    for (int a = -5; a < 5; ++a)
        for (int b = -5; b < 5; ++b)
        {
            rtvar choose_mat = lvgm::rand01();
            rtvec center(a + 0.9 * lvgm::rand01(), 0.2, b + 0.9*lvgm::rand01());
            if ((center - rtvec(4, 0.2, 0)).normal()>0.9)
            {
                if (choose_mat < 0.55)
                    list[size++] = new moving_sphere(center, center + rtvec(0, 0.5*lvgm::rand01(), 0), 0., 1., 0.2,
                        new lambertian(rtvec(lvgm::rand01()*lvgm::rand01(), lvgm::rand01()*lvgm::rand01(), lvgm::rand01()*lvgm::rand01())));

                else if (choose_mat < 0.85)
                    list[size++] = new sphere(center, 0.2,
                        new metal(rtvec(0.5*(1 + lvgm::rand01()), 0.5*(1 + lvgm::rand01()), 0.5*(1 + lvgm::rand01())), 0.5*lvgm::rand01()));

                else
                    list[size++] = new sphere(center, 0.2,
                        new dielectric(1.5));
            }
        }

    list[size++] = new sphere(rtvec(0, 1, 0), 1.0, new dielectric(1.5));
    list[size++] = new moving_sphere(rtvec(-4.5, 1, 0.65), rtvec(-4,1,0.15), 0., 1., 1.0, 
        new lambertian(rtvec(0.4, 0.2, 0.1)));
    list[size++] = new sphere(rtvec(4, 1, 0), 1.0, new metal(rtvec(0.7, 0.6, 0.5), 0.));

    return new intersections(list, size);
}

void build_12_1()
{
    stds ofstream file("graph1-2.ppm");
    size_t W = 200, H = 100, sample = 100;

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

        intersect* world = random_sphere();

        rtvec lookfrom(13, 2, 3);
        rtvec lookat(0, 0, 0);
        float dist_to_focus = 10.0;
        float aperture = 0.0;
        camera cma(lookfrom, lookat, rtvec(0, 1, 0), 20, rtvar(W) / rtvar(H), aperture, 0.7*dist_to_focus, 0., 1.);

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

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

int main()
{
    build_12_1();
}

 

 

有點傷心且不知所措，但其實還是學了很多相交的知識的~

 

感謝您的閱讀，生活愉快~