Saturday, May 27, 2017

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!

No comments:

Post a Comment