【Ray Tracing The Next Week 超詳解】 光線追蹤2-2
Chapter 2：Bounding Volume Hierarchies
今天我們來講層次包圍盒，乍一看比較難，篇幅也多，但是咱們一步一步來，相信大家應該都能聽懂
BVH 和 Perlin textures是這本書中最難的兩章，為什么把BVH放在第二章講呢，據(jù)說，層次包圍盒的渲染效率是比較高的，考慮到大家的渲染時間開銷，所以先講這個
光線 - 物體相交是光線跟蹤器中的主要時間瓶頸，運行時間與物體數(shù)量成線性關(guān)系。但它是對同一模型的重復搜索，所以我們可以采用二分搜索的速度進行對數(shù)級復雜度搜索。因為我們在同一模型上發(fā)送數(shù)百萬到數(shù)十億的光線，我們可以對模型進行分類，然后每個光線交叉點可以是次線性搜索。兩個最常見的兩類是
1）劃分空間
2）劃分對象
后者通常更容易編碼，并且運行速度與大多數(shù)模型一樣快。 關(guān)鍵思想是找到一個完全包圍（邊界）所有對象的體積。例如，你計算得出了10個對象包圍盒。任何未與包圍盒相交的射線肯定不會和十個對象相交。如果射線擊中了包圍盒，那么它可能會擊中十個物體中的一個。所以偽代碼：
if (ray hits bounding object) return whether ray hits bounded objects else return false
最關(guān)鍵的一件事情是我們?nèi)绾螌⒏鱾€物體劃分到每個子集中，單個子集為一個包圍盒
引用書上一張圖
藍色和紅色邊界包含在紫色邊界中，但它們可能重疊，并且它們不是有序的，它們只是單純地被包含在內(nèi)部。
右邊是左圖的樹結(jié)構(gòu)
對應的檢測偽代碼是：
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
好了，我們下面就來實現(xiàn)上面的偽代碼框架
我們需要將場景中的物體進行劃分，且包圍盒需要非常緊湊，以及考慮光線與包圍盒相交的方法，計算量盡量少。 在大多數(shù)模型的實踐中，軸對齊的盒子（axis - aligned bounding box，即AABB）比較好一些,我們只需要知道光線是否擊中盒體，而無需操心撞擊點的任何信息
有一種常用的“slab”的方法，它是基于n個維度的AABB，就是取n個坐標軸上的區(qū)間表示，稱為“slabs”
一維空間的一個區(qū)間，比如：x∈【3,5】，它是一條線段
二維空間的一個區(qū)間，比如：x∈【3,5】，y∈【3,5】,它是一塊矩形區(qū)域
我們來確定光線與區(qū)間的相交信息
>假設(shè)，上述圖中的這種情況，我們的光線和 x = x0, x = x1 相交于 t0 和 t1
回顧視線方程：p(t) = a + t *b，若為x區(qū)間相交，則方程寫為 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）
綜上，我們發(fā)現(xiàn)：
如果 兩個區(qū)間的左端點最大值 小于 右端點的最小值
光線一定和區(qū)域有交點
反之
光線和區(qū)域相離
>如果是三維的,那么同理，步驟如下：
計算Q1（tx0,tx1）
計算Q2（ty0,ty1）
計算Q3（tz0,tz1）
Q1、Q2和Q3是否有交集
代碼簡單描寫為：
_min 和 _max 分別指的是包圍盒中三個維度區(qū)間的左端點集合和右端點集合
它們是aabb的兩個數(shù)據(jù)成員
（這個代碼只用來理解用）
后來我們觀察發(fā)現(xiàn)，一般情況下滿足 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; }
同時，我們的相交類也要改一下
我們從現(xiàn)在開始增加各個子類實現(xiàn)
關(guān)于sphere，球體的三個維度的左端點集合和右端點集合分別為
aabb sphere::getbox()const { return aabb(_heart - rtvec(_radius, _radius, _radius), _heart + rtvec(_radius, _radius, _radius)); }
對于moving_sphere我們需要綜合開始時刻和結(jié)束時刻兩個球的盒體的邊界
aabb _surrounding_box(aabb box1, aabb box2);
但是，出于某種考慮，我覺得把它放在aabb盒體類中作為靜態(tài)成員函數(shù)比較好
/// 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)); }
現(xiàn)在我們開始著手，劃分物體，并解決“光線是否擊中了當前盒體”這個開篇的問題
首先，我們需要創(chuàng)建像開篇那張圖中的一顆盒體范圍樹
樹節(jié)點定義：
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; }
構(gòu)造函數(shù)中那兩個時間實在不知道有什么用（=.=）
之后我們就需要寫hit函數(shù)了，其實很好寫
樹結(jié)構(gòu)，遍歷左子樹遍歷右子樹，返回離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; } }
構(gòu)造函數(shù)設(shè)計：
1）隨機選擇一個軸
2）使用庫qsort對物體進行排序
3）在每個子樹中放一半物體
并且特判了兩種情況（物體個數(shù)為1 或者 2）
如果我只有一個元素，我在每個子樹中復制它。兩個物體的話，一邊一個。
明確檢查三個元素并且只跟隨一個遞歸可能會有所幫助，但我認為整個方法將在以后進行優(yōu)化。即：
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); }
比較函數(shù)以_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; }
整個方法在之后可能會優(yōu)化，但是目前確實不咋好：
測試代碼
#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(); }
有點傷心且不知所措，但其實還是學了很多相交的知識的~
感謝您的閱讀，生活愉快~
浙公網(wǎng)安備 33010602011771號