来自 技术 2019-01-18 00:00 的文章

【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 objectselse return false

最关键的一件事情是我们如何将各个物体划分到每个子集中,单个子集为一个包围盒

引用书上一张图

蓝色和红色边界包含在紫色边界中,但它们可能重叠,并且它们不是有序的,它们只是单纯地被包含在内部。

右边是左图的树结构

对应的检测伪代码是:

if (hits purple) hit0 = hits blue enclosed objects hit1 = hits red enclosed objectsif (hit0 or hit1) return true and info of closer hitreturn 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 classclass 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 classinline 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();}main.cpp

有点伤心且不知所措,但其实还是学了很多相交的知识的~

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