Not really happy with the representation of the interval's length. Please send suggestions!

I'm having a surprisingly hard time grokking the code, will continue to try to figure it out, but is "moving average" conceptually comparable to the (quad) stack blur method where we use an accumulator so that we only have to do one addition, subtraction and division per data element? https://observablehq.com/@jobleonard/mario-klingemans-stackblur

What I did here is only to trim the code from [blur](https://github.com/d3/d3-contour/blob/master/src/blur.js), which is arguably a bit more readable. The algorithm is that for each line you start by summing r times the first value (sr = w * source[k0]), then iterate by adding to sr the new value (at i) and substracting the oldest value (at i - w). Then same for columns except the indices are reversed (i + nj becomes j + ni). In other words, and if I understand correctly, 1 iteration corresponds to the first part of your write-up, up until "The stackblur algorithm is taking this summing trick and applying it to itself one more time". (We don't have this extra trick, nor the third level of "quadratic stack blur".) Other key differences are that this runs only on 1 general purpose "channel" (of numbers not RGBA); uses Float32Array instead of Uint8Array; and needs to run 3 iterations to have something that looks good, instead of just 1 for stackblur. Because of all this (??) StackBlur seems to be faster on 2D color image data. On grayscale (1 channel), I have mixed results, depending on what you take into account (conversion from RGBA data to the Float32Array then back to RGBA is longer than 1 blurring iteration). So if we want to make this blur on par with stackblur I need to check if we can integrate the sum-of-sum trick into that loop, and how we can pass info like "work on channel 1 of 4".

I've just added a chart that shows how one pixel is blurred https://observablehq.com/@fil/moving-average-blur#PROFILE

I now have a variant (https://observablehq.com/d/63ded3c5ba9dcd1f) that takes the sum of sums approach (stack-blur). It was really difficult to get right (& still needs a bit of polish on the edges). Plus, the code is even more difficult to read… for no notable speed gain (maybe 10%). I'm leaving this for the moment.

Happy to report that with a bit a polishing we get again 30% more speed. Overall we're at 2.1ms (sum of sums) against 2.9ms (streamlined d3-contour) vs 6.4ms (d3-contour original code). Overall 3 times faster.

Nice! And somewhat reassuring that you also struggled, since I *still* haven't figured out where the off-by-one error is in my quad. stack blur code. Tangent: I'm kind of wondering whether the accuracy of the sum-of-sums and sum-of-sums-of-sums approaches turns out to be worse or better. Especially in the last approach the summed value can get very large, which can lead to rounding errors if the values being summed differ too much in order of magnitude. On the other hand the summing happens in steps, so the values that are being summed over should not differ too much in said orders of magnitude.

Yes that's why I divide by the blocks’ count before adding to the dr2 "integral". A very nice way to look at the result is to return any of s1, s2, (sr2 - 2 * sr1 + sr0) or dr instead of dr2, and the "profile" graph shows how the "bell shape" is constructed from the difference of two blocky rectangles, which become triangles (first integration) then "bell" (second integration). I've added a menu (profile variable) that allows to explore this. (This is how I was debugging, in fact, except I was changing the returned value in the code ^^).

You can also check this algorithm, proposed by Ivan Kutskir: http://blog.ivank.net/fastest-gaussian-blur.html A slightly modified version to handle Float arrays worked well for me with very good performance.

Nice suggestion, thank you! I've pasted the code in https://observablehq.com/d/d85e3cc1db985503 (removed Math.round so that it works with floats). At first glance the results seem to be on par with my proposal in terms of speed.