How to build a BVH – part 9b: Massive

In this series we explore the wonderful world of the bounding volume hierarchy: the mysterious data structure that enables real-time ray tracing on recent GPUs, in complex, animated scenes. Despite the broad adoption of ray tracing in games, in-depth understanding of the underlying technology seems to be reserved to a select few. In this article we show how a proper BVH is implemented without support from DXR / RTX / OptiX / Vulkan, in easy-to-read ‘sane C++’. Whether this is actually useful for you, or just interesting as background information, depends on your needs, obviously.

This series consist of ten articles. For each article, full source code is available on Github, in easy-to-read ‘Sane C++’. Contents:

1. Basics: data structure, simple construction and ray traversal.
2. Faster rays: the Surface Area Heuristic and ordered traversal.
3. Faster construction: binned BVH building.
4. BVHs for animated models: refitting and rebuilding.
5. BVHs for animated scenes: the top-level BVH.
6. TLAS & BLAS part 2 – all together now.
7. Consolidating, part 1: ray tracing textured and shaded meshes.
8. Consolidating, part 2: Whitted-style ray tracing.
9. OpenCL: GPU ray tracing, part 1.
10. OpenCL: GPU ray tracing, part 2 (SERIES FINAL, this article).
11. P.S.: Concurrency (under construction)

This is the end. Last time we ported the ray tracing code to the GPU, or at least: BLAS traversal and basic visualization. All that is left to do is some cleaning up, and TLAS traversal. With that in place, we will have the same power on CPU and GPU, in terms of animation: we will be able to specify a transform (rotate, translate, scale) for each instance, on top of the capability to rebuild or refit meshes on-the-fly.

In this final article we implement TLAS traversal and the recursive Whitted-style ray tracer on the GPU. For that a few details need to be sorted out. The final product deserves a fitting final demo, for which I will take some inspiration from the work of two of my students at Utrecht University.

Preparations

To keep the article 9a code intact we start with a new file, titled raytracer.cl, which is initially a copy of kernels.cl. The new file will be loaded from ::Init() on the CPU side:

The OpenCL source file currently includes template/common.h. We add cl/tools.h, which currently just contains some functions for generating random numbers.

We can now move most of the code from raytracer.cl to tools.h. I propose to move the various structs (Intersection, Ray, Tri, TriEx, BVHNode, TLASNode and BVHInstance), as well as the color format conversion function RGB32FtoRGB8 and the intersection functions IntersectTri and IntersectAABB. And finally, also BVHIntersect, for BLAS traversal. All that remains in raytracer.cl is then:

• float3 Trace(...), which will replicate the functionality of the Trace function we used on the CPU, and
• __kernel void render(...), the OpenCL kernel, which we call from the C++ code.

The reorganized code lets us deal with high-level things in raytracer.cl, while the low-level details are hidden away in tools.h.

To match the theme of the article, a new mesh has been added to the assets folder: dragon.obj.

This version of the Stanford Dragon consists of 19332 triangles, and it replaces the 1024 triangle Utah teapot. Some changes to the Mesh class are needed to accommodate the larger number of polygons:

Similar changes are needed in the implementation of the Mesh constructor.

Here, 11050 is approximately the number of vertices in the dragon mesh.

With these changes, the code still works, but there is a small issue:

The white shape that obscures the sky is a massive dragon. In the CPU code we would resolve this by assigning a scaling matrix to the instance of the dragon mesh. Or rather: we would apply the inverse of this matrix to the rays that intersects the dragon.

Time to add TLAS traversal.

GPU TLAS

The concept of the TLAS (a.k.a. top level BVH) was introduced in article 5 and article 6. The idea is to build and traverse a BVH over BVHs, where each of those BVHs can have a 4×4 matrix transform, for translation, rotation, and scaling. Using the same BVH multiple times but with a different matrix enables instancing.

We already have the necessary data on the GPU. There is the TLASNode struct:

And an array of BVHInstance objects, indexed by the BLAS field of the TLASNode:

The strange dummy array at the end is there to ensure that the size and layout of the OpenCL struct matches the CPU struct, at least to the extend that fields transform and invTransform are in the correct location. More on that in a second, because I made some serious mistakes here that need fixing…

We can now port the TLAS::Intersect(...) method. The process is similar to porting the BVH::Intersect method, so I paste the result here without going over every step in detail.

We call two functions from TLASIntersect that we have not defined yet.

The first one is particularly interesting: We call IntersectAABB, with a TLASNode* rather than a BVHNode*. To my surprise, this just works. The TLASNode struct has, for the purpose of AABB intersection, the same layout as a BVHNode. And OpenCL doesn’t complain about the pointer type mismatch… Great, for now, but definitely something to keep in mind: Next time this ‘flexibility’ may not be what we expect.

The second function requires more work. Here is the code for InstanceIntersect, based on BVHInstance::Intersect from the CPU code:

This function in turn calls TransformRay, which is defined as follows:

This implements the TransformVector and TransformPosition matrix functions we used in the CPU code: In OpenCL, we have to implement these ourselves.

Note how the float16 data is addressed here: To make a float3 out of cells 0, 1 and 2, we write invTransform->s012. A same technique is available for all OpenCL vector types. For instance, to access x of ray->D we can simply write ray->D.x, but we can also get a float2 based on the fields: e.g., ray->D.xy, but also: ray->D.xx, which would duplicate the first field, or ray->D.zyx, which would flip the three available fields and yields a float3. Any combination of x, y and z works. The technique is called swizzling. For a float16, x, y and z are obviously insufficient. In this case the fields are labeled 0..9 and A..F (or a..f, if you prefer that): this is the hexadecimal index of the 16 matrix cells.

Debugging

There is a problem with all this new code: it doesn’t work. When that happens, we debug, of course. But in OpenCL that is easier said than done: Breakpoints do not exist, and would be somewhat pointless anyway, considering that the code runs on (potentially) thousands of GPU threads. We do have some tools however. We can for example use printf from a kernel:

Besides that, we can visualize assumptions. For example:

In the latter case, Trace will return a red pixel if we passed it an non-existent skyPixels or triData buffer. Other colors can be used to return information on different problems. Since the visualization is per-pixel, this can yield more detailed information than the printf.

With these limited tools I found several issues. The first is in struct TLASNode, which is very similar to BVHNode, and should thus be ‘just fine’:

If this has you scratching your head: look no further, because the problem is not here, it is in the CPU version of this struct.

Now, one would assume that the two structs have the same size. Sadly, they don’t. The float3 type takes 16 bytes in OpenCL, and therefore, the CPU-side version mimics this, to ensure data compatibility. We thus have a CPU TLASNode that is substantially larger than 32 bytes.

The problem is quickly resolved by changing the CPU-side struct:

That is very ugly, but: Now the two 16-byte structs overlap the float3 fields in memory, which causes fields leftRight and BLAS to be exactly in the right location.

After this fix, we still get no useful results out of TLASIntersect. It turns out that there is another problematic struct; this time it is BVHInstance. The struct contains two float16 fields, which have the same size as a mat4 in the CPU code. The problem is: OpenCL dictates that any field must be positioned in a struct such that its offset is a multiple of its size. In this case, the size of a float16 is 64 bytes; we have only 12 bytes before the first float16 and therefore OpenCL will inject 52 unused bytes between field idx and the first transform…

Fixing this is again easier on the CPU side. We modify the BVHInstance struct to:

In the new layout, transform and invTransform are properly aligned: transform now has an offset of 0, and invTransform is at offset 64. The array of dummy integers at the end ensures that the size of the whole struct is now a multiple of 64 bytes: This way, in an array of BVHInstances object other than the first also have proper offsets for the float16 fields.

With these fixes, we finally have the full TLAS & BLAS functionality on the GPU, albeit without any fancy shading:

Restoring Light

The hard part is now done. Adding some shading is a matter of porting the body of the Trace function. This should not have too many surprises:

Note that the size of the brick texture is now hardcoded. We’ll want to change that at a later stage, to make the renderer more versatile. The light source is also still very much hardcoded, as it was in the CPU version. The three constants for the light source can be defined at global scope, somewhere above the Trace function:

To convert pixels read from the texture to float3 we also need a port of the RGB8toRGB32F function. It’s straightforward, and should go to tools.cl:

And finally, transforming the normal vector using the matrix for the intersected instance now uses a new function TransformVector, which in hindsight would come in handy in the TransformRay function as well. It is defined as:

And with that we have properly textured and illuminated dragons.

Mirror Mirror

The next challenge for the OpenCL port of the recursive ray tracer is: Recursion, which we need for mirrors, as well as glass and water (dielectrics). There is however a problem with recursion in OpenCL: It is forbidden. Yet, in the original ::Trace function on the CPU, it was rather instrumental to our needs:

The highlighted line (line 12) calls the Trace function itself with a just constructed secondary ray, which starts at the intersection point, and extends into the direction calculated using the specular reflection formula on line 8.

Luckily, we have ways to bypass recursion. We already did this to traverse the BVH: Instead of calling BVHIntersect for the left and right child nodes, we switch processing to the nearest child node, while we push the far node on a stack. And, whenever we would normally return, we now pop work from the stack, until the stack is empty.

We can deploy the same scheme here. Except, it will be simpler: Until we add glass or water, the path of a ray does not split up. It ends whenever we hit the skydome, but also when we hit a diffuse object. Only when we hit a reflective dragon we bounce, but in that case we can just restart the Trace function after replacing the original ray by a reflected ray. Like this:

Basically, this loop makes a ray bounce around, until it does not hit a mirror. Of course we need to be careful with infinite ‘recursion’ (in this case: an infinite loop). So, we keep track of the ‘bounce depth’ and once it exceeds a certain value, we terminate the loop. In that case the last surface we hit must have been a mirror; we thus return white, i.e. (float3)(1, 1, 1), which is probably better than any other color.

The scheme is quite effective, and brings back the shiny dragons:

MASSIVE

One advantage of ray tracing is how it handles massive scenes. For simple scenes, traditional rasterization is typically a lot faster, but at some point, the tables turn. This has everything to do with the BVH: once build (that’s a rather important condition!), a BVH over N triangles requires on average log_2 N steps to get from the root of the tree to a leaf node. For 1,000 triangles we get to a single leaf in 10 steps; for 1,000,000 we get there in 20 steps: Twice the cost, for a thousand times the complexity. Even a billion triangles is doable; this only take three times more work than one thousand.

To put that to the test we will create a scene, in which every vertex of the dragon model is replaced by a dragon. The dragon has 11,042 vertices, and 19,332 polygons; the scene thus will have a bit more than 213 million triangles. You’ll have to believe me when I say that adding another zero to that number will hardly matter.

The scene is constructed in the ::Init() function:

To make that work, the BVHInstance array in the application class definition in massive.h is changed to a pointer:

The ::Init() function needs access to vertex positions and normals. These were stored in local arrays (P and N) in the Mesh constructor. We add pointers to this data to the Mesh class, and make a small change to the Mesh constructor:

Tiny but important detail: Indices in obj files start at 1; we thus either make the arrays 1 element larger, or subtract 1 from every index that we use – hence the 11043 figure.

We are almost ready to show the scene: The scene is constructed, the (single) BLAS for the 11,042 dragon instances is built, and we have a TLAS – which constructs in about half a second on my machine. Half a second: That means the TLAS will have to be static, for now. We thus build it and copy it once to the GPU, in the ::Init() function. The tlasData->CopyToDevice() and instData->CopyToDevice() commands can now also be moved to ::Init(); we will render a static scene.

If you typed along, you now have a ‘dragon of dragons’, and it runs at a respectable speed.

Camera

Such a scene deserves a proper view.

For that, we will move the camera along a smooth path. The Catmull-Rom spline is ideal for this: it is smooth, and moves through a series of points.

Given P_0 .. P_3, we can find a point on the segment P_1 .. P_2 with the following code:

By varying t from 0..1, this function allows us to smoothly travel the curve. Once we reach t=1 (or t=0, if we move backwards), we simply select the next set of four points on the curve; this way, the smooth path can consist of any number of points.

Let’s take this one step further: If we move not just the camera over a smooth curve, but also a point at which the camera is aimed, we get a pretty versatile camera system. Consider the following set of points:

Each pair in this array contains one camera position (on the left) and a camera target (on the right). I obtained these points by manually navigating the scene, and recording 20 points to a file. The code for that is left in the program; feel free to replace the spline with a new path.

The camera can now be updated in the ::Tick() function, with a little help from the new SplineCam function:

And with that, we have reached our final destination.

Bandwidth

There’s one thing that is rather silly about the current demo, and it has to do with the template used. It is however a generic thing, so let me explain.

To render an image, the CPU code creates a buffer in VRAM. The OpenCL kernel fills this buffer with pixels. Next, we copy the buffer back to the CPU. And then what? Well, under hood the template maintains a single texture, the size of which matches the window size. The rendered pixels are copied to this texture once per frame, after ::Tick() completes. OpenGL is then used to draw a screen filling quad using this texture. This is how a frame ends up on the monitor. We thus have:

• Pixels rendered on the GPU using OpenCL
• Pixels transferred to CPU RAM
• Pixels copied from CPU RAM to a texture in VRAM

That feels like a rather long route with unnecessary detours, and it is. Luckily, there is a shortcut, and the name of the shortcut is GL/CL interop. This is a technique where OpenGL makes a resource (in this case: a texture) available to OpenCL, so OpenCL can render directly to it. And this saves time, obviously.

To enable the shortcut, we replace one line in ::Init(). The line that currently reads

is to be replaced by

This is inevitably highly template-specific. The first line creates the buffer using an OpenGL texture identifier, obtained from the render target of the template. The second line is used to tell the template to stop sending frames to the OpenGL texture – we are using OpenCL for that now.

A few more changes are needed. In ::Tick(), we delete the two lines that ‘obtain the rendered result’. And, in the kernel code, the first argument is now different:

And finally, plotting to the target changes. We used a unsigned int array before, but now we have a floating point OpenGL texture to work with. We use the following line to write to it:

After this change the demo works again, this time without wasting time sending pixel data back and forth.

Future Work

Where do we go from here? Well, as a wise man once said, you can be “Ray Tracing – the Rest of Your Life”…

Faster BVH building – The BVH for the dragon is constructed in just a few milliseconds, on a single core. For multiple meshes multiple cores can be used – until they run out. Take it further by optimizing the BVH builder: Some SIMD magic should bring down the time significantly, and beyond that, well, there’s ways to make it very fast. An important trick is to build your BVH while the GPU is tracing rays: the only thing that needs to happen between frames is sending the updated data. And even that isn’t strictly necessary.

Faster traversal – We trace hundreds of millions of rays per second now on the GPU – without RTX. But, the code is far from optimal… There are BVH layouts that are better suited for the GPU, and there’s 4-wide and 8-wide BVHs, which further improve performance. On the CPU packet traversal makes a huge difference, especially for primary rays. By the way, we never implemented proper shadow rays, so that should be on the top of your list.

Better rendering – A few reflecting dragons is nice, but what about more interesting materials, normal maps, depth of field, motion blur, soft shadows, indirect light? This is a field in which you can explore endlessly. Start with basic path tracing, and quickly your code will be too slow again, forcing you to look into filtering techniques and ReSTIR for efficiently sampling many lights. A life is not going to be enough, I’m afraid.

Engine flexibility – The code of this tutorial is short and hopefully clear, but definitely not robust. Handling multiple meshes and textures and being able to delete them and add them at runtime will require a lot of administration – especially if you still want to render on the GPU. And while you’re at it: A fixed-size skydome is not how things should be, obviously. Do mind compile times, once you start cleaning things up: you may have noticed that the tutorial code compiles in a flash, and that’s not just a lucky accident.

OptiX-CL – NVIDIA provides RTX support from CUDA in the form of a library called OptiX. It would be nice if the functionality of this library became available outside of CUDA, and thus for all devices of all vendors. If you see a challenge here: Let me know! Perhaps we can join forces. Especially the now obsolete OptiX Prime API would be very interesting (and relatively simple) to replicate.

Lightouse 2 – Speaking of OptiX: My Lighthouse 2 project uses OptiX to quickly render complex scenes on the GPU, but it is NVIDIA only, essentially. A proper render core using OpenCL or even using the CPU would be great to have. One day, perhaps.

That’s just my ideas though. Let me know what your plans are, but whatever you do, keep tracing.

And with that, allow me to sign off for once as if it’s 1999 on Flipcode:

Love you all,

Jacco Bikker, a.k.a. “The Phantom”.

Closing Remarks

The complete code for this project, in ten convenient incremental stages, can be found on GitHub (for Windows / Visual Studio).

Enjoyed the article? Spotted an error?

Ping me on Twitter (@j_bikker), comment on this article, or send an email to bikker.j@gmail.com.
You can also find me on Github: https://github.com/jbikker.

Or read some more articles: https://jacco.ompf2.com. Some favorites:

• Wavefront path tracing
• Probability Theory for Physically Based Rendering – part 1
• Probability Theory for Physically Based Rendering – part 2
• OPT#3: SIMD (part 1 of 2) – also check the other optimization blog posts!