The breakthrough, or "Jesus Christ, will this backstory ever end?"

Breakthrough! Finally managed to come up with a merging algorithm that performs well enough to be practical! But first, some more back story.

I've been trying different approaches to generating and merging these visibility curves. First of all, I've transitioned to using density curves instead, which are easier to work with. Density changes more linearly than visibility, as it simply increases along the view rays out from the camera when they hit smoke volumes. In addition, in general it looks much better when interpolating density instead of visibility. At the end, the visibility of a certain density threshold can be easily calculating using exp(-density). My visualizations still work with visibility as the exponential values are more intuitive to look at and debug though.

As the algorithm I had developed turned out to be such a bad fit for GPU hardware, I shifted my focus to finding an algorithm that could accomplish the same thing while taking advantage of GPU hardware. One of the most powerful features of the GPU is the ability to sample textures with hardware accelerated filtering, allowing you to get expensive filtering of textures essentially for free.

This lead me to consider some kind of binary search based algorithm. By integrating through a smoke volume, it is trivial to generate a 3D texture, where each pixel on the screen gets a 1D array of density values at different depths. Such a 3D texture is essentially what the Frostbite Engine presentation used directly, and is easy to generate by simply integrating through a smoke volume. The end result is that we have a texture storing the total density along the view ray as a function of depth for each pixel. Now we need to figure out where this density function reaches the density thresholds that we want to store the depths of.

By using hardware linear filtering to interpolate between these density values, we essentially turn this discrete function into a continuous one. As total density along the view ray can only ever increase, we can simply essentially do a binary search on this continuous function to narrow down the depth range in which we know the density reaches the threshold value.

1. Initialize minDepth and maxDepth to the full depth range of the smoke volume.
2. Loop N times:
    1. Calculate centerDepth = (minDepth + maxDepth)/2.
    2. Sample the density at that depth into centerDensity.
    3. If centerDensity is less than targetDensity, then maxDepth = centerDepth.
    4. If centerDensity is greater than targetDensity, then minDepth = centerDepth.

With a sufficient number of iterations, this algorithm will converge to a very narrow depth range that can be used as an approximation for when the function reaches certain density values. Each iteration will halve the size of the depth range, so the algorithm converges quickly even for very high depth complexity. By simply running one of these binary searches with a fixed number of iterations for each density threshold we're interested in, we can quickly generate a usable density function.

This does not solve the problem of merging an existing density curve with a new smoke volume though. However, this algorithm can be expanded to handle that too! In essence, what we want to do is figure out when the sum of the two curves reaches the density thresholds we care about. However, the existing curve stores depth as a function of density, while the integration results store density as a function of depth. It turns out that we can use this to our advantage!

1. Initialize minDepth and maxDepth to the depth range of the smoke volume AND the density curve.
2. Loop N times:
    1. Calculate centerDepth = (minDepth + maxDepth)/2.
    2. Sample from the integration curve the density at that depth into centerDensity.
    3. If centerDenstity is greater than targetDensity, then minDepth = centerDepth.
    4. else
        1. Calculate missingDensity = targetDensity - centerDensity.
        2. Sample the depth at which the existing curve reaches missingDensity into existingDepth.
        3. centerDepth and existingDepth now form a depth range in which the two functions will exactly sum up to targetDensity somewhere.
        4. Update minDepth and maxDepth based on this new depth range.

This algorithm only requires two texture samples and a tiny amount of math per iteration and can be executed repeatedly. In addition, the algorithm can actually converge FASTER than the original one, as both the minDepth and the maxDepth can be narrowed down at the same time in one iteration. If centerDepth accidentally lands right on top of the of the point where the threshold is, the depth sampled from the existing function will be exactly the same depth, meaning that the depth range will be updated to minDepth=maxDepth=centerDepth.

Although the idea was good, in practice it turned out to be pretty slow anyway. A lot of iterations, something like 15-20, were required to get decent results as the depth values needed to be very precise. With 32 thresholds and 20 iterations per threshold each requiring two texture fetches, I ended up with over 1000 texture fetches per pixel, which was simply too slow. I tried to greatly narrow down the depth range that needed to be searched for each threshold, but even if I managed to massively reduce the number of average iterations it was still too slow. Regardless, it was still a solid 3x faster than my initial approach.

I guess at this point we've essentially caught up with where I currently am today, which leads me back to... Breakthrough! I managed to write a very clever merging algorithm that required a minimal number of dynamic loops and texture fetches to very efficiently merge an existing density curve with the integrated density of a new smoke volume.

 - My original implementation loaded in both curves into arrays, which the GPU could not access with dynamic indices. Hence, merging the two arrays into a single array became extremely slow as very inefficient O(n^2) algorithms were needed to work around this limitation. The new implementation uses texture fetches to do dynamic reads as there is no limitation on the texture coordinates used to fetch data from a texture which allows me to use an unrollable loop, and also avoids needing to hold both curves in memory at all times. only the current and previous values of each of the two curves (4 values in total) need to be stored in memory now.

 - After merging the two curves into a single list, I need to loop over it to figure out where the curve reaches certain thresholds. This was originally done with a linear search through the entire merged curve per threshold, an O(n^2) operation. This could technically be done during the merge operation, again allowing me to discard old results and not having to hold the entire merged array in memory. However, this would require dynamic array indexing/writing which is not supported for normal arrays on GPUs... buuuut you can actually do dynamic indexing with shared memory! Hence, I allocate enough shared memory to store the output curve for each shader invocation and can safely write to the shared memory. There's still a dynamic loop needed for each iteration, but it's still massively faster.

With these optimizations, the algorithm went from requiring several O(n^2) operations to being a single unrolled loop, which is obviously O(n). This made the algorithm around one magnitude faster than the original algorithm, and compared to the binary search based algorithm it was still 3-4x faster!

The implementation, or "I hate GPUs."

The next step in the process was to

1. Come up with an algorithm that could calculate the depth values for when a smoke volume reaches certain visibility levels to calculate a visibility curve.

2. Come up with an algorithm that could take a pre-existing visibility curve and "add" another smoke volume to it, outputting a new visibility curve. The purpose of this algorithm is to allow me to compress multiple smoke clouds into a single visibility curve, so that lighting only has to be done once for the final curve.

3. Come up with a good GPU implementation of these algorithms.



I started out by writing a small visualization program to help me develop and test the initial algorithms. My first approach to the problem was to first integrate/raytrace through the smoke volume and generate a visibility curve, then merge that curve with the existing visibility curve for that pixel. This can then be done repeatedly to accumulate all smoke volumes into a single visibility curve.



Generating a visibility curve from a smoke volume

The solution to algorithm 1 that I came up with was:

1. Integrate through the smoke volume.
2. Whenever the visibility reaches the next threshold/milestone, calculate the exact depth and store that.

Simple enough. Here's an example of a generated visibility curve from a simple smoke volume that starts halfway into the scene.

Here's a second curve from a constant-density fog. For this curve, one can simply calculate the exact point the visibility drops to a certain value mathematically instead of integrating.

Merging two visibility curves into a single visibility curve

For algorithm 2, the algorithm for merging two visibility functions, I went with an exact mathematical solution. Given the input curves A and B:

1. For each point in curve A, find the corresponding visibility at that depth in curve B and multiply it together with its original visibility to calculate a new visibility for that point, storing the output in curve C. In the resulting curve, each point will NOT have a set visibility, so both depth and visibility will vary for these points.

2. Do the same for curve B, finding the corresponding visibility of each point in A and multiplying it together, storing the result in curve D.

3. Merge curve C and D into a single curve, curve E, sorting the nodes by visibility into a single curve.

4. Loop through all nodes in E to generate a new fixed-visibility curve. Whenever the visibility drops below the next threshold in visibility, store that depth into the next node.



The following is an example of a merge operation. Given the above two curves with blue and red nodes, the result of multiplying each curve by the point on the other curve and merging the two into one curve results in the following:

This curve has twice as many nodes as it should, and the nodes do not correspond to the fixed visibility values that we need. Hence, we continue with step 4 in the algorithm, compressing the curve to the correct number of nodes, where each node has a predetermined visibility, just like the original input curves.

This merging algorithm can then be performed again to add a third volume to the visibility curve.

The GPU implementation

The first GPU implementation revealed a number of problems that I had overlooked. The biggest one was how to deal with the case where the visibility never reaches all predefined values. For example, if the first visibility threshold is 0.95 and the visibility only ever drops to 0.96, the first threshold would never be reached, meaning that the smoke would simply never show up. Essentially, it meant that the visibility along a given curve was always rounded up to the closest threshold. This caused a kind of banding as the fog reached certain visibility levels.

The solution here was to adjust the visibility thresholds so that they went from 1.0 down to lowest visibility reached for each individual pixel. This redistributed the visibility levels dynamically based on the demands of each pixel so that all visibility levels were always reached, which lead to a huge increase in quality. When merging two visibility curves, the lowest visibility of the result is simply the lowest visibility of each of the curves multiplied together. This completely eliminated this issue.

There were however severe problems with the performance of the implementations. Generating a visibility curve during integration of a smoke volume requires dynamic writes and dynamic nested loops. This was a bad fit for GPUs.

The merge operation was even more inefficient. Step 1 and 2 each turned into O(n^2) operations, as each point in curve A needs a linear search through curve B and vice versa. In addition, the simple merging of two sorted lists turned into yet another nightmare, as that requires dynamic reads from the two arrays. Finally, the fourth step suffered from the same problem as the integration did, requiring dynamic writes and dynamic nested loops. 

The future was looking grim...

The algorithm, or "I just swapped the X and Y axes and called it a day."

With the background and the current state-of-the-art algorithm established in the previous post, it's time to take a look at my ultra sleek revolutionary new algorithm. Spoiler alert: it's pretty simple.

First of all, let's take a look at a simple scenario from the perspective of a single pixel in the scene. For each pixel, we're interested in being able to store a representation of the smoke's density along the view vector of that pixel. In other words, the entire problem boils down to creating a nice way of storing the mathematical function for density of the smoke over distance from the camera for each pixel.

In our little test scenario, we have a single, thick smoke cloud at a fairly high distance from the camera. The particular pixel we're looking at therefore has a zero density for a long distance, followed by a rapid increase in density. Instead of plotting the density directly, I'll be plotting the visibility of the smoke at a certain distance. Essentially, the visibility at a certain depth from the camera depends on how much smoke was between that point and the camera; the more smoke, the lower visibility. This is a bit more intuitive than plotting density, and will help with the visualization I want to make.

Here's the visibility curve. The x-axis plots depth, while the y-axis plots visibility.

As you can see, the visibility starts at 1.0 (100%) and remains at that level for half the view range. Then comes a rapid drop caused by the mentioned smoke cloud, causing the visibility to quickly approach zero, but never actually reach there (the visibility is equal to exp(-distance*density), meaning the value never reaches zero no matter how much smoke we have).

The approach that the algorithm mentioned in the previous blog post takes is that it simply samples this function at certain fixed depths, computes the visibility there and essentially interpolates the rest.

The red line is a pretty bad approximation of the original green curve. The sample points are simply spaced too far away. Even worse, this representation would be even worse if the depth complexity of the scene increased. In this case, the samples would need to be distributed over a bigger depth range, reducing the accuracy of the representation even further. Of course, this example only uses 8 samples, but it's clear that the sample count would need to be increased vastly to be able to accurately represent this function. Although the samples does not have to be evenly spaced, modifying the distribution of the samples can merely help in specific cases, and will be a trade-off regardless (for example having more samples closer to the camera and fewer farther away).

Something that I mentioned in the previous blog post was also the inefficiencies in lighting this representation. The lighting would be computed by looping over different points and computing lighting at those points. However, a vast majority of these points either have no smoke (meaning they won't reflect light) or have zero visibility (so any reflected light would be invisible). A lighting algorithm running on a GPU would not be able to efficiently branch to avoid unnecessary work, meaning that a huge amount of work will go wasted, while not enough lighting samples will be taken where it really matters (high density areas).

Looking at this, I made an observation on the similarities between this problem and the problem of order-independent transparency. That problems also essentially boils down to the problem of efficiently computing a visibility function for each pixel so that blending can be done by computing a weighted sum using the visibility function. An example of this is Adaptive OIT, which can efficiently compute a step curve which is an extremely good representation of a visibility function that is completely independent of scene depth complexity. This got me thinking on a different approach.

The froxel algorithm essentially works by quantizing the depth of the function into discrete values while keeping the visibility/density at a high precision. What if we quantized the visibility/density instead and kept the depth at a high precision? In other words, instead of storing the visibility at certain predefined depths, we store the high precision depth values of certain predefined visibility values! How would such a scheme do at representing the above curve?

Well, well, well, would you look at that? That is actually an extremely good representation of the original function with the same number of samples! There are no samples wasted on the vast empty space in front of the smoke cloud and no samples wasted behind it either. All the samples are put to excellent use!

This kind of representation has a large number of advantages compared to quantizing depth. In the previous algorithm, aside from the problems with depth complexity if the smoke cloud was translated towards or away from the camera, the result will wobble and flicker as the smoke cloud travels over the quantization points. This can result in inconsistent sorting of smoke and cause lots of unwanted lighting artifacts, which the Frostbite paper needed heavy temporal filtering to fix. This quantization is completely unaffected by both the depth complexity of the scene and translation of the smoke clouds. However, the greatest advantage of this scheme stems from the fact that this kind of representation also gives optimal lighting sample points. By properly choosing the predefined visibility values to quantize to, simply calculating the lighting at the sample points (and possibly points interpolated between these points), the lighting samples will be perfectly distributed to have the same contribution to the final scene. This is a huge advantage that can be abused to get away with much lower sample counts that are traditionally needed to light smoke, as every single sample will have maximum contribution to the final scene!

The algorithm seems to have lots of potential! Of course, there's just one minor detail left... Implementing the sucker on modern GPUs! Oh, what could possibly go wrong? Stay tuned for the next episode: The implementation, or "I hate GPUs."!

Background, or "Oh god I need to make a blog?"

So apparently I completely missed the fact that I need to write a blog for this project. I'm getting close to completing the project, so the first few blog posts will essentially be recaps of the progress so far.


I'll dive directly into the actual issue I'm investigating: Rendering smoke!

Smoke rendering is a very complex problem in real-time 3D graphics, because there are a lot of interacting parts. As it deals with transparency, the ordering becomes non-trivial and an order-independent rendering technique is ideal to avoid the requirement of sorting. Further complexity arises when lighting is supposed to be applied to the smoke, as the smoke can self-shadow. In addition, multiple separate smoke "systems" with different properties and scale (say a cloud, some cigarette smoke, fog on an early morning and atmospheric-scale scattering effects) can all interact and require merging and sorting together for accurate effects, as the different smoke clouds cannot be handled separately.

My project will NOT be about how to simulate smoke motion, and will not involve complicated self-shadowing. although the algorithm described IS compatible with such advanced techniques. Instead, the project will focus solely on the issue of sorting, merging and lighting multiple intersecting smoke clouds in an efficient way on modern GPUs.

Smoke rendering, like any kind of volumetric rendering, requires integration across the smoke volumes per pixel. This is essentialy ray-tracing, which is commonly avoided for real-time graphics due to its steep performance cost. Still, GPUs nowadays are fast enough to take hundreds or even thousands of samples for each pixel on the screen and still maintain real-time performance, but even with thousands of samples it is still often not enough to accurately represent effects of vastly different scale.

Physically-based & Unified Volumetric Rendering in Frostbite is a very interesting presentation that details an implementation of the most common modern approach to this problem. In essence, the view frustum is split up into "froxels" (a frustum-aligned voxel), which each contain the smoke parameters of a certain part of the frustum's volume. With the froxel buffer set up, different smoke clouds can be ray-traced independently and injected into the froxel buffer. Once all smoke clouds have been integrated into the froxel buffer, the buffer itself can the be lit by lights and shaded independently of the different smoke systems that were used to build the buffer. This means that a cloud, some cigarette smoke, some fog and the smoke from a large fire can all be integrated into a single data structure, which can then be lit independently of these smoke systems. This decouples smoke "sorting"/merging from smoke lighting, allowing for much more scalable performance. It's essentially the smoke version of deferred shading/lighting, which provides a similar advantage for opaque rendering where it decouples lighting from the objects and triangles that are being rendered.

Froxels have some big problems though. They use a lot of memory as they need to store a lot of volumetric data. In practice, this needs to be mitigated by significantly dropping the resolution of the smoke rendering, which makes it hard to catch detailed smoke effects like thin cigarette smoke. In addition, the limited number of layers along the depth of the camera can cause issues as well as the spacing of these depth layers depend on the depth complexity of the scene. For a large-scale outdoor scene, the layers will have to be spaced so far away from each other that the algorithm no longer provides a meaningful way of sorting smoke clouds along depth, instead simply merging the two into the same froxel layer. This can again increase blurriness and make it difficult to mix smoke effects on different scales with adequate quality.

The algorithm that I'll be investigating in this project provides a related approach, but there are some key differences that can provide significant advantages over the depth layer/"froxel" approach. My next blog post will detail the basic idea behind my new algorithm.