Note on "Ray Tracing The Next Week" - Part 2

Note on "Ray Tracing The Next Week" - Part 2

Ayano Kagurazaka Lv3
  2024-02-10 2024-02-10 Created   2025-01-15 16:22:46 2025-01-15 16:22:46 Updated      2k Words  12 Mins  

From linear to sub-linear

If we take a look at the hit function in HittableList class:

1
2
3
4
5
6
7
8
9
10
11
12
13
bool HittableList::hit(const Ray &r, Interval interval, HitRecord &record) const {
    HitRecord temp_rec;
    bool if_hit = false;
    auto closest_t = interval.max;
    for (const auto& i : objects) {
        if (i->hit(r, Interval(interval.min, closest_t), temp_rec)) {
            if_hit = true;
            closest_t = temp_rec.t;
            record = temp_rec;
        }
    }
    return if_hit;
}

We can see it got linear time complexity. It’s not hard to imagine that if we have a lot of objects in the scene, the performance will be terrible. So we need to find a way to improve it. One way to do it is to use a bounding volume hierarchy (BVH).

Bounding Volume Hierarchy

Imagine we have many objects in a scene, we can use a box to wrap all the objects up, and two sub-boxes to wrap half of the objects up. Then we can do the same in the sub-box, and so on. In this case, we have a tree structure. And finding a node in the tree will be O(nlog⁡(n))O(n\log(n)) time complexity, which is much better than O(n)O(n).

One question is how can we warp the objects up? The bounding should be fast and compact. One way to do it is aligning the bounding box with the axis, or make an axis-aligned bounding box(AABB). In R2\mathbb{R}^2 we know if we want to have a bound on something, we take intersection of two intervals. In R3\mathbb{R}^3 we can take intersection of two boxes. So we can use this to warp the objects up.

bound

To detect whether the ray hits the box, we can use the ray’s equation:

R=A+tb\mathbf{R}=\mathbf{A} + t\mathbf{b}

If we only consider the xx axis, we can have:

Rx=Ax+tbxR_x = A_x + t b_x

and planes that is orthogonal to yz plane has this equation:

x=x0x=x1x = x_0 \\ x = x_1

We can then simply plugin the ray’s equation into the plane’s equation and solve for tt:

Similar thing goes for t1t_1. Since one could be either greater or small than each other, we make t0t_0 the max, and t1t_1 the min. Then we can have the interval of tt that the ray hits the box. We can do the same for yy and zz axis. We now define a AABB class to represent the box:

icon-padding

MathUtil.h

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
class AABB {
public:
	Interval x, y, z;

	AABB() = default;

	AABB(const Point3 &a, const Point3 &b);

	AABB(const Interval& x, const Interval& y, const Interval& z);

	Interval axis(int i) const;

	[[nodiscard]] bool hit(const Ray &r, Interval ray_int) const;

};

icon-padding

MathUtil.cpp

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
AABB::AABB(const Point3 &a, const Point3 &b) {
	x = Interval(fmin(a.x, b.x), fmax(a.x, b.x));
	y = Interval(fmin(a.y, b.y), fmax(a.y, b.y));
	z = Interval(fmin(a.z, b.z), fmax(a.z, b.z));
}

AABB::AABB(const Interval &x, const Interval &y, const Interval &z) : x(x), y(y), z(z) {

}

Interval AABB::axis(int i) const {
	switch (i) {
		case 0:
			return x;
		case 1:
			return y;
		case 2:
			return z;
	}
}

bool AABB::hit(const Ray &r, Interval ray_int) const {
	for (int i = 0; i < 3; i++) {
		auto invD = 1.0 / r.dir()[i];
		auto orig = r.pos()[i];

		auto t0 = (axis(i).min - orig) * invD;
		auto t1 = (axis(i).max - orig) * invD;

		if (invD < 0.0) {
			std::swap(t0, t1);
		}

		if (t0 > ray_int.min) ray_int.min = t0;
		if (t1 < ray_int.max) ray_int.max = t1;
		
		if (ray_int.max < ray_int.min) {
			return false;
		}
	}
	return true;
}

Now we can use the AABB class to wrap one stationary sphere up:

1
2
3
4
Sphere::Sphere(double radius, Vec3 position, std::shared_ptr<IMaterial> mat) : radius(radius), position(std::move(position)), material(std::move(mat)) {
	auto rvec = Vec3{radius, radius, radius};
	bbox = AABB(this->position - rvec, this->position + rvec);
}

For moving sphere, we want to wrap the sphere up in the whole time interval. To make our life easier, we define other constructors for Interval and AABB:

icon-padding

MathUtil.h

1
2
3
4
5
6
7
8
9
10
11
12
13
14
class Interval {
public:
	// ...
	Interval(const Interval& first, const Interval& second);
    // ...

};

class AABB {
public:
    // ...
	AABB(const AABB& up, const AABB& down);
    // ...
};

icon-padding

MathUtil.cpp

1
2
3
4
5
6
7
8
9
10
Interval::Interval(const Interval &first, const Interval &second) {
	min = fmin(first.min, second.min);
	max = fmax(first.max, second.max);
}

AABB::AABB(const AABB &up, const AABB &down) {
	x = Interval(up.x, down.x);
	y = Interval(up.y, down.y);
	z = Interval(up.z, down.z);
}

And the AABB of a list of hittable objects is just an AABB that wraps all the AABBs up:

icon-padding

HittableList.cpp

1
2
3
4
5
6
7
8
9
10
11
12
13
14
HittableList::HittableList(const std::shared_ptr<IHittable>& obj) {
    add(obj);
	bbox = obj->boundingBox();
}

void HittableList::add(const std::shared_ptr<IHittable>& obj) {
    objects.push_back(obj);
	bbox = AABB(bbox, obj->boundingBox());
}

void HittableList::clear() {
    objects.clear();
	bbox = AABB();
}

Building a tree

BVHs are also hittables. They are containers, but they can also respond to hits(that’s what we want). Since we will build a binary tree structure, we will define a BVHNode class:

icon-padding

MathUtil.h

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
class BVHNode : public IHittable {
public:
	BVHNode() = default;

	BVHNode(const HittableList& list);

	BVHNode(const std::vector<std::shared_ptr<IHittable>>& objects, size_t start, size_t end);
	
	bool hit(const Ray &r, Interval interval, HitRecord &record) const override;
	
	AABB boundingBox() const override;
private:
	std::shared_ptr<IHittable> left;
	std::shared_ptr<IHittable> right;
	AABB bbox;
};

To divide the whole scene into tree, we will follow these steps:

  1. Randomly choose one axis to divide the scene into two parts.
  2. Sort the objects based on the axis index.
  3. Recursively divide the scene into two parts.

This diagram shows how the tree is built:

bvhTree

and this is the corresponding bounding box:

To sort the objects, we need to define the comparison function:

icon-padding

MathUtil.h

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
class BVHNode : public IHittable {
public:
	// ...
private:
	std::shared_ptr<IHittable> left;
	std::shared_ptr<IHittable> right;

	static bool compare(const std::shared_ptr<IHittable>& a, const std::shared_ptr<IHittable>& b, int axis);

	static bool compareX(const std::shared_ptr<IHittable>& a, const std::shared_ptr<IHittable>& b);

	static bool compareY(const std::shared_ptr<IHittable>& a, const std::shared_ptr<IHittable>& b);

	static bool compareZ(const std::shared_ptr<IHittable>& a, const std::shared_ptr<IHittable>& b);

	AABB bbox;
};

icon-padding

MathUtil.cpp

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
bool BVHNode::compare(const std::shared_ptr<IHittable> &a, const std::shared_ptr<IHittable> &b, int axis) {
	return a->boundingBox().axis(axis).min < b->boundingBox().axis(axis).min;
}

bool BVHNode::compareX(const std::shared_ptr<IHittable> &a, const std::shared_ptr<IHittable> &b) {
	return compare(a, b, 0);
}

bool BVHNode::compareY(const std::shared_ptr<IHittable> &a, const std::shared_ptr<IHittable> &b) {
	return compare(a, b, 1);
}

bool BVHNode::compareZ(const std::shared_ptr<IHittable> &a, const std::shared_ptr<IHittable> &b) {
	return compare(a, b, 2);
}

And to divide the scene based on randomized axis, we need to define a randomInt function:

icon-padding

MathUtil.h

1
2
3
inline int randomInt(int min, int max) {
	return static_cast<int>(randomDouble(min, max + 1));
}

Then we can define the BVHNode class:

icon-padding

MathUtil.cpp

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
bool BVHNode::hit(const Ray &r, Interval interval, HitRecord &record) const {
	if (!bbox.hit(r, interval)) {
		return false;
	}
	bool hit_left = left->hit(r, interval, record);
	bool hit_right = right->hit(r, Interval(interval.min, hit_left ? record.t : interval.max), record);
	return hit_left || hit_right;
}

AABB BVHNode::boundingBox() const {
	return bbox;
}

bool BVHNode::compare(const std::shared_ptr<IHittable> &a, const std::shared_ptr<IHittable> &b, int axis) {
	return a->boundingBox().axis(axis).min < b->boundingBox().axis(axis).min;
}

bool BVHNode::compareX(const std::shared_ptr<IHittable> &a, const std::shared_ptr<IHittable> &b) {
	return compare(a, b, 0);
}

bool BVHNode::compareY(const std::shared_ptr<IHittable> &a, const std::shared_ptr<IHittable> &b) {
	return compare(a, b, 1);
}

bool BVHNode::compareZ(const std::shared_ptr<IHittable> &a, const std::shared_ptr<IHittable> &b) {
	return compare(a, b, 2);
}

BVHNode::BVHNode(const HittableList& list) : BVHNode(list.objects, 0, list.objects.size()) {

}

BVHNode::BVHNode(const std::vector<std::shared_ptr<IHittable>> &objects, size_t start, size_t end) {
	auto modifiable_objects = objects;
	auto axis = randomInt(0, 2);
	auto comparator = (axis == 0) ? compareX : (axis == 1) ? compareY : compareZ;
	auto object_span = end - start;

	if (object_span == 1) {
		left = right = objects[start];
	}
	else if (object_span == 2) {
		if (comparator(objects[start], objects[start + 1])) {
			left = objects[start];
			right = objects[start + 1];
		}
		else {
			left = objects[start + 1];
			right = objects[start];
		}
	}
	else {
		std::sort(modifiable_objects.begin() + start, modifiable_objects.begin() + end, comparator);
		auto mid = start + object_span / 2;
		left = std::make_shared<BVHNode>(modifiable_objects, start, mid);
		right = std::make_shared<BVHNode>(modifiable_objects, mid, end);
	}

	bbox = AABB(left->boundingBox(), right->boundingBox());
}

Now let’s apply that in our main function:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
#include "ImageUtil.h"
#include "GraphicObjects.h"
#include "Material.h"
#include "Camera.h"

int main() {

	auto camera = Camera(1920, 16.0 / 9.0, 30, {0, 0, 0}, {-13, 2, 3}, 0.6);
	auto world = HittableList();
	auto ground_material = std::make_shared<Metal>(Metal(Color(0.69, 0.72, 0.85), 0.4));
	auto left_ball_material = std::make_shared<Lambertian>(Lambertian(Color(0.357, 0.816, 0.98)));
	auto center_ball_material = std::make_shared<Metal>(Metal(Color(0.965, 0.671, 0.729), 0.4));
	auto right_ball_material = std::make_shared<Dielectric>(Dielectric(1.5, Color(0.8, 0.8, 0.8)));
	world.add(std::make_shared<Sphere>(Sphere(1000, {0, -1000, -1.0}, ground_material)));
	world.add(std::make_shared<Sphere>(Sphere(1, {0, 1, 0}, center_ball_material)));
	world.add(std::make_shared<Sphere>(Sphere(1, {4, 1, 0}, right_ball_material)));
	world.add(std::make_shared<Sphere>(Sphere(1, {-4, 1, 0}, left_ball_material)));
	camera.setSampleCount(100);
	camera.setShutterSpeed(1.0/24.0);
	int obj = 0;
	for(int i = -15; i < 15; i += 2) {
		for (int j = -15; j < 15; j += 2) {
			obj++;
			auto coord = Vec3((i + randomDouble(-1, 1)), 0.2,(j + randomDouble(-1, 1)));
			auto displacement = Vec3{0, randomDouble(0, 0), 0};
			auto material = static_cast<int>(3.0 * randomDouble());
			if ((coord - Vec3{0, 1, 0}).length() > 0.9) {
				Vec3 color;
				std::shared_ptr<IMaterial> sphere_mat;
				switch (material) {
					case 0:
						color = Vec3::random() * Vec3::random();
						sphere_mat = std::make_shared<Lambertian>(color);
						world.add(std::make_shared<Sphere>(0.2, coord, coord + displacement, sphere_mat));
						break;
					case 1:
						color = Vec3::random() * Vec3::random();
						sphere_mat = std::make_shared<Metal>(color, randomDouble(0.2, 0.5));
						world.add(std::make_shared<Sphere>(0.2, coord, coord + displacement, sphere_mat));
						break;
					case 2:
						color = Vec3::random(0.7, 1);
						sphere_mat = std::make_shared<Dielectric>(randomDouble(1, 2), color);
						world.add(std::make_shared<Sphere>(0.2, coord, coord + displacement, sphere_mat));
						break;
					default:
						break;
				}

			}
		}
	}
	world = HittableList(std::make_shared<BVHNode>(world));
	camera.setRenderDepth(50);
	camera.setRenderThreadCount(12);
	camera.setChunkDimension(64);
	camera.Render(world, "out", "test.ppm");
	spdlog::info("total object: {}", obj + 4);
	return 0;
}

Notice we are converting the HittableList to a BVHNode before rendering, which the latter uses a tree structure to store the objects and searches node in logarithmic time complexity rather than linear time complexity, but the result is not so different:

bvhResult

Recommended reading~
Note on Note on "Ray tracing in one week" - Part 3 Note on Note on "Ray Tracing The Next Week" - Part 5 Note on Note on "Ray tracing in one week" - Part 5
Recommended reading~
Note on Note on "Ray tracing in one week" - Part 3 Note on Note on "Ray Tracing The Next Week" - Part 5
Comments
On this page
Note on "Ray Tracing The Next Week" - Part 2
  1. From linear to sub-linear
  2. Bounding Volume Hierarchy
  3. Building a tree