Chapter 2 - Bounding Volume Hierarchy (BVH)

The construction of an effective and well-optimized Bounding Volume Hierarchy (BVH) plays a pivotal role in determining the efficiency of these calculations. Choosing the right BVH construction algorithm can greatly impact the speed and accuracy of solar shading mask calculations.

Numerous techniques have been developed to construct BVHs efficiently, and the choice of technique can significantly affect the structure and traversal behavior of the resulting BVH. Techniques vary from classic bottom-up approaches like the Surface Area Heuristic (SAH) to top-down techniques like Binary Space Partitioning (BSP).

Each technique has its strengths and weaknesses, and their impact becomes even more pronounced when applied to solar shading mask calculations. For instance, scenes with highly detailed geometry might benefit from BVHs constructed with the SAH algorithm (Cost Functions), as it efficiently reduces the number of unnecessary ray-object intersection tests. On the other hand, scenes with varying object densities could benefit from median splits and axis-cycling, ensuring that the splitting axes are evenly distributed and can help maintain overall balance during the BVH construction.

Furthermore, the trade-off between construction time and traversal efficiency cannot be overlooked. While some techniques may yield BVHs with faster traversal times, they might require longer construction times, but this would be acceptable for the computation of solar shading masks, since the BVH is constructed only once and then used for many ray-object intersection tests.

A visual definition of a BVH structure using Axis-Aligned Bounding-Boxes (from here)

In this context, N1 would be the bounding volume for the entire object or scene, N2 and N3 would be those respectively containing N4 and N5, and so on, until we reach the leaves of the tree, which are directly the bounding boxes of the primitives.

In the code snippet below, written by Luca Berti, a recursive top-down construction algorithm was proposed for further improvements. We can see that he adopts a recursive build function, which is called to create each node of the tree. Here is the proposed implementation:

This function is responsible for constructing the BVH tree from the primitives. It’s called recursively and each time it either creates a leaf node if there’s only one primitive left, or an internal node with two child nodes. The primitives are split by choosing a cutting dimension and sorting them by their centroids along this dimension, and then the data is divided into two equally sized parts, for each of which a new node is created.

The cutting dimension is cycled between 0, 1, 2 (representing the x, y, and z axes in a 3D space) by using (cut_dimension+1)%nDim in the recursive calls (Axis Cycling). This is the main "Divide and Conquer" idea behind this top-down construction algorithm.

It then sorts the primitives by their centroids along the cutting dimension, using the std::nth_element function, which partially sorts the primitives so that the element at the mid index will be in the place it would be in a fully sorted array, and all elements before it are less than or equal to the elements after it. The comparison function :

is used to sort the elements based on their centroids along the cutting dimension (Median Cut).

Finally, the data is divided into two equally sized parts when calculating the midpoint of the primitives' indexes.

Other splitting algorithms can be used, such as the Surface Area Heuristic (SAH) or the Middle Split Heuristic (MSH), which are listed and explained in the Spatial Splits and Clustering section.

Instead of starting with all scene primitives in one cluster and recursively splitting them, bottom-up construction algorithms start with each primitive in its own cluster and recursively merge the closest pairs. This is done either until the desired number of clusters is reached, or each cluster contains a maximum number of primitives. The clusters are then used as the primitives for the next level of the tree. This process is repeated until the root node is reached.

But such a method wouldn’t be the best for our purpose, first of all, looking at its inefficiency in finding pairs, we observe that every time we need to merge two clusters, we have to determine which two are the closest. This operation can be computationally expensive and can result in O(n^2) time complexity if not optimized. For a large number of primitives, this can slow down the construction process considerably. Most important, such techniques suffer from a strong lack in flexibility: the algorithm involves the creation and deletion of nodes dynamically, and the management of an array of active nodes. This can lead to fragmentation and additional memory overhead, which is unacceptable due to the need of scaling up the algorithm to be performed on GPUs, where memory is very limited. The bottom-up approach generally lacks the flexibility to adapt the BVH structure based on specific traversal use-cases. For instance, in solar shading masks computations, certain areas might need more detailed BVH structures (higher resolution) than others. Bottom-up methods don’t inherently allow for this granularity during construction, as opposed to the top-down approach, where the BVH structure can be adapted to our needs at each step of the construction process.

Below is an example of a bottom-up construction algorithm, found in [real-time-collision-detection] by Christer Ericson:

Introduced by Walter et al., bottom-up construction by agglomerative clustering proposes to start with all scene primitives considered as individual clusters and recursively merges the closest pairs (the distance function being for example the surface area of a bounding box enclosing both clusters). In general, these trees tend to have lower global costs, but the construction is more time-consuming.

These can be defined thanks to various algorithms presented on NVIDIA’s website.

The major inconvenience when using bottom-up algorithms is that the upper nodes are poorly locally optimized and thus the research for the closest neighbor can be very costly. To prevent this, Gu et al proposed to recursively perform spatial median splits based on Morton codes until each subtree contains less than a chosen number of clusters. The clusters are merged using agglomerative clustering. Using this at all levels in the BVH, even the top level nodes' split will be locally optimized.

Meister and Bittner proposed a GPU-based algorithm using k-means clustering (found in [real-time-collision-detection]): scene primitives are subdivided into k clusters using k-means clustering. When done recursively, a k-ary BVH is built, which can also be converted to a binary tree by constructing intermediate levels using agglomerative clustering.

Introduced by Meister and Bittner, the key observation is that the distance functions have a non-decreasing property, meaning that once we found two mutually corresponding nearest neighbors, we can immediately merge their clusters since no other closer one will be found. The clusters are kept sorted along the Morton Curve, finding the nearest cluster by searching both sides of the sorted cluster array, testing a predefined number of clusters. Since it does not rely on distance matrices, it is GPU-friendly, and only a small number of iterations are needed to build the whole tree.

The hierarchical nature of the BVH prevents a straightforward parallelization of the construction algorithm. But now, the BVH construction can be reduced to sorting scene primitives along the Morton curve (the order is given by Morton codes of fixed length, 32 or 64 bits), and using optimized sorting algorithms such as the radix sort, it can be done in 2n-1 time. The Morton code implicitly encodes a BVH constructor by spatial median splits.

One straightforward approach is to choose the axis with the longest extent of the bounding volume as the separating axis. This can help effectively divide the scene along its largest dimension, potentially leading to more balanced partitions.

Another method involves cycling through the three axes (X, Y, Z) and selecting the next axis in a cyclic manner for each spatial split. This approach ensures that the splitting axes are evenly distributed and can help maintain overall balance in the BVH construction. This is the approach proposed by Luca Berti, presented in the original code of this project, like seen during the call to the recursive build function:

The 2nd value representing the cutting dimension is cycled between 0, 1 and 2, representing the x, y and z axes of our 3 dimensional euclidean space, by using (cut_dimension+1)%nDim in the recursive calls. At each call, it is incremented by 1, enabling a different splitting axis to be used at each level of the tree. After choosing the splitting axis, the median value along that axis is computed and used as the splitting position, also know as a median cut, discussed right below.

The median cut strategy involves computing the median value along a specific axis and using it as the splitting position. This method aims to divide the scene into two halves containing an equal number of objects, which can help achieve good load balancing. This is implemented in the following line of the recursive build method, when calling the std::nth_element function:

Here is an example of a construction using the median cut strategy, but on a single splitting axis:

And such spacial divisions can lead to a similar tree as the following:

We will combine this method with the Axis Cycling one, to ensure that the splitting axes are evenly distributed and can help maintain overall balance during the BVH construction.

The quality of a particular BVH can be estimated in terms of the expected number of operations needed for finding the nearest intersection with a given ray. It can be estimated thanks to the recurrence equation according to Daniel Meister et al. [survey-BVHs-for-ray-tracing]:

Let’s look at a natural way of structuring the tree by mapping its nodes in a breadth-first level-by-level manner:

This way, we always know that a parent’s children can be found at positions \(2i+1\) and \(2i+2\) in the array, usually inducing wasted memory unless dealing with a complete tree.

When preordering them in traversal order, the left child will always follow its parent, and only one link is needed to point to the right child.

Flattening the tree in this way allows us to store the tree in a single array, with each node containing a pointer to its right child and a flag indicating whether it has a left child or not. This way, we can easily traverse the tree by following the right child pointers and using the left child flags to determine whether we should follow the left child or not and avoid the need for a stack and storage of 2 pointers per node (only one is necessary). This method is also cache-friendly since the nodes are stored in a linear array.

When using modern architecture, execution time is mostly limited by cache issues when fetching data from memory. One possible way of adopting a cache-friendlier solution would be by merging the sets of binary tree nodes into a 'tri-node' containing the parent and its children, preventing it from needing internal links. Below we can see an example representing a complete 4-level binary tree with 14 internal links with a 2-level tri-node tree storing only 4 internal links. Even better, this representation can also be combined with other optimizing structures seen before.

Flattening a tri-node tree is similar to flattening a binary tree, except that we need to store the parent’s index in the array as well as the left and right child flags. The right child pointer is replaced by a flag indicating whether the parent has a right child or not, the left and parent’s one are replaced in the same manner. The root node is a special case, since it has no parent, signified by a special flag. Three new structures (GPUNode, GPURay and GPUTree) were introduced, storing only critical information for it to be of small enough size to be copied-by-value to the GPU.

In our own GPU implementation, which had to be able to reproduce exactly the same BVHs as the original one presented by Luca Berti, we had no other choice than copying the informations containend in each BVH node to the GPU, and force it to retain the exact same structure. In order to do so, a suitable structure for the rays, the nodes and the global BVH had to be found, meaning excluding every unsupported data structure, such as std::vector or Eigen::MatrixXd. The whole BVH had to be flattened into a single static array, preventing different threads to have to access shared memory, which would have been very costly. In order to do so, the following structures where chosen (available in the Bvh_GPU.cuh file):

And implementing only the methods necessary to perform the basic traversal on the axis-aligned bounding boxes, and not directly on the primitives, which could cause memory overhead.

Looking at the sizes of the differents structures and classes, we can compute the average size of a BVH node (supposing a 64 bits architecture):

Total size for BVH Node: 8 + 8 + 8 + 4 + 4 + 4 + 3 x (4 + 8 + 3 x 8) = 144 bytes

Having to transfer \(144 \times N_{nodes}\) could cause memory overhead, contrary to the structure implemented for the GPU:

Which totals out at 72 bytes per node, which is half the size of the original structure, and thus much more efficient in terms of memory usage. This enables us also to be able to create pointers on the GPU to perform the stack-less traversal without the need of preordering the nodes in a specific order such as described in the Preordering Algorithm section.

When looking at the storage needed for the computation of ShadingMasks, we can pass the whole BVH structure and make a copy of it directly on the GPU’s shared memory. This way, we can avoid the need to transfer the BVH structure from the CPU to the GPU constantly, which can be a very expensive operation. Using such a method may cause problems depending on the size of the BVH structure, since the GPU’s shared memory is limited. However, we can use the its structure’s size as a parameter to determine whether or not we should use this method. If the BVH structure is too big, we can implement smaller structures to hold the BVH, preordering the nodes in flattened 1D array’s only containing useful information (and not all methods and attributes of the BVH structure). This way, we can reduce the size of the BVH structure and make it fit in the GPU’s shared memory. And since the traversal is performed nRays * nElements times (more than 5000 rays per element), we can compute the array’s once by indexing the nodes, its children and its parent.

Depending on the number of bounding boxes representing the acceleration structure, we can use the same stack-less algorithm as used on the CPU before, only adding atomic operations to count the number of intersections between a certain ray and leaf nodes. If a leaf is intersection, we append the firstPrimOffset to the results list, enabling the CPU to access the primitives and perform the intersection tests. This way, we avoid the need to transfer twice the amount of data.

This method can be coupled with different parallelization methods, since it serves a more generall purpose, namely the traversal of the BVH structure by a single ray. Leveraging the GPU’s parallelization capabilities, we can write a kernel that will be executed by all the available threads, each one performing the traversal of a single ray.

In this implementation, each thread is assigned a single ray and will perform the traversal of the BVH structure, storing the results in a local array. Once the traversal is done, the results are stored in the global results array, which will be accessed by the CPU to perform the intersection tests on the primitives.

Finally, the last parallelization occurs when the wrapper method calls the GPUraySearch method, which will be executed on all available GPUs in the cluster, each one maxing out their number of threads.

In order to take further this work, one could enlarge the last method to take into account different GPU models, each one having a specific amount of memory and threads available. This way, we could optimize the number of rays per device, and thus the number of threads per block, to maximize the number of rays processed per second.

Further optimizations can be made by assigning different parts of the shading mask matrix to different GPUs or threads, which would enable the processing of very large meshes, unable to fit into on single GPU’s memory. Problems at the boundaries (buildings from zone A casting shadows on zone B) could be solved by using a buffer zone, which would be computed by the CPU and then sent to the GPU to be processed, or simply by limiting the origin of the rays to the zone’s boundaries but not the rays' directions (which would specifically cause problems in this case).

When dealing with the computation of shading masks, view factors or radiative transport, we use static geometry to realistically represent the scene. Only few topological changes have to be taken into account, hence the decision of also optimizing the build for the BVH tree’s quality in order to reduce traversal operations. Even if the construction speed is important, we are not developing a real-time application, but rather trying to compute physically realistic results. We can build the BVH once and reuse it for multiple ray tracing operations without the need to update or rebuild the BVH. This approach can significantly improve performance, as constructing the BVH is a computationally expensive operation.

Even when taking into account the changes occurring due to the seasonality of the chosen districts and cities (french cities are subdued to changing weather conditions, leaves are falling and trees do not cast as big of a shadow in winter than in summer).

With its sole focus on accelerating ray tracing on GPUs, NVIDIA’s OptiX API offers a suite of capabilities that often outshines traditional hand-crafted traversal algorithms. We will discuss the inner workings of OptiX, its advantages, and how it can be harnessed to achieve better results.

Spatial data structures exploit the spatial locality of scene primitives. But this isn’t the only way of leveraging spatial locality. To further accelerate the whole process, we could map rays to interior nodes deeper in the tree during the traversal, skipping top-level nodes. A major caveat of such methods is that there is no guarantee that the found intersection corresponds to the closest one. But when computing shading masks, the lack of distance consideration is not a drawback. Instead, we solely focus on determining whether an object is present along the path of the ray.

Another way to optimize the ray generation would be to exploit the graphics card’s instancing of objects, enabling it to create multiple copies of one object in record time. Benthin and Wald decided that, instead of tracing the rays sequentially, they would generate bounding frusta of coherent rays simultaneously harnessing the potential of a SIMD unit (as many rays in one frustum as the SIMD unit is wide).

This could be taken further, by assigning parts of a matrix to a specific block in the GPU, leveraging the constant memory and launching the frustum of rays in the respective direction defined by the block-assigned resulting matrix. This way, the rays are processed in a more coherent manner, and the GPU’s constant memory is used to its full potential. Moreover, the frustum could be instantiated directly on the GPU, and the identical rays could be transformed and translated through random values, generated by the mersene twister algorithm that can be implemented on a CUDA kernel, and therefore be naturally processed in parallel. This would result in a more efficient memory transfer, since the rays shouldn’t be transferred back to the CPU, but only the resulting intersected leaves.

(image found here)