31:T1e0a,{ const width = 240 + 16 * (30 + 20); const height = 64 + 5 * 178; const ctx = DOM.context2d(width, height); yield ctx.canvas; ctx.canvas.style.maxWidth = "100%"; await visibility(); const deltas = [ [8, 128], [16, 64], [32, 32], [64, 16], [128, 8] ]; const graph = new PMJButterflyGraph(4); console.log({ graph }); function lerp(a, b, t) { return a * (1 - t) + b * t; } function easedLine(x0, y0, x1, y1, t0 = 1 / 16, t1 = 1 - t0) { const dx = x1 - x0; const dy = y1 - y0; ctx.beginPath(); ctx.moveTo(x0, y0); // Thanks for reminding me bezierCurveTo exists, Fabian :) ctx.bezierCurveTo(x0, lerp(y0, y1, t0), x1, lerp(y0, y1, t1), x1, y1); ctx.stroke(); } function drawBase() { ctx.fillStyle = "#FFF"; ctx.fillRect(0, 0, width, height); ctx.strokeStyle = "black"; ctx.lineWidth = 2; for (let i = 0; i < 5; i++) { const yi = 50 + i * 178; const x0 = 50; ctx.fillStyle = "#DDD"; ctx.fillRect(x0, yi, 128, 128); const xMax = 16 >> i; const yMax = 16 >> (4 - i); const [dx, dy] = deltas[i]; ctx.fillStyle = "#F00"; for (let x = 0; x < xMax; x++) { for (let y = 0; y < yMax; y++) { // tiles ctx.strokeRect(x0 + x * dx, yi + y * dy, dx, dy); const j = x + y * xMax; const xj = x0 + 160 + j * 50; const yn = yi + 64 - 15; // nodes if (graph.nodes[i * 16 + j]) { // fill occupied nodes ctx.fillRect(xj, yn, 30, 30); } ctx.strokeRect(xj, yn, 30, 30); // connections if (i < 4) { const jnl = (x >> 1) + y * xMax; const jnr = (x >> 1) + y * xMax + (xMax >> 1); const xjnl = x0 + 160 + jnl * 50 + 15; const xjnr = x0 + 160 + jnr * 50 + 15; easedLine(xj + 15, yn + 30, xjnl, yn + 178); easedLine(xj + 15, yn + 30, xjnr, yn + 178); } } } // points ctx.lineWidth = 2; ctx.strokeStyle = "#000"; ctx.fillStyle = "#FFF"; for (let [xp, yp] of graph.points) { xp = ((xp - 1) * 128) / 16; yp = ((yp - 1) * 128) / 16; ctx.fillRect(60 + xp, 10 + yp + yi, 4, 4); ctx.strokeRect(60 + xp, 10 + yp + yi, 4, 4); } } for (const [x, y] of graph.points) { drawTakenPaths(x*8 + 50, y*8 + 50); } } const splitDown = ([xNode, yNode]) => [ xNode >> 1, xNode >> 1, yNode << 1, (yNode << 1) + 1 ]; function splitPoints(points, split) { const nextPoints = []; for (const point of points) { let [xn1, xn2, yn1, yn2] = split(point); nextPoints.push([xn1, yn1], [xn2, yn2]); } return nextPoints; } function colorTiles(points, j) { const stride = 16 >> j; const [dx, dy] = deltas[j]; for (const [x, y] of points) { ctx.fillRect(51 + x * dx, 51 + j * 178 + y * dy, dx - 2, dy - 2); ctx.fillRect(210 + (x + y * stride) * 50, 100 + j * 178, 30, 30); } } function drawConnections(points, nextPoints, yl, yln, stride, nextStride) { for (let i = 0; i < points.length; i++) { let [x, y] = points[i]; let [xn1, yn1] = nextPoints[i * 2]; let [xn2, yn2] = nextPoints[i * 2 + 1]; x = (x + y * stride) * 50 + 225; xn1 = (xn1 + yn1 * nextStride) * 50 + 225; xn2 = (xn2 + yn2 * nextStride) * 50 + 225; easedLine(x, yl, xn1, yln); easedLine(x, yl, xn2, yln); } } function drawNodeHover(x, y) { const x0 = x - 210; if (x0 < 16 * 50 && x0 >= 0 && x0 % 50 < 30) { const y0 = y - 100; if (y0 % 178 >= 0 && y0 % 178 < 30) { // find which node we overlap with const i = (x0 / 50) | 0; const j = (y0 / 178) | 0; // node mouse cursor is on ctx.fillStyle = "#00BFFF"; ctx.fillRect(211 + i * 50, 100 + j * 178, 28, 28); // tile matching the current node const stride = 16 >> j; const xNode = i % stride; const yNode = (i / stride) | 0; const [dx, dy] = deltas[j]; ctx.fillRect( 51 + xNode * dx, 51 + j * 178 + yNode * dy, dx - 2, dy - 2 ); for (let l = 1, points = [[xNode, yNode]]; j + l < 5; l++) { ctx.fillStyle = `#00BFFF${(1 + (16 >> l)).toString(16)}${( 1 + (16 >> l) ).toString(16)}`; const nextPoints = splitPoints(points, splitDown); colorTiles(nextPoints, j + l); const stride = 16 >> (j - 1 + l); const nextStride = 16 >> (j + l); const yl = 130 + (j - 1 + l) * 178; const yln = 100 + (j + l) * 178; ctx.lineWidth = 4.75 - 0.75 * l; drawConnections(points, nextPoints, yl, yln, stride, nextStride); points = nextPoints; } } } } function drawTakenPaths(x, y) { if (x >= 50 && x < 178) { const x0 = x - 50; const y0 = (y - 50) % 178; if (y0 >= 0 && y0 < 128) { for (let i = 0; i < 5; i++) { // red tiles const stride = 16 >> i; const xi = x0 >> (3 + i); const yi = y0 >> (3 + 4 - i); const xn = (xi + yi * stride) * 50 + 210; const yn = 100 + i * 178; // connections ctx.lineWidth = 4; if (i < 4) { const xj = x0 >> (4 + i); const yj = y0 >> (6 - i); const xnj = (xj + yj * (8 >> i)) * 50 + 210; const ynj = yn + 178; easedLine(xn + 15, yn + 30, xnj + 15, ynj); } ctx.stroke(); } } } } yield ctx.canvas; while (true) { graph.clear(); for (let i = 0; i < 16; i++) { await visibility(); let x; { const offset = graph.size; const j = Math.floor(Math.random() * (graph.stride - offset)); x = graph.picks[offset + j]; graph.picks[j + offset] = graph.picks[offset]; graph.picks[offset] = x; } // draw top pick drawBase(); drawNodeHover(225 + x * 50, 115); await Promises.tick(300); await visibility(); // starting from position x, tumble down the graph // to our final point y. let y = 0; const { nodes, layers, stride } = graph; // set the top layer to occupied nodes[x] = true; for (let j = 1; j < layers; j++) { const layerOffset = j * stride; const layerStride = stride >> j; const xj = x >> j; y <<= 1; const yc1 = y * layerStride; const yc2 = (y + 1) * layerStride; const left = nodes[layerOffset + xj + yc1]; const right = nodes[layerOffset + xj + yc2]; // if left is occupied, right muse be unoccupied // therefore we can pick the right. If left is // unoccupied and right is also unoccupied it's //a 50/50 coin flip which side we pick let pick; if (left || (!right && Math.random() < 0.5)) { // pick right side and mark it occupied y += 1; pick = xj + yc2; nodes[layerOffset + xj + yc2] = true; } else { // pick left side and mark it occupied pick = xj + yc1; nodes[layerOffset + xj + yc1] = true; } drawBase(); drawNodeHover(225 + pick * 50, 115 + j * 178); await Promises.tick(200); await visibility(); } const point = [x, y]; graph.points.push(point); drawBase(); await Promises.tick(200); } await Promises.tick(5000); } }32:T735,This notebook builds on the work in: * [Progressive Multi-Jittered Sample Sequences](https://graphics.pixar.com/library/ProgressiveMultiJitteredSampling/) by Christensen et al. * [Efficient Generation of Points that Satisfy Two-Dimensional Elementary Intervals](http://jcgt.org/published/0008/01/04/) by Matt Pharr. This notebook explains how to create a simple data-structure that lets us generate a *"progressive multi-jittered sequence with strata along all base-2 elementary intervals"* (or `pmj02` sequence for short) of *n* points in ${tex`O(n \log{n})`}. This is an improvement on the method described in the second paper listed above, which is ${tex`O(n \sqrt{n})`}, but I don't work in graphics so perhaps this approach already well known. In the off chance that it is something novel I hope people will get some use out of it. ***Addendum:*** after publishing this notebook Andrew Helmer reached out to me on twitter to point out that he and his co-authors describe a method for generating these sequences in ${tex`O(n)`} in the supplementary material of [*Stochastic Generation of (t, s) Sample Sequences*](https://diglib.eg.org/handle/10.2312/sr20211287), with only the sample points as memory overhead. C++ implementation [here](https://github.com/Andrew-Helmer/stochastic-generation/blob/main/sampling/ssobol.cpp#L101). **If you wish to generate these sequences, use their method!** I think this notebook still is nice for developing some intuitive sense of how these sequences work though. ***Addendum 2:*** It turns out that I independently discovered a different formulation of the "Netshuffle" algorithm described by Abdalla G. M. Ahmed and Peter Wonka in [*Optimizing Dyadic Nets*](https://abdallagafar.com/publications/dyadic-nets/), which was published in August 2021, just four months prior to me sharing this notebook.33:T416,{ const ctx = DOM.context2d(5 * 148 + 100, 50 + 3 * 178); ctx.canvas.style.maxWidth = "100%"; ctx.fillStyle = "#DDD"; ctx.strokeStyle = "black"; ctx.lineWidth = 2; const deltas = [ [8, 128], [16, 64], [32, 32], [64, 16], [128, 8] ]; for (let i = 0; i < 5; i++) { const x = 50 + 148 * i; const y = 50; ctx.fillRect(x, y, 128, 128); ctx.fillRect(x, y + 178, 128, 128); ctx.fillRect(x, y + 2 * 178, 128, 128); const [dx, dy] = deltas[i]; for (let x0 = 0; x0 < 128; x0 += dx) { for (let y0 = 0; y0 < 128; y0 += dy) { ctx.strokeRect(x + x0, y + y0, dx, dy); ctx.strokeRect(x + x0, y + y0 + 178, dx, dy); ctx.strokeRect(x + x0, y + y0 + 2 * 178, dx, dy); } } } ctx.fillStyle = "#FF000088"; ctx.fillRect(198, 50, 32, 64); ctx.fillRect(346, 50, 32, 64); ctx.fillRect(346, 50 + 178, 64, 32); ctx.fillRect(494, 50 + 178, 64, 32); ctx.fillRect(494, 50 + 2 * 178, 128, 16); ctx.fillRect(642, 50 + 2 * 178, 128, 16); return ctx.canvas; }34:T1bfb,{ const width = 240 + 16 * (30 + 20); const height = 64 + 5 * 178; const dpr = window.devicePixelRatio || 1; const ctx = DOM.context2d(width, height, dpr); ctx.canvas.style.maxWidth = "100%"; const deltas = [ [8, 128], [16, 64], [32, 32], [64, 16], [128, 8] ]; function lerp(a, b, t) { return a * (1 - t) + b * t; } function easedLine(x0, y0, x1, y1, t0 = 1 / 16, t1 = 1 - t0) { const dx = x1 - x0; const dy = y1 - y0; ctx.beginPath(); ctx.moveTo(x0, y0); // Thanks for reminding me bezierCurveTo exists, Fabian :) ctx.bezierCurveTo(x0, lerp(y0, y1, t0), x1, lerp(y0, y1, t1), x1, y1); ctx.stroke(); } function drawBase() { ctx.fillStyle = "#FFF"; ctx.fillRect(0, 0, width, height); ctx.fillStyle = "#DDD"; ctx.strokeStyle = "black"; ctx.lineWidth = 2; for (let i = 0; i < 5; i++) { const yi = 50 + i * 178; const x0 = 50; ctx.fillRect(x0, yi, 128, 128); const xMax = 16 >> i; const yMax = 16 >> (4 - i); const [dx, dy] = deltas[i]; for (let x = 0; x < xMax; x++) { for (let y = 0; y < yMax; y++) { // tiles ctx.strokeRect(x0 + x * dx, yi + y * dy, dx, dy); const j = x + y * xMax; const xj = x0 + 160 + j * 50; const yn = yi + 64 - 15; // nodes ctx.strokeRect(xj, yn, 30, 30); // connections if (i < 4) { const jnl = (x >> 1) + y * xMax; const jnr = (x >> 1) + y * xMax + (xMax >> 1); const xjnl = x0 + 160 + jnl * 50 + 15; const xjnr = x0 + 160 + jnr * 50 + 15; easedLine(xj + 15, yn + 30, xjnl, yn + 178); easedLine(xj + 15, yn + 30, xjnr, yn + 178); } } } } } function drawTileHover(x, y) { if (x >= 50 && x < 178) { const x0 = x - 50; const y0 = (y - 50) % 178; if (y0 >= 0 && y0 < 128) { for (let i = 0; i < 5; i++) { // red tiles const stride = 16 >> i; const xi = x0 >> (3 + i); const yi = y0 >> (3 + 4 - i); const xn = (xi + yi * stride) * 50 + 210; const yn = 100 + i * 178; ctx.fillStyle = "#F00"; ctx.fillRect(xn, yn, 30, 30); const [dx, dy] = deltas[i]; const xt = xi * dx + 50; const yt = yi * dy + 50 + i * 178; ctx.fillRect(xt, yt, dx, dy); // connections ctx.lineWidth = 4; if (i < 4) { const xj = x0 >> (4 + i); const yj = y0 >> (6 - i); const xnj = (xj + yj * (8 >> i)) * 50 + 210; const ynj = yn + 178; easedLine(xn + 15, yn + 30, xnj + 15, ynj); } ctx.stroke(); // mouse cursor ctx.fillStyle = "#FFF"; ctx.strokeStyle = "#000"; ctx.lineWidth = 1; ctx.fillRect(x - 4, y0 - 5 + 50 + i * 178, 8, 8); ctx.strokeRect(x - 4, y0 - 5 + 50 + i * 178, 8, 8); } } } } const splitUp = ([xNode, yNode]) => [ xNode << 1, (xNode << 1) + 1, yNode >> 1, yNode >> 1 ]; const splitDown = ([xNode, yNode]) => [ xNode >> 1, xNode >> 1, yNode << 1, (yNode << 1) + 1 ]; function splitPoints(points, split) { const nextPoints = []; for (const point of points) { let [xn1, xn2, yn1, yn2] = split(point); nextPoints.push([xn1, yn1], [xn2, yn2]); } return nextPoints; } function colorTiles(points, j) { const stride = 16 >> j; const [dx, dy] = deltas[j]; for (const [x, y] of points) { ctx.fillRect(51 + x * dx, 51 + j * 178 + y * dy, dx - 2, dy - 2); ctx.fillRect(210 + (x + y * stride) * 50, 100 + j * 178, 30, 30); } } function drawConnections(points, nextPoints, yl, yln, stride, nextStride) { for (let i = 0; i < points.length; i++) { let [x, y] = points[i]; let [xn1, yn1] = nextPoints[i * 2]; let [xn2, yn2] = nextPoints[i * 2 + 1]; x = (x + y * stride) * 50 + 225; xn1 = (xn1 + yn1 * nextStride) * 50 + 225; xn2 = (xn2 + yn2 * nextStride) * 50 + 225; easedLine(x, yl, xn1, yln); easedLine(x, yl, xn2, yln); } } function drawNodeHover(x, y) { const x0 = x - 210; if (x0 < 16 * 50 && x0 >= 0 && x0 % 50 < 30) { const y0 = y - 100; if (y0 % 178 >= 0 && y0 % 178 < 30) { // find which node we overlap with const i = (x0 / 50) | 0; const j = (y0 / 178) | 0; // node mouse cursor is on ctx.fillStyle = "#F00"; ctx.fillRect(211 + i * 50, 100 + j * 178, 28, 28); // tile matching the current node const stride = 16 >> j; const xNode = i % stride; const yNode = (i / stride) | 0; const [dx, dy] = deltas[j]; ctx.fillRect( 51 + xNode * dx, 51 + j * 178 + yNode * dy, dx - 2, dy - 2 ); for (let l = 1, points = [[xNode, yNode]]; j - l >= 0; l++) { ctx.fillStyle = `#F00${(2 + (16 >> l)).toString(16)}`; const nextPoints = splitPoints(points, splitUp); colorTiles(nextPoints, j - l); const stride = 16 >> (j + 1 - l); const nextStride = 16 >> (j - l); const yl = 100 + (j + 1 - l) * 178; const yln = 130 + (j - l) * 178; ctx.lineWidth = 4.75 - 0.75 * l; drawConnections(points, nextPoints, yl, yln, stride, nextStride); points = nextPoints; } for (let l = 1, points = [[xNode, yNode]]; j + l < 5; l++) { ctx.fillStyle = `#F00${(2 + (16 >> l)).toString(16)}`; const nextPoints = splitPoints(points, splitDown); colorTiles(nextPoints, j + l); const stride = 16 >> (j - 1 + l); const nextStride = 16 >> (j + l); const yl = 130 + (j - 1 + l) * 178; const yln = 100 + (j + l) * 178; ctx.lineWidth = 4.75 - 0.75 * l; drawConnections(points, nextPoints, yl, yln, stride, nextStride); points = nextPoints; } } } } function hoverDraw(e) { const box = ctx.canvas.getBoundingClientRect(); const x = (((e.clientX - box.left) / box.width) * ctx.canvas.width) / dpr; const y = (((e.clientY - box.top) / box.height) * ctx.canvas.height) / dpr; drawBase(); drawTileHover(x, y); drawNodeHover(x, y); // mouse cursor ctx.fillStyle = "#FFF"; ctx.strokeStyle = "#000"; ctx.lineWidth = 1; ctx.fillRect(x - 4, y - 4, 8, 8); ctx.strokeRect(x - 4, y - 4, 8, 8); } ctx.canvas.addEventListener("mousemove", hoverDraw); // initiate a point to give a hint to the user that they can hover over the image // yield early so box is not zero yield html`

${ctx.canvas}

(hover or tap to change highlight)

`; const box = ctx.canvas.getBoundingClientRect(); hoverDraw({ clientX: (dpr * (275 + 5 * 50 + box.left) * box.width) / ctx.canvas.width, clientY: (dpr * (115 + 1 * 178 + box.top) * box.height) / ctx.canvas.height }); }35:T6be,Now assume that we have a graph that meets all the restrictions of a `pmj02` sequence. Then any tile is unoccupied if, and only if, at least one of its children is unoccupied, and at least one of its parents is unoccupied. This can be shown with a simple proof by contradiction: suppose a tile was unoccupied but both of its children were. In that case, both of its children must be overlapping with the tile's partner, which would break one of the restrictions of a `pmj02` sequence. Therefore this cannot be true. Due to the symmetry in our graph we can apply the same argument to show that given an unoccupied tile at least one of its parents must be unoccupied. And by induction this also immediately shows that *if* we pick an unoccupied tile, there *has* to be a path of unoccupied tiles from top to bottom going through it. Putting all of this together the algorithm becomes quite simple: we can pick *any* random unoccupied tile on *any* layer. Then we repeatedly pick an unoccupied parent until we reach the top layer. Starting from our original random tile we also pick a random unoccupied child until we reach the top and bottom layer. The tile at the top layer gives us our x-position, the tile at the bottom our y-position. Because the entire path is guaranteed to be unoccupied we also know that it must be a valid one. All that's left is now marking the tiles in this path as occupied, and we are ready to repeat the whole process. This is kind of similar to that of a [pachinko tree](https://observablehq.com/@jobleonard/pachinko-tree), but applied to our special graph structure. In the visualization below we simplify this by always picking an unoccupied tile in the top layer and moving down from there:36:T1ac2,{ const width = 240 + 16 * (30 + 20); const height = 64 + 5 * 178; const ctx = DOM.context2d(width, height); ctx.canvas.style.maxWidth = "100%"; yield ctx.canvas; await visibility(); const deltas = [ [8, 128], [16, 64], [32, 32], [64, 16], [128, 8] ]; const graph = new PMJButterflyGraph(4); console.log({ graph }); function lerp(a, b, t) { return a * (1 - t) + b * t; } function easedLine(x0, y0, x1, y1, t0 = 1 / 16, t1 = 1 - t0) { const dx = x1 - x0; const dy = y1 - y0; ctx.beginPath(); ctx.moveTo(x0, y0); // Thanks for reminding me bezierCurveTo exists, Fabian :) ctx.bezierCurveTo(x0, lerp(y0, y1, t0), x1, lerp(y0, y1, t1), x1, y1); ctx.stroke(); } function drawBase() { ctx.fillStyle = "#FFF"; ctx.fillRect(0, 0, width, height); ctx.strokeStyle = "black"; ctx.lineWidth = 2; for (let i = 0; i < 5; i++) { const yi = 50 + i * 178; const x0 = 50; ctx.fillStyle = "#DDD"; ctx.fillRect(x0, yi, 128, 128); const xMax = 16 >> i; const yMax = 16 >> (4 - i); const [dx, dy] = deltas[i]; ctx.fillStyle = "#F00"; for (let x = 0; x < xMax; x++) { for (let y = 0; y < yMax; y++) { // tiles ctx.strokeRect(x0 + x * dx, yi + y * dy, dx, dy); const j = x + y * xMax; const xj = x0 + 160 + j * 50; const yn = yi + 64 - 15; // nodes if (graph.nodes[i * 16 + j]) { // fill occupied nodes ctx.fillRect(xj, yn, 30, 30); } ctx.strokeRect(xj, yn, 30, 30); // connections if (i < 4) { const jnl = (x >> 1) + y * xMax; const jnr = (x >> 1) + y * xMax + (xMax >> 1); const xjnl = x0 + 160 + jnl * 50 + 15; const xjnr = x0 + 160 + jnr * 50 + 15; easedLine(xj + 15, yn + 30, xjnl, yn + 178); easedLine(xj + 15, yn + 30, xjnr, yn + 178); } } } // points ctx.lineWidth = 2; ctx.strokeStyle = "#000"; ctx.fillStyle = "#FFF"; for (let [xp, yp] of graph.points) { xp = ((xp - 1) * 128) / 16; yp = ((yp - 1) * 128) / 16; ctx.fillRect(60 + xp, 10 + yp + yi, 4, 4); ctx.strokeRect(60 + xp, 10 + yp + yi, 4, 4); } } } const splitDown = ([xNode, yNode]) => [ xNode >> 1, xNode >> 1, yNode << 1, (yNode << 1) + 1 ]; function splitPoints(points, split) { const nextPoints = []; for (const point of points) { let [xn1, xn2, yn1, yn2] = split(point); nextPoints.push([xn1, yn1], [xn2, yn2]); } return nextPoints; } function colorTiles(points, j) { const stride = 16 >> j; const [dx, dy] = deltas[j]; for (const [x, y] of points) { ctx.fillRect(51 + x * dx, 51 + j * 178 + y * dy, dx - 2, dy - 2); ctx.fillRect(210 + (x + y * stride) * 50, 100 + j * 178, 30, 30); } } function drawConnections(points, nextPoints, yl, yln, stride, nextStride) { for (let i = 0; i < points.length; i++) { let [x, y] = points[i]; let [xn1, yn1] = nextPoints[i * 2]; let [xn2, yn2] = nextPoints[i * 2 + 1]; x = (x + y * stride) * 50 + 225; xn1 = (xn1 + yn1 * nextStride) * 50 + 225; xn2 = (xn2 + yn2 * nextStride) * 50 + 225; easedLine(x, yl, xn1, yln); easedLine(x, yl, xn2, yln); } } function drawNodeHover(x, y) { const x0 = x - 210; if (x0 < 16 * 50 && x0 >= 0 && x0 % 50 < 30) { const y0 = y - 100; if (y0 % 178 >= 0 && y0 % 178 < 30) { // find which node we overlap with const i = (x0 / 50) | 0; const j = (y0 / 178) | 0; // node mouse cursor is on ctx.fillStyle = "#00BFFF"; ctx.fillRect(211 + i * 50, 100 + j * 178, 28, 28); // tile matching the current node const stride = 16 >> j; const xNode = i % stride; const yNode = (i / stride) | 0; const [dx, dy] = deltas[j]; ctx.fillRect( 51 + xNode * dx, 51 + j * 178 + yNode * dy, dx - 2, dy - 2 ); for (let l = 1, points = [[xNode, yNode]]; j + l < 5; l++) { ctx.fillStyle = `#00BFFF${(1 + (16 >> l)).toString(16)}${( 1 + (16 >> l) ).toString(16)}`; const nextPoints = splitPoints(points, splitDown); colorTiles(nextPoints, j + l); const stride = 16 >> (j - 1 + l); const nextStride = 16 >> (j + l); const yl = 130 + (j - 1 + l) * 178; const yln = 100 + (j + l) * 178; ctx.lineWidth = 4.75 - 0.75 * l; drawConnections(points, nextPoints, yl, yln, stride, nextStride); points = nextPoints; } } } } while (true) { graph.clear(); for (let i = 0; i < 16; i++) { await visibility(); let x; { const offset = graph.size; const j = Math.floor(Math.random() * (graph.stride - offset)); x = graph.picks[offset + j]; graph.picks[j + offset] = graph.picks[offset]; graph.picks[offset] = x; } // draw top pick drawBase(); drawNodeHover(225 + x * 50, 115); yield Promises.tick(1000, ctx.canvas); await visibility(); // starting from position x, tumble down the graph // to our final point y. let y = 0; const { nodes, layers, stride } = graph; // set the top layer to occupied nodes[x] = true; for (let j = 1; j < layers; j++) { const layerOffset = j * stride; const layerStride = stride >> j; const xj = x >> j; y <<= 1; const yc1 = y * layerStride; const yc2 = (y + 1) * layerStride; const left = nodes[layerOffset + xj + yc1]; const right = nodes[layerOffset + xj + yc2]; // if left is occupied, right muse be unoccupied // therefore we can pick the right. If left is // unoccupied and right is also unoccupied it's //a 50/50 coin flip which side we pick let pick; if (left || (!right && Math.random() < 0.5)) { // pick right side and mark it occupied y += 1; pick = xj + yc2; nodes[layerOffset + xj + yc2] = true; } else { // pick left side and mark it occupied pick = xj + yc1; nodes[layerOffset + xj + yc1] = true; } drawBase(); drawNodeHover(225 + pick * 50, 115 + j * 178); yield Promises.tick(1000, ctx.canvas); await visibility(); } const point = [x, y]; graph.points.push(point); drawBase(); yield Promises.tick(1000, ctx.canvas); } yield Promises.tick(5000, ctx.canvas); } }37:T571,## Future Work The approach described above can be optimized in a number of ways. First, it is trivial to represent the graph in sequential array form (see appendix below), and each node (tile) in the graph requires only one bit of information to store whether or not it is occupied, keeping the graph representation extremely compact. To pick a parent or child at random we only need a single bit of randomness. Finally, there is room for parallelization: after picking an unoccupied tile the ascent up the parent nodes to determine the x position and the descent down the children nodes to determine the y position are completely independent of each other. Of these optimizations only the sequential array form is used in the implementation below. One way to improve the noise quality a bit would be to use a separate pachinko tree (specifically, a quad-tree version) to pick tiles from middle tile layer (the one with the square tiles) and then work our way up and down. This would ensure that the coarse-grained structure of the order in which the tiles are picked is evenly distributed. I will probably do this in a future notebook once I worked out some other ideas regarding pachinko quadtrees (specifically, how to exploit their inherent structure for fast nearest neighbor search, as that might make them very practical for best-candidate based approaches to blue noise generation).38:Tdbb,viewof PMJButterflyGraph = { class PMJButterflyGraph { constructor(power) { this.layers = power + 1; this.stride = 2 ** power; // This could be optimized to be one bit per node, but // booleans are clearer for demonstration purposes. this.nodes = [false]; for (let i = 1; i < this.stride * this.layers; i++) { this.nodes[i] = false; } this.picks = [0]; for (let i = 1; i < this.stride; i++) this.picks[i] = i; this.points = []; } clear() { this.nodes.fill(false); this.points.length = 0; } get size() { return this.points.length; } fill() { this.clear(); while (this.points.length < this.stride) this.pick(); return this; } pick() { // we can't generate more points than the stride length if (this.points.length === this.stride) return; // to keep the code brief we just pick an unoccupied tile // at from the top layer and descend - in other words: // pick a random x-coordinate of those still available. // the values in in the range [0, i) have already been // picked, and in [i, this.picks.length) are still // available. After picking a point, we swap it to the // front. Note that we do not need to move the whole // array, since we do not care about the order of un- // picked values (only take values from it at random // positions). let x; { const offset = this.size; const j = Math.floor(Math.random() * (this.stride - offset)); x = this.picks[offset + j]; this.picks[j + offset] = this.picks[offset]; this.picks[offset] = x; } // starting from position x, tumble down the graph // to our final point y. let y = 0; const { nodes, layers, stride } = this; // set the top layer to occupied nodes[x] = true; for (let i = 1; i < layers; i++) { // see appendix for an explanatio nof how to // translate between (x,y) ranges and indices // of the node array. const layerOffset = i * stride; const layerStride = stride >> i; const xi = x >> i; y <<= 1; const yc1 = y * layerStride; const yc2 = (y + 1) * layerStride; const left = nodes[layerOffset + xi + yc1]; const right = nodes[layerOffset + xi + yc2]; // if left is occupied, right must be unoccupied // therefore we can pick the right. If left is // unoccupied and right is also unoccupied it's //a 50/50 coin flip which side we pick // left | right | out // F | F | random pick // F | T | F // T | F | T // T | T | invalid if (left || (!right && Math.random() < 0.5)) { // pick right side and mark it occupied y += 1; nodes[layerOffset + xi + yc2] = true; } else { // pick left side and mark it occupied nodes[layerOffset + xi + yc1] = true; } // Note that we could also do this branchlessly // using bitwise operators: // y += left ^ ( !(left|right) & randomBit ) // nodes[layerOffset + xi + y * layerStride] = true } const point = [x, y]; this.points.push(point); return point; } } const view = md`~~~javascript ${PMJButterflyGraph.toString()} ~~~`; view.value = PMJButterflyGraph; return view; }39:T6c8,{ const width = 240 + 16 * (30 + 20); const height = 64 + 5 * 178; const dpr = window.devicePixelRatio || 1; const ctx = DOM.context2d(width, height, dpr); ctx.canvas.style.maxWidth = "100%"; const deltas = [ [8, 128], [16, 64], [32, 32], [64, 16], [128, 8] ]; function lerp(a, b, t) { return a * (1 - t) + b * t; } function easedLine(x0, y0, x1, y1, t0 = 1 / 16, t1 = 1 - t0) { const dx = x1 - x0; const dy = y1 - y0; ctx.beginPath(); ctx.moveTo(x0, y0); // Thanks for reminding me bezierCurveTo exists, Fabian :) ctx.bezierCurveTo(x0, lerp(y0, y1, t0), x1, lerp(y0, y1, t1), x1, y1); ctx.stroke(); } function drawBase() { ctx.fillStyle = "#FFF"; ctx.fillRect(0, 0, width, height); ctx.fillStyle = "#DDD"; ctx.strokeStyle = "black"; ctx.lineWidth = 2; for (let i = 0; i < 5; i++) { const yi = 50 + i * 178; const x0 = 50; ctx.fillRect(x0, yi, 128, 128); const xMax = 16 >> i; const yMax = 16 >> (4 - i); const [dx, dy] = deltas[i]; for (let x = 0; x < xMax; x++) { for (let y = 0; y < yMax; y++) { // tiles ctx.strokeRect(x0 + x * dx, yi + y * dy, dx, dy); const j = x + y * xMax; const xj = x0 + 160 + j * 50; const yn = yi + 64 - 15; // nodes ctx.strokeRect(xj, yn, 30, 30); // connections if (i < 4) { const jnl = (x >> 1) + y * xMax; const jnr = (x >> 1) + y * xMax + (xMax >> 1); const xjnl = x0 + 160 + jnl * 50 + 15; const xjnr = x0 + 160 + jnr * 50 + 15; } } } } } drawBase(); return ctx.canvas; }3a:Tc81,{ const width = 240 + 16 * (30 + 20); const height = 64 + 5 * 178; const dpr = window.devicePixelRatio || 1; const ctx = DOM.context2d(width, height, dpr); ctx.font = "bold 16px sans"; ctx.canvas.style.maxWidth = "100%"; const deltas = [ [8, 128], [16, 64], [32, 32], [64, 16], [128, 8] ]; function drawBase() { ctx.fillStyle = "#FFF"; ctx.fillRect(0, 0, width, height); ctx.strokeStyle = "#000"; ctx.lineWidth = 2; for (let z = 0; z < 5; z++) { ctx.fillStyle = "#DDD"; const yz = 50 + z * 178; const x0 = 50; ctx.fillRect(x0, yz, 128, 128); const xMax = 16 >> z; const yMax = 16 >> (4 - z); const [dx, dy] = deltas[z]; for (let x = 0; x < xMax; x++) { for (let y = 0; y < yMax; y++) { // tiles ctx.strokeRect(x0 + x * dx, yz + y * dy, dx, dy); const j = x + y * xMax; const xj = x0 + 160 + j * 50; const yn = yz + 64 - 15; // nodes ctx.strokeRect(xj, yn, 30, 30); // connections if (z < 4) { const jnl = (x >> 1) + y * xMax; const jnr = (x >> 1) + y * xMax + (xMax >> 1); const xjnl = x0 + 160 + jnl * 50 + 15; const xjnr = x0 + 160 + jnr * 50 + 15; } } } ctx.fillStyle = "#000"; } } function drawTileHover(x, y) { if (x >= 50 && x < 178) { const x0 = x - 50; const y0 = (y - 50) % 178; if (y0 >= 0 && y0 < 128) { for (let i = 0; i < 5; i++) { // red tiles const stride = 16 >> i; const xi = x0 >> (3 + i); const yi = y0 >> (3 + 4 - i); const xn = (xi + yi * stride) * 50 + 210; const yn = 100 + i * 178; ctx.fillStyle = "#F00"; ctx.fillRect(xn, yn, 30, 30); const [dx, dy] = deltas[i]; const xt = xi * dx + 50; const yt = yi * dy + 50 + i * 178; ctx.fillRect(xt, yt, dx, dy); // mouse cursor ctx.fillStyle = "#FFF"; ctx.strokeStyle = "#000"; ctx.lineWidth = 1; ctx.fillRect(x - 4, y0 - 5 + 50 + i * 178, 8, 8); ctx.strokeRect(x - 4, y0 - 5 + 50 + i * 178, 8, 8); } } } } yield ctx.canvas; while (true) { for (let z = 0; z < 5; z++) { const [dx, dy] = deltas[z]; for (let y = 5; y < 128; y += dy) { for (let x = 4; x < 128; x += dx) { await visibility(); await Promises.tick(400); drawBase(); drawTileHover(50 + x, 50 + y); // fade out other rows than current one ctx.fillStyle = "#FFFB"; ctx.fillRect(0, 0, width, 40 + z * 178); ctx.fillRect(0, 40 + (z + 1) * 178, width, height); ctx.fillStyle = "#000"; const i = (x - 4) / dx, j = (y - 5) / dy; ctx.fillText( ` z = ${z}, ${ z === 4 ? `idx = j = ${j}`: !z ? `idx = i = ${i}` : `idx = i + j × ${16 >> z} = ${i} + ${j} × ${16 >> z} = ${i + j * (16>>z)}` }`, 205, 75 + 178 * z ); } } } } }3b:T5fb,That means that if we have a pair of column/row indices *(i, j)* on layer *z*, the column/row indices of the child nodes will given by: i_z+1 = i_z / 2 (same for both coordinates) j_z+1 = j_z × 2 and (j_z × 2) + 1 The column/row indices of the parents are similarly given by: i_z-1 = i_z × 2 and (i_z × 2) + 1 j_z-1 = j_z / 2 (same for both coordinates) Now let's look at this from the perspective of *(x, y)* coordinates in our grid. The grid itself needs to have side lengths equal to the number of samples. We can translate the above equations into the following optimized bits of code: idx = (x >> z) + ((y >> (zMax - z)) << (zMax - z)) Here `z` is the specific layer that we are looking at (starting from zero), and `zMax` is the bottom layer, which is log₂(samples), so 4 in our example. Since we only use integers and powers of two, we use bit-shifting operators for a bit of increased performance. That is per layer. If we compress that into one big flat array, we get: idx = (x >> z) + ((y >> (zMax - z)) << (zMax - z)) + (z << zMax)3c:Td54,{ const width = 240 + 16 * (30 + 20); const height = 64 + 5 * 178; const dpr = window.devicePixelRatio || 1; const ctx = DOM.context2d(width, height, dpr); ctx.font = "bold 16px sans"; ctx.canvas.style.maxWidth = "100%"; const deltas = [ [8, 128], [16, 64], [32, 32], [64, 16], [128, 8] ]; function drawBase() { ctx.fillStyle = "#FFF"; ctx.fillRect(0, 0, width, height); ctx.strokeStyle = "#000"; ctx.lineWidth = 2; for (let i = 0; i < 5; i++) { ctx.fillStyle = "#DDD"; const yi = 50 + i * 178; const x0 = 50; ctx.fillRect(x0, yi, 128, 128); const xMax = 16 >> i; const yMax = 16 >> (4 - i); const [dx, dy] = deltas[i]; for (let x = 0; x < xMax; x++) { for (let y = 0; y < yMax; y++) { // tiles ctx.strokeRect(x0 + x * dx, yi + y * dy, dx, dy); const j = x + y * xMax; const xj = x0 + 160 + j * 50; const yn = yi + 64 - 15; // nodes ctx.strokeRect(xj, yn, 30, 30); // connections if (i < 4) { const jnl = (x >> 1) + y * xMax; const jnr = (x >> 1) + y * xMax + (xMax >> 1); const xjnl = x0 + 160 + jnl * 50 + 15; const xjnr = x0 + 160 + jnr * 50 + 15; } } } } } function drawTileHover(x, y) { if (x >= 50 && x < 178) { const x0 = x - 50; const y0 = (y - 50) % 178; if (y0 >= 0 && y0 < 128) { for (let i = 0; i < 5; i++) { // red tiles const stride = 16 >> i; const xi = x0 >> (3 + i); const yi = y0 >> (3 + 4 - i); const xn = (xi + yi * stride) * 50 + 210; const yn = 100 + i * 178; ctx.fillStyle = "#F00"; ctx.fillRect(xn, yn, 30, 30); const [dx, dy] = deltas[i]; const xt = xi * dx + 50; const yt = yi * dy + 50 + i * 178; ctx.fillRect(xt, yt, dx, dy); // mouse cursor ctx.fillStyle = "#FFF"; ctx.strokeStyle = "#000"; ctx.lineWidth = 1; ctx.fillRect(x - 4, y0 - 5 + 50 + i * 178, 8, 8); ctx.strokeRect(x - 4, y0 - 5 + 50 + i * 178, 8, 8); } } } } function drawLabels(x, y) { ctx.fillStyle = "#000"; x >>= 3; y >>= 3; for (let i = 0; i < 5; i++) { ctx.fillText( ` z = ${i}, ${ i === 4 ? `idx = y = ${y}` : !i ? `idx = x = ${x}` : `idx = (x >> ${i}) + ((y >> ${4 - i}) << ${ 4 - i }) = (${x} >> ${i}) + ((${y} >> ${4 - i}) << ${4 - i}) = ${ (x >> i) + ((y >> (4 - i)) << (4 - i)) }` }`, 205, 75 + 178 * i ); } } yield ctx.canvas; while (true) { for (let i = 0; i < 5; i++) { const [dx, dy] = deltas[i]; for (let y = 5; y < 128; y += dy) { for (let x = 4; x < 128; x += dx) { await visibility(); await Promises.tick(1000); drawBase(); drawTileHover(50 + x, 50 + y); drawLabels(x, y); // fade out other rows than current one ctx.fillStyle = "#FFFB"; ctx.fillRect(0, 0, width, 40 + i * 178); ctx.fillRect(0, 40 + (i + 1) * 178, width, height); } } } } }3d:[] 2f:{"title":"Fast Progressive Multi-Jittered Sample Sequences","cells":[{"id":691,"value":"# Fast Progressive Multi-Jittered Sample Sequences","pinned":false,"hidden":false,"mode":"md"},{"id":812,"value":"$31","pinned":false,"hidden":false,"mode":"ojs"},{"id":0,"value":"$32","pinned":false,"hidden":false,"mode":"md"},{"id":106,"value":"viewof enableSequences = Inputs.toggle({\n label: \"Generate sequences\",\n value: false\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":27,"value":"sequences = {\n const graph = new PMJButterflyGraph(12);\n const div = 3;\n const width = graph.stride >> div;\n const height = graph.stride >> div;\n const ctx = DOM.context2d(width * 2 + 50, height);\n ctx.canvas.style.maxWidth = \"100%\";\n yield html`

${ctx.canvas}

Left: pmj02-sequence, right: white noise`;\n do {\n ctx.fillStyle = \"#FFF1\";\n ctx.fillRect(0, 0, width * 2 + 50, height);\n ctx.fillStyle = \"#0008\";\n graph.clear();\n graph.fill();\n const { random } = Math;\n for (const [x, y] of graph.points) {\n ctx.fillRect(x >> div, y >> div, 2, 2);\n ctx.fillRect(\n (width + 50 + random() * width) | 0,\n (random() * height) | 0,\n 2,\n 2\n );\n }\n await visibility();\n await Promises.delay(0);\n } while (enableSequences);\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":821,"value":"What is a `pm02` sequence exactly? It is easiest to explain via a minimal example. Imagine that we wish to generate 16 random points in the unit square:","pinned":false,"hidden":false,"mode":"md"},{"id":149,"value":"viewof random16 = {\n restart2;\n const ctx = DOM.context2d(1 * 148 + 100, 100 + 128);\n ctx.canvas.style.maxWidth = \"100%\";\n ctx.fillStyle = \"#DDD\";\n ctx.strokeStyle = \"#000\";\n ctx.lineWidth = 2;\n\n const x = 50;\n const y = 50;\n ctx.fillRect(x, y, 128, 128);\n ctx.strokeRect(x, y, 128, 128);\n ctx.fillStyle = \"#FFF\";\n ctx.strokeStyle = \"#000\";\n ctx.canvas.value = Array.from({ length: 16 }, () => [\n ((Math.random() * 124) | 0) / 128,\n ((Math.random() * 124) | 0) / 128\n ]);\n\n for (const [rx, ry] of ctx.canvas.value) {\n const dx = (rx * 128) | 0;\n const dy = (ry * 128) | 0;\n ctx.strokeRect(3 + x + dx, 3 + y + dy, 2, 2);\n ctx.fillRect(3 + x + dx, 3 + y + dy, 2, 2);\n }\n\n return ctx.canvas;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":146,"value":"However, just picking any sixteen random points won't do: we wish to restrict their placement such that no two points are in the same base-2 elementary interval. With \"elementary base-2 interval\" we really just mean all the grids we can create where the number of rows times the number of columns is the number of samples (all of which will be a power of two). In our example we have sixteen samples, and we can create five grids that satisfy this requirement, as seen below:","pinned":false,"hidden":false,"mode":"md"},{"id":135,"value":"{\n const ctx = DOM.context2d(5 * 148 + 100, 100 + 128);\n ctx.canvas.style.maxWidth = \"100%\";\n ctx.fillStyle = \"#DDD\";\n ctx.strokeStyle = \"black\";\n ctx.lineWidth = 2;\n const deltas = [\n [8, 128],\n [16, 64],\n [32, 32],\n [64, 16],\n [128, 8]\n ];\n for (let i = 0; i < 5; i++) {\n const x = 50 + 148 * i;\n const y = 50;\n ctx.fillRect(x, y, 128, 128);\n const [dx, dy] = deltas[i];\n for (let x0 = 0; x0 < 128; x0 += dx) {\n for (let y0 = 0; y0 < 128; y0 += dy) {\n ctx.strokeRect(x + x0, y + y0, dx, dy);\n }\n }\n }\n return ctx.canvas;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":158,"value":"From left to right, we have sixteen columns by one row, eight columns by two rows, four by four, two by eight, one by sixteen. Our goal is to place points at random with the restriction that if we were to place all these grids over the unit square, no two points are allowed to be in the same rectangle of any of the grids. Something like this:","pinned":false,"hidden":false,"mode":"md"},{"id":365,"value":"{\n restart3;\n const ctx = DOM.context2d(5 * 148 + 100, 100 + 128);\n ctx.canvas.style.maxWidth = \"100%\";\n ctx.lineWidth = 2;\n const deltas = [\n [8, 128],\n [16, 64],\n [32, 32],\n [64, 16],\n [128, 8]\n ];\n ctx.strokeStyle = \"#000\";\n\n const graph = new PMJButterflyGraph(4).fill();\n for (let i = 0; i < 5; i++) {\n const x = 50 + 148 * i;\n const y = 50;\n ctx.fillStyle = \"#DDD\";\n ctx.fillRect(x, y, 128, 128);\n const [dx, dy] = deltas[i];\n for (let x0 = 0; x0 < 128; x0 += dx) {\n for (let y0 = 0; y0 < 128; y0 += dy) {\n ctx.strokeRect(x + x0, y + y0, dx, dy);\n }\n }\n for (const [rx, ry] of graph.points) {\n const dx = (3 + rx * 8) | 0;\n const dy = (3 + ry * 8) | 0;\n ctx.fillStyle = \"#FFF\";\n ctx.strokeRect(x + dx, y + dy, 2, 2);\n ctx.fillRect(x + dx, y + dy, 2, 2);\n }\n }\n return ctx.canvas;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":457,"value":"viewof restart3 = html``","pinned":false,"hidden":false,"mode":"ojs"},{"id":376,"value":"It is almost certain that the sixteen points picked completely at random in our first example do not meet those requirements:","pinned":false,"hidden":false,"mode":"md"},{"id":330,"value":"{\n const ctx = DOM.context2d(5 * 148 + 100, 100 + 128);\n ctx.canvas.style.maxWidth = \"100%\";\n ctx.lineWidth = 2;\n const deltas = [\n [8, 128],\n [16, 64],\n [32, 32],\n [64, 16],\n [128, 8]\n ];\n ctx.strokeStyle = \"#000\";\n for (let i = 0; i < 5; i++) {\n const x = 50 + 148 * i;\n const y = 50;\n ctx.fillStyle = \"#DDD\";\n ctx.fillRect(x, y, 128, 128);\n const [dx, dy] = deltas[i];\n for (let x0 = 0; x0 < 128; x0 += dx) {\n for (let y0 = 0; y0 < 128; y0 += dy) {\n ctx.strokeRect(x + x0, y + y0, dx, dy);\n }\n }\n for (const [rx, ry] of random16) {\n const dx = (rx * 128) | 0;\n const dy = (ry * 128) | 0;\n ctx.fillStyle = \"#FFF\";\n ctx.strokeRect(3 + x + dx, 3 + y + dy, 2, 2);\n ctx.fillRect(3 + x + dx, 3 + y + dy, 2, 2);\n }\n }\n return ctx.canvas;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":346,"value":"viewof restart2 = html``","pinned":false,"hidden":false,"mode":"ojs"},{"id":311,"value":"For our example set of sixteen samples it would be doable to solve this by brute force (meaning: start with an empty grid, generate new points and keep trying until we find one that fits into an empty spot, repeating this for as long as it takes to find sixteen points that fit). However, if we were to increase our sample size this quickly becomes unbearably slow. For example, for 2²⁰ points (roughly a million), we would have to check 21 grids (a single row of million columns, two rows of half a million, etc, all the way to a single column of a million rows). Once the grids fill up the odds of randomly picking a spot that is empty in all 21 layers becomes vanishingly small. Luckily, a bit of maths can help us out.\n\nLook at the marked tiles below and notice how they perfectly overlap:","pinned":false,"hidden":false,"mode":"md"},{"id":226,"value":"{\n const ctx = DOM.context2d(5 * 148 + 100, 100 + 128);\n ctx.canvas.style.maxWidth = \"100%\";\n ctx.fillStyle = \"#DDD\";\n ctx.strokeStyle = \"black\";\n ctx.lineWidth = 2;\n const deltas = [\n [8, 128],\n [16, 64],\n [32, 32],\n [64, 16],\n [128, 8]\n ];\n for (let i = 0; i < 5; i++) {\n const x = 50 + 148 * i;\n const y = 50;\n ctx.fillRect(x, y, 128, 128);\n const [dx, dy] = deltas[i];\n for (let x0 = 0; x0 < 128; x0 += dx) {\n for (let y0 = 0; y0 < 128; y0 += dy) {\n ctx.strokeRect(x + x0, y + y0, dx, dy);\n }\n }\n }\n ctx.fillStyle = \"#FF000088\";\n ctx.fillRect(50, 50, 16, 128);\n ctx.fillRect(198, 50, 16, 128);\n return ctx.canvas;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":241,"value":"In fact, we can find overlapping pairs between each adjacent grid layer:","pinned":false,"hidden":false,"mode":"md"},{"id":171,"value":"$33","pinned":false,"hidden":false,"mode":"ojs"},{"id":527,"value":"Because of our geometric constraints, each tile has two overlapping tiles above it, and overlaps with two tiles below it (unless it is the top or bottom layer, for obvious reasons). We can think of this as a [graph structure](https://en.wikipedia.org/wiki/Graph_%28discrete_mathematics%29): each tile is a node in the graph, connected to the overlapping tiles directly above and below it. The two overlapping tiles above it are its parents, the two overlapping tiles below it its children.\n\nFurthermore, each tile will have a \"sibling\" with whom it shares the same two parents, and a \"partner\" with whom it will have the same two children. In the image below, the red tile is the sibling of the blue tile, and partnered with the green tile:","pinned":false,"hidden":false,"mode":"md"},{"id":1112,"value":"{\n const ctx = DOM.context2d(5 * 148 + 100, 228);\n ctx.canvas.style.maxWidth = \"100%\";\n ctx.fillStyle = \"#DDD\";\n for (let i = 0; i < 5; i++) {\n const x = 50 + 148 * i;\n const y = 50;\n ctx.fillRect(x, y, 128, 128);\n }\n\n ctx.fillStyle = \"#4444FF88\";\n ctx.fillRect(198, 50, 32, 64);\n ctx.fillRect(346, 82, 32, 32);\n ctx.fillStyle = \"#FF000088\";\n ctx.fillRect(346, 50, 32, 32);\n ctx.fillStyle = \"#00CC0088\";\n ctx.fillRect(378, 50, 32, 32);\n ctx.fillRect(494, 50, 64, 32);\n\n ctx.strokeStyle = \"black\";\n ctx.lineWidth = 2;\n const deltas = [\n [8, 128],\n [16, 64],\n [32, 32],\n [64, 16],\n [128, 8]\n ];\n for (let i = 0; i < 5; i++) {\n const x = 50 + 148 * i;\n const y = 50;\n const [dx, dy] = deltas[i];\n for (let x0 = 0; x0 < 128; x0 += dx) {\n for (let y0 = 0; y0 < 128; y0 += dy) {\n ctx.strokeRect(x + x0, y + y0, dx, dy);\n }\n }\n }\n return ctx.canvas;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1108,"value":"The resulting graph structure produced by this is effectively a [butterfly diagram](https://en.wikipedia.org/wiki/Butterfly_diagram):","pinned":false,"hidden":false,"mode":"md"},{"id":464,"value":"$34","pinned":false,"hidden":false,"mode":"ojs"},{"id":246,"value":"$35","pinned":false,"hidden":false,"mode":"md"},{"id":743,"value":"$36","pinned":false,"hidden":false,"mode":"ojs"},{"id":746,"value":"$37","pinned":false,"hidden":false,"mode":"md"},{"id":12,"value":"$38","pinned":false,"hidden":false,"mode":"ojs"},{"id":914,"value":"## Appendix: Array Representation Of A Butterfly Diagram\n\nSo I first published this notebook saying the *\"array representation is trivial\"*, and then I remembered how much I despise the sentence *\"x is trivial and left as an exercise to the reader\"*. This is my attempt to make up a little for repeating that sin myself.\n\nLet's take the previous example with `n` = 16 samples as our concrete starting point again. We saw before that this produces a stack of five grids (specifically: given `n` samples, where `n` needs to be a power of two, we end up with an `n` × `n` grid, and log₂(`n`) + 1 layers).\n\nThe easiest way to represent this stack in an array would be to treat each layer as a row of sixteen elements, and arrange these rows from top to bottom in one big array:","pinned":false,"hidden":false,"mode":"md"},{"id":970,"value":"$39","pinned":false,"hidden":false,"mode":"ojs"},{"id":975,"value":"The next question then becomes how to map each tile in a layer to their row. The simplest way to do that is *idx = i + j × stride*, where *i* the column index, *j* the row index, and *stride* the number of columns in that layer:","pinned":false,"hidden":false,"mode":"md"},{"id":986,"value":"$3a","pinned":false,"hidden":false,"mode":"ojs"},{"id":1030,"value":"$3b","pinned":false,"hidden":false,"mode":"md"},{"id":1024,"value":"$3c","pinned":false,"hidden":false,"mode":"ojs"}],"files":"$W3d","stdlib":"1","titling":"auto"}