Repository: d3/d3-array
Pull request #151 (2020-06-13 22:11:10)

d3.blur
fixes https://github.com/d3/d3-array/issues/56

see https://observablehq.com/@fil/moving-average-blur

- [x] remove iterations from the public API (keep it only for debugging)
- [x] <del>stackblur algo (a bit faster, but more difficult to understand)</del>
- [x] consistent ES6
- [x] verify that it can be used by d3-contour, in particular:
  - [x] allow different radius on x/y (https://github.com/d3/d3-contour/issues/27)
  - [x] the shape around a single dot should be circular(ish)
- [x] add fractional blurring (radius = .5)


curran: Amazing!
Fil: I'm pushing a new version which reduces the problem of the extremities being "dampened" by duplicating (*radius* times) the value at each extremity.

Bonus, this approach makes it even faster (we're now under 50% of the original time).

before:
![before1](https://user-images.githubusercontent.com/7001/84588821-c11d3e00-ae2a-11ea-9cd4-d5b63d3b506f.png)
after:
![after1](https://user-images.githubusercontent.com/7001/84588823-c2e70180-ae2a-11ea-9d4d-7beb3ab0d240.png)



On the random image and a large radius, it creates a visible artifact: the values on the edges tends to impose themselves. Not strictly better or worse than the original, only different. (The solution when this is a problem is to compose first a small blur radius, then a larger one.)

before:
<img width="397" alt="before" src="https://user-images.githubusercontent.com/7001/84588901-764ff600-ae2b-11ea-9012-fb9ffadd9d99.png">
after:
<img width="398" alt="after" src="https://user-images.githubusercontent.com/7001/84588902-77812300-ae2b-11ea-936b-b923174e15ae.png">

curran: Fascinating. Yeah I thought there was something going on with the extremities being "dampened", with the Gaussian blur approach in general. It's great that your implementation accounts for this!

Some interesting related reading that might touch upon this issue:

 * [Discrete Gaussian](https://en.wikipedia.org/wiki/Gaussian_function#Discrete_Gaussian)
 * [The discrete Gaussian kernel](https://en.wikipedia.org/wiki/Scale_space_implementation#The_discrete_Gaussian_kernel)
Fil: Another use case for this type of blurring: "hue" blurring on an image. We first read each pixel's value in the Lab model (or HCL, HSL…), then blur a and b, and reassemble the image:
https://observablehq.com/d/b729443cd180e56f
Fil: Finally I don't think the rgba thing is in the scope of this pull-request. It creates a bit of complications (stride, channel, dtype) "just" for specific optimizations, and it's possible to do without.
Fil: I've added a notion that the shape created around a dot (or single pixel) should be circular(ish), so that d3-contours are more "round" than "square". Unfortunately it seems that the stackblur's distribution profile results in shapes that are a bit squarish:

simple sum * 3 iterations
![sum blur](https://user-images.githubusercontent.com/7001/85152409-42266c00-b255-11ea-8b79-3dcef4eacdf7.png)

sum of sums
![sumsum (stack) blur](https://user-images.githubusercontent.com/7001/85152399-405ca880-b255-11ea-8972-c38b04e34dad.png)

So I think I have to either go back on that part, offer the two possibilities, or implement it separately in d3-contours.
Fil: I went back to the original algorithm (to have "circular" shapes for d3-contour), and introduced non-integer radius, which help blur "just a little", for example with blur().radius(0.1).

I feel this PR is now complete.
curran: Wow, this is awesome! Thank you @Fil for your work on this.
Fil: Another demo with the "cloud contours" notebook: https://observablehq.com/d/c9d89106eafefb4f
curran: Wow that is so cool!
Fil: I find myself using this code more often than I thought. For example, when doing image FFT, to display the contours of a magnitude chart (which is too spiky if you take the raw values).
curran: Nice. This PR is super cool. Can't wait for it to be merged!

Is there any more work needed here to complete the PR? Looks ready to me.
curran: If this could improve `d3-contour`, _and_ provide new generic blur functionality independent of `d3-contour`, could be a big win for D3 on the whole.
mbostock: This is a somewhat superficial reaction, but I wonder what this would look like without the fluent-style API. So, instead of:

```js
const blurred = d3.blur()
      .radius(radius)
      .width(detail.width)
      .value(d => d.b)
  (pixels);
```

You might say:

```js
const blurred = d3.blur(pixels, radius, detail.width, d => d.b);
```

The other thing that feels weird is that d3-array methods are, generally, 1-dimensional methods. Here d3.blur is two-dimensional, expressed in terms of *x* and *y* (or horizontal and vertical). This makes me think that this method should perhaps live in a different module that’s more related to 2D images. Alternatively, would it be possible to have a slightly lower-level API with an optional [*stride*](https://en.wikipedia.org/wiki/Stride_of_an_array) argument instead (where *stride* = 1 for “horizontal” and *stride* = *width*? for vertical), similar to [ndarray](https://github.com/scijs/ndarray)?
mbostock: Another consideration is whether d3.blur should happen in-place, and eliminate the accessor. Would that be possible or does it require a copy?

```js
const blurred = Float64Array.from(pixels, d => d.b);
d3.blur(blurred, radius, detail.width);
```
Fil: It's also 1-d, with a no-op vertical blur in that case. Internally there are two copies, that we swap at each iteration, so I'm not sure it can be done "in place", though maybe by tracing the parity of swaps it might be possible to save 1 copy allocation. In practice though I think I prefer this "immutable" version.

(As I needed it on several experiments I was just at this instant publishing it on github and npm (as array-blur) — but happy to revert that if it can be integrated into D3 in some way.)

curran: Thrilled to see that this was published in https://github.com/Fil/array-blur !
cdaringe: traversing thru the observablehq forums lead the thru the d3 docs and back to here. was looking for EMWA or movingAverage for 1d data, and found array-blur instead (which is actually more appropriate for my use case anyway :)). great work fil! are the outstanding blockers around whether or not such functionality is relevant/welcome here, or rather outstanding qualms with the interface? just curious :)
Fil: Thank you, Christopher!

While it may never be included in d3-array, you can use this function by calling the [array-blur](https://github.com/Fil/array-blur) module. I've packaged it just like the d3 modules, which means you can even attach it to the d3 symbol if you like, with:
```js
const d3 = require("d3@7", "array-blur@1");
const blur = d3.blur().radius(1.3).width(w);
blur(pixels); // blurred image
```

Another possibility is that it would exist as a stand-alone module d3-blur (_i.e._ a renaming of array-blur); this would keep this 2D notion separate from d3-array, and allow d3-contour to use this faster and better blurring method.

I've made this notebook to summarize various use cases, in particular how to use it with Observable Plot (by passing _blur().radius(r)_ to Plot.mapY):
https://observablehq.com/@fil/array-blur-d3-plot
Fil: We'll need two documentation notebooks:
* https://observablehq.com/d3/d3-blur (we can show the profile of the kernel in 1d and 2d, also how to use it on a grayscale data matrix)
* https://observablehq.com/d3/d3-blurimage (showing image blur and laplacian)

Patch:
diff --git a/README.md b/README.md
index 1bcab86..f6c9adf 100644
--- a/README.md
+++ b/README.md
@@ -683,6 +683,44 @@ Returns an array of arrays, where the *i*th array contains the *i*th element fro
 d3.zip([1, 2], [3, 4]); // returns [[1, 3], [2, 4]]
 ```
 
+#### Blur
+
+<a name="blur" href="#blur">#</a> d3.<b>blur</b>(*data*, *radius*) · [Source](https://github.com/d3/d3-array/blob/main/src/blur.js), [Examples](https://observablehq.com/@d3/d3-blur)
+
+Blurs an array of *data* in-place by applying three iterations of a moving average transform, for a fast approximation of a gaussian kernel of the given *radius*, a non-negative number, and returns the array.
+
+```js
+const randomWalk = d3.cumsum({length: 1000}, () => Math.random() - 0.5);
+blur(randomWalk, 5);
+```
+
+Copy the data if you don’t want to smooth it in-place:
+```js
+const smoothed = blur(randomWalk.slice(), 5);
+```
+
+<a name="blur2" href="#blur2">#</a> d3.<b>blur2</b>({*data*, *width*[, *height*]}, *rx*[, *ry*]) · [Source](https://github.com/d3/d3-array/blob/main/src/blur.js), [Examples](https://observablehq.com/@d3/d3-blur)
+
+Blurs a matrix of the given *width* and *height* in-place, by applying an horizontal blur of radius *rx* and a vertical blur or radius *ry* (which defaults to *rx*). The matrix *data* is stored in a flat array, used to determine the *height* if it is not specified. Returns the blurred {data, width, height}.
+
+```js
+data = [
+  1, 0, 0,
+  0, 0, 0,
+  0, 0, 1
+];
+blur2({data, width: 3}, 1);
+```
+
+<a name="blurImage" href="#blurImage">#</a> d3.<b>blurImage</b>(*imageData*, *rx*[, *ry*]) · [Source](https://github.com/d3/d3-array/blob/main/src/blur.js), [Examples](https://observablehq.com/@d3/d3-blurimage)
+
+Blurs an [ImageData](https://developer.mozilla.org/en-US/docs/Web/API/ImageData) structure in-place, blurring each of the RGBA layers independently by applying an horizontal blur of radius *rx* and a vertical blur or radius *ry* (which defaults to *rx*). Returns the blurred ImageData.
+
+```js
+const imData = context.getImageData(0, 0, width, height);
+blurImage(imData, 5);
+```
+
 ### Iterables
 
 These are equivalent to built-in array methods, but work with any iterable including Map, Set, and Generator.
diff --git a/src/blur.js b/src/blur.js
new file mode 100644
index 0000000..7021cbf
--- /dev/null
+++ b/src/blur.js
@@ -0,0 +1,115 @@
+export function blur(values, r) {
+  if (!((r = +r) >= 0)) throw new RangeError("invalid r");
+  let length = values.length;
+  if (!((length = Math.floor(length)) >= 0)) throw new RangeError("invalid length");
+  if (!length || !r) return values;
+  const blur = blurf(r);
+  const temp = values.slice();
+  blur(values, temp, 0, length, 1);
+  blur(temp, values, 0, length, 1);
+  blur(values, temp, 0, length, 1);
+  return values;
+}
+
+export const blur2 = Blur2(blurf);
+
+export const blurImage = Blur2(blurfImage);
+
+function Blur2(blur) {
+  return function(data, rx, ry = rx) {
+    if (!((rx = +rx) >= 0)) throw new RangeError("invalid rx");
+    if (!((ry = +ry) >= 0)) throw new RangeError("invalid ry");
+    let {data: values, width, height} = data;
+    if (!((width = Math.floor(width)) >= 0)) throw new RangeError("invalid width");
+    if (!((height = Math.floor(height !== undefined ? height : values.length / width)) >= 0)) throw new RangeError("invalid height");
+    if (!width || !height || (!rx && !ry)) return data;
+    const blurx = rx && blur(rx);
+    const blury = ry && blur(ry);
+    const temp = values.slice();
+    if (blurx && blury) {
+      blurh(blurx, temp, values, width, height);
+      blurh(blurx, values, temp, width, height);
+      blurh(blurx, temp, values, width, height);
+      blurv(blury, values, temp, width, height);
+      blurv(blury, temp, values, width, height);
+      blurv(blury, values, temp, width, height);
+    } else if (blurx) {
+      blurh(blurx, values, temp, width, height);
+      blurh(blurx, temp, values, width, height);
+      blurh(blurx, values, temp, width, height);
+    } else if (blury) {
+      blurv(blury, values, temp, width, height);
+      blurv(blury, temp, values, width, height);
+      blurv(blury, values, temp, width, height);
+    }
+    return data;
+  };
+}
+
+function blurh(blur, T, S, w, h) {
+  for (let y = 0, n = w * h; y < n;) {
+    blur(T, S, y, y += w, 1);
+  }
+}
+
+function blurv(blur, T, S, w, h) {
+  for (let x = 0, n = w * h; x < w; ++x) {
+    blur(T, S, x, x + n, w);
+  }
+}
+
+function blurfImage(radius) {
+  const blur = blurf(radius);
+  return (T, S, start, stop, step) => {
+    start <<= 2, stop <<= 2, step <<= 2;
+    blur(T, S, start + 0, stop + 0, step);
+    blur(T, S, start + 1, stop + 1, step);
+    blur(T, S, start + 2, stop + 2, step);
+    blur(T, S, start + 3, stop + 3, step);
+  };
+}
+
+// Given a target array T, a source array S, sets each value T[i] to the average
+// of {S[i - r], …, S[i], …, S[i + r]}, where r = ⌊radius⌋, start <= i < stop,
+// for each i, i + step, i + 2 * step, etc., and where S[j] is clamped between
+// S[start] (inclusive) and S[stop] (exclusive). If the given radius is not an
+// integer, S[i - r - 1] and S[i + r + 1] are added to the sum, each weighted
+// according to r - ⌊radius⌋.
+function blurf(radius) {
+  const radius0 = Math.floor(radius);
+  if (radius0 === radius) return bluri(radius);
+  const t = radius - radius0;
+  const w = 2 * radius + 1;
+  return (T, S, start, stop, step) => { // stop must be aligned!
+    if (!((stop -= step) >= start)) return; // inclusive stop
+    let sum = radius0 * S[start];
+    const s0 = step * radius0;
+    const s1 = s0 + step;
+    for (let i = start, j = start + s0; i < j; i += step) {
+      sum += S[Math.min(stop, i)];
+    }
+    for (let i = start, j = stop; i <= j; i += step) {
+      sum += S[Math.min(stop, i + s0)];
+      T[i] = (sum + t * (S[Math.max(start, i - s1)] + S[Math.min(stop, i + s1)])) / w;
+      sum -= S[Math.max(start, i - s0)];
+    }
+  };
+}
+
+// Like blurf, but optimized for integer radius.
+function bluri(radius) {
+  const w = 2 * radius + 1;
+  return (T, S, start, stop, step) => { // stop must be aligned!
+    if (!((stop -= step) >= start)) return; // inclusive stop
+    let sum = radius * S[start];
+    const s = step * radius;
+    for (let i = start, j = start + s; i < j; i += step) {
+      sum += S[Math.min(stop, i)];
+    }
+    for (let i = start, j = stop; i <= j; i += step) {
+      sum += S[Math.min(stop, i + s)];
+      T[i] = sum / w;
+      sum -= S[Math.max(start, i - s)];
+    }
+  };
+}
diff --git a/src/index.js b/src/index.js
index 5fb58e8..29fc334 100644
--- a/src/index.js
+++ b/src/index.js
@@ -1,6 +1,7 @@
 export {default as bisect, bisectRight, bisectLeft, bisectCenter} from "./bisect.js";
 export {default as ascending} from "./ascending.js";
 export {default as bisector} from "./bisector.js";
+export {blur, blur2, blurImage} from "./blur.js";
 export {default as count} from "./count.js";
 export {default as cross} from "./cross.js";
 export {default as cumsum} from "./cumsum.js";


Test patch:
diff --git a/test/blur-test.js b/test/blur-test.js
new file mode 100644
index 0000000..014757a
--- /dev/null
+++ b/test/blur-test.js
@@ -0,0 +1,202 @@
+import assert from "assert";
+import {blur, blur2} from "../src/index.js";
+
+it("blur(values, r) returns values", () => {
+  const V = [0, 0, 0, 0, 0, 0, 27, 0, 0, 0, 0, 0, 0, 0];
+  assert.strictEqual(blur(V, 1), V);
+  assert.deepStrictEqual(V, [0, 0, 0, 1, 3, 6, 7, 6, 3, 1, 0, 0, 0, 0]);
+});
+
+it("blur(values, r) observes the expected integer radius r", () => {
+  assert.deepStrictEqual(blur([0, 0, 0, 0, 0, 0, 27, 0, 0, 0, 0, 0, 0, 0], 0.0).map(round), [0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 27.00, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000]);
+  assert.deepStrictEqual(blur([0, 0, 0, 0, 0, 0, 27, 0, 0, 0, 0, 0, 0, 0], 1.0).map(round), [0.000, 0.000, 0.000, 1.000, 3.000, 6.000, 7.000, 6.000, 3.000, 1.000, 0.000, 0.000, 0.000, 0.000]);
+  assert.deepStrictEqual(blur([0, 0, 0, 0, 0, 0, 27, 0, 0, 0, 0, 0, 0, 0], 2.0).map(round), [0.216, 0.648, 1.296, 2.160, 3.240, 3.888, 4.104, 3.888, 3.240, 2.160, 1.296, 0.648, 0.216, 0.000]);
+  assert.deepStrictEqual(blur([0, 0, 0, 0, 0, 0, 27, 0, 0, 0, 0, 0, 0, 0], 3.0).map(round), [1.023, 1.338, 1.732, 2.204, 2.598, 2.834, 2.913, 2.834, 2.598, 2.204, 1.653, 1.181, 0.787, 0.472]);
+});
+
+it("blur(values, r) observes the expected fractional radius r", () => {
+  assert.deepStrictEqual(blur([0, 0, 0, 0, 0, 0, 27, 0, 0, 0, 0, 0, 0, 0], 0.5).map(round), [0.000, 0.000, 0.000, 0.422, 2.531, 6.328, 8.438, 6.328, 2.531, 0.422, 0.000, 0.000, 0.000, 0.000]);
+  assert.deepStrictEqual(blur([0, 0, 0, 0, 0, 0, 27, 0, 0, 0, 0, 0, 0, 0], 1.5).map(round), [0.053, 0.316, 0.949, 2.004, 3.322, 4.430, 4.852, 4.430, 3.322, 2.004, 0.949, 0.316, 0.053, 0.000]);
+  assert.deepStrictEqual(blur([0, 0, 0, 0, 0, 0, 27, 0, 0, 0, 0, 0, 0, 0], 2.5).map(round), [0.672, 1.078, 1.609, 2.234, 2.813, 3.188, 3.313, 3.188, 2.813, 2.234, 1.594, 1.031, 0.594, 0.281]);
+  assert.deepStrictEqual(blur([0, 0, 0, 0, 0, 0, 27, 0, 0, 0, 0, 0, 0, 0], 3.5).map(round), [1.266, 1.503, 1.780, 2.057, 2.294, 2.452, 2.505, 2.452, 2.294, 2.030, 1.701, 1.371, 1.081, 0.844]);
+});
+
+it("blur(values, r) repeats starting values before the window", () => {
+  assert.deepStrictEqual(blur([27, 0, 0, 0, 0, 0, 0, 0], 0.0).map(round), [27.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000]);
+  assert.deepStrictEqual(blur([27, 0, 0, 0, 0, 0, 0, 0], 1.0).map(round), [13.000, 9.000, 4.000, 1.000, 0.000, 0.000, 0.000, 0.000]);
+  assert.deepStrictEqual(blur([27, 0, 0, 0, 0, 0, 0, 0], 2.0).map(round), [11.016, 9.072, 6.696, 4.104, 2.160, 0.864, 0.216, 0.000]);
+  assert.deepStrictEqual(blur([27, 0, 0, 0, 0, 0, 0, 0], 3.0).map(round), [10.233, 8.974, 7.478, 5.825, 4.093, 2.676, 1.574, 0.787]);
+});
+
+it("blur(values, r) approximately preserves total value", () => {
+  assert.strictEqual(blur([0, 0, 0, 0, 0, 0, 27, 0, 0, 0, 0, 0, 0, 0], 0.0).reduce((p, v) => p + v), 27);
+  assert.strictEqual(blur([0, 0, 0, 0, 0, 0, 27, 0, 0, 0, 0, 0, 0, 0], 0.5).reduce((p, v) => p + v), 27);
+  assert.strictEqual(blur([0, 0, 0, 0, 0, 0, 27, 0, 0, 0, 0, 0, 0, 0], 1.0).reduce((p, v) => p + v), 27);
+  assert.strictEqual(blur([0, 0, 0, 0, 0, 0, 27, 0, 0, 0, 0, 0, 0, 0], 1.5).reduce((p, v) => p + v), 27);
+  assert.strictEqual(blur([0, 0, 0, 0, 0, 0, 27, 0, 0, 0, 0, 0, 0, 0], 2.0).reduce((p, v) => p + v), 27.000000000000004);
+  assert.strictEqual(blur([0, 0, 0, 0, 0, 0, 27, 0, 0, 0, 0, 0, 0, 0], 2.5).reduce((p, v) => p + v), 26.640625);
+  assert.strictEqual(blur([0, 0, 0, 0, 0, 0, 27, 0, 0, 0, 0, 0, 0, 0], 3.0).reduce((p, v) => p + v), 26.370262390670547);
+});
+
+const unit = {
+  width: 11,
+  height: 11,
+  data: [
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+    0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0,
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0
+  ]
+};
+
+it("blur2(data, r) modifies in-place", () => {
+  const copy = copy2(unit);
+  assert.strictEqual(blur2(copy, 1), copy);
+});
+
+it("data.height is redundant for blur2", () => {
+  const copy = copy2(unit);
+  delete copy.height;
+  assert.deepStrictEqual(blur2(copy, 1).data.map(round), [
+    0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000,
+    0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000,
+    0.000, 0.000, 0.001, 0.004, 0.008, 0.010, 0.008, 0.004, 0.001, 0.000, 0.000,
+    0.000, 0.000, 0.004, 0.012, 0.025, 0.029, 0.025, 0.012, 0.004, 0.000, 0.000,
+    0.000, 0.000, 0.008, 0.025, 0.049, 0.058, 0.049, 0.025, 0.008, 0.000, 0.000,
+    0.000, 0.000, 0.010, 0.029, 0.058, 0.067, 0.058, 0.029, 0.010, 0.000, 0.000,
+    0.000, 0.000, 0.008, 0.025, 0.049, 0.058, 0.049, 0.025, 0.008, 0.000, 0.000,
+    0.000, 0.000, 0.004, 0.012, 0.025, 0.029, 0.025, 0.012, 0.004, 0.000, 0.000,
+    0.000, 0.000, 0.001, 0.004, 0.008, 0.010, 0.008, 0.004, 0.001, 0.000, 0.000,
+    0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000,
+    0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000
+  ]);
+});
+
+it("blur2(data, r) observes the expected integer radius r", () => {
+  assert.deepStrictEqual(blur2(copy2(unit), 0), unit);
+  assert.deepStrictEqual(blur2(copy2(unit), 1).data.map(round), [
+    0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000,
+    0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000,
+    0.000, 0.000, 0.001, 0.004, 0.008, 0.010, 0.008, 0.004, 0.001, 0.000, 0.000,
+    0.000, 0.000, 0.004, 0.012, 0.025, 0.029, 0.025, 0.012, 0.004, 0.000, 0.000,
+    0.000, 0.000, 0.008, 0.025, 0.049, 0.058, 0.049, 0.025, 0.008, 0.000, 0.000,
+    0.000, 0.000, 0.010, 0.029, 0.058, 0.067, 0.058, 0.029, 0.010, 0.000, 0.000,
+    0.000, 0.000, 0.008, 0.025, 0.049, 0.058, 0.049, 0.025, 0.008, 0.000, 0.000,
+    0.000, 0.000, 0.004, 0.012, 0.025, 0.029, 0.025, 0.012, 0.004, 0.000, 0.000,
+    0.000, 0.000, 0.001, 0.004, 0.008, 0.010, 0.008, 0.004, 0.001, 0.000, 0.000,
+    0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000,
+    0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000
+  ]);
+  assert.deepStrictEqual(blur2(copy2(unit), 2).data.map(round), [
+    0.001, 0.001, 0.002, 0.003, 0.003, 0.004, 0.003, 0.003, 0.002, 0.001, 0.001,
+    0.001, 0.002, 0.004, 0.006, 0.007, 0.007, 0.007, 0.006, 0.004, 0.002, 0.001,
+    0.002, 0.004, 0.006, 0.010, 0.012, 0.012, 0.012, 0.010, 0.006, 0.004, 0.002,
+    0.003, 0.006, 0.010, 0.014, 0.017, 0.018, 0.017, 0.014, 0.010, 0.006, 0.003,
+    0.003, 0.007, 0.012, 0.017, 0.021, 0.022, 0.021, 0.017, 0.012, 0.007, 0.003,
+    0.004, 0.007, 0.012, 0.018, 0.022, 0.023, 0.022, 0.018, 0.012, 0.007, 0.004,
+    0.003, 0.007, 0.012, 0.017, 0.021, 0.022, 0.021, 0.017, 0.012, 0.007, 0.003,
+    0.003, 0.006, 0.010, 0.014, 0.017, 0.018, 0.017, 0.014, 0.010, 0.006, 0.003,
+    0.002, 0.004, 0.006, 0.010, 0.012, 0.012, 0.012, 0.010, 0.006, 0.004, 0.002,
+    0.001, 0.002, 0.004, 0.006, 0.007, 0.007, 0.007, 0.006, 0.004, 0.002, 0.001,
+    0.001, 0.001, 0.002, 0.003, 0.003, 0.004, 0.003, 0.003, 0.002, 0.001, 0.001
+  ]);
+});
+
+it("blur2(data, rx, 0) does horizontal blurring", () => {
+  assert.deepStrictEqual(blur2(copy2(unit), 0, 0).data, [
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+    0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0,
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0
+  ]);
+  assert.deepStrictEqual(blur2(copy2(unit), 1, 0).data.map(round), [
+    0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000,
+    0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000,
+    0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000,
+    0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000,
+    0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000,
+    0.000, 0.000, 0.037, 0.111, 0.222, 0.259, 0.222, 0.111, 0.037, 0.000, 0.000,
+    0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000,
+    0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000,
+    0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000,
+    0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000,
+    0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000
+  ]);
+  assert.deepStrictEqual(blur2(copy2(unit), 2, 0).data.map(round), [
+    0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000,
+    0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000,
+    0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000,
+    0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000,
+    0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000,
+    0.024, 0.048, 0.080, 0.120, 0.144, 0.152, 0.144, 0.120, 0.080, 0.048, 0.024,
+    0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000,
+    0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000,
+    0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000,
+    0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000,
+    0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000
+  ]);
+});
+
+it("blur2(data, 0, ry) does vertical blurring", () => {
+  assert.deepStrictEqual(blur2(copy2(unit), 0, 0).data, [
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+    0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0,
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0
+  ]);
+  assert.deepStrictEqual(blur2(copy2(unit), 0, 1).data.map(round), [
+    0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000,
+    0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000,
+    0.000, 0.000, 0.000, 0.000, 0.000, 0.037, 0.000, 0.000, 0.000, 0.000, 0.000,
+    0.000, 0.000, 0.000, 0.000, 0.000, 0.111, 0.000, 0.000, 0.000, 0.000, 0.000,
+    0.000, 0.000, 0.000, 0.000, 0.000, 0.222, 0.000, 0.000, 0.000, 0.000, 0.000,
+    0.000, 0.000, 0.000, 0.000, 0.000, 0.259, 0.000, 0.000, 0.000, 0.000, 0.000,
+    0.000, 0.000, 0.000, 0.000, 0.000, 0.222, 0.000, 0.000, 0.000, 0.000, 0.000,
+    0.000, 0.000, 0.000, 0.000, 0.000, 0.111, 0.000, 0.000, 0.000, 0.000, 0.000,
+    0.000, 0.000, 0.000, 0.000, 0.000, 0.037, 0.000, 0.000, 0.000, 0.000, 0.000,
+    0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000,
+    0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000
+  ]);
+  assert.deepStrictEqual(blur2(copy2(unit), 0, 2).data.map(round), [
+    0.000, 0.000, 0.000, 0.000, 0.000, 0.024, 0.000, 0.000, 0.000, 0.000, 0.000,
+    0.000, 0.000, 0.000, 0.000, 0.000, 0.048, 0.000, 0.000, 0.000, 0.000, 0.000,
+    0.000, 0.000, 0.000, 0.000, 0.000, 0.080, 0.000, 0.000, 0.000, 0.000, 0.000,
+    0.000, 0.000, 0.000, 0.000, 0.000, 0.120, 0.000, 0.000, 0.000, 0.000, 0.000,
+    0.000, 0.000, 0.000, 0.000, 0.000, 0.144, 0.000, 0.000, 0.000, 0.000, 0.000,
+    0.000, 0.000, 0.000, 0.000, 0.000, 0.152, 0.000, 0.000, 0.000, 0.000, 0.000,
+    0.000, 0.000, 0.000, 0.000, 0.000, 0.144, 0.000, 0.000, 0.000, 0.000, 0.000,
+    0.000, 0.000, 0.000, 0.000, 0.000, 0.120, 0.000, 0.000, 0.000, 0.000, 0.000,
+    0.000, 0.000, 0.000, 0.000, 0.000, 0.080, 0.000, 0.000, 0.000, 0.000, 0.000,
+    0.000, 0.000, 0.000, 0.000, 0.000, 0.048, 0.000, 0.000, 0.000, 0.000, 0.000,
+    0.000, 0.000, 0.000, 0.000, 0.000, 0.024, 0.000, 0.000, 0.000, 0.000, 0.000
+  ]);
+});
+
+function copy2({data, width, height}) {
+  return {data: data.slice(), width, height};
+}
+
+function round(x) {
+  return Math.round(x * 1000) / 1000 || 0;
+}