6 ms·

Implementing Mario's Stack Blur 15 times in C++ (with tests and benchmarks)

111 points by wlll 3y ago

bhouston 3y ago
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. [1] https://en.wikipedia.org/wiki/Separable_filter https://en.wikipedia.org/wiki/Separable_filter
- mrob 3y ago
  It's already implemented this way: https://melatonin.dev/blog/implementing-marios-stack-blur-15-times-in-cpp/#horizontal-and-vertical-stack-blur-passes https://melatonin.dev/blog/implementing-marios-stack-blur-15...
dxhdr 3y ago
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.
- evrimoztamur 3y ago
  From what I understand that’s what stack blur does, approximating the normal distribution.
  - Jasper_ 3y ago
    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: https://developer.apple.com/documentation/accelerate/vimage/pixelbuffer/3951789-tentconvolve https://developer.apple.com/documentation/accelerate/vimage/...
    - iainmerrick 3y ago
      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?
      - a_e_k 3y ago
        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.
    - pornel 3y ago
      Box blur applied twice gives the same result. It will quickly converge to Gaussian.
- sudara 3y ago
  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/blurring_an_image?language=objc https://developer.apple.com/documentation/accelerate/blurrin...
BryanLegend 3y ago
I was so disappointed that this had nothing to do with the new Mario Wonder game. :(
- sudara 3y ago
  ha!
Tempest1981 3y ago
How were the images in this web page created? Nicely done.
- sudara 3y ago
  Thanks! It's Figma, which I have some rudimentary proficiency in, really easy to make boxes and move em around.
lionkor 3y ago
instead of a deque or circular buffer, wouldnt it make more sense to use a span?
- sudara 3y ago
  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?
  - lionkor 3y ago
    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
shove 3y ago
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.
abbeyj 3y ago
> 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.
- sudara 3y ago
  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.
  - deleted 3y ago
    [deleted]
  - klyrs 3y ago
    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.
sebastianmestre 3y ago
In my experience, 3 passes of box blur looks good enough for small kernels (less than 50px). Is this faster than that?
- danbruc 3y ago
  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. This is the box filter inner loop. sum += data[right++] - data[left++] This is the stack blur inner loop. lsum += data[middle] - data[left++] rsum += data[right++] - data[middle++] sum += rsum - lsum
bitwize 3y ago
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? Well, regardless, today I learned.
LoganDark 3y ago
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 https://github.com/LoganDark/stackblur-iter Please fix :) AMA
- sudara 3y ago
  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!
  - LoganDark 3y ago
    > 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 :)
    - sudara 3y ago
      > 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!
corysama 3y ago
Probably would have been much easier to do 15 times in https://halide-lang.org/ https://halide-lang.org/ 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.
- sudara 3y ago
  OMG would have loved to do that instead.
TylerGlaiel 3y ago
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)
- sudara 3y ago
  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.
deleted 3y ago
[deleted]
anonymoushn 3y ago
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. 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/ 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.

111 points by wlll 3y ago

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. [1] https://en.wikipedia.org/wiki/Separable_filter https://en.wikipedia.org/wiki/Separable_filter

mrob 3y ago
It's already implemented this way: https://melatonin.dev/blog/implementing-marios-stack-blur-15-times-in-cpp/#horizontal-and-vertical-stack-blur-passes https://melatonin.dev/blog/implementing-marios-stack-blur-15...

dxhdr 3y ago

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.

evrimoztamur 3y ago
From what I understand that’s what stack blur does, approximating the normal distribution.
- Jasper_ 3y ago
  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: https://developer.apple.com/documentation/accelerate/vimage/pixelbuffer/3951789-tentconvolve https://developer.apple.com/documentation/accelerate/vimage/...
  - iainmerrick 3y ago
    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?
    - a_e_k 3y ago
      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.
  - pornel 3y ago
    Box blur applied twice gives the same result. It will quickly converge to Gaussian.
sudara 3y ago
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/blurring_an_image?language=objc https://developer.apple.com/documentation/accelerate/blurrin...

BryanLegend 3y ago

I was so disappointed that this had nothing to do with the new Mario Wonder game. :(

sudara 3y ago
ha!

Tempest1981 3y ago

How were the images in this web page created? Nicely done.

sudara 3y ago
Thanks! It's Figma, which I have some rudimentary proficiency in, really easy to make boxes and move em around.

lionkor 3y ago

instead of a deque or circular buffer, wouldnt it make more sense to use a span?

sudara 3y ago
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?
- lionkor 3y ago
  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

shove 3y ago

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.

abbeyj 3y ago

> 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.

sudara 3y ago
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.
- deleted 3y ago
  [deleted]
- klyrs 3y ago
  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.

sebastianmestre 3y ago

In my experience, 3 passes of box blur looks good enough for small kernels (less than 50px). Is this faster than that?

danbruc 3y ago
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. This is the box filter inner loop. sum += data[right++] - data[left++] This is the stack blur inner loop. lsum += data[middle] - data[left++] rsum += data[right++] - data[middle++] sum += rsum - lsum

bitwize 3y ago

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? Well, regardless, today I learned.

LoganDark 3y ago

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 https://github.com/LoganDark/stackblur-iter Please fix :) AMA

sudara 3y ago
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!
- LoganDark 3y ago
  > 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 :)
  - sudara 3y ago
    > 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!

corysama 3y ago

Probably would have been much easier to do 15 times in https://halide-lang.org/ https://halide-lang.org/ 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.

sudara 3y ago
OMG would have loved to do that instead.

TylerGlaiel 3y ago

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)

sudara 3y ago
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.

deleted 3y ago

[deleted]

anonymoushn 3y ago

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. 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/ 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.