2D Radial blur is a separable filter [1], which means you can calculate it via first doing a 1D blur in one direction, and then using that intermediate result and doing a 1D blur in the other direction and it is equivalent to doing a 2D radial blur.
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.
One thing the article didn't mention is that a gaussian blur can be approximated with multiple passes of a box blur. Not sure how that would relate performance-wise to the discussed stack blur algorithm. The code would probably be a bit simpler, anyway.
Stack blur is a fast, discrete tent filter (weights are 1, 2, 3, 2, 1, normalize, it makes a triangle/tent shape). You can find tent filter blurs all over, including in Apple's own image processing library:
Does that actually look good? I would have thought it would have really obnoxious diagonal artifacts.
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?
The convolution of a box filter with a box filter is a tent filter. So yes, a tent filter blur can be done with two box filter blurs.
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.
I'm interested in this idea. I think I got confused at some point and mistakenly thought box blur was a 2D kernel and so it wouldn't perform great. vImage does contain a box blur but I haven't checked its performance (I did check the tent blur and it was so-so...) https://developer.apple.com/documentation/accelerate/blurrin...
Probably! The dequeue implementation was just for fun, as I miss working in higher level languages :) Home cooked circular buffers is usually how one usually sees real-time audio code written... I haven't had much hands on experience with std::span, but seeing as it's a non-owning view, maybe it would be perfect?
I must confess I didnt read the whole article yet, but yeah, if you are looking for a view onto some data, spans are awesome and surprisingly easy to use
Mario is absolutely legendary in the creative coding space. Everything he does has that “how the hell did he do that?!” vibe, from the Flash/ActionScript stuff back in the day to some of the only actually interesting ML image generation work I’ve seen.
> Each time the [std::deque] slides, every value in the queue has to change position in memory
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.
Thanks for the proofread! I had so much text juggling between this and the README that it was guaranteed some things would fall through the cracks! I updated the things you mentioned about `stackSum` and thanks for the catch on the `sumIn` definition.
> 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.
std::deque will never beat your hand-rolled ring buffer. It's a container, it will spend time managing memory every time you walk off the end of a chunk. If it's implemented right*, it will hold onto a chunk which fell off, and will only allocate twice. If it's implemented wrong*, it will allocate after `chunksize` pushes, and free every `chunksize` pops (which might be handled properly by your allocator). But it's still going to need to shuffle those chunks around, adding unnecessary overhead because you're using a ring in a way that is a perfect fit for your application: you use the value you're popping every time you push.
* for your purpose! Therein lies the challenge of writing standard libraries... choices must be made.
> std::deque will never beat your hand-rolled ring buffer.
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.)
After some thought, I figured out a way to implement a chunked array of the sort specified by the c++ standard library, which does not require reallocation if the queue length stays below the chunksize. It involves storing cyclic offsets for each chunk.
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.
Will be essentially the same performance. A box filter has to add one pixel entering the kernel and remove one pixel exiting the kernel. This one has one pixel entering and leaving each half of the kernel and then combining the two halfs. So three box blurs have six additions and once this filter is also six additions. Maybe the box filter could be somewhat slower because it has to increment two pointers three times while stack blur has to increment three pointers once and also incurs the loop overhead only once. For a fixed radius you could however use only one pointer and realize the offsets between the pointers with an appropriate addressing mode if available on the target architecture.
Every time I see an unadorned "Mario" I think of the video game plumber. So I was a bit stoked to see some sort of algorithm for a special effect for Super Mario Bros. Wonder that escaped my notice, maybe a way to blur the player sprite for invincibility power-up effects by stacking multiple copies of the sprite of varying opacity on top of one another?
Oh hey, I'm credited in the article! However, that rustrepo website you linked is a really bad mirror of the actual GitHub repo: https://github.com/LoganDark/stackblur-iter
Fixed! Sorry about that. The rustrepo stuff confused me!
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!
> My question is — was solving the edge bleed worth the (assumed) performance tradeoff?
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 :)
> If you take a look at the video in my README, you'll see how bad it is.
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!
It seems like it should still be a win to use libdivide or something like libdivide even if the divisor varies within a row - you get to reuse those divisors for each row!
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.
Biggest place for improvement performance-wise I can see from this is that you aren't taking into account cache coherency here. Horizontal blurs with this method are going to be ripping fast, but vertical blurs are going to constantly cache miss. There's a number of ways to potentially fix this, the fact that you only need to store 3 sums for one line means you could do a few (64) columns at a time and just store that many distinct sets of sums (sorta like doing 64 columns in parallel, just making sure to use the data thats in the cache while its still in the cache)
Cache locality and specifically the vertical pass was top of my mind when trying to come up with good ways to vectorize. In the end (at least in my vector implementations) the difference between the passes weren't too large. But most of the them ended up having to do things like first convert the incoming row/col to its own float vector.
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.
It seems like for the 1d case you should be able to compute the double-prefix sum (the prefix sum of the prefix sum of the input row) and then add 2 elements of that array and subtract 2 other elements to get the output pixel.
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.
40 comments
[ 1.8 ms ] story [ 35.3 ms ] threadCould 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.
[1] https://en.wikipedia.org/wiki/Separable_filter
https://melatonin.dev/blog/implementing-marios-stack-blur-15...
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.
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.
This is the box filter inner loop.
This is the stack blur inner loop.Well, regardless, today I learned.
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. :)
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.
Here 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.