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

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.