Implementing Mario's Stack Blur 15 times in C++ (with tests and benchmarks)
melatonin.dev
melatonin.dev
I don't think that this is true. `std::deque` has guaranteed O(1) complexity for insertions or deletions at the beginning or end of the container. It couldn't provide this if it was moving every value in the queue.
> How do we alter the stackSum when the queue moves?
You never quite answer this question. How are `sumOut` and `sumIn` related to `stackSum`? This is available in the code but I feel like there should at least be a one-line summary available without having to click through to the code.
> sumOut is what we need to add to the stack.
This is supposed to be `sumIn`?
The examples that you provide use a radius of 2, thus having a `sizeOfStack` of (radius+1)*2 = 9 and requiring a division by 9 for each pixel. Have you considered using a radius of 1 or a radius of 3 instead? This would require division by (1+1)*2 = 4 or (3+1)*2 = 16 respectively which could then be done with a bitshift instead of a division. You'd have to template the code on `radius` instead of passing it in as a runtime parameter so that the compiler could lower the divisions to bitshifts.
> It couldn't provide this if it was moving every value in the queue.
I actually don't remember anymore what std::deque does under the hood, I did look into it, but the only thing I remember is that it was quite slow!
> You'd have to template the code on `radius` instead of passing it in as a runtime parameter so that the compiler could lower the divisions to bitshifts.
Yes, I really like this idea. Especially because radii only really vary between 1-48px for most drop shadow needs. It would be nice to have a handful of the common radii be ripping fast.
* for your purpose! Therein lies the challenge of writing standard libraries... choices must be made.
std::deque may never beat a hand-rolled ring buffer of fixed size that maintains cursors, but it will beat a hand-rolled linked-list implementation.
The typical implementation will generally exhibit better cache locality than a linked list, and will outperform it for random access as required by the specification (amortized constant, vs linear).
(The typical implementation per cppreference, although I don't think this is formally required by the spec, is a sequence of fixed-length arrays.)
So what I wrote isn't quite correct about needing to fiddle with memory regardless. But there is still overhead -- even if there's well-optimized code for the single chunk case, the container needs to check that condition.
Could one modify the stack blur algorithm to also operate in a separable fashion? It seems to me that one can combine Stack Blur with the separable concept to get the best of both worlds.
Of course traversing the data in a way not aligned with its in-memory data layout may cause issues, so one would have to figure out an efficient way to do that. I can think of a few ways of doing that using prefetching.
https://melatonin.dev/blog/implementing-marios-stack-blur-15...
The idea behind Halide is that scheduling memory access patterns is critical to performance. But, access patterns being interwoven into arithmetic algorithms makes them difficult to modify separately.
So, in Halide you specify the arithmetic and the schedule separately so you can rapidly iterate on either.
One main issue I never resolved is in the middle of the main loop, data has to be converted and written back to the source image and the incoming pixels have to be converted and loaded in. Even when doing all rows or cols in bulk (which was always faster somehow than doing batches of 32/64), that seemed pretty brutal.
I also wondered whether it might be more efficient to rotate the entire image before and after the vertical pass, but in my implementations at least, there wasn't a huge difference in the pass timings.
Well, regardless, today I learned.
https://developer.apple.com/documentation/accelerate/vimage/...
Unless it can be tuned to behave like a repeated box blur (and therefore approximate a true Gaussian). A tent shape would just be two box blurs, right?
And repeated box filters asymptotically approach a Gaussian filter, though in practice it converges very quickly. A single box filter can be thought of as a piece-wise constant approximation of a Gaussian filter. Two passes of box filters produce a piece-wise linear approximation (the tent). And three passes of box filters produce a piece-wise quadratic approximation. Four passes for a cubic approximation, etc. Usually, just the three box filters are used since they are close enough.
I'd need to look at the math in detail here, but I suspect that it's ultimately just cleverly interleaving the running sums for the two box filter to perform them in a single pass.
If so, it might be an interesting challenge to see if the idea can be extended to compute the higher quality three box filter approximation in a single pass.
This is the box filter inner loop.
sum += data[right++] - data[left++]
lsum += data[middle] - data[left++]
rsum += data[right++] - data[middle++]
sum += rsum - lsumHere are notes about how to compute the prefix sum, which can be adapted to compute the double-prefix sum: https://en.algorithmica.org/hpc/algorithms/prefix/
Since you access members of the double-prefix sum array in a fixed pattern you can also use SIMD here, i.e. load x4, add, add, sub, divide by a constant = (mul, shift), shuffle, store to compute 8 pixels on AVX2. This assumes you have a radius > 15 and are using u32 elements, for <= 15 you can use u16 elements in the double-prefix sum array, but I'm not sure if this is a win because I don't know how to do implement the divide step. It's cheap to do only part of the row and pick up the double-prefix sum computation later, so length of the buffer you work on should be chosen to fit in L1.
This technique should work with padding: just pad the array before computing the double-prefix sum, so the full double-prefix sum array has n + 2r + 2 entries, but as above your work buffer should be smaller than this if the image is very wide, and the whole array doesn't need to exist at the same time. It should also work without padding: you need to do libdivide or fastmod-style computations with ~r different fixed divisors at the edge of the image. This is still worth it, you'll be doing dozens or thousands of divisions with each edgy divisor. An ugly corner case appears here when 2r > n. Sorry.
Once your 1d case is much, much faster, I would expect transposing and then transposing back to be a clear win over not transposing at all, because SIMD gather instructions are dreadful and if you want to index individual bytes with a large stride they don't even exist.
I'm interested in comparing this to the approach used in TFA, but the github repository seems to be missing a bunch of header files...
As an aside, some other comments mentioned that stack blur is "a box blur of a box blur" and the "prefix sum of a prefix sum" idea in this comment corresponds to that "box blur of a box blur" idea. If you wanted, it wouldn't be that much more expensive to add 1 more level of prefix summing (= 1 more level of box blur) to get a kernel that looks much more like a gaussian.
Please fix :)
AMA
My question is — was solving the edge bleed worth the (assumed) performance tradeoff? There were a couple vendor APIs that seemed to have smarter edge bleed options, but they performed worse. I never got around to actually visually comparing the images, though...
Thanks for all the great comments in your repo, they were quite helpful when trying to figure out how Stack Blur worked!
Absolutely. Edge bleed is a distracting visual artifact! If you take a look at the video in my README, you'll see how bad it is. It's slightly less bad if you implement the sRGB transfer functions properly, but it's still bad.
As for the performance tradeoff - it's probably not as big as you think. However, I can say there is a huge performance tradeoff in using real division, and the only reason I can't use libdivide is the fact that my denominator changes. That's some 30% of the runtime - replacing it with a multiplication (as in libdivide) would probably shave at least 20% off the runtime of my algorithm.
> Thanks for all the great comments in your repo, they were quite helpful when trying to figure out how Stack Blur worked!
You're welcome! I will admit, I mainly wanted to prove that I actually came up with a better version of the algorithm, and that my library was therefore the best :)
Ahh, the video threw me off originally, I think because of the FPS glitches on the full width — but I see it now! Big bands of blue on the left and right edges. I'll have to cook up some examples to reproduce. Also curious what it means (if anything) in the drop shadow context.
> I can say there is a huge performance tradeoff in using real division
Makes sense. It would also throw a wrench in my vector implementations, where the expectation is to perform the same operation efficiently across groups of pixels.
> and that my library was therefore the best :)
Haha, well it was for me! I don't know Rust, but it helped me figure out my first C++ implementations!
Nothing if the edges of your bitmap (up to the blur's radius) are transparent. :)