See Monte Carlo ray tracing technology series Monte Carlo ray tracing technology

This part is the most difficult and complex part of the ray tracker we are studying so far. I put it in Chapter 2 so that the code can run faster, and because it refactored some hitable s, when I added rectangles and boxes, we didn't have to go back and refactor them.

The intersection of light and object is the main time bottleneck of ray tracker, and time is linear with the number of objects. But it is the first mock exam for the same model, so we should be able to use the speed of two point search to search for logarithmic complexity. Because we send millions to billions of rays on the same model, we can do a simulation of the sorting model, and then each ray intersection can be a sublinear search. The two most common sorting methods are 1) space partition and 2) object partition. For most models, the latter is usually easier to write and run just as fast.

The key idea of bounding volume on a set of primitives is to find a body that completely surrounds (surrounds) all objects. For example, suppose you calculated a bounding sphere for 10 objects. Any ray that misses the boundary sphere will surely miss all ten objects. If the ray hits the boundary sphere, it may hit one of ten objects. So the form of boundary code is always:

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

The key point is that we divide the objects into subsets. We don't care about the screen or the volume. Any object is in only one boundary body, but the boundary bodies can overlap.

To make things sublinear, we need to make the boundary volume hierarchical. For example, if we divide a group of objects into two groups, red and blue, and use a rectangular boundary volume, we will get:

Note that the blue and red boundary volumes are contained in the purple boundary body, but they may overlap, and they have no order - they are all inside. So the tree shown on the right has no concept of sorting the left and right child elements; they are just inside. The code is:

If (hit purple) hit0 = hit blue closed object hit1 = hit red closed object If (hit0 or hit1) Return true and close up hit information Return false

To do all this, we need a way to get good segmentation instead of bad segmentation, and a way to intersect Ray with boundary volume. The intersection of light boundary bodies needs to be calculated quickly, and the boundary volume needs to be very compact. In most model practices, axis aligned boxes work better than other alternatives, but this design choice always needs to be remembered if unusual model types are encountered.

From now on, we'll call it axis aligned bounding rectangle parallelepiped (actually, that's what they need to be called, if accurate) axis aligned bounding box, or AABBs. Any way you want to use crossover with AABB rays is good. What we need to know is whether we hit it or not; we don't need to hit the midpoint, the normal, or anything we need to show the object.

Most people use the "slab" method. This is based on an observation that an n-dimensional AABB is just the intersection of n-axis alignment intervals, commonly known as "slab". The spacing is just a point between two endpoints, such as X, making 3 < x < 5, or more succinctly x in (3,5). In 2D, two overlapping intervals form a 2D AABB (rectangle):

To get the light to an interval, we first need to determine whether the light reaches the boundary. For example, also in 2D, this is the ray parameters t0 and t1. (if the light is parallel to the plane, these are undefined.)

In 3D, these boundaries are planes. The equations for the plane are x=x0 and x=x1. Where does the ray hit which plane? Recall that a ray can be seen as a function of a given t return position p(t):

p(t) = A + tB

This equation applies to all three x/y/z coordinates. For example, x(t)=Ax+t*Bx. X-ray strikes plane x=x0 at t meeting the following equation:

x0 = Ax + t0* Bx

So t at the impact point is:

t0 = (x0 - Ax) / Bx

We get a similar expression t1 = (x1 ax) / Bx for x1.

The key observation of turning one-dimensional mathematics into hit test is that for hit, t interval needs to overlap. For example, in 2D, green and blue overlaps occur only on hits:

What "overlapping of t intervals in the board? Hope in the code is similar to:

compute (tx0, tx1) compute (ty0, ty1) return overlap?( (tx0, tx1), (ty0, ty1))

This is very simple. In fact, the 3D version can also work. That's why people like slab method:

compute (tx0, tx1) compute (ty0, ty1) compute (tz0, tz1) return overlap?( (tx0, tx1), (ty0, ty1), (tz0, tz1))

There are a few caveats that make it less beautiful than it initially looked. First, assume that the ray moves in the negative x direction. The interval (tx0, tx1) calculated above may be the opposite, for example (7, 3). Second, the boundary here can give us infinite. If the origin of the light is on a flat boundary, we can get a NaN. In AABB of various ray trackers, there are many ways to deal with these problems (there are also some vectorization problems, such as SIMD, which we will not discuss here). If you want to learn more about speed vectorization, INGO Ward's paper is a good starting point.) As far as our purpose is concerned, as long as we do it fairly fast, it is unlikely to become a major bottleneck, so let's do the simplest, which is usually the fastest! First let's look at the calculation interval:

tx0 = (x0 - Ax) / Bx tx1 = (x1 - Ax) / Bx

The trouble is that a fully effective ray will have Bx=0, resulting in dividing by 0. Some of the rays are in the plate, while others are not. In addition, under IEEE floating-point, zero will have a sign. The good news for Bx=0 is that if it is not between x0 and x1, tx0 and tx1 are both + infty or - infty. So, using min and max should get the right answer:

tx0 = min((x0 - Ax) / Bx, (x1 - Ax) / Bx); tx1 = max((x0 - Ax) / Bx, (x1 - Ax) / Bx);

If we do this, the rest of the trouble is if Bx=0, x0 AX = 0, or X1 AX = 0, then we get a NaN. In this case, we may accept a hit or miss answer, but we will discuss it later.

Now, let's look at overlapping functions. Suppose we can assume that the interval is not reversed (so the first value is less than the second value in the interval), and we want to return true in this case. The Boolean overlaps of overlapping intervals (F, f) of calculation intervals (D, d) and (E, e) are as follows:

bool overlap(d, D, e, E, f, F) f = max(d, e) F = min(D, E) return (f < F)

If there is any NaN running there, the comparison will return false, so we need to make sure that our bounding box is a little filled if we are concerned about the edge wipe (we should probably be concerned, because in the raytracer everything eventually happens). When all three-dimensional spaces are in a cycle and are transferred in the interval tmin, tmax, we get:

#ifndef __AABB_H__ #define __AABB_H__ #include "vec3.h" #include "ray.h" inline float ffmin(float a, float b) { return a < b ? a : b; } inline float ffmax(float a, float b) { return a > b ? a : b; } class aabb { public: aabb(){} aabb(const vec3&a, const vec3&b) { _min = a;_max = b; } vec3 min()const { return _min; } vec3 max()const { return _max; } bool hit(const ray& r,float tmin,float tmax)const { for (int a = 0;a < 3;a++) { float t0 = ffmin((_min[a]-r.origin()[a])/r.direction()[a], (_max[a]-r.origin()[a])/r.direction()[a]); float t1 = ffmax((_min[a]-r.origin()[a])/r.direction()[a], (_max[a] - r.origin()[a]) / r.direction()[a]); tmin = ffmax(t0, tmin); tmax = ffmin(t1, tmax); if (tmax <= tmin)return false; } } vec3 _min; vec3 _max; }; #endif

Note that the built-in fmax() is replaced by ffmax(), which is much faster because it doesn't worry about NaNs and other exceptions. In reviewing this intersection method, Andrew Kensler of Pixar did some experiments and proposed this version of code, which is very effective on many compilers. I regard it as my starting method:

Now we need to add a function to calculate the bounding boxes of all hitable. Then we'll build a hierarchy of boxes on all primitives, and individual primitives (such as spheres) will be on the leaves. This function returns a Boolean value because not all elements have bounding boxes (such as infinite planes). In addition, the object moves, so the frame interval takes time 1 and time 2, and the bounding box binds the object that moves through that interval.

bool sphere::bounding_box(float t0, float t1, aabb& box)const { box = aabb(center - vec3(radius, radius, radius), center + vec3(radius, radius, radius)); return true; }

For the moving sphere, we can take the box of the sphere at t0 and the box of the sphere at t1, and then calculate the boxes of the two boxes:

aabb surrounding_box(aabb box0, aabb box1) { vec3 small(fmin(box0.rmin().x(),box1.rmin().x()), fmin(box0.rmin().y(),box1.rmin().y()), fmin(box0.rmin().z(), box1.rmin().z())); vec3 big(fmax(box0.rmax().x(), box1.rmax().x()), fmax(box0.rmax().y(), box1.rmax().y()), fmax(box0.rmax().z(), box1.rmax().z())); } bool moving_sphere::bounding_box(float t0, float t1, aabb&box)const { vec3 v_radius(radius, radius, radius); box = surrounding_box(aabb(center0-v_radius,center0+v_radius),aabb(center1 - v_radius, center1 + v_radius)); return true; }

BVH will also be a hitable - just like the hitable on the list. It's actually a container, but it can respond to "this ray hit you?" A design problem is whether we have two classes, one for the tree and one for the nodes in the tree, or whether we have only one class and the root is only the node we point to. If possible, I'm a fan of class design. We design such a class:

class bvh_node :public hitable { public: bvh_node() {} bvh_node(hitable **l, int n, float time0,float time1); virtual bool hit(const ray&r, float t_min, float t_max, hit_record&rec)const; virtual bool bounding_box(float t0, float t1, aabb&box)const; hitable *left; hitable *right; aabb box; }; bool bvh_node::bounding_box(float t0, float t1, aabb& b)const { b = box; return true; }

Note that the child pointer points to the generic hitables. They can be other bvh_notes, or spheres, or any other hitable.

The hit function is very simple: check whether the frame of the node is hit, and if so, check the child nodes and collate any details:

bool bvh_node::hit(const ray&r, float t_min, float t_max, hit_record &rec)const { if (box.hit(r, t_min, t_max)) { hit_record left_rec, right_rec; bool hit_left = left->hit(r, t_min, t_max, left_rec); bool hit_right = right->hit(r, t_min, t_max, right_rec); if (hit_left && hit_right) { if (left_rec.t < right_rec.t) rec = left_rec; else rec = right_rec; return true; } else if (hit_left) { rec = left_rec; return true; } else if (hit_right) { rec = right_rec; return true; } else return false; } else return false; }

The most complex part of any efficiency structure, including bvh, is building it. We do this in the constructor. One of the cool things about bvh is that as long as bvh_ The list of objects in the node is divided into two sub lists, and the hit function can work. If the segmentation is done well, the bounding box of the two sub objects is smaller than that of the parent object, and the effect is the best, but this is for speed rather than correctness. I'll select the middle position and split the list along an axis at each node. I want to be simple:

1) Random axis selection

2) Sorting elements using library qsort

3) Half for each subtree

I use the old C sort qsort instead of C + + sort because I need a different comparison operator depending on the axis. q sort takes a comparison function instead of the less than operator. I passed in a pointer to the pointer - this is just the "pointer array" of C, because the pointer in C can also just point to the first element of the array. When the incoming list is two elements, I put an element in each subtree and end the recursion.

The traversal algorithm should be smooth, without checking the null pointer, so if I have only one element, I will copy it in each subtree. It might be helpful to explicitly examine three elements and follow a recursion, but I think the whole method will be optimized later. This will result in:

bvh_node::bvh_node(hitable **l, int n, float time0, float time1) { int axis = int(3 * myRandom()); if (axis == 0) qsort(l, n, sizeof(hitable*), box_x_compare); else if (axis == 1) qsort(l, n, sizeof(hitable *), box_y_compare); else qsort(l, n, sizeof(hitable *), box_z_compare); if (n == 1) { left = right = l[0]; } else if (n == 2) { left = l[0]; right = l[1]; } else { left = new bvh_node(l,n/2,time0,time1); right = new bvh_node(l+n/2,n-n/2,time0,time1); } aabb box_left, box_right; if (!left->bounding_box(time0, time1, box_left) || !right->bounding_box(time0, time1, box_right)); box = surrounding_box(box_left, box_right); }

Check if there is a bounding box in case you send something similar to an infinite plane without a bounding box. We don't have these elements, so we don't need to check them before we add them. The compare function must accept the null pointer you converted. This is the old c function, remind me why C + + was invented. I really messed it up, so I got all the pointers right. If you like this part, you will become a person who can implement the system in the future!

int box_x_compare(const void *a, const void *b) { aabb box_left, box_right; hitable *ah = *(hitable**)a; hitable *bh = *(hitable**)b; if (!ah->bounding_box(0, 0, box_left) || !bh->bounding_box(0, 0, box_right)) return 0;//Output error message if (box_left.rmin().x() - box_right.rmin().x() < 0.0)return -1; else return 1; } int box_y_compare(const void *a, const void *b) { aabb box_left, box_right; hitable *ah = *(hitable**)a; hitable *bh = *(hitable**)b; //if (!ah->bounding_box(0, 0, box_left) || !bh->bounding_box(0, 0, box_right)) // return 0; / / output error information if (box_left.rmin().y() - box_right.rmin().y() < 0.0)return -1; else return 1; } int box_z_compare(const void *a, const void *b) { aabb box_left, box_right; hitable *ah = *(hitable**)a; hitable *bh = *(hitable**)b; //if (!ah->bounding_box(0, 0, box_left) || !bh->bounding_box(0, 0, box_right)) // return 0; / / output error information if (box_left.rmin().z() - box_right.rmin().z() < 0.0)return -1; else return 1; }

After that, let's test it in the next section.