2c:T7ee,{ const height = Math.min(576, width); const ctx = DOM.context2d(width, height); ctx.canvas.style.maxWidth = "100%"; const newSet = htl.html``; yield htl.html`${ctx.canvas}
${newSet}`; const colorramp = makeRamp(64, "#warm"); const tilesize = 64; const log = Math.log2(tilesize); const valuePositions = Array.from({length: tilesize}, () => new Uint16Array(tilesize*2)); const posIndices = new Uint16Array(tilesize); const total = tilesize*tilesize; const coprime = findCoprimeApprox(tilesize, 2**0.5) const draw = async () => { ctx.fillStyle = "#DDDC"; ctx.fillRect(32, 32, tilesize*8, tilesize*8); ctx.lineWidth = 1; ctx.strokeStyle = "#000"; ctx.strokeRect(31.5, 31.5, tilesize*8 + 1, tilesize*8 + 1); await Promises.delay(20); const {points, paletteOrder} = (await makeTile({ tilesize, graphStructure: "random", shuffleColors: "random", colorscale: "warm", coprimeFactor: ((1 + 5**0.5)/2), graphCoPrimeFactor: 2**0.5, showTile: "H-curve pinwheel reverse", visibility })).value; for (let y = 0; y < tilesize; y++) { let row = y * tilesize; for (let x = 0; x < tilesize; x++) { const v = points[row++]; const positions = valuePositions[v]; let posIdx = posIndices[v]; positions[posIdx++] = x; positions[posIdx++] = y; posIndices[v] = posIdx; } } for (let n = 0; n < tilesize; n++) { // const idx = ( (j * coprime) % tilesize) * 2; const positions = valuePositions[(n*coprime)%tilesize]; ctx.fillStyle = colorramp[paletteOrder[n]]; for (let j = 0; j < tilesize; j++) { const k = positions[j*2]; const l = positions[j*2+1]; const x = 32; const y = 32; ctx.fillRect(x + k * 8, y + l * 8, 8, 8); } await Promises.delay(20); await visibility(); } }; draw(); newSet.addEventListener("click", draw); }2d:T630,{ const w = 5 * 148 + 100; const h = 60 + 128; const ctx = DOM.context2d(w, h); ctx.canvas.style.maxWidth = "100%"; const newSet = htl.html``; const fold = htl.html``; yield htl.html`${fold}
${ctx.canvas}
${newSet}`; const deltas = [ [8, 128], [16, 64], [32, 32], [64, 16], [128, 8] ]; const graph = new PMJButterflyGraph(4).fill(0); const draw = offset => { ctx.fillStyle = "#FFF"; ctx.fillRect(0, 0, w, h); for (let i = 0; i < 5; i++) { const x = 50 + (148 * offset * i); const y = 30; const [dx, dy] = deltas[i]; ctx.fillStyle = "#E0E0E080"; ctx.fillRect(x, y, 128, 128); ctx.fillStyle = "#333"; ctx.strokeStyle = "#FFF"; ctx.lineWidth = 2.5; for (const [rx, ry] of graph.points) { const dx = (2 + rx * 8) | 0; const dy = (2 + ry * 8) | 0; ctx.strokeRect(x + dx - 0.5, y + dy -0.5, 6, 6); ctx.fillRect(x + dx, y + dy, 5, 5); } ctx.strokeStyle = "#000"; ctx.lineWidth = 1; for (let x0 = 0; x0 < 128; x0 += dx) { for (let y0 = 0; y0 < 128; y0 += dy) { ctx.strokeRect(x + x0 + 0.5, y + y0 + 0.5, dx, dy); } } } }; draw(1); let foldValue = 1; fold.addEventListener("input", e => { if (fold.value !== foldValue) draw(foldValue = fold.value); }); newSet.addEventListener("click", () => { graph.fill(0); draw(foldValue); }); }2e:T4fe,This is an example of a so-called "progressive multi-jittered sequence" (`pmj02`). In a [previous notebook](https://observablehq.com/@jobleonard/fast-pmj02-sequences) I described a way to create these `pmj02` sequences in `O(nlogn)` time. Later I found out that I had 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/), published in August 2021, four months prior to me publishing my notebook. Of course, they are domain experts, and I am a random programmer who was just messing around with some interesting ideas for a month or two, so luckily I didn't have to feel the academic sting of being scooped and being unable to publish a paper. I admit I was a tiny bit disappointed anyway, since I had gotten really excited about the idea of being the first to discover this algorithm. But a while ago I realized something kinda neat, and I think Ahmed and Wonka might have missed that. By using the symmetries and the graph structure inherent to the problem, we can relatively easily create a special *set* of `pmj02` sequences, where none of the sequences overlap, and where the set as a whole covers the entire grid:2f:Tac4,{ const w = 5 * 148 + 100; const h = 60 + 128; const ctx = DOM.context2d(w, h); ctx.canvas.style.maxWidth = "100%"; const newSet = htl.html``; const fold = htl.html``; yield htl.html`${fold}
${ctx.canvas}
${newSet}`; const colorramp = makeRamp(16, "#warm").map(c => c + "CC"); const tilesize = 16; const log = Math.log2(tilesize); const deltas = []; const upscale = 8; for (let i = tilesize, j = 1; i > 0; i = i >>> 1) { deltas.push([i * upscale, j * upscale]); j *= 2 } const valuePositions = Array.from({length: tilesize}, () => new Uint16Array(tilesize*2)); const posIndices = new Uint16Array(tilesize); let paletteOrder; const generateTile = async () => { const tile = (await makeTile({ tilesize, graphStructure: "random", shuffleColors: "random", colorscale: "warm", coprimeFactor: ((1 + 5**0.5)/2), graphCoPrimeFactor: 2**0.5, showTile: "H-curve pinwheel reverse", visibility })).value; paletteOrder = tile.paletteOrder; for (let y = 0; y < tilesize; y++) { let row = y * tilesize; for (let x = 0; x < tilesize; x++) { const v = tile.points[row++]; const positions = valuePositions[v]; let posIdx = posIndices[v]; positions[posIdx++] = x; positions[posIdx++] = y; posIndices[v] = posIdx; } } }; await generateTile(); const draw = offset => { ctx.fillStyle = "#FFF"; ctx.fillRect(0, 0, w, h); // ctx.fillStyle = "#EEE8"; // for (let i = 0; i < 5; i++) { // const x = 50 + 148 * offset * i; // const y = 30; // ctx.fillRect(x, y, 128, 128); // } for (let i = 0; i < 5; i++) { const x = 50 + 148 * offset * i; const y = 30; ctx.strokeStyle = "#0008"; ctx.lineWidth = 1.5; ctx.strokeRect(x, y, 128, 128) for (let n = 0; n < tilesize; n++) { const positions = valuePositions[n]; ctx.fillStyle = colorramp[paletteOrder[n]]; for (let k = 0; k < tilesize; k++) { const j = positions[(k%tilesize)*2]; ctx.fillRect(x + 1 + j * 8, y + 1 + k * 8, 7, 7); } } 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 + 0.5, y + y0 + 0.5, dx, dy); } } } }; draw(1); let foldValue = 1; fold.addEventListener("input", e => { if (fold.value !== foldValue) draw(foldValue = fold.value); }); newSet.addEventListener("click", async () => { await generateTile(); draw(foldValue); }); }30:T9b4,{ const text1 = `

Aside: why do I think they missed it?

Ahmed and Wonka provided a set of graphs for their paper that can be used to create point sets with a good blue noise distribution. If I use those graphs, then the first sequence indeed has a nice even distribution associated with blue noise. However, the others sequences are more unevenly distributed again, and look just like "regular" pmj02 sequences:

`, text2 = `

The blue dots are the points generated by the blue noise graph provided by A&W, the red dots the other sets generated by the same graph

Having said that, I don't even know if it's really possible to create a tile where all sequences come reasonably close to having this blue noise property.

`; yield html`${text1}${text2}`; let showTile = 0; const tilesize = 256; const tile1 = await makeTile({ tilesize, graphStructure: "blue", showTile: "row", visibility }); const tile2 = await makeTile({ tilesize, graphStructure: "blue", showTile: "row", visibility }); const {ctx: ctx1, imgData, points, paletteOrder} = tile1.value; const {ctx: ctx2} = tile2.value; const valuePositions = Array.from({length: tilesize}, () => new Uint16Array(tilesize*2)); const posIndices = new Uint16Array(tilesize); for (let y = 0; y < tilesize; y++) { let row = y * tilesize; for (let x = 0; x < tilesize; x++) { const v = points[row++]; const positions = valuePositions[v]; let posIdx = posIndices[v]; positions[posIdx++] = x; positions[posIdx++] = y; posIndices[v] = posIdx; } } yield html`${text1}${tile1}${tile2}${text2}`; let j = 0; const coprime = findCoprimeApprox(tilesize); const render = (ctx, positions, fill1, fill2) => { ctx.fillStyle = "#000"; ctx.fillRect(0, 0, tilesize, tilesize); ctx.fillStyle = fill1; for (let i = 0; i < tilesize; i++) { ctx.fillRect(positions[i<<1]-1, positions[(i<<1)+1]-1, 3, 3); } ctx.fillStyle = fill2; for (let i = 0; i < tilesize; i++) { ctx.fillRect(positions[i<<1], positions[(i<<1)+1], 1, 1); } }; // blue noise points render(ctx1, valuePositions[0], "#BBBBEE80", "#BBBBEE") // all other points while (true) { j = (j + coprime) % tilesize; render(ctx2, valuePositions[j], "#FFFF0080", "#FFFF0FF"); await visibility(); await Promises.delay(1000); } }31:Tc81,If you would put the grid to its left on top of it, the two blue tiles would overlap with the red tile. Let's call those tiles the "parent tiles". The red tile also has a "sibling" that is also covered by the same parents, which is marked blue as well (for the colorblind: the sibling and parent tiles are marked with a hollow rectangle, as is the red tile). If you would put the red tile's grid on top of the grid on its right, the red tile would overlap with the two green tiles. We can call those its "children". There also is a "partner" that covers the same two children (for the colorblind: the partner and children are marked with a filled rectangle, as is the red tile). This is true for every tile in every grid: it always has two parents and a sibling (except the leftmost grid), and a partner and two children (except the rightmost grid). If you want to verify this, the image is interactive and you can hover to highlight another tile. Now remember that the rule for `pmj02` sequences was that each tile must overlap with *exactly* one point. With that in mind, let's look at the illustration again. If a point overlaps with the red tile, it *must* overlap with one of the parent tiles (we are assuming there are no points on edges between tiles, so a discrete grid). Now we want to place another point in the sibling tile. We are not allowed to break our rule of one-point-per-tile. So that means we *have* to place it in such a way that it overlaps with the other parent. The same is true for the red tile's partner and children: the point in the red tile will overlap with one of the two children, and any valid point overlapping the green partner tile *must* overlap with the other child. Imagine starting with an empty grid, and placing a first point, marking one tile on each layer as occupied. What we've just discovered via these exclusion principles is that doing so also halves the areas where a point is allowed to be placed for the sibling and partner tiles on each grid. It's basically like a geometric version of sudoku! The reverse must also be true: if a tile still has some area left that does not overlap with another occupied tile, then this must also be true for at least one of its parents, and at least one of its children. The parent and children areas must also partially overlap. And this repeats: it must also be true for at least one of the parent's parents, and one of that child's children, and so on. This means that there has to be an area representing an unoccupied path down the entire "family tree". This leads to the core insight that allows us generate these sequences efficiently: **pick any unoccupied node anywhere in the graph, and there will be at least one possible path of connected unoccupied nodes (and therefore overlapping tiles) up to the top layer, and down to the bottom layer. And since the node at top layer determines the exact x-position, and the node at the bottom layer the exact y-position, that path must therefore determine a valid point in our sequence!** This is easier to see if we treat all tiles as nodes, and connect all parents to their children with vertices. The result ends up looking like a butterfly diagram:32:T1b8b,{ 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 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]; for (let x = 0; x < xMax; x++) { for (let y = 0; y < yMax; y++) { const j = x + y * xMax; const xj = x0 + 160 + j * 50; const yn = yi + 64 - 15; // parent pair highlights if (!(x&1) && i < 4) { ctx.fillStyle = "#D0D0E8"; ctx.fillRect(xj - 4, yn - 9, 90, 46); } // fill occupied nodes ctx.strokeRect(xj, yn, 30, 30); // tiles ctx.strokeRect(x0 + x * dx, yi + y * dy, dx, dy); // nodes ctx.fillStyle = "#FFF"; 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(ctx, xj + 15, yn + 30, xjnl, yn + 178); easedLine(ctx, 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(ctx, 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(ctx, x, yl, xn1, yln); easedLine(ctx, 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 ctx.canvas; 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 }); }33:T1d10,{ 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); 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 // parent pair highlights if (!(x&1) && i < 4) { ctx.fillStyle = "#D0D0E8"; ctx.fillRect(xj - 4, yn - 9, 90, 46); } // fill occupied nodes ctx.fillStyle = graph.nodes[i * 16 + j] ? "#F00" : "#FFF"; 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(ctx, xj + 15, yn + 30, xjnl, yn + 178); easedLine(ctx, 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(ctx, x, yl, xn1, yln); easedLine(ctx, 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) { ctx.strokeStyle = "#F44"; 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 = 5; 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(ctx, xn + 15, yn + 30, xnj + 15, ynj); } ctx.stroke(); } ctx.strokeStyle = "#000"; } } } yield ctx.canvas; 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(400/((i>>2)+1)); 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((300 - 20*j)/((i>>2)+1)); await visibility(); } const point = [x, y]; graph.points.push(point); drawBase(); await Promises.tick(200); } restart1; }34:T5dd,This was the approach I discovered in my previous notebook, which as mentioned turned out to be more or less the same as the netshuffle algorithm. ## How To Build These Sets, Step Two: The Netshuffle Algorithm I said that the netshuffle algorithm is essentially the same algorithm, but there is one difference: it represents the graph as a *number*, and can do the whole process of tumbling down from x to y with some bitshifting and bitmasking inside a loop. How do we do that? Well, remember how we saw earlier that whenever a path connects one parent to one of their children, the other parent must connect with the other? This is a consequence of the fact that we only pick empty paths: when one parent has already picked a path connecting with one child, then the other parent *must* pick the other child. This symmetry onto the graph that lets us encode it as a binary string. Let's use a simplified example as our base: we make a sequence of points where the left parent always picks the leftmost child, and the right parent always picks the rightmost child. Remember that the parent pairs are highlighted light blue to make it easier to spot them. Furthermore, and I will explain why in a bit, we'll label each node and each connecting line. Each node will be labeled by its row index in binary (so 0 to 15, but written as 0000 to 1111). Each outgoing connection to the left child will be labeled as a `0`, and each outgoing connection the right child as `1`. The result looks like this:35:T16bb,{ 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); 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 = "#FFF"; 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; if (i < 4) { // parent pair highlights if (!(x&1)) { ctx.fillStyle = "#D0D0E8"; ctx.fillRect(xj - 4, yn - 9, 90, 46); ctx.fillStyle = "#FFF"; } // connections 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(ctx, xj + 15, yn + 30, xjnl, yn + 178); easedLine(ctx, xj + 15, yn + 30, xjnr, yn + 178); } // nodes if (graph.nodes[i * 16 + j]) { // fill occupied nodes ctx.fillStyle = "#F00"; ctx.fillRect(xj, yn, 30, 30); } ctx.fillRect(xj, yn, 30, 30); ctx.strokeRect(xj, yn, 30, 30); } } // 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); } // draw bits ctx.fillStyle = "#000"; ctx.strokeStyle = "#FFF"; ctx.lineWidth = 4; for (let i = 0; i < 5; i++) { const yi = 50 + i * 178; const x0 = 50; const xMax = 16 >> i; const yMax = 16 >> (4 - i); const [dx, dy] = deltas[i]; for (let y = 0; y < yMax; y++) { for (let x = 0;x < xMax; x++) { // tiles const j = x + y * xMax; const xj = x0 + 160 + j * 50; const yn = yi + 64 - 15; const k = x + y * xMax // bit labels ctx.font = "12px monospace"; ctx.fillStyle = "#000"; ctx.strokeText("1", xj + 21, yn + 33); ctx.strokeText("0", xj + 4, yn + 33); ctx.fillStyle = "#000"; ctx.fillText("1", xj + 21, yn + 33); ctx.fillText("0", xj + 4, yn + 33); const gray = k.toString(2).padStart(4, "0"); const black = (k >>(4-i)).toString(2).padStart(i, "0"); ctx.strokeText(gray, xj + 3, yn + 5) ctx.fillStyle = "#888"; ctx.fillText(gray, xj + 3, yn + 5); if (i) { ctx.fillStyle = "#000"; ctx.fillText(black, xj + 3, yn + 5); } if (i < 4) { ctx.font = "bold 12px monospace"; ctx.fillText((k&1) ? " 1" : " 0", xj + 3, yn + 5); } } } } } 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(ctx, x, yl, xn1, yln); easedLine(ctx, x, yl, xn2, yln); } } function drawTakenPaths(x, y) { if (x >= 50 && x < 178) { const x0 = x - 50; const y0 = (y - 50) % 178; if (y0 >= 0 && y0 < 128) { ctx.strokeStyle = "#F44"; 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 = 5; 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(ctx, xn + 15, yn + 30, xnj + 15, ynj); } ctx.stroke(); } ctx.strokeStyle = "#000"; } } } graph.clear(); for (let i = 0; i < 16; i++) { graph.points.push([i, reverseNbits(i%16, 4)]); } drawBase(); }36:Tc98,Except for the bit string labels on the bottom layer, each of them have their least significant bit (last bit on the right) bolded. Also, some parts of the labels are in gray. Why? Well, let's explain the least significant bit first. Each parent pair has a left and a right parent. Do you notice that the left parent always has a 0-bit and the right parent a 1-bit? That's just a consequence of odd and even indices. But now try picking any `x`-position (so any node at the top), and trace its red path down to the connected `y`-position at the bottom. E.g. if we look at the x-node for `0101`, the path down goes through these nodes:

If you look at the vertices in the illustration above, then you'll see that the last bit of the node label predicts where the outgoing red path! Odd node indices go to the right, even indices to the left. There are more patterns hiding in there. If you string the last bits of each of the first four labels together, then those four bits are the same as the bitstring of the bottom node! Also, the italiced parts of the node label match the path bit string created so far. And finally: each `x` connects to a `y` that is its bit string in reverse, and this matches the bit string created by the taken path! This is true for our `pmj02` sequence that follows the rule that every left (or even) parent connects with the left child, and ever right (or odd) parent connects with the right child. What about all the other sequences? Well, via our earlier "mutual exclusion" principle, we can swap left and right for an individual parent pair and the result is still a valid sequence.. Putting all these insights together lets us encode the generation of these sequences in a different way. We take the above "base" graph where the left parent always turns left, and the right parent always turns right. Then for each parent pair we have a random bit that determines if that parent pair flips this output. We can then do this "bit reversal" trick to the x position to get the y position, but just have to weave in these random flips where we turn a different direction. This ~~simplifies~~ shortens the code to this: ~~~ javascript const totalPoints = 16, // can be any power of two tpLog = Math.log2(totalPoints), points = [], randomGraph = [], mask = totalPoints - 1, offsetMask = mask << tpLog; for (let i = 0; i < totalPoints / 2 * tpLog; i++) { randomGraph[i] = Math.random() < 0.5 ? 0 : 1; } for (let x = 0; x < totalPoints; x++) { let y = x; for (let row = 0; row < tpLog; row++) { const odd = y & 1, yMod = y & (mask >> row), yOffset = y & (offsetMask >> row), // flip direction (and therefore y-bit?) randomBit = randomGraph[row * totalPoints / 2 + (y >> 1)], turnRight = odd ^ randomBit; y = yOffset + (yMod >> 1) + (turnRight ? ((totalPoints / 2) >> row) : 0); } points.push([x, y]); } ~~~ I'll admit that it's quite dense, but don't worry too much about it. Look at the graph below! I added the "flip" bits between them so you can see it in practice:37:T1aae,viewof parentpairgraph = { const points = [] const totalPoints = 16; const tpLog = Math.log2(totalPoints) const randomGraph = new Uint8Array(totalPoints / 2 * tpLog); const randomize = () => { points.length = 0; for (let i = 0; i < randomGraph.length; i++) randomGraph[i] = Math.random() < 0.5 ? 0 : 1; const mask = totalPoints - 1; const offsetMask = mask << tpLog; for (let x = 0; x < totalPoints; x++) { let y = x; for (let row = 0; row < tpLog; row++) { const odd = y & 1; const yMod = y & (mask >> row); const yOffset = y & (offsetMask >> row); const randomBit = randomGraph[row * 8 + (y >> 1)]; const turnRight = odd ^ randomBit; y = yOffset + (yMod >> 1) + (turnRight ? (8 >> row) : 0); } points.push([x, y]); } } // Drawing part const width = 240 + 16 * (30 + 20); const height = 64 + 5 * 178; const ctx = DOM.context2d(width, height); 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 = "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 = "#FFF"; for (let x = 0; x < xMax; x++) { for (let y = 0; y < yMax; y++) { const j = x + y * xMax; const xj = x0 + 160 + j * 50; const yn = yi + 64 - 15; 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; if (i < 4) { // parent pair highlights if (!(x&1)) { ctx.fillStyle = "#D0D0E8"; ctx.fillRect(xj - 4, yn - 9, 90, 46); ctx.fillStyle = "#FFF"; } // connections easedLine(ctx, xj + 15, yn + 30, xjnl, yn + 178); easedLine(ctx, xj + 15, yn + 30, xjnr, yn + 178); } ctx.fillRect(xj, yn, 30, 30); ctx.strokeRect(xj, yn, 30, 30); // tiles ctx.strokeRect(x0 + x * dx, yi + y * dy, dx, dy); } } // points ctx.lineWidth = 2; ctx.strokeStyle = "#000"; ctx.fillStyle = "#FFF"; for (let [xp, yp] of 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 points) { drawTakenPaths(x*8 + 50, y*8 + 50); } // draw bits ctx.fillStyle = "#000"; ctx.strokeStyle = "#FFF"; ctx.lineWidth = 4; for (let i = 0; i < 5; i++) { const yi = 50 + i * 178; const x0 = 50; const xMax = 16 >> i; const yMax = 16 >> (4 - i); const [dx, dy] = deltas[i]; for (let y = 0; y < yMax; y++) { for (let x = 0;x < xMax; x++) { // tiles const j = x + y * xMax; const xj = x0 + 160 + j * 50; const yn = yi + 64 - 15; const k = x + y * xMax // bit labels ctx.font = "12px monospace"; ctx.fillStyle = "#000"; ctx.strokeText("1", xj + 21, yn + 33); ctx.strokeText("0", xj + 4, yn + 33); ctx.fillText("1", xj + 21, yn + 33); ctx.fillText("0", xj + 4, yn + 33); if (i < 4 && (k&1)) { ctx.strokeText(randomGraph[(k>>>1) + i * 8], xj - 13, yn + 20); ctx.fillText(randomGraph[(k>>>1) + i * 8], xj - 13, yn + 20); } const gray = k.toString(2).padStart(4, "0"); const black = (k >>(4-i)).toString(2).padStart(i, "0"); ctx.strokeText(gray, xj + 3, yn + 5) ctx.fillStyle = "#888"; ctx.fillText(gray, xj + 3, yn + 5); if (i) { ctx.fillStyle = "#000"; ctx.fillText(black, xj + 3, yn + 5); } if (i < 4) { ctx.font = "bold 12px monospace"; ctx.fillText((k&1) ? " 1" : " 0", xj + 3, yn + 5); } } } } } 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(ctx, x, yl, xn1, yln); easedLine(ctx, x, yl, xn2, yln); } } function drawTakenPaths(x, y) { if (x >= 50 && x < 178) { const x0 = x - 50; const y0 = (y - 50) % 178; if (y0 >= 0 && y0 < 128) { ctx.strokeStyle = "#F44"; 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 = 5; 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(ctx, xn + 15, yn + 30, xnj + 15, ynj); } ctx.stroke(); } ctx.strokeStyle = "#000"; } } } /** for (let i = 0; i < tilesize; i++) { let j = 0; for (let k = 0; k < log; k++) { const randomBit = graph[k*tilesize/2 + (j>>1)]; j = (j <<1) + (((i>>k)^randomBit^(n>>(k)))&1) } ctx.fillRect(i, j, 1, 1); } */ ctx.canvas.value = () => { randomize(); drawBase(); } ctx.canvas.style.maxWidth = "100%"; yield ctx.canvas; await visibility(); randomize(); drawBase(); }38:T4b8,This formulation should essentially be equivalent to the netshuffle algorithm by Ahmed and Wonka. Modeling the problem this way gives us some neat powers too. For example, we have just determined all possible `pmj02` sequences of any given tile size! Because it has to be equivalent to all the possible combinations of whether each parent pair is swapped or not. This is equal to stringing together the total number of bits in the "random" graph as a binary number. So in the above 16 by 16 grid we have 32 bits for our random graph, or roughly 4 billion possible configurations (the number is much smaller if we don't count patterns that are rotations, translations or mirror images of each other; there are a lot of configurations that are connected to each other like that). ## How To Build These Sets, Step Three: Each Path Is A Number Too Now we're where I wanted to go with this notebook: how to go from one sequence to tiling an entire grid with non-overlapping sequences. In order to do that we have to match every `x` to every `y` in the grid. For any given `x` that means traversing down every possible path exactly once, which is like connecting one node at the top to every node at the bottom:39:T1f40,{ const points = [] const totalPoints = 16; const tpLog = Math.log2(totalPoints) const randomGraph = new Uint8Array(totalPoints / 2 * tpLog); const randomize = () => { points.length = 0; for (let i = 0; i < randomGraph.length; i++) randomGraph[i] = Math.random() < 0.5 ? 0 : 1; const mask = totalPoints - 1; const offsetMask = mask << tpLog; for (let x = 0; x < totalPoints; x++) { let y = x; for (let row = 0; row < tpLog; row++) { const odd = y & 1; const yMod = y & (mask >> row); const yOffset = y & (offsetMask >> row); const randomBit = randomGraph[row * 8 + (y >> 1)]; const turnRight = odd ^ randomBit; y = yOffset + (yMod >> 1) + (turnRight ? (8 >> row) : 0); } points.push([x, y]); } } // Drawing part const width = 240 + 16 * (30 + 20); const height = 64 + 5 * 178; const ctx = DOM.context2d(width, height); 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 = "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 = "#FFF"; for (let x = 0; x < xMax; x++) { for (let y = 0; y < yMax; y++) { const j = x + y * xMax; const xj = x0 + 160 + j * 50; const yn = yi + 64 - 15; 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; if (i < 4) { // parent pair highlights if (!(x&1)) { ctx.fillStyle = "#D0D0E8"; ctx.fillRect(xj - 4, yn - 9, 90, 46); ctx.fillStyle = "#FFF"; } // connections easedLine(ctx, xj + 15, yn + 30, xjnl, yn + 178); easedLine(ctx, xj + 15, yn + 30, xjnr, yn + 178); } ctx.fillRect(xj, yn, 30, 30); ctx.strokeRect(xj, yn, 30, 30); // tiles ctx.strokeRect(x0 + x * dx, yi + y * dy, dx, dy); } } // points ctx.lineWidth = 2; ctx.strokeStyle = "#000"; ctx.fillStyle = "#FFF"; for (let [xp, yp] of 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); } } drawNodeHover(225, 115); for (const [x, y] of points) { drawTakenPaths(x*8 + 50, y*8 + 50); } // draw bits ctx.fillStyle = "#000"; ctx.strokeStyle = "#FFF"; ctx.lineWidth = 4; for (let i = 0; i < 5; i++) { const yi = 50 + i * 178; const x0 = 50; const xMax = 16 >> i; const yMax = 16 >> (4 - i); const [dx, dy] = deltas[i]; for (let y = 0; y < yMax; y++) { for (let x = 0;x < xMax; x++) { // tiles const j = x + y * xMax; const xj = x0 + 160 + j * 50; const yn = yi + 64 - 15; const k = x + y * xMax // bit labels ctx.font = "12px monospace"; ctx.fillStyle = "#000"; ctx.strokeText("1", xj + 21, yn + 33); ctx.strokeText("0", xj + 4, yn + 33); ctx.fillText("1", xj + 21, yn + 33); ctx.fillText("0", xj + 4, yn + 33); if (i < 4 && (k&1)) { ctx.strokeText(randomGraph[(k>>>1) + i * 8], xj - 13, yn + 20); ctx.fillText(randomGraph[(k>>>1) + i * 8], xj - 13, yn + 20); } const gray = k.toString(2).padStart(4, "0"); const black = (k >>(4-i)).toString(2).padStart(i, "0"); ctx.strokeText(gray, xj + 3, yn + 5) ctx.fillStyle = "#888"; ctx.fillText(gray, xj + 3, yn + 5); if (i) { ctx.fillStyle = "#000"; ctx.fillText(black, xj + 3, yn + 5); } if (i < 4) { ctx.font = "bold 12px monospace"; ctx.fillText((k&1) ? " 1" : " 0", xj + 3, yn + 5); } } } } } 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(ctx, x, yl, xn1, yln); easedLine(ctx, x, yl, xn2, yln); } } function drawTakenPaths(x, y) { if (x >= 50 && x < 178) { const x0 = x - 50; const y0 = (y - 50) % 178; if (y0 >= 0 && y0 < 128) { ctx.strokeStyle = "#F44"; 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 = 5; 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(ctx, xn + 15, yn + 30, xnj + 15, ynj); } ctx.stroke(); } ctx.strokeStyle = "#000"; } } } 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); 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); points = nextPoints; } } } } ctx.canvas.value = () => { randomize(); drawBase(); } ctx.canvas.style.maxWidth = "100%"; yield ctx.canvas; await visibility(); randomize(); drawBase(); }3a:T9ec,But notice something about each node in our graph here: there always is a red and a black vertex going out of it. So we can think of our vertices as having a bit value too: red as `0`, and black as `1`. So far, we've only taken the red-red-red-red path from `x` to `y`, or the `0000` "path value". But why not visit every other binary path value? One thing we can immediately see is that when we visit every path value for every `x`, we must have visited each (x, y) pair exactly once, meaning we visited every point in the grid exactly once! That was one of our goals. The only thing left to do is figure out a way to do that so that it produces a set of valid non-overlapping `pmj02` sequences. That turns out to be fairly simple: we know that if we flip a parent pair we still produce a `pmj02` sequence. So we have to find a way to flip enough parent pairs to produce a second sequence that does not overlap with the first. That turns out to be the case when each point in a given sequence uses the same path value (and the same randomization graph, but that was implied). In short: as long as we use one randomization graph for all paths, and one specific path for each `pmj02` sequence, then none of the resulting sequences will with each other, and the set of all of these sequences will visit each point in the grid exactly once. So all we have to do to make this work is extend the earlier algorithm as follows: ~~~ javascript const totalPoints = 16, tpLog = Math.log2(totalPoints), points = [], randomGraph = [], mask = totalPoints - 1, offsetMask = mask << tpLog; for (let i = 0; i < totalPoints / 2 * tpLog; i++) { randomGraph[i] = Math.random() < 0.5 ? 0 : 1; } // value essentially iterates over every // flip bit combination that we can take through // our butterfly diagram for a given random graph // without passing a vertex more than once. for (let value = 0; value < totalPoints; value++) { // ⬅️ for (let x = 0; x < totalPoints; x++) { let y = x; for (let row = 0; row < tpLog; row++) { const odd = y & 1, yMod = y & (mask >> row), yOffset = y & (offsetMask >>> row), randomBit = randomGraph[row * (totalPoints / 2) + (y >> 1)], valueBit = (value >>> row) & 1; // ⬅️ turnRight = odd ^ randomBit ^ valueBit; // ⬅️ y = yOffset + (yMod >> 1) + (turnRight ? ((totalPoints / 2) >> row) : 0); } points.push([x, y, value]); // ⬅️ } } ~~~ And we can cover our grid:3b:T190a,{ restart6 const totalPoints = 16; const tpLog = Math.log2(totalPoints) const randomGraph = new Uint8Array(totalPoints / 2 * tpLog); for (let i = 0; i < randomGraph.length; i++) randomGraph[i] = Math.random() < 0.5 ? 0 : 1; // Drawing part 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] ]; ctx.fillStyle = "#FFF"; ctx.fillRect(0, 0, width, height); 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]; 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; 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(ctx, xj + 15, yn + 30, xjnl, yn + 178); easedLine(ctx, xj + 15, yn + 30, xjnr, yn + 178); } } } } function drawBase() { ctx.strokeStyle = "black"; ctx.lineWidth = 2; // draw bits ctx.fillStyle = "#000"; ctx.strokeStyle = "#FFF"; ctx.lineWidth = 4; for (let i = 0; i < 5; i++) { const yi = 50 + i * 178; const x0 = 50; const xMax = 16 >> i; const yMax = 16 >> (4 - i); const [dx, dy] = deltas[i]; for (let y = 0; y < yMax; y++) { for (let x = 0;x < xMax; x++) { // tiles const j = x + y * xMax; const xj = x0 + 160 + j * 50; const yn = yi + 64 - 15; const k = x + y * xMax // bit labels ctx.font = "10px monospace"; ctx.fillStyle = "#000"; ctx.strokeText("1", xj + 21, yn + 33); ctx.strokeText("0", xj + 4, yn + 33); ctx.fillText("1", xj + 21, yn + 33); ctx.fillText("0", xj + 4, yn + 33); if (i < 4 && (k&1)) { ctx.strokeText(randomGraph[(k>>>1) + i * 8], xj - 13, yn + 20); ctx.fillText(randomGraph[(k>>>1) + i * 8], xj - 13, yn + 20); } const gray = k.toString(2).padStart(4, "0"); const black = (k >>(4-i)).toString(2).padStart(i, "0"); ctx.fillStyle = "#888"; ctx.strokeText(gray, xj + 3, yn + 5) ctx.fillText(gray, xj + 3, yn + 5); if (i) { ctx.fillStyle = "#000"; ctx.fillText(black, xj + 3, yn + 5); } if (i < 4) { ctx.font = "bold 10px monospace"; ctx.fillText((k&1) ? " 1" : " 0", xj + 3, yn + 5); } } } } } 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(ctx, x, yl, xn1, yln); easedLine(ctx, x, yl, xn2, yln); } } function drawTakenPaths(x, y, style = "#F44") { if (x >= 50 && x < 178) { const x0 = x - 50; const y0 = (y - 50) % 178; if (y0 >= 0 && y0 < 128) { ctx.strokeStyle = style; 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 = 5; 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(ctx, xn + 15, yn + 30, xnj + 15, ynj); } ctx.stroke(); } ctx.strokeStyle = "#000"; } } } // points const coprime = findCoprimeApprox(totalPoints); const palette = makeRamp(totalPoints, "#viridis"); const mask = totalPoints - 1; const offsetMask = mask << tpLog; for (let i = 0; i < totalPoints; i++) { const points = [] for (let x = 0; x < totalPoints; x++) { let y = x; for (let row = 0; row < tpLog; row++) { const odd = y & 1; const yMod = y & (mask >> row); const yOffset = y & (offsetMask >> row); const pathBit = (i >>> row) & 1 const randomBit = randomGraph[row * 8 + (y >> 1)]; const turnRight = odd ^ randomBit ^ pathBit; y = yOffset + (yMod >> 1) + (turnRight ? (8 >> row) : 0); } points.push([x, y]); } ctx.lineWidth = 2; ctx.strokeStyle = "#000"; ctx.fillStyle = palette[i]; for (let [xp, yp] of points) { xp = ((xp - 1) * 128) / 16; yp = ((yp - 1) * 128) / 16; for (let j = 0; j < 5; j++) { const yj = 50 + j * 178; const x0 = 50; ctx.fillRect(59 + xp, 9 + yp + yj, 6, 6); } } ctx.beginPath(); for (const [x, y] of points) { drawTakenPaths(x*8 + 50, y*8 + 50, palette[i]); } drawBase(); yield ctx.canvas; await Promises.delay(100) } yield ctx.canvas; }3c:T839,

Aside: how one might use value(x, y, seed) for shuffling.

This is a condensed explanation, see ["Constant Time Stateless Shuffling and Grouping"](https://www.ea.com/seed/news/constant-time-stateless-shuffling) by Alan Wolfe & William Donnelly for more details.

When we shuffle an array we keep all original values, but give each of them a new index. That means we assign each index exactly once too, otherwise we would copy multiple values to the same position in the array. We can think of these indices as forming their own "permutation array" that we apply to our input array:

  
  shuffled[i] = input[indices[i]]

In the textures we produce each row/column/stratified tile contains each value exactly once. So for arrays of the same size, they could be used as shuffle indices too!

  
  shuffled[i] = input[value(i, pass, seed)]

We treat `x` as the index, and `y` as the `n`th pass of our shuffle. And that's essentially it.

What is really powerful about this, is that this lets you determine the value at *any* index of the shuffled array at *any* (valid) pass, without shuffling any other point in the array! And it lets you keep the original order of `input` for other purposes (also nice when dealing with immutable data).

This has trade-offs: we add a call to `value(i, pass, seed)` for every look-up. We also impose an extreme limitation of how randomized our shuffle really is. The latter however might be desirable, since there are still many ways we can shape this noise.

Another aspect is that `pass` being the y-coordinate implies that for an array of length `l` (which also is the length of the sides of our texture), we can only shuffle `l` times before we run out of unique shuffles. There are ways to work around that too, but let's save that for another discussion. And of course, by default our textures are highly regular, so this assumes the textures themselves have been shuffled as well. We will come back to that later in the notebook.

3d:T772,Let's again start with looking at what happens with a graph with only zero bits (i.e. a "zero seed"), when generation points for the zero value (i.e. zero path). This limits our `(x,y)` pairs to the simplest `pmj02` sequence. We discovered before that `y` is equal to the reversed bitstring of `x` in that case:

  y = rbits(x)
  x = rbits(y)

What about other `(x,y)` values? Well, we could model each `(x,y)` pair as `xor`'ing some unknown value `q` with our bit string:

              y  =  q  ^  rbits(x)
  y  ^  rbits(x) =  q
        rbits(x) =  y  ^  q
              x  =  rbits(y  ^  q)
              x  =  rbits(y) ^ rbits(q)

`q` must clearly be determined by the random graph and the `value` variable in our earlier algorithm. If the random graph is zero, then `q` turns out to be `rbits(value)`:

              y  = rbits(value)  ^  rbits(x)
        rbits(y) = value  ^  x
  x  ^  rbits(y) = value

So that answers `value(x, y, 0)`. Now we just need to extend this to non-zero seed values. Recall that in the absence of a random graph `turnRight` depends on the path taken through the graph, which depends on value. The random graph only adds a random bit to `xor` with to the mix:

  pathBit = (value >> n) & 1
  randomBit = randomGraph[row * (totalPoints / 2) + (y >> 1)],
  turnRight = pathBit ^ odd ^ randomBit

This gives us enough information to rewrite the earlier function: - we tumble from `x` down the graph to `y` like before - we know the effective path taken through the graph, because we know which `y` we end up with - we *also* know the direction `y` "wants" to take, in the absence of a value, based on `x` and the random graph - therefore, the `value` can be "extracted" by `xor`ing the bit of the direction the graph would have gone in its absence with the bit of the effective path That results in the following function:3e:T47f,{ const totalPoints = 256; const {log2} = Math; const ctx = DOM.context2d(totalPoints, totalPoints, 1); ctx.canvas.style.width = ctx.canvas.style.height = dprAdjusted; ctx.canvas.className = "pixelated"; yield ctx.canvas; await visibility(); // co-prime random graph const randomGraph = []; const rgLength = log2(totalPoints) * totalPoints / 2; for (let i = 0; i < rgLength; i++) { randomGraph[i] = i*2 > rgLength ? 0 : 1; } const coprimes = [ findCoprimeApprox(rgLength), // findCoprimeApprox(rgLength, 2**0.5) ]; const temp = randomGraph.slice(); for (const coprime of coprimes) { for (let i = 0, j = 0; i < rgLength; i++) { temp[i] = randomGraph[j]; j = (j + coprime) % rgLength } for (let i = 0; i < rgLength; i++) { randomGraph[i] = temp[i]; } } const gen = makeValueGen(totalPoints, randomGraph); const palette = makeRamp(totalPoints, "#grayscale"); for (let y = 0; y < totalPoints; y++) { for (let x = 0; x < totalPoints; x++) { const v = gen(x, y); ctx.fillStyle = palette[v]; ctx.fillRect(x, y, 1, 1); } } return ctx.canvas }3f:Ta1c,{ const totalPoints = 1024; const tileSizes = [ [1024, 1], [512, 2], [256, 4], [128, 8], [64, 16], [32, 32], [16, 64], [8, 128], [4, 256], [2, 512], [1, 1024] ]; const passed = Array.from({length: 11}, _ => false); const makeSummary = () => { let summary = "\n\n"; for(let i = 0; i < 10; i++) { const [dx, dy] = tileSizes[i]; const check = passed[i]; summary += `- ${check ? "☑" : "☐"} ${dx} × ${dy}\n`; } summary += `- ${passed[10] ? "☑" : "☐"} total values\n` return summary; }; yield md`***Initiating test, creating generator***${makeSummary()}`; await visibility(); const gen = makeValueGen(totalPoints); const points = new Uint16Array(totalPoints*totalPoints); const values = new Map(); const insert = v => { let count = values.get(v) | 0; values.set(v, count + 1); }; yield md`***Initiating test, generating points***${makeSummary()}`; for (let y = 0; y < totalPoints; y++) { for (let x = 0; x < totalPoints; x++) { const p = gen(x, y); points[x + y*totalPoints] = p } } const tileCheck = new Uint16Array(totalPoints); for (let i = 0; i < 10; i++) { const [dx, dy] = tileSizes[i]; yield md`***Testing tile ${dx} × ${dy}***${makeSummary()}`; for (let x0 = 0; x0 < totalPoints/dx; x0 += dx) { for(let y0 = 0; y0 < totalPoints/dy; y0 += dy) { tileCheck.fill(0); for (let x = 0; x < dx; x++) { for (let y = 0; y < dy; y++) { const p = points[x0 + x + (y0 + y) * totalPoints]; if (tileCheck[p]++) { console.log(points) throw new Error(`tile [${dx}, ${dy}] ${x0} ${y0} at ${x0+x},${y0+y} contains ${p} more than once!`); } } } } } passed[i] = true; } yield md`***Testing total values***${makeSummary()}`; for (let y = 0; y < totalPoints; y++) { for (let x = 0; x < totalPoints; x++) { insert(points[x + y*totalPoints]); } } const finalCount = [...values]; // sort by key finalCount.sort((a, b) => a[0] - b[0]); // verify that each point occurs exactly totalPoints times for (let i = 0; i < totalPoints; i++) { const p = finalCount[i]; if (!p[0] === i) { throw new Error(`value ${i} expected, found ${p[0]} at $[{p.toString()}]`); } if (!p[1] === totalPoints) { throw new Error(`value ${i} did not occur ${totalPoints} times, but ${p[1]}`); } } passed[10] = true; yield md`***All tests passed!***${makeSummary()}`; }40:T463,Well, that certainly makes it look more scrambled, but we also still see (or more like "feel") patterns in the output. If you don't, just shuffle a few times. You'll also notice that the results sometimes are lumpier and "more textured" than others, just like with regular white noise. That's combining the bad aspects of both worlds: the lumpy bits of white noise, and the patterns of non-random output. However, there's a fairly simple way to improve that, courtesy of Alan Wolfe and Marc Reynolds. The trick is to shuffle the array using a highly irrational number, like the golden ratio or the square root of 2 or 3. We don't use the irrational number directly, but look for an integer approximation that is coprime to the length of the array that we wish to shuffle. This way we ensure that we project from each index to an index in the array that has not been visited yet. Alan Wolfe explains this trick in more detail [here](https://blog.demofox.org/2020/07/26/irrational-numbers/) and [here](https://blog.demofox.org/2024/05/19/a-low-discrepancy-shuffle-iterator-random-access-inversion/), among other blog posts.41:T41b,{ const ctx = DOM.context2d(512, 512); ctx.fillStyle = "#000" ctx.fillRect(0, 0, 256, 512); ctx.strokeStyle = "#000"; ctx.font = "bold 16px monospace"; ctx.textAlign = "center" ctx.lineWidth = 4; const palette = makeRamp(16, "#grayscale"); for (let i = 0; i < 16; i++) { const ri = reverseNbits(i, 4); ctx.fillStyle = palette[i]; ctx.fillRect(0, i*512/16, 256, 512/16); ctx.fillStyle = palette[ri]; ctx.fillRect(256, i*512/16, 256, 512/16); ctx.fillStyle = "#FFF"; ctx.strokeText(`${i.toString(2).padStart(4, "0")} (${i.toString().padStart(2, " ")})`, 128, i*512/16 + 24) ctx.fillText(`${i.toString(2).padStart(4, "0")} (${i.toString().padStart(2, " ")})`, 128, i*512/16 + 24) ctx.strokeText(`${ri.toString(2).padStart(4, "0")} (${ri.toString().padStart(2, " ")})`, 256+128, i*512/16 + 24) ctx.fillText(`${ri.toString(2).padStart(4, "0")} (${ri.toString().padStart(2, " ")})`, 256+128, i*512/16 + 24) } ctx.canvas.style.width = ctx.canvas.style.height = dprAdjusted; return ctx.canvas; }42:T429,When we go up one step at a time normally, the least significant bit changes the most often (odd, even, odd), the second-least significant bit half as often, the third least significant bit half as often as that, and so on. The most significant bit only changes once: halfway through. Reversing the bit string flips that around, and because the most significant bit always encodes the biggest change in value, the second most significant bit the second most, and so on, the result is that we "maximally alternate" our numbers will still visiting each number in the range. Which is great if we want to generate something noisy-looking! As you might expect based on the above explanation this works best when applied to a path that follows a smooth gradient. This is what our space-filling curves essentially do: they project a one-dimensional gradient into a 2D space. So one could expect that reversing the values of that output might also result in an even spread of values in 2D. Let's first look at the output without a randomized curve to see if that holds up:43:T191c,async function makeTile({ tilesize: n = 256, graphStructure = "ascending", shuffleColors = "", colorscale = "grayscale", coprimeFactor = 2/(1+5**0.5), graphCoPrimeFactor = 2**0.5, showTile = "Hilbert", coprimeRow = 0, coprimeCol = 0, visibility = null } = {}){ if (visibility) await visibility(); const ctx = DOM.context2d(n, n, 1); ctx.canvas.style.width = ctx.canvas.style.height = dprAdjusted; ctx.canvas.className = "pixelated"; ctx.fillStyle = "red"; ctx.fillRect(0, 0, n, n); const imgData = ctx.getImageData(0, 0, n, n); const points = makeNxNpoints(n, graphStructure, graphCoPrimeFactor); const paletteOrder = new Uint32Array(n); const sqrtn = Math.sqrt(n); const {log2, random} = Math; if (shuffleColors === "random" || shuffleColors === "random rows/columns") { for (let i = 0; i < n; i++) { const j = Math.random() * (i+1) | 0; paletteOrder[i] = paletteOrder[j]; paletteOrder[j] = i; } if (shuffleColors === "random rows/columns") { const project = points.slice(); let rowShuffle = [], colShuffle = []; for (let i = 0; i < n; i++) { rowShuffle.push(i); colShuffle.push(i); } for (let i = 0; i < n; i++) { const j = random() * (i+1) | 0; const k = random() * (i+1) | 0; let t = rowShuffle[j]; rowShuffle[j] = rowShuffle[i]; rowShuffle[i] = t; t = colShuffle[k]; colShuffle[k] = colShuffle[i]; colShuffle[i] = t; } for (let row = 0; row < n; row++) { const ri = row * n; const rp = rowShuffle[row] * n; for (let col = 0; col < n; col++) { points[col + ri] = project[colShuffle[col] + rp]; } } } } else if (showTile === "row") { if (shuffleColors === "coprime") { const step = findCoprimeApprox(n, coprimeFactor); for (let i = 0, j = random()*n|0; i < n; i++) { paletteOrder[j] = i; j = (j + step) % n; } } else { for (let i = 0; i < n; i++) paletteOrder[i] = i; } } else if ((showTile || showTile === 0) && Math.sqrt(n)%1 === 0) { const logn = log2(n); const p0 = n*n - sqrtn*n; for (let y = 0; y < sqrtn; y++) { for (let x = 0; x < sqrtn; x++) { const p = points[p0 + x + y * n]; switch(showTile) { case "Bayer": { let v = 0, yv = y, xv = x; yv ^= xv; for (let i = 0; i < (log2(sqrtn)); ++i) { v = (v << 2) | (xv & 1) + ((yv & 1) << 1); xv >>= 1; yv >>= 1; } paletteOrder[p] = v >>> 0; break; } case "Bayer reverse": { let v = 0, yv = y, xv = x; yv ^= xv; for (let i = 0; i < (log2(sqrtn)); ++i) { v = (v << 2) | (xv & 1) + ((yv & 1) << 1); xv >>= 1; yv >>= 1; } paletteOrder[p] = reverseNbits(v >>> 0, logn); break; } case "Hilbert": paletteOrder[p] = hilbert(x, y); break; case "Hilbert reverse": paletteOrder[p] = reverseNbits(hilbert(x, y), logn); break; case "H-curve": paletteOrder[p] = hcurve[sqrtn].gridToLine[x + y * sqrtn]; break case "H-curve reverse": paletteOrder[p] = reverseNbits(hcurve[sqrtn].gridToLine[x + y * sqrtn], logn); break case "H-curve pinwheel": paletteOrder[p] = pinwheel[sqrtn].gridToLine[x + y * sqrtn]; break; case "H-curve pinwheel reverse": paletteOrder[p] = reverseNbits(pinwheel[sqrtn].gridToLine[x + y * sqrtn], logn); break; case "Morton": { let c = 0; for (let i = sqrtn-1; i >= 0; i--) { c = (c<<1) + ((x >> i) &1); c = (c<<1) + ((y >> i) &1); } paletteOrder[p] = c; break; } case "Morton reverse": { let c = 0; for (let i = sqrtn-1; i >= 0; i--) { c = (c<<1) + ((x >> i) &1); c = (c<<1) + ((y >> i) &1); } paletteOrder[p] = reverseNbits(c, logn); break; } // all values in range 0-n case (Number.isInteger(showTile) ? showTile : !showTile): { const c = (showTile >= 0 && showTile < n) ? showTile : 0; for (let i = 0; i < n; i++) { paletteOrder[i] = (i === c) ? paletteOrder.length-1 : 0; } break; } case "row reverse": case "reverse": paletteOrder[p] = reverseNbits(x + y * sqrtn, logn); break; default: case "gradient": paletteOrder[p] = x + y * sqrtn; break; } } } if (shuffleColors === "coprime") { let step = findCoprimeApprox(n, coprimeFactor); let j = Math.random() * n | 0; let coprimeshuffle = []; for (let i = 0; i < n; i++) { coprimeshuffle[j] = i; j = (j + step) % n } for (let i = 0; i < n; i++) { paletteOrder[i] = coprimeshuffle[paletteOrder[i]]; } } } else { for (let i = 0; i < n; i++) { paletteOrder[points[i]] = i; } } if (coprimeRow > 0 || coprimeCol > 0) { const project = points.slice(); const rowStep = findCoprimeApprox(n, coprimeRow); const colStep = findCoprimeApprox(n, coprimeCol); let rowShuffle = [], colShuffle = []; for (let i = 0, j = rowStep, k = colStep; i < n; i++) { rowShuffle.push(j); colShuffle.push(k); j = (j + rowStep) % n; k = (k + colStep) % n; } for (let i = 0; i < points.length; i++) { points[i] = project[(colShuffle[i % n] + rowShuffle[(i/n)|0] * n) % points.length]; } } const colorramp = makeRamp(n, colorscale); for (let i = 0, {data} = imgData; i < points.length; i++) { const color = paletteOrder[points[i]]; // const color = points[i] === 0 ? 0 : 0xFF; data[i * 4] = colorramp[color*3] data[i * 4 + 1] = colorramp[color*3+1] data[i * 4 + 2] = colorramp[color*3+2] } ctx.putImageData(imgData, 0, 0); ctx.canvas.value = {ctx, canvas: ctx.canvas, imgData, points, paletteOrder}; return ctx.canvas; } 44:T72f,makeNxNpoints = { let randomGraph = new Uint8Array(0); return function makeNxNpoints(n, graphStructure = "ascending", graphCoPrimPow = 1, debug = false) { const logn = Math.log2(n); if ((logn|0) !== logn) throw new Error("only works with n == whole power of two"); // we only need logn layers of n/2 bits to encode the graph, // but I'm lazy and use a byte per bit instead. if (randomGraph.length < logn * n/2) { randomGraph = new Uint8Array(logn * n/2); } else { randomGraph.fill(0); } switch(graphStructure) { case "random": for (let i = 0, {random} = Math; i < randomGraph.length; i++) { randomGraph[i] = random() < 0.5 ? 0 : 1; } break; case "coprime": { const step = findCoprimeApprox((n/2)**graphCoPrimPow); for (let i = 0, j = 0; i < randomGraph.length; i++) { randomGraph[j] ^= i&1; j = (j + step) % randomGraph.length; } break; } case "ascending": break; default: case "blue": randomGraph.set(blueGraphs[n]) break; } const points = new Int32Array(n*n); points.fill(-1); const mask = n - 1; const offsetMask = mask * n; for (let path = 0, odd; path < n; path++) { for(let x = 0; x < n; x++) { let y = x; for (let row = 0; row < logn; row++) { const pathBit = (path >> row) & 1 odd = (y & 1); const yMod = y & (mask >> row); const yOffset = y & (offsetMask >> row); const randomBit = randomGraph[row * (n >> 1) + (y >> 1)]; const turnRight = odd ^ randomBit ^ pathBit; y = yOffset + (yMod >> 1) + (turnRight ? ((n >> 1) >> row) : 0); } points[x + y * n] = path; } } return points; } }45:T6a6,function makeRamp(n, rampName){ let returnString = rampName[0] === "#"; if (returnString) { rampName = rampName.substring(1, rampName.length); } const palette = new Uint8Array(n * 3); const strColors = [] if (rampName === "grayscale") { if (returnString) { for (let i = 0; i < n; i++) { const c = ((i * 256 / n) | 0).toString(16).padStart(2, "0"); strColors[i] = `#${c}${c}${c}`; } } else { for (let i = 0; i < n; i++) { palette[i*3] = i * 256 / n; palette[i*3+1] = i * 256 / n; palette[i*3+2] = i * 256 / n; } return palette; } } else { const rampFn = ({ spectral: d3.interpolateSpectral, cividis: d3.interpolateCividis, viridis: d3.interpolateViridis, magma: d3.interpolateMagma, inferno: d3.interpolateInferno, plasma: d3.interpolatePlasma, jet: d3.interpolateTurbo, turbo: d3.interpolateTurbo, sinebow: d3.interpolateSinebow, warm: d3.interpolateWarm })[rampName] || d3.interpolateViridis; for (let i = 0; i < n; i++){ let c = rampFn(i / n), colors; if (c.startsWith("#")) { // already "#RRGGBB" format colors = [c.slice(1, 3), c.slice(3, 5), c.slice(5, 7)].map(str => +("0x"+str)); } else { // "rgb(R, G, B)", convert to "#RRGGBB" format c = c.substring(4, c.length - 1); colors = c.split(", ").map(str => +str); c = colors.reduce((a, v) => a+v.toString(16).padStart(2, "0"), "#"); } strColors.push(c); palette[i*3] = colors[0] palette[i*3+1] = colors[1] palette[i*3+2] = colors[2] } } return returnString ? strColors : palette; }46:Tdc5,viewof PMJButterflyGraph = { class PMJButterflyGraph { constructor(power) { const layers = this.layers = power + 1; const stride = this.stride = 2 ** power; // Instead of booleans use counters for each node this.nodes = new Uint32Array(stride * layers); const picks = this.picks = new Uint32Array(stride); for (let i = 0; i < stride; i++) picks[i] = i; this.points = []; } clear() { this.nodes.fill(0); const {stride, picks} = this; for (let i = 0; i < stride; i++) picks[i] = i; this.points.length = 0; } get size() { return this.points.length; } fill() { this.clear(); while (this.points.length < this.stride) this.pick(); return this; } pick(leftLeaning = true) { // 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 explanation of 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]++; } else { // pick left side and mark it occupied nodes[layerOffset + xi + yc1]++; } // 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; }47:T1fa5,function makeParentSiblingPartnerChildFigure() { const dpr = window.devicePixelRatio || 1; const ctx = DOM.context2d(5 * 148 + 100, 228, dpr); ctx.canvas.style.maxWidth = "100%"; ctx.font = "bold 12px monospace"; ctx.textBaseline = "middle"; function determineLayer(x, y) { if (y > 30 && y < 158) { let i = 5; while (i-- > 0) { if (x > 50 + 148 * i && x < 50 + 128 + 148 * i) return i; } return i; } } function hoverDraw(mx, my) { const layer = determineLayer(mx, my); const parentFill = "#AAAAFF"; const partnerFill = "#66DD33"; const tileFill = "#DD6666"; const parentStroke = "#CCCCFF"; const childStroke = "#99EE66"; const tileStroke = "#FF9999"; const parentGlyph = (x, y, stroke = parentStroke) => { ctx.lineWidth = 4; ctx.strokeStyle = stroke; ctx.strokeRect(x+1, y-4, 4, 4); ctx.lineWidth = 2; ctx.strokeStyle = "#000"; ctx.strokeRect(x+1, y-4, 4, 4); } const childGlyph = (x, y, stroke = childStroke) => { ctx.lineWidth = 3.5; ctx.strokeStyle = stroke; ctx.strokeRect(x+1, y-4, 4, 4); ctx.fillRect(x+.5, y-4.5, 5, 5); } if (layer > -1) { ctx.fillStyle = "#FFF"; ctx.fillRect(0, 0, 5*148+100, 228); ctx.fillStyle = "#DDD"; for (let i = 0; i < 5; i++) { const x = 50 + 148 * i; const y = 30; ctx.fillRect(x, y, 128, 128); } const xOffset = 50 + layer * 148; const xi = mx - xOffset; const yi = my - 30; switch(layer) { case 0: { // no parents, only partner and children const tile = xi >> 3; const partner = (tile >> 1) << 1; ctx.fillStyle = partnerFill ctx.fillRect(xOffset + partner*8, 30, 16, 128); ctx.fillRect(xOffset + 148 + partner*8, 30, 16, 128); // current tile ctx.fillStyle = tileFill ctx.fillRect(xOffset + tile*8, 30, 8, 128); // symbols ctx.fillStyle = "#000"; ctx.strokeStyle = "#FFF"; const x0 = xOffset + 1 + partner*8; const x1 = xOffset + 148 + 5 + partner*8 const y0 = 64; const y1 = 128; childGlyph(x0, y0); childGlyph(x0 + 8, y1); childGlyph(x1, y0); childGlyph(x1, y1); break; } case 1: { const xtile = xi >> 4; const ytile = yi > 64 ? 94 : 30; const sibling = yi > 64 ? 30 : 94; const partner = (xtile >> 1) << 1; ctx.fillStyle = partnerFill ctx.fillRect(xOffset + partner*16, ytile, 32, 64); ctx.fillRect(xOffset + 148 + partner*16, ytile, 32, 64); ctx.fillStyle = tileFill ctx.fillRect(xOffset + xtile*16, ytile, 16, 64); ctx.fillStyle = parentFill ctx.fillRect(xOffset + xtile*16, sibling, 16, 64); ctx.fillRect(xOffset - 148 + xtile*16, 30, 16, 128); ctx.fillStyle = "#000"; ctx.strokeStyle = "#FFF"; parentGlyph(xOffset + xtile*16 + 5, 64); parentGlyph(xOffset + xtile*16 + 5, 128); parentGlyph(xOffset + xtile*16 + 1 - 148, 64); parentGlyph(xOffset + xtile*16 + 9 - 148, 128); childGlyph(xOffset + partner*16 + 5, ytile+16); childGlyph(xOffset + partner*16 + 21, ytile+48); childGlyph(xOffset + 148 + 12 + partner*16, ytile+16); childGlyph(xOffset + 148 + 12 + partner*16, ytile+48); break; } case 2: { const xtile = xi >> 5; const ytile = yi >> 5; const partner = (xtile + ((xtile&1) ? -1 : 1))*32; const sibling = 30 + (ytile + ((ytile&1) ? -1 : 1))*32; ctx.fillStyle = tileFill ctx.fillRect(xOffset + xtile*32, 30 + ytile*32, 32, 32); // partner + children ctx.fillStyle = partnerFill ctx.fillRect(xOffset + partner, 30 + ytile*32, 32, 32); ctx.fillRect(xOffset + 148 + partner, 30 +ytile*32, 32, 32); ctx.fillRect(xOffset + 148 + xtile*32, 30 +ytile*32, 32, 32); // parents + sibling ctx.fillStyle = parentFill ctx.fillRect(xOffset + xtile*32, sibling, 32, 32); ctx.fillRect(xOffset - 148 + xtile*32, sibling, 32, 32); ctx.fillRect(xOffset - 148 + xtile*32, 30 + ytile*32, 32, 32); ctx.fillStyle = "#000"; ctx.strokeStyle = "#FFF"; const y0 = ((ytile>>1)<<1); parentGlyph(xOffset + xtile*32 + 12, y0*32 + 48); parentGlyph(xOffset + xtile*32 + 12, y0*32 + 80); parentGlyph(xOffset - 148 + 4 + xtile*32, y0*32 + 48); parentGlyph(xOffset - 148 + 20 + xtile*32, y0*32 + 80); childGlyph(xOffset + (xtile>>1)*64 + 12, ytile*32 + 40); childGlyph(xOffset + (xtile>>1)*64 + 44, ytile*32 + 56); childGlyph(xOffset + (xtile>>1)*64 + 176, ytile*32 + 40); childGlyph(xOffset + (xtile>>1)*64 + 176, ytile*32 + 56); break; } case 3: { const xtile = xi >> 6; const ytile = yi >> 4; const partner = (xtile + ((xtile&1) ? -1 : 1))*64; const sibling = 30 + (ytile + ((ytile&1) ? -1 : 1))*16; ctx.fillStyle = tileFill ctx.fillRect(xOffset + xtile*64, 30 + ytile*16, 64, 16); // partner + children ctx.fillStyle = partnerFill ctx.fillRect(xOffset + partner, 30 + ytile*16, 64, 16); ctx.fillRect(xOffset + 148 + partner, 30 +ytile*16, 64, 16); ctx.fillRect(xOffset + 148 + xtile*64, 30 +ytile*16, 64, 16); // parents + sibling ctx.fillStyle = parentFill ctx.fillRect(xOffset + xtile*64, sibling, 64, 16); ctx.fillRect(xOffset - 148 + xtile*64, sibling, 64, 16); ctx.fillRect(xOffset - 148 + xtile*64, 30 + ytile*16, 64, 16); ctx.fillStyle = "#000"; ctx.strokeStyle = "#FFF"; parentGlyph(xOffset + xtile*64 + 12, (ytile>>1)*32 + 40); parentGlyph(xOffset + xtile*64 + 44, (ytile>>1)*32 + 56); parentGlyph(xOffset + xtile*64 + 12 - 148, (ytile>>1)*32 + 40); parentGlyph(xOffset + xtile*64 + 44 - 148, (ytile>>1)*32 + 56); childGlyph(xOffset + 28, ytile * 16 + 36); childGlyph(xOffset + 64 + 28, ytile * 16 + 44); childGlyph(xOffset + 28 + 148, ytile * 16 + 36); childGlyph(xOffset + 64 + 28 + 148, ytile * 16 + 44); break; } case 4: { const tile = yi >> 3; const sibling = (tile >> 1) << 1; // parents + sibling ctx.fillStyle = parentFill ctx.fillRect(xOffset, 30 + sibling*8, 128, 16); ctx.fillRect(xOffset - 148, 30 + sibling*8, 128, 16); ctx.fillStyle = tileFill ctx.fillRect(xOffset, 30 + tile*8, 128, 8); ctx.fillStyle = "#000"; ctx.strokeStyle = "#FFF"; parentGlyph(xOffset + 28, (tile>>1) * 16 + 36); parentGlyph(xOffset + 64 + 28, (tile>>1) * 16 + 44); parentGlyph(xOffset + 28 - 148, (tile>>1) * 16 + 36); parentGlyph(xOffset + 64 + 28 - 148, (tile>>1) * 16 + 44); break; } default: break; } } 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 = 30; 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); } } } } hoverDraw(400, 95); ctx.canvas.addEventListener("mousemove", e => { const box = ctx.canvas.getBoundingClientRect(); const mx = (((e.clientX - box.left) / box.width) * ctx.canvas.width) / dpr; const my = (((e.clientY - box.top) / box.height) * ctx.canvas.height) / dpr; hoverDraw(mx, my); }); return ctx.canvas; }48:[{"id":"51fc56299a74f5682ac4ae8685c9f20f26b9ad1870ccb7d79b5b98b9cf1e29cba2a2ba7344f66def42ed72ce89e5f979330dff8212ee785fddf13a6c16dc5e92","href":"https://static.observableusercontent.com/files/51fc56299a74f5682ac4ae8685c9f20f26b9ad1870ccb7d79b5b98b9cf1e29cba2a2ba7344f66def42ed72ce89e5f979330dff8212ee785fddf13a6c16dc5e92","name":"4096.png","type":"image/png","lastModified":1775724809755,"size":3181},{"id":"01928d88593ae92fe75ebcbfcbe4dce0c9ee95c88ba8303b4f2ec22d21cf5578e297b1344f72fa755c20d8983fc69868c42cd3e33c1aac04280591345cac935a","href":"https://static.observableusercontent.com/files/01928d88593ae92fe75ebcbfcbe4dce0c9ee95c88ba8303b4f2ec22d21cf5578e297b1344f72fa755c20d8983fc69868c42cd3e33c1aac04280591345cac935a","name":"0128.png","type":"image/png","lastModified":1775724809753,"size":156},{"id":"169e95fdd3280027a69249b7f99dc3d7792c696e3773b913ecf2ff435dfb4564d2aac537c6f3653c60bf40b365cf2560816e04d165f94bfec265d3feae4722ca","href":"https://static.observableusercontent.com/files/169e95fdd3280027a69249b7f99dc3d7792c696e3773b913ecf2ff435dfb4564d2aac537c6f3653c60bf40b365cf2560816e04d165f94bfec265d3feae4722ca","name":"0016.png","type":"image/png","lastModified":1775724809749,"size":97},{"id":"41305ffc02ce417dbf36c60ed17a68ba9869309d710dd72d259b532cdc8a901f0d6cb62996a9a960efb97c8b8bb7ea3b11998545ddc8499d5f431546e7fcbf82","href":"https://static.observableusercontent.com/files/41305ffc02ce417dbf36c60ed17a68ba9869309d710dd72d259b532cdc8a901f0d6cb62996a9a960efb97c8b8bb7ea3b11998545ddc8499d5f431546e7fcbf82","name":"0064.png","type":"image/png","lastModified":1775724809753,"size":121},{"id":"a809aeb4ccb4f688be2eb3d13987e442e341cd71d64df4592d21e9949e13e499ecf773c1381871cfb5003b8d3f0b0cfeaf39d1870771c74b503e37526c1e635d","href":"https://static.observableusercontent.com/files/a809aeb4ccb4f688be2eb3d13987e442e341cd71d64df4592d21e9949e13e499ecf773c1381871cfb5003b8d3f0b0cfeaf39d1870771c74b503e37526c1e635d","name":"0512.png","type":"image/png","lastModified":1775724809753,"size":393},{"id":"525d868a98db0a91e7dedb9366b7b7112e56bd56b923a0058be798243eadd32e9cc473ebba53c9fbb213e6af00fc1cf153b1eaceedd191186fe00cbb853be644","href":"https://static.observableusercontent.com/files/525d868a98db0a91e7dedb9366b7b7112e56bd56b923a0058be798243eadd32e9cc473ebba53c9fbb213e6af00fc1cf153b1eaceedd191186fe00cbb853be644","name":"0032.png","type":"image/png","lastModified":1775724809757,"size":106},{"id":"7cdf373d71c6a61f14f398937228a3b0640b1131a2dd4c0d4c761485712dac43dd4842f5304ef4e2be708d620209fb124ac03cfeb5727e574918514ba2dceb57","href":"https://static.observableusercontent.com/files/7cdf373d71c6a61f14f398937228a3b0640b1131a2dd4c0d4c761485712dac43dd4842f5304ef4e2be708d620209fb124ac03cfeb5727e574918514ba2dceb57","name":"2048.png","type":"image/png","lastModified":1775724809753,"size":1516},{"id":"4ce95b4f5a0a647e3d1dbebe6cffec7238172a606881d770a1eb85c271cf8383317de44853e406acca65e0d8129fa9fe6b7c6e513901880c0bdcbd5554bc8c71","href":"https://static.observableusercontent.com/files/4ce95b4f5a0a647e3d1dbebe6cffec7238172a606881d770a1eb85c271cf8383317de44853e406acca65e0d8129fa9fe6b7c6e513901880c0bdcbd5554bc8c71","name":"0256.png","type":"image/png","lastModified":1775724809754,"size":229},{"id":"709ba2fbf2cbb81421fdd0bcd27533405358ea53a69274716170fc8dce092c06fc4d053d59707052023f7c0c9384ecbbf1a048460b9267b7cdfa098b67408756","href":"https://static.observableusercontent.com/files/709ba2fbf2cbb81421fdd0bcd27533405358ea53a69274716170fc8dce092c06fc4d053d59707052023f7c0c9384ecbbf1a048460b9267b7cdfa098b67408756","name":"1024.png","type":"image/png","lastModified":1775724809761,"size":746}] 2a:{"title":"Covering the grid with progressive multi-jittered sample sequences","cells":[{"id":1929,"value":"# Covering the grid with progressive multi-jittered sample sequences","pinned":false,"hidden":false,"mode":"md"},{"id":1931,"value":"$2c","pinned":false,"hidden":false,"mode":"ojs"},{"id":1780,"value":"Suppose we want to place sixteen points in a rectangle, and in such a way that if we divide that rectangle into 16×1, 8×2, 4×4, 2×8, and 1×16 grids, each tile in those grids is covered by with exactly one point:","pinned":false,"hidden":false,"mode":"md"},{"id":15,"value":"$2d","pinned":false,"hidden":false,"mode":"ojs"},{"id":0,"value":"$2e","pinned":false,"hidden":false,"mode":"md"},{"id":1791,"value":"$2f","pinned":false,"hidden":false,"mode":"ojs"},{"id":4324,"value":"$30","pinned":false,"hidden":false,"mode":"ojs"},{"id":1789,"value":"The resulting set of points has some very interesting qualities. By construction, each row contains exactly one point of each sequence. Each column contains exactly one point of each sequence. In fact, each stratified tile meets this requirement! It's kind of like a sudoku with somewhat different constraining rules.\n\nThis might be useful in various ways. For example, as a noise or pseudo-random number generator, where the noise is \"shaped\" by these guarantees. The stratification ensures a certain evenness of the output, and you could use that in a game to guarantee that a player can't get unlucky with random dice rolls too long in a row. It might also have some applications in 3D graphics - I'm not a 3D graphics person, but I understand `pmj02` sequences have some applications there. Similarly in numerical integration, apparently, but I'll leave that to the experts.\n","pinned":false,"hidden":false,"mode":"md"},{"id":6201,"value":"## How To Build These Sets, Step One: How To Build One Sequence\n\nLet's first quickly go over the explanation of the previous notebook, since we need that basis for the rest of the algorithm. Look at the red tile in the image below.","pinned":false,"hidden":false,"mode":"md"},{"id":6210,"value":"makeParentSiblingPartnerChildFigure()","pinned":false,"hidden":false,"mode":"ojs"},{"id":1918,"value":"$31","pinned":false,"hidden":false,"mode":"md"},{"id":1920,"value":"$32","pinned":false,"hidden":false,"mode":"ojs"},{"id":6307,"value":"*(The light blue highlight boxes mark all nodes representing a pair of parents)*\n\nYou can hover over a tile to highlight all of the tiles that overlap with that position, as well as the path this creates through the graph. Hover over a node to highlight its matching tile, all of the tiles that partially overlap with that tile, and all of the nodes associated with those tiles.\n\nNow as mentioned, any empty node always has at least one valid path to an empty node at the top and bottom of the graph. So to generate points, all we that have to do is pick one of those paths at random to find a valid (x, y) position for a new point. We add that point to the sequence, mark all nodes in the chosen path as occupied, then repeat until all points are generated and the entire graph is occupied. The simplest way to do that is to pick a random empty x-position and \"tumble down\" to an empty y-position.","pinned":false,"hidden":false,"mode":"md"},{"id":39,"value":"$33","pinned":false,"hidden":false,"mode":"ojs"},{"id":90,"value":"viewof restart1 = html``","pinned":false,"hidden":false,"mode":"ojs"},{"id":2052,"value":"$34","pinned":false,"hidden":false,"mode":"md"},{"id":2057,"value":"$35","pinned":false,"hidden":false,"mode":"ojs"},{"id":2293,"value":"$36","pinned":false,"hidden":false,"mode":"md"},{"id":2381,"value":"$37","pinned":false,"hidden":false,"mode":"ojs"},{"id":2489,"value":"{\n const button = html``\n console.log(button)\n button.onclick = () => parentpairgraph();\n return button;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":2494,"value":"$38","pinned":false,"hidden":false,"mode":"md"},{"id":6415,"value":"$39","pinned":false,"hidden":false,"mode":"ojs"},{"id":2765,"value":"$3a","pinned":false,"hidden":false,"mode":"md"},{"id":2526,"value":"$3b","pinned":false,"hidden":false,"mode":"ojs"},{"id":2557,"value":"viewof restart6 = html``","pinned":false,"hidden":false,"mode":"ojs"},{"id":2596,"value":"And there we have it! A way to generate `n` sets of `pmj02(n)` sequences that cover an `n×n` grid in. Each sequence takes `O(nlogn)` time, so that suggests that our algorithm scales in `O(n²logn)` time. Except that's a bit misleading, since that n² is really just the surface area of the resulting grid, which equals the total number of points. Which is what we should set `n` to, meaning the algorithm is actually `O(nlogn)`.\n\nYou might have noticed that the basic output looks highly regular, even with a randomized graph","pinned":false,"hidden":false,"mode":"md"},{"id":3909,"value":"makeTile({\n tilesize: 256, graphStructure: \"ascending\",\n shuffleColors: \"\",\n colorscale: \"grayscale\",\n showTile: \"row\",\n visibility\n});","pinned":false,"hidden":false,"mode":"ojs"},{"id":3906,"value":"makeTile({\n tilesize: 256,\n graphStructure: \"blue\",\n shuffleColors: \"\",\n colorscale: \"grayscale\",\n showTile: \"row\",\n visibility\n});","pinned":false,"hidden":false,"mode":"ojs"},{"id":3495,"value":"However, there are ways to shuffle the output in a way that makes it look like random noise while still keeping all the previous stratified qualities:","pinned":false,"hidden":false,"mode":"md"},{"id":3496,"value":"makeTile({\n tilesize: 256,\n graphStructure: \"random\",\n shuffleColors: \"coprime\",\n colorscale: \"grayscale\",\n coprimeFactor: 1/2**0.5,\n showTile: \"H-curve pinwheel reverse\",\n visibility\n});","pinned":false,"hidden":false,"mode":"ojs"},{"id":5426,"value":"It's a bit of a long side-track though, so I'll postpone that discussion for later. There is a different problem that I want to solve first.\n\nThe algorithm so far generates `(x, y)` points by iterating over the `x`-coordinates to determine the `y` value, using a randomization graph (the seed), and a value. In math terms our function is `y(x, seed, value)`.\n\nA more practical formulation would be `value(x, y, seed)`. This would open up other possible uses. For example, you could use it for so-called [stateless shuffle iterators](https://www.ea.com/seed/news/constant-time-stateless-shuffling). Also, a `value(x, y, seed)` formulation should theoretically be much easier to implement as shader, but that is not really my specialty.","pinned":false,"hidden":false,"mode":"md"},{"id":5051,"value":"$3c","pinned":false,"hidden":false,"mode":"md"},{"id":5424,"value":"$3d","pinned":false,"hidden":false,"mode":"md"},{"id":3818,"value":"function makeValueGen(totalPoints, randomGraph) {\n const tpLog = Math.log2(totalPoints),\n mask = totalPoints - 1,\n offsetMask = mask << tpLog\n\n if (!randomGraph || randomGraph.length < totalPoints) {\n randomGraph = []; \n for (let i = 0; i < totalPoints / 2 * tpLog; i++) {\n randomGraph[i] = Math.random() < 0.5 ? 0 : 1;\n }\n };\n return function value(x, y) {\n const y0 = y;\n let ry = reverseNbits(y, tpLog);\n const effPath = x ^ ry;\n let randomBit, odd,\n yOffset, yMod,\n turnRight,\n value = 0;\n y = x;\n for (let row = 0; row < tpLog; row++) {\n odd = y & 1;\n yMod = y & (mask >> row);\n yOffset = y & (offsetMask >>> row);\n randomBit = randomGraph[row * totalPoints / 2 + (y >>> 1)];\n turnRight = odd ^ ((effPath >>> row) & 1);\n y = yOffset + (yMod >> 1) + (turnRight ? ((totalPoints / 2) >> row) : 0);\n value = (value << 1) ^ turnRight ^ randomBit ^ odd;\n }\n return value;\n };\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":3881,"value":"Look... I did say *relatively* easy at the start of the notebook, ok? Anyway, let's test it out with the zero-graph:","pinned":false,"hidden":false,"mode":"md"},{"id":3823,"value":"{\n const totalPoints = 256;\n const ctx = DOM.context2d(totalPoints, totalPoints, 1);\n ctx.canvas.style.width = ctx.canvas.style.height = dprAdjusted;\n ctx.canvas.className = \"pixelated\";\n yield ctx.canvas;\n await visibility();\n // blank random graph\n const randomGraph = [];\n for (let i = 0; i < 8 * 128; i++) {\n randomGraph[i] = 0;\n }\n const gen = makeValueGen(totalPoints, randomGraph);\n const palette = makeRamp(totalPoints, \"#grayscale\");\n for (let y = 0; y < totalPoints; y++) {\n for (let x = 0; x < totalPoints; x++) {\n const v = gen(x, y);\n ctx.fillStyle = palette[v];\n ctx.fillRect(x, y, 1, 1);\n }\n }\n\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":3884,"value":"Well, that's our previous result, just rotated 90 degrees (or maybe flipped along the diagonal axis? I can't tell, honestly). Now with a random graph:","pinned":false,"hidden":false,"mode":"md"},{"id":3896,"value":"$3e","pinned":false,"hidden":false,"mode":"ojs"},{"id":4741,"value":"Let's do a quick test to see if our generator produces properly stratified tiles.","pinned":false,"hidden":false,"mode":"md"},{"id":4733,"value":"$3f","pinned":false,"hidden":false,"mode":"ojs"},{"id":3927,"value":"With that out of the way, let's take a look at...\n\n## Scrambling the generated point set\n\nSo we established earlier that our noise generator is way too regular. We need to find a way to scramble the output, preferably in a way that keeps our stratified qualities.\n\nOne first thing to notice is that since a row or column uses every value exactly once, we could shuffle an array of the same values and use that to project each value in a point set to a new one. Because this is a reversible permutation this doesn't affect the stratified qualities.","pinned":false,"hidden":false,"mode":"md"},{"id":3931,"value":"restart7, makeTile({\n tilesize: 256,\n graphStructure: \"random\",\n shuffleColors: \"random\",\n colorscale: \"grayscale\",\n coprimeFactor: 0,\n graphCoPrimeFactor: 0,\n showTile: \"row\",\n visibility\n});","pinned":false,"hidden":false,"mode":"ojs"},{"id":4759,"value":"viewof restart7 = Inputs.button(\"Shuffle!\")","pinned":false,"hidden":false,"mode":"ojs"},{"id":2647,"value":"$40","pinned":false,"hidden":false,"mode":"md"},{"id":718,"value":"// See the following blog posts for more details:\n// https://blog.demofox.org/2020/07/26/irrational-numbers/\n// https://blog.demofox.org/2024/05/19/a-low-discrepancy-shuffle-iterator-random-access-inversion/\nfunction findCoprimeApprox(n, phi = 2/(1 + 5**0.5)) {\n // initiate by multiplying n by some number (default: golden ratio) phi and rounding to odd\n let coprime = n * phi | 1 >>> 0;\n if (coprime%1) coprime--\n while (gcd(n, coprime) !== 1 && coprime > 1) {\n coprime -= 2;\n }\n return coprime;\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":711,"value":"function gcd(a, b) {\n b = b|0;\n a = a|0\n if (!(a > 0 && b > 0)) return 0;\n let t = 0;\n // requires that neither parameter is negative\n while (b !== 0) {\n t = a;\n a = b;\n b = t % b;\n }\n return a;\n};","pinned":true,"hidden":false,"mode":"ojs"},{"id":4028,"value":"makeTile({\n tilesize: 256,\n graphStructure: \"random\",\n shuffleColors: \"coprime\",\n colorscale: \"grayscale\",\n coprimeFactor: 2/(1 + 5**0.5),\n showTile: \"row\",\n visibility\n});","pinned":true,"hidden":false,"mode":"ojs"},{"id":4039,"value":"That's a lot better! However, if we zoom in there are still some subtly discernible patterns, in the form of \"sawtooth\" gradients.\n\nHere's another thing we can try: when our texture has sides that are a power of four we are guaranteed to also have square stratified tiles (e.g. a 256×256 grid is guaranteed to be divisible into 16×16 tiles, each containing 16×16 points. Each of those tiles contains the total 256 values exactly once). Instead of randomly shuffling the value projection, we could also map all the values in such a way that one of the tiles, say the bottom left produces a gradient, and see what effect that has on the rest of the texture:","pinned":false,"hidden":false,"mode":"md"},{"id":2668,"value":"makeTile({\n tilesize: 256,\n graphStructure: \"ascending\",\n showTile: \"gradient\",\n visibility\n});","pinned":true,"hidden":false,"mode":"ojs"},{"id":4105,"value":"This is highly regular again, but it's a starting point for some other strategies. For example, what if we apply our coprime-shuffle trick to this tile instead of to a row or column?","pinned":false,"hidden":false,"mode":"md"},{"id":4108,"value":"makeTile({\n tilesize: 256,\n graphStructure: \"random\",\n shuffleColors: \"coprime\",\n coprimeFactor: 2/(1+5**0.5),\n showTile: \"gradient\",\n visibility\n});","pinned":true,"hidden":false,"mode":"ojs"},{"id":2680,"value":"Not much better, and now we see a clear strip pattern in the bottom left tile.\n\nOk, what about mapping that tile using the Hilbert curve instead, or something like the [H-curve](https://observablehq.com/@jobleonard/a-simple-algorithm-to-generate-h-curves), or the [HI and pinwheel variations](https://observablehq.com/@jobleonard/the-hi-curve?collection=%40jobleonard%2Fspace-filling-curves):","pinned":false,"hidden":false,"mode":"md"},{"id":2672,"value":"makeTile({\n tilesize: 256,\n graphStructure: \"random\",\n showTile: \"Hilbert\",\n visibility\n});","pinned":true,"hidden":false,"mode":"ojs"},{"id":3975,"value":"makeTile({\n tilesize: 256,\n graphStructure: \"random\",\n showTile: \"H-curve pinwheel\",\n visibility\n});","pinned":true,"hidden":false,"mode":"ojs"},{"id":4102,"value":"Cool, but not what we're looking for.\n\nWhat if we apply our coprime shuffling to a line of values projected by these space-filling curves instead?","pinned":false,"hidden":false,"mode":"md"},{"id":4064,"value":"makeTile({\n tilesize: 256,\n graphStructure: \"random\",\n shuffleColors: \"coprime\",\n coprimeFactor: 2/(1+5**0.5),\n showTile: \"Hilbert\",\n visibility\n});","pinned":true,"hidden":false,"mode":"ojs"},{"id":4123,"value":"Getting there, but I can still spot vertical \"edges\" in the noise.\n\nOk, how about reversing the bits of our space-filling curves? You see, reversing a bit-string is a bit like the 1D equivalent of applying [a Bayer index matrix](https://en.wikipedia.org/wiki/Ordered_dithering#Threshold_map), which is probably easier to intuit if I just show you:","pinned":false,"hidden":false,"mode":"md"},{"id":4767,"value":"$41","pinned":false,"hidden":false,"mode":"ojs"},{"id":4794,"value":"$42","pinned":false,"hidden":false,"mode":"md"},{"id":4798,"value":"makeTile({\n tilesize: 256,\n graphStructure: \"ascending\",\n showTile: \"Hilbert reverse\",\n visibility\n});","pinned":true,"hidden":false,"mode":"ojs"},{"id":4802,"value":"makeTile({\n tilesize: 256,\n graphStructure: \"ascending\",\n showTile: \"H-curve pinwheel reverse\",\n visibility\n});","pinned":true,"hidden":false,"mode":"ojs"},{"id":4807,"value":"Well, it sort of does, but also produces regular patterns again. Let's try applying a randomized graph:","pinned":false,"hidden":false,"mode":"md"},{"id":4125,"value":"makeTile({\n tilesize: 256,\n graphStructure: \"random\",\n showTile: \"Hilbert reverse\",\n visibility\n});","pinned":true,"hidden":false,"mode":"ojs"},{"id":4127,"value":"makeTile({\n tilesize: 256,\n graphStructure: \"random\",\n showTile: \"H-curve reverse\",\n visibility\n});","pinned":true,"hidden":false,"mode":"ojs"},{"id":4249,"value":"Still blocky patterns. Let's try is combining all tricks together: first we assign values to a subtile using a space-filling curves with the reversed bits. Then we use the \"irrational coprime\" trick to shuffle those values to another position in that subtile:","pinned":false,"hidden":false,"mode":"md"},{"id":4132,"value":"makeTile({\n tilesize: 256,\n graphStructure: \"random\",\n shuffleColors: \"coprime\",\n coprimeFactor: 1/2**0.5,\n showTile: \"H-curve pinwheel reverse\",\n visibility\n});","pinned":true,"hidden":false,"mode":"ojs"},{"id":4251,"value":"There we go! Something noisy looking that still follows our stratified constraints! Compare our results to white noise and how lumpy that output is:","pinned":false,"hidden":false,"mode":"md"},{"id":4254,"value":"{\n const size = 256;\n const ctx = DOM.context2d(size, size);\n ctx.canvas.style.width = ctx.canvas.style.height = dprAdjusted;\n ctx.canvas.className = \"pixelated\";\n const palette = makeRamp(size, \"#grayscale\");\n const {random} = Math;\n for (let x = 0; x < size; x++) {\n for (let y = 0; y < size; y++) {\n ctx.fillStyle = palette[random()*palette.length | 0];\n ctx.fillRect(x, y, 1, 1);\n }\n }\n return ctx.canvas;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":5430,"value":"## TODO: write aside for how to combine this with `value(x, y, seed)`","pinned":false,"hidden":false,"mode":"md"},{"id":5593,"value":"## TODO: try \"recursive\" scrambling\n\nSo far, a tile of `n²` points only has `n` distinct values. Could we actually turn that into `n²` distinct values?\n\n- if `n` is a power of 4 (meaning there is a stratification with square tiles with sides of length `n**(1/2)`) we could try generating a tile of `n**(1/2)` distinct points per `pmj02` sequence and use that assign each point in that sequence a unique value\n- How do we generate an `n**(1/2)` sized tile? Why, do the same thing there but `n**(1/4)` sized tiles, all the way down to 2x2 sized tiles, which will be distributed randomly\n- being able to generate textures with `n²` unique points would also let us use them for assigning an increasing value to each sequence","pinned":false,"hidden":false,"mode":"md"},{"id":4321,"value":"## TODO: something resembling a conclusion\n\n*Many thanks to Fabian Iwand and Alan Wolfe for their feedback on the draft version of this notebook - it is much better thanks to their proofreading and tips.*\n\n

","pinned":false,"hidden":false,"mode":"md"},{"id":5523,"value":"dprAdjusted = `${2**(Math.log2(Math.min(512, width))|0)}px`","pinned":false,"hidden":false,"mode":"ojs"},{"id":104,"value":"$43","pinned":true,"hidden":false,"mode":"ojs"},{"id":343,"value":"$44","pinned":true,"hidden":false,"mode":"ojs"},{"id":392,"value":"$45","pinned":true,"hidden":false,"mode":"ojs"},{"id":292,"value":"// import {hilbert} from \"@mourner/hilbert-curve-packing\"\nhilbert = (x, y, s) => {\n if (!s) s = 256;\n let d = 0;\n for (s = s*s >> 1; s > 0; s >>= 1) {\n const rx = +((x & s) > 0);\n const ry = +((y & s) > 0);\n d += s * s * ((3 * rx) ^ ry);\n\n if (ry === 0) {\n if (rx === 1) {\n x = s - 1 - x;\n y = s - 1 - y;\n }\n let t = x;\n x = y;\n y = t;\n }\n }\n return d;\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":786,"value":"pinwheel = ({\n 4: pinwheelIndices(0),\n 8: pinwheelIndices(1),\n 16: pinwheelIndices(2),\n 32: pinwheelIndices(3),\n 64: pinwheelIndices(4),\n 128: pinwheelIndices(5),\n 256: pinwheelIndices(6)\n})","pinned":true,"hidden":false,"mode":"ojs"},{"id":4827,"value":"hcurveIndices(5)","pinned":true,"hidden":false,"mode":"ojs"},{"id":821,"value":"hcurve = ({\n 4: hcurveIndices(0),\n 8: hcurveIndices(1),\n 16: hcurveIndices(2),\n 32: hcurveIndices(3),\n 64: hcurveIndices(4),\n 128: hcurveIndices(5),\n 256: hcurveIndices(6)\n})","pinned":true,"hidden":false,"mode":"ojs"},{"id":784,"value":"import {pinwheelIndices} from \"@jobleonard/the-hi-curve\"","pinned":true,"hidden":false,"mode":"ojs"},{"id":818,"value":"import {hcurveIndices} from \"@jobleonard/a-simple-algorithm-to-generate-h-curves\"","pinned":true,"hidden":false,"mode":"ojs"},{"id":11,"value":"$46","pinned":true,"hidden":false,"mode":"ojs"},{"id":109,"value":"/* \n * Actually, this set of fallbacks is probably overkill\n * (does Observable even work on IE8?) but it won't\n * harm anyone so I didn't bother to simplify it\n */\npixelatedCSS = html``","pinned":true,"hidden":false,"mode":"ojs"},{"id":2069,"value":"function reverseNbits(v, n) {\n switch(n) {\n case 32:\n v = (v >>> 16) + (v << 16);\n case 16:\n v = ((v & 0xFF00FF00) >>> 8) + ((v & 0xFF00FF) << 8);\n case 8:\n v = ((v & 0xF0F0F0F0) >>> 4) + ((v & 0xF0F0F0F) << 4);\n case 4:\n v = ((v & 0b11001100110011001100110011001100) >>> 2) + ((v & 0b110011001100110011001100110011) << 2);\n case 2:\n v = ((v & 0b10101010101010101010101010101010) >>> 1) + ((v & 0b1010101010101010101010101010101) << 1);\n return v;\n default: {\n let rv = 0;\n for (let i = 0; i < n; i++) {\n rv += rv + ((v >>> i)&1);\n }\n return rv;\n }\n }\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":1854,"value":" function blueImgToGraph(img){\n const {width: w, height: h} = img;\n const ctx = DOM.context2d(w, h);\n ctx.drawImage(img, 0, 0, w, h);\n const imageData = ctx.getImageData(0, 0, w, h);\n const graph = [];\n for (let i = 0; i < imageData.data.length; i += 4) {\n graph.push(imageData.data[i] ? 0 : 1);\n }\n return graph\n }\n","pinned":true,"hidden":false,"mode":"ojs"},{"id":1848,"value":"blueGraphs = ({\n 16: blueImgToGraph(await FileAttachment(\"0016.png\").image()),\n 32: blueImgToGraph(await FileAttachment(\"0032.png\").image()),\n 64: blueImgToGraph(await FileAttachment(\"0064.png\").image()),\n 128: blueImgToGraph(await FileAttachment(\"0128.png\").image()),\n 256: blueImgToGraph(await FileAttachment(\"0256.png\").image()),\n 512: blueImgToGraph(await FileAttachment(\"0512.png\").image()),\n 1024: blueImgToGraph(await FileAttachment(\"1024.png\").image()),\n 2048: blueImgToGraph(await FileAttachment(\"2048.png\").image()),\n 4096: blueImgToGraph(await FileAttachment(\"4096.png\").image())\n}) ","pinned":true,"hidden":false,"mode":"ojs"},{"id":1905,"value":"$47","pinned":false,"hidden":false,"mode":"ojs"},{"id":6331,"value":"function lerp(a, b, t) {\n return a * (1 - t) + b * t;\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":6330,"value":"function easedLine(ctx, x0, y0, x1, y1, t0 = 1 / 4, t1 = 1 - t0) {\n const dx = x1 - x0;\n const dy = y1 - y0;\n ctx.beginPath();\n ctx.moveTo(x0, y0);\n // Thanks for reminding me bezierCurveTo exists, Fabian :)\n ctx.bezierCurveTo(x0, lerp(y0, y1, t0), x1, lerp(y0, y1, t1), x1, y1);\n ctx.stroke();\n}","pinned":true,"hidden":false,"mode":"ojs"}],"files":"$W48","stdlib":"1","titling":"auto"}