31:T57b,{ const w = 240; const h = 150; const { ctx, imgData } = getCtx(w, h, Math.min(800, width)); // Background let gradient = ctx.createLinearGradient(0, 0, w, 0); gradient.addColorStop(0, "black"); gradient.addColorStop(1, "#444"); ctx.fillStyle = gradient; ctx.fillRect(0, 0, w, h); ctx.font = `${((h * 2) / 9) | 0}px serif`; ctx.textBaseline = "middle"; gradient = ctx.createLinearGradient(w / 8, 0, (w * 7) / 8, 0); gradient.addColorStop(0, "#FFF"); gradient.addColorStop(1, "#888"); ctx.fillStyle = gradient; ctx.textAlign = "left"; ctx.fillText("Pseudo", w / 8, h / 4); ctx.fillText("Blue", (w - ctx.measureText("Noise").width) / 2, h / 2); ctx.textAlign = "right"; ctx.fillText("Noise", (w * 7) / 8, (h * 3) / 4); const title = ctx.getImageData(0, 0, w, h); while (true) { const indices = noise(w, h); for (let y = 0; y < h; y++) { for (let x = 0; x < w; x++) { const idx = (x + y * w) * 4; const treshold = (indices[idx >> 2] * 256) / indices.length; const pixel = title.data[idx] < treshold ? 0 : 255; imgData.data[idx] = pixel; imgData.data[idx + 1] = pixel; imgData.data[idx + 2] = pixel; imgData.data[idx + 3] = 255; // make sure image is not transparent } } ctx.putImageData(imgData, 0, 0); await visibility(); yield Promises.tick(100, ctx.canvas); } }32:T5a8,md`Blue noise is a noise color[[1](https://en.wikipedia.org/wiki/Colors_of_noise)] that is particularly well-suited for image dithering purposes (see also [[2](https://bartwronski.com/2016/10/30/dithering-part-three-real-world-2d-quantization-dithering/)][[3](https://momentsingraphics.de/BlueNoise.html)]). Unfortunately it's also computationally intensive to generate. For example, Bart Wronski mentions that generating a *"64×64 wrapping blue noise-like sequence"* took *"a couple hours on an old MacBook."* This was in 2016, but unless dramatic progress has been made in the algorithms to generate blue noise this will not be significantly faster today. What if instead of trying to generate *proper* blue noise, we look at the features that make it useful for image dithering? Then we can try to come up with a cheaper noise function that also (mostly) has these features. ## Understanding the problem we are trying to solve Let's start at the beginning: we are talking about using blue noise for the purpose of image dithering. To quote Wikipedia, image dithering *"is used in computer graphics to create the illusion of color depth in images on systems with a limited color palette"*[[4](https://en.wikipedia.org/wiki/Dither#Digital_photography_and_image_processing)]. For example, say we have only black and white pixels (or printer ink) available to us as output, and want to represent this grayscale ramp as accurately as possible:`33:T93b,md`## How the Bayer index matrix spreads out pixels so evenly ${render(width >> 2, 15, (x, y) => toU32Gray(get8x8threshold(x, y) < (x * 256) / (width >> 2) ? 0 : 255), Math.min(width, 800) )} Let's start with taking a deeper look at the logic behind classic ordered dithering. This is based on Bayer index matrices, which were introduced in 1973 by Bryce Bayer[[7](https://web.archive.org/web/20130512190753/http://white.stanford.edu/~brian/psy221/reader/Bayer.1973.pdf)] (hence the name). Imagine we have a 2 × 2 grid arrangement, so a total of four pixels. What is the most evenly spread out order in which we could flip the pixels of this grid from white to black? Bayer concluded it must be the following: ${render(16, 4, (x, y) => toU32Gray(255 - [ [1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1], [1,0,0,1,4,0,1,4,0,1,4,4,1,4,4,1], [1,0,0,1,0,0,1,0,4,1,0,4,1,4,4,1], [1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1] ][y][x] * 63))} If we flip one pixel to black it will be in one of the corners. For the sake of evenly distributing points doesn't actually matter which corner for now (we will get back to this, however). But when we flip the *second* pixel, and if we want that to be as evenly spread in our grid as possible, then we *must* flip the pixel in the opposite corner. Seems pretty obvious, right? The order in which we fill the third and last pixels doesn't seem to matter again, but like the first pixel there is a caveat: it is true that when we look at the 2 × 2 grid *in isolation* this doesn't matter. However, if we are tiling a space with 2 × 2 grids, then the most even results will be achieved if we always follow the same order in each tile (in this notebook we'll start in the top-left and end in the bottom right). What Bayer realized was that we can apply this logic recursively: if we have, say, a 8 × 8 grid of pixels, we can divide it into four 4 × 4 sub-grids. We first place a pixel in the top-left sub-grid, then the bottom-right sub-grid, and so on. But then the question is: wher to place those pixels in each 4 × 4 sub-grid? Well, we can divide each of the sub-grids into four 2 × 2 sub-grids. Then we place the first pixel in the top left 2 × 2, and so on. And then we are back at our original grid size so we know how to place those pixels. So to determine the pixel position we just keep dividing all the way down:`34:T609,Did you notice that each sub-grid has a sequence of numbers increasing in steps of four, but that they start a step apart? And that those numbers are again divided in this "top left, bottom right, top right, bottom left" pattern? ```text full grid: 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, top-left: 0, 4, 8, 12 bottom-right: 1, 5, 9, 13 top-right: 2, 6, 10, 14 bottom-left: 3, 7, 11, 15 ``` That is the ingenious part of doing this division recursively: at the big scale we "spread out" our available indices as evenly as possible into four grids, then for each sub-grid we divide the new sub-set of numbers evenly again, ensuring that this even spread is retained all the way down to the smallest sub-grid. So to recap: Bayer's approach divides a sequence of numbers, e.g. *[0, 1, 2, ..., K]*, into four new sequences, *[0, 4, 8, ..., K-3], [1, 5, 9, ..., K-2], [2, 6, 10, ..., K-1], [3, 7, 11, ..., K]*, passing those four sequences to the sub-grids, first top-left, then bottom-right, then top-right, and finally bottom-left. Each of these sub-grids then repeats this process, dividing each of their sequences into four smaller sequences again and multiplying the step size by four as a result as well. It then repeats this until we only have one number to place on a 1 × 1 grid, giving us the final position of that index in the grid as a whole:35:T8f0,{ // Like the above split function, but with a max recursion depth for the data visualization function split(stride, indices, step) { if (step-- <= 0) return indices; let topLeft = [], topRight = [], bottomLeft = [], bottomRight = []; // Divide indices across four sections in Bayer-order for (let i = 0; i < indices.length; i++) { const idx = indices[i]; switch (i % 4) { case 0: topLeft.push(idx); continue; case 1: bottomRight.push(idx); continue; case 2: topRight.push(idx); continue; case 3: bottomLeft.push(idx); continue; } } // Recursively subdivide indices in each quadrant const stride2 = stride / 2; topLeft = split(stride2, topLeft, step); topRight = split(stride2, topRight, step); bottomLeft = split(stride2, bottomLeft, step); bottomRight = split(stride2, bottomRight, step); // Merge back indices. We merge each quadrant row by row const output = []; // First the top two quadrants for (let y = 0; y < stride2; y++) { const line = y * stride2; for (let x = 0; x < stride2; x++) { output.push(topLeft[x + line]); } for (let x = 0; x < stride2; x++) { output.push(topRight[x + line]); } } // Then the bottom quadrants for (let y = 0; y < stride2; y++) { const line = y * stride2; for (let x = 0; x < stride2; x++) { output.push(bottomLeft[x + line]); } for (let x = 0; x < stride2; x++) { output.push(bottomRight[x + line]); } } return output; } const length = 4 ** 4; const rootIndices = []; for (let i = 0; i < length; i++) rootIndices[i] = i; const finalIndices = split(2 ** 4, rootIndices, maxBayerRecursionSteps); const { ctx, imgData } = getCtx(16, 16, 512); for (let y = 0; y < 16; y++) { for (let x = 0; x < 16; x++) { const idx = x + y * 16; const v = (finalIndices[idx] * 256) / length; imgData.data[idx * 4] = v; imgData.data[idx * 4 + 1] = v; imgData.data[idx * 4 + 2] = v; imgData.data[idx * 4 + 3] = 255; } } ctx.putImageData(imgData, 0, 0); return ctx.canvas; }36:T753,// This is a *very* inefficient approach, // but it makes it easier to follow the steps (I hope) function makeBayerIndex(recursionSteps) { function split(stride, indices) { if (stride === 1) return indices; let topLeft = [], topRight = [], bottomLeft = [], bottomRight = []; // Divide indices across four sections in Bayer-order for (let i = 0; i < indices.length; i++) { const idx = indices[i]; switch (i % 4) { case 0: topLeft.push(idx); continue; case 1: bottomRight.push(idx); continue; case 2: topRight.push(idx); continue; case 3: bottomLeft.push(idx); continue; } } // Recursively subdivide indices in each quadrant const stride2 = stride / 2; topLeft = split(stride2, topLeft); topRight = split(stride2, topRight); bottomLeft = split(stride2, bottomLeft); bottomRight = split(stride2, bottomRight); // Merge back indices. We merge each quadrant row by row const output = []; // First the top two quadrants for (let y = 0; y < stride2; y++) { const line = y * stride2; for (let x = 0; x < stride2; x++) { output.push(topLeft[x + line]); } for (let x = 0; x < stride2; x++) { output.push(topRight[x + line]); } } // Then the bottom quadrants for (let y = 0; y < stride2; y++) { const line = y * stride2; for (let x = 0; x < stride2; x++) { output.push(bottomLeft[x + line]); } for (let x = 0; x < stride2; x++) { output.push(bottomRight[x + line]); } } return output; } const length = 4 ** recursionSteps; const rootIndices = []; for (let i = 0; i < length; i++) rootIndices[i] = i; return split(2 ** recursionSteps, rootIndices); }37:Ta11,## Arbitrary Bayer indices What if I don't want to have an index matrix with sides that are a power of two? What if I don't want them to be *equal* sides? Can we still make this work? The answer is "sort of", but we will not get a perfectly even and regular grid out of it. But that is acceptable, since we are not aiming for that when creating noise anyway! To make it easier to follow how this works we will simplify the problem to a 1D "matrix" for now (so instead of looking at a 2D surface area, we look at a 1D line). Our recursive function then simplifies to dividing up our indices in a left and right section in alternating fashion. So let's say we have an eight pixel line. In that case we start with these linear indices: [0, 1, 2, 3, 4, 5, 6, 7] In the first call, we split these into two sets of indices: [[0, 2, 4, 6], [1, 3, 5, 7]] These two sets of indices then are passed to two recursive calls, then each get further split into two sets of alternating indices again: [[[0, 4], [2, 6]], [[1, 5], [3, 7]]] After all recursion is done and our function is finished, our final ordering becomes: [0, 4, 2, 6, 1, 5, 3, 7] The problem is that when our line length (or later side length) is not a neat power of two, we cannot repeatedly divide it into two without coming across an odd number at some point. Which means that at that point we will have an uneven left and right half. For example, if we wish to divide ten pixels: [0, 1, 2, 3, 4, 5, 6, 7, 8, 9] Then after the first recursive step we have: [[0, 2, 4, 6, 8], [1, 3, 5, 7, 9]] Now how do we divide these sets of five indices into two subsets? Obviously, one set should get three indices, and the other two. We also want our algorithm to alternate between the two as evenly as possible. The trick to achieve that is using fractions to decide which sub-division gets the next index. Each half gets a counter and a fraction proportional to their relative size. Then we repeatedly add the fractions to their respective counters, one by one. Whenever a counter has increased by a whole number we assign an index to the half it belongs to. This will obviously happen more often for the larger fraction. We repeat this until we run out of indices to assign. If you work out the math, you will see that this ensures both halves get the right number of indices, and in an alternating fashion as well (this is essentially the logic used in Bresenham's line algorithm[[8](https://en.wikipedia.org/wiki/Bresenham%27s_line_algorithm)], but repurposed for dividing our indices).38:T7a7,function split1D(length) { const rootIndices = []; function splitHalf(indices) { // split in half const left = Math.floor(indices.length / 2); const right = indices.length - left; let leftIndices = [], rightIndices = []; // We use Bresenham line algorithm-style error accumulators // (in other words: fractions) to decide when to divide an // index to the left half and when to the right half. let lNumerator = 0, rNumerator = 0; // We make the denominator the *smallest* of the two halves. // See below for an explanation why const denominator = right < left ? right : left; let i = 0; for (let i = 0; i < indices.length; ) { // grow our fractions. Because the denominator is the // smallest of the two, this is effectively ensuring // that the smaller fraction is equal to "one" every // loop, and the larger one... well, whatever the // proportion of large half : small half is lNumerator += left; rNumerator += right; // Now for both fractions we pass on indices to their // respective half for as long as they are "bigger than one", // which is guaranteed to be true at least once due to // using the smallest of the numerators as the denominator). do { lNumerator -= denominator; leftIndices.push(indices[i++]); } while (lNumerator >= denominator); do { rNumerator -= denominator; rightIndices.push(indices[i++]); } while (rNumerator >= denominator); } // And of course we repeat this proces recursively on both indices if (leftIndices.length > 2) leftIndices = splitHalf(leftIndices); if (rightIndices.length > 2) rightIndices = splitHalf(rightIndices); // copy divided indices back from buffer return [...leftIndices, ...rightIndices]; } for (let i = 0; i < length; i++) rootIndices[i] = i; return splitHalf(rootIndices); }39:T7b9,{ const w = Math.min(width, 1024); const h = (w / 8) | 0; const { ctx, imgData } = getCtx(w, h, 1024); // Like the above split function, but with a max recursion depth for the data visualization function split1DMax(length, maxRecursion = length) { const rootIndices = []; function splitHalf(indices, depth) { if (depth >= maxRecursion) return indices; // split in half const left = Math.floor(indices.length / 2); const right = indices.length - left; let leftIndices = [], rightIndices = []; let lNumerator = 0, rNumerator = 0; const denominator = right < left ? right : left; let i = 0; for (let i = 0; i < indices.length; ) { lNumerator += left; rNumerator += right; do { lNumerator -= denominator; leftIndices.push(indices[i++]); } while (lNumerator >= denominator); do { rNumerator -= denominator; rightIndices.push(indices[i++]); } while (rNumerator >= denominator); } // And of course we repeat this proces recursively on both indices if (leftIndices.length > 2) leftIndices = splitHalf(leftIndices, depth + 1); if (rightIndices.length > 2) rightIndices = splitHalf(rightIndices, depth + 1); return [...leftIndices, ...rightIndices]; } for (let i = 0; i < length; i++) rootIndices[i] = i; return splitHalf(rootIndices, 0); } const finalIndices = split1DMax(bayerLength1D, maxBayer1DRecursionSteps); for (let y = 0; y < h; y++) { for (let x = 0; x < w; x++) { const v = ((finalIndices[((x * bayerLength1D) / w) | 0] * 255) / (bayerLength1D - 1)) | 0; const idx = (x + y * w) * 4; imgData.data[idx] = v; imgData.data[idx + 1] = v; imgData.data[idx + 2] = v; imgData.data[idx + 3] = 255; } } ctx.putImageData(imgData, 0, 0); return ctx.canvas; }3a:T7cd,{ const w = Math.min(width, 1024); const h = (w / 8) | 0; const { ctx, imgData } = getCtx(w, h, 1024); // Like the above split function, but with a max recursion depth for the data visualization function split1DMax(length, maxRecursion = length) { const rootIndices = []; function splitHalf(indices, depth) { if (depth >= maxRecursion) return indices; // split in half const left = Math.floor(indices.length / 2); const right = indices.length - left; let leftIndices = [], rightIndices = []; let lNumerator = 0, rNumerator = 0; const denominator = right < left ? right : left; let i = 0; for (let i = 0; i < indices.length; ) { lNumerator += left; rNumerator += right; do { lNumerator -= denominator; leftIndices.push(indices[i++]); } while (lNumerator >= denominator); do { rNumerator -= denominator; rightIndices.push(indices[i++]); } while (rNumerator >= denominator); } // And of course we repeat this proces recursively on both indices if (leftIndices.length > 2) leftIndices = splitHalf(leftIndices, depth + 1); if (rightIndices.length > 2) rightIndices = splitHalf(rightIndices, depth + 1); return [...leftIndices, ...rightIndices]; } for (let i = 0; i < length; i++) rootIndices[i] = i; return splitHalf(rootIndices, 0); } const finalIndices = split1DMax(bayerLength1D, maxBayer1DRecursionSteps); for (let y = 0; y < h; y++) { for (let x = 0; x < w; x++) { const v = finalIndices[((x * bayerLength1D) / w) | 0] / (bayerLength1D - 1) < (h - y) / h ? 0 : 255; const idx = (x + y * w) * 4; imgData.data[idx] = v; imgData.data[idx + 1] = v; imgData.data[idx + 2] = v; imgData.data[idx + 3] = 255; } } ctx.putImageData(imgData, 0, 0); return ctx.canvas; }3b:T84c,function split2D(width, height) { const { random } = Math; // Instead of creating new arrays all the time I ping-pong // between two arrays of large enough size, and use a length // plus offset to divide indices into smaller parts function splitHalf(w, h, source, target, offset) { const left = w >> 1; const right = w - left; let lNumerator = 0, rNumerator = 0; const denominator = right < left ? right : left; let i = offset, j = offset + left * h, k = offset; const end = offset + w * h; while (k < end) { lNumerator += left; rNumerator += right; while (lNumerator >= denominator) { target[i++] = source[k++]; lNumerator -= denominator; } while (rNumerator >= denominator) { target[j++] = source[k++]; rNumerator -= denominator; } } return left; } function divide(x, y, w, h, source, target, offset) { if (w === 1 && h === 1) { indices[x + y * width] = source[offset]; } else if (w > h) { const left = splitHalf(w, h, source, target, offset); const right = w - left; // note that we swap the target and source in the recursive call divide(x, y, left, h, target, source, offset); divide(x + left, y, right, h, target, source, offset + left * h); } else { const top = splitHalf(h, w, source, target, offset); const bottom = h - top; divide(x, y, w, top, target, source, offset); divide(x, y + top, w, bottom, target, source, offset + top * w); } } const indices = new Int32Array(width * height); // Use a single arraybuffer for better cache locality // (probably premature optimization, but also harmless) // See also: https://devdocs.io/javascript/global_objects/typedarray/subarray const buffer = new Int32Array(width * height * 2); const source0 = buffer.subarray(0, width * height); const target0 = buffer.subarray(width * height, width * height * 2); for (let i = 0; i < width * height; i++) source0[i] = i; divide(0, 0, width, height, source0, target0, 0); return indices; }3c:Ta1a,function split2D_v2(width, height) { const { random } = Math; const flip = () => random() < 0.5; // Instead of creating new arrays all the time I ping-pong // between two arrays of large enough size, and use a length // plus offset to divide indices into smaller parts function splitHalf(w, h, source, target, offset) { const left = w >> 1; const right = w - left; // Note that we set the denominator to the bigger half, // as this increases randomization of the partitions. let lNumerator = 0, rNumerator = 0; const denominator = right < left ? left : right; let i = offset, j = offset + left * h, k = offset; const end = offset + w * h; while (k < end) { lNumerator += left; rNumerator += right; if (flip()) { if (lNumerator >= denominator) { target[i++] = source[k++]; lNumerator -= denominator; } if (rNumerator >= denominator) { target[j++] = source[k++]; rNumerator -= denominator; } } else { if (rNumerator >= denominator) { target[j++] = source[k++]; rNumerator -= denominator; } if (lNumerator >= denominator) { target[i++] = source[k++]; lNumerator -= denominator; } } } return left; } function divide(x, y, w, h, source, target, offset) { if (w === 1 && h === 1) { indices[x + y * width] = source[offset]; } else if (w > h) { const left = splitHalf(w, h, source, target, offset); const right = w - left; // note that we swap the target and source in the recursive call divide(x, y, left, h, target, source, offset); divide(x + left, y, right, h, target, source, offset + left * h); } else { const top = splitHalf(h, w, source, target, offset); const bottom = h - top; divide(x, y, w, top, target, source, offset); divide(x, y + top, w, bottom, target, source, offset + top * w); } } const indices = new Int32Array(width * height); // Use a single arraybuffer for better cache locality // (probably premature optimization, but also harmless) // See also: https://devdocs.io/javascript/global_objects/typedarray/subarray const buffer = new Int32Array(width * height * 2); const source0 = buffer.subarray(0, width * height); const target0 = buffer.subarray(width * height, width * height * 2); for (let i = 0; i < width * height; i++) source0[i] = i; divide(0, 0, width, height, source0, target0, 0); return indices; }3d:Ta6d,function split2D_v3(width, height) { const { random } = Math; const flip = () => random() < 0.5; // Instead of creating new arrays all the time I ping-pong // between two arrays of large enough size, and use a length // plus offset to divide indices into smaller parts function splitHalf(w, h, source, target, offset) { let left = w >> 1; let right = w - left; // randomize which side we round to if (flip()) [left, right] = [right, left]; // Note that we set the denominator to the bigger half, // as this increases randomization of the partitions. let lNumerator = 0, rNumerator = 0; const denominator = right < left ? left : right; let i = offset, j = offset + left * h, k = offset; const end = offset + w * h; while (k < end) { lNumerator += left; rNumerator += right; if (flip()) { if (lNumerator >= denominator) { target[i++] = source[k++]; lNumerator -= denominator; } if (rNumerator >= denominator) { target[j++] = source[k++]; rNumerator -= denominator; } } else { if (rNumerator >= denominator) { target[j++] = source[k++]; rNumerator -= denominator; } if (lNumerator >= denominator) { target[i++] = source[k++]; lNumerator -= denominator; } } } return left; } function divide(x, y, w, h, source, target, offset) { if (w === 1 && h === 1) { indices[x + y * width] = source[offset]; } else if (w > h) { const left = splitHalf(w, h, source, target, offset); const right = w - left; // note that we swap the target and source in the recursive call divide(x, y, left, h, target, source, offset); divide(x + left, y, right, h, target, source, offset + left * h); } else { const top = splitHalf(h, w, source, target, offset); const bottom = h - top; divide(x, y, w, top, target, source, offset); divide(x, y + top, w, bottom, target, source, offset + top * w); } } const indices = new Int32Array(width * height); // Use a single arraybuffer for better cache locality // (probably premature optimization, but also harmless) // See also: https://devdocs.io/javascript/global_objects/typedarray/subarray const buffer = new Int32Array(width * height * 2); const source0 = buffer.subarray(0, width * height); const target0 = buffer.subarray(width * height, width * height * 2); for (let i = 0; i < width * height; i++) source0[i] = i; divide(0, 0, width, height, source0, target0, 0); return indices; }3e:Tbe7,function split2D_v4(width, height) { const { random, round } = Math; const flip = () => random() < 0.5; const indices = new Int32Array(width * height); // Instead of creating new arrays all the time I ping-pong // between two arrays of large enough size, and use a length // plus offset to divide indices into smaller parts // Use a single arraybuffer for better cache locality // (probably premature optimization, but also harmless) // See also: https://devdocs.io/javascript/global_objects/typedarray/subarray const buffer = new Int32Array(width * height * 2); const source0 = buffer.subarray(0, width * height); const target0 = buffer.subarray(width * height, width * height * 2); for (let i = 0; i < width * height; i++) source0[i] = i; function splitHalf(w, h, source, target, offset) { let left, right; if (w > 3) { left = 3*w/7 + random() * (w/7 + 1) | 0; right = w - left; } else { left = w >> 1; right = w - left; } // randomize which side we round to if (flip()) [left, right] = [right, left]; // Note that we set the denominator to the bigger half, // as this increases randomization of the partitions. let lNumerator = 0, rNumerator = 0; const denominator = right < left ? left : right; let i = offset, j = offset + left * h, k = offset; const end = offset + w * h; while (k < end) { lNumerator += left; rNumerator += right; // randomize tie-breaks if (flip()) { if (lNumerator >= denominator) { target[i++] = source[k++]; lNumerator -= denominator; } if (rNumerator >= denominator) { target[j++] = source[k++]; rNumerator -= denominator; } } else { if (rNumerator >= denominator) { target[j++] = source[k++]; rNumerator -= denominator; } if (lNumerator >= denominator) { target[i++] = source[k++]; lNumerator -= denominator; } } } return left; } function splitQuart(w, h, source, target, offset) { } function divide(x, y, w, h, source, target, offset) { if (w === 1 && h === 1) { indices[x + y * width] = source[offset]; } else if (w > h) { const left = splitHalf(w, h, source, target, offset); const right = w - left; // note that we swap the target and source in the recursive call divide(x, y, left, h, target, source, offset); divide(x + left, y, right, h, target, source, offset + left * h); } else { // if (h >= 2 * w) { // note that we "rotate" splitHalf to work on heigh instead // of width simply by swapping w and h when calling it. const top = splitHalf(h, w, source, target, offset); const bottom = h - top; divide(x, y, w, top, target, source, offset); divide(x, y + top, w, bottom, target, source, offset + top * w); } } divide(0, 0, width, height, source0, target0, 0); return indices; }3f:T45a,### References [1] [*Colors of noise*](https://en.wikipedia.org/wiki/Colors_of_noise), Wikipedia [2] [*Dithering in games pt. 2*](https://bartwronski.com/2016/10/30/dithering-part-two-golden-ratio-sequence-blue-noise-and-highpass-and-remap/) by Bart Wronski, 2016 [3] [*Free blue noise textures*](https://momentsingraphics.de/BlueNoise.html) by Christoph Peters, 2016 [4] [*Dither - Digital photography and image processing*](https://en.wikipedia.org/wiki/Dither#Digital_photography_and_image_processing), Wikipedia [5] [*Ordered dithering*](https://en.wikipedia.org/wiki/Ordered_dithering), Wikipedia [6] [*Next Generation Post Processing in Call of Duty: Advanced Warfare*](http://www.iryoku.com/next-generation-post-processing-in-call-of-duty-advanced-warfare) by Jorge Jimenez, 2014 [7] [*An optimum method for two-level rendition of continuous-tone pictures*](https://web.archive.org/web/20130512190753/http://white.stanford.edu/~brian/psy221/reader/Bayer.1973.pdf) (PDF link), Bryce Bayer, 1973 [8] [*Bresenham's line algorithm*](https://en.wikipedia.org/wiki/Bresenham%27s_line_algorithm), Wikipedia40:T633,function split1DNoise(length) { const rootIndices = []; const { random } = Math; const flip = () => random() < 0.5; function splitHalf(indices) { // split in half let left = Math.floor(indices.length / 2); let right = indices.length - left; // randomize which side we round to if (flip()) [left, right] = [right, left]; let leftIndices = [], rightIndices = []; let lNumerator = 0, rNumerator = 0; const denominator = right < left ? left : right; let i = 0; for (let i = 0; i < indices.length; ) { lNumerator += left; rNumerator += right; if (flip()) { if (lNumerator >= denominator) { lNumerator -= denominator; leftIndices.push(indices[i++]); } if (rNumerator >= denominator) { rNumerator -= denominator; rightIndices.push(indices[i++]); } } else { if (rNumerator >= denominator) { rNumerator -= denominator; rightIndices.push(indices[i++]); } if (lNumerator >= denominator) { lNumerator -= denominator; leftIndices.push(indices[i++]); } } } // And of course we repeat this proces recursively on both indices if (leftIndices.length > 2) leftIndices = splitHalf(leftIndices); if (rightIndices.length > 2) rightIndices = splitHalf(rightIndices); // copy divided indices back from buffer return [...leftIndices, ...rightIndices]; } for (let i = 0; i < length; i++) rootIndices[i] = i; return splitHalf(rootIndices); }41:[{"id":"e0fdbd2d089b5287c64fe8f27578d3b4a1ff458c36fac47f2c82d7c3d0a9eea24712d253278e4ada0ee840db7a52beae2c3be8048152a974e8266628616cefc0","href":"https://static.observableusercontent.com/files/e0fdbd2d089b5287c64fe8f27578d3b4a1ff458c36fac47f2c82d7c3d0a9eea24712d253278e4ada0ee840db7a52beae2c3be8048152a974e8266628616cefc0","name":"image@2.png","type":"image/png","lastModified":1637098635837.93,"size":9332},{"id":"633f14c280d45f8f551f8ac6e8b9b2544ecf860178e4ab07ebdd547ba26b0f8f5ef7727be5c3b00dd10c17bbe32bbe855b136778aa9a2dfbfc987c246ad9be54","href":"https://static.observableusercontent.com/files/633f14c280d45f8f551f8ac6e8b9b2544ecf860178e4ab07ebdd547ba26b0f8f5ef7727be5c3b00dd10c17bbe32bbe855b136778aa9a2dfbfc987c246ad9be54","name":"aliasing.png","type":"image/png","lastModified":1637099833461.866,"size":2285},{"id":"182978e5e1610b41699b87c0bcd828a102b32e0c5ffada303a25c419f53d04fd7df151c6f102a07f6ac49aaa20676c86a2b1cd9d4024b562f54d9c4f8d46c65c","href":"https://static.observableusercontent.com/files/182978e5e1610b41699b87c0bcd828a102b32e0c5ffada303a25c419f53d04fd7df151c6f102a07f6ac49aaa20676c86a2b1cd9d4024b562f54d9c4f8d46c65c","name":"Bayer Matrix 1.gif","type":"image/gif","lastModified":1637108573347.809,"size":86408},{"id":"5b384a6278acfdec4bf5e2ffe01775eb15f0419f44b1a79d791ff477a69849edd3924d0e3125af269a31fde9898d02c03ed59b1da00a6dd7e512fa2d17721f2a","href":"https://static.observableusercontent.com/files/5b384a6278acfdec4bf5e2ffe01775eb15f0419f44b1a79d791ff477a69849edd3924d0e3125af269a31fde9898d02c03ed59b1da00a6dd7e512fa2d17721f2a","name":"Bayer Matrix 3.png","type":"image/png","lastModified":1637108660121.149,"size":44244},{"id":"bce57816c29e94a1c7dc414d6904b31e9dba0b2c0e8c5992c64bf646b06a338eabd5dd8de03e82c029b608e58fa67e9c202d5fa0d5fa054130f61499679c6b54","href":"https://static.observableusercontent.com/files/bce57816c29e94a1c7dc414d6904b31e9dba0b2c0e8c5992c64bf646b06a338eabd5dd8de03e82c029b608e58fa67e9c202d5fa0d5fa054130f61499679c6b54","name":"Bayer Matrix 2.gif","type":"image/gif","lastModified":1637108648799.997,"size":86670}] 2f:{"title":"Pseudo-Blue Noise","cells":[{"id":3827,"value":"html``","pinned":false,"hidden":false,"mode":"ojs"},{"id":3835,"value":"$31","pinned":false,"hidden":false,"mode":"ojs"},{"id":0,"value":"$32","pinned":false,"hidden":false,"mode":"ojs"},{"id":1610,"value":"label(\n \"Grayscale (for the nitpickers: a naive sRGB color ramp, so technically not a perceptually linear change)\",\n render(width, 150, (x, y) => toU32Gray((255 - (x * 256) / width) | 0))\n)","pinned":false,"hidden":false,"mode":"ojs"},{"id":1643,"value":" We basically wish to round our gray levels for every pixel to black and white in such a way that the sum of the resulting black and white pixels create an illusion of gray. A very simple method would be to flip a coin for each pixel to decide whether to turn it black or white, but bias the coin by the input color. The result looks like this:","pinned":false,"hidden":false,"mode":"md"},{"id":479,"value":"label(\n \"Random dithering (white noise)\",\n render(width, 150, (x, y) => toU32Gray(Math.random() * width < x ? 0 : 255))\n)","pinned":false,"hidden":false,"mode":"ojs"},{"id":1660,"value":"This sort of works, but as you can see the result does not look very \"even\". True randomness is lumpier than human intuition expects it to be! So we want something that is a bit smoother in its transition. One approach is to use ordered dithering[[5](https://en.wikipedia.org/wiki/Ordered_dithering)], which often uses Bayer ordering:","pinned":false,"hidden":false,"mode":"md"},{"id":1603,"value":"label(\n \"8 × 8 Bayer Dithering\",\n render(width, 150, (x, y) =>\n toU32Gray(get8x8threshold(x, y) < (256 * x) / width ? 0 : 255)\n )\n)","pinned":false,"hidden":false,"mode":"ojs"},{"id":1705,"value":"This produces a smoother gradient, but it introduces a new problem: the clearly visible patterns. This can be somewhat mitigated by using a different ordering, for example Interleaved Gradient Noise[[6](http://www.iryoku.com/next-generation-post-processing-in-call-of-duty-advanced-warfare)]:","pinned":false,"hidden":false,"mode":"md"},{"id":490,"value":"label(\n \"Interleaved Gradient Noise Dithering\",\n render(width, 150, (x, y) =>\n toU32Gray(interleavedGradientNoise(x, y) < x / width ? 0 : 255)\n )\n)","pinned":false,"hidden":false,"mode":"ojs"},{"id":1716,"value":"Interleaved Gradient Noise does not have the rectangular patterns of Bayer dithering. That said, it still has repeating patterns, and in both cases these the patterns can result in aliasing effects when we zoom in or out:","pinned":false,"hidden":false,"mode":"md"},{"id":1780,"value":"md`*The effect of 80% browser zoom on the above two forms of dithering*\n${await FileAttachment(\"aliasing.png\").image()}`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1742,"value":"What we ideally want is something that has no obvious repeating patterns, like the \"random\" dithering, but also is \"evenly spread out\", like ordered dithering. This apparent contradiction is the problem that blue noise solves in the context of image dithering:","pinned":false,"hidden":false,"mode":"md"},{"id":1746,"value":"md`*Blue noise dithering and grayscale ramp for comparison
([source: Christoph Peters' \"Free Blue Noise Textures\" blog post](https://momentsingraphics.de/BlueNoise.html))*\n${await FileAttachment(\"image@2.png\").image()}`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1834,"value":"md`Blue noise is noise with the low frequency components - the \"lumpiness\" in white noise - removed. Sadly, generating it is rather slow. Yes, technically speaking [Bridson's algorithm](https://observablehq.com/@mbostock/bridsons-algorithm) is ${tex`\\mathcal{O}(N)`}, but ${tex`N`} in this case is the number of samples one wishes to generate evenly in a space, and it assumes there is enough empty space between the samples (it also has relatively high overhead per sample). It does not work when trying to create a random value for each pixel in a grid.\n\nThis is where pseudo-blue noise comes in:`","pinned":false,"hidden":false,"mode":"ojs"},{"id":483,"value":"{\n const indices = noise(width, 150);\n return label(\n \"Pseudo-blue noise\",\n render(width, 150, (x, y) =>\n toU32Gray(indices[x + y * width] < x * 150 ? 0 : 255)\n )\n );\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1851,"value":"returnFromAppendixA = md`While not quite as smoothly distributed and visually pleasing as proper blue noise, it is a lot better than white noise. The algorithm to generate it assigns a value to each pixel, and is ${tex`\\mathcal{O}(Nlog(N))`}, where ${tex`N`} is the total number of pixels. It also has much lower constant overhead per item than Bridson's algorithm. *Update:* Fil managed to come up with a point-based shader version, meaning it runs in ${tex`\\mathcal{O}(N)`}, [check out his notebook](https://observablehq.com/@fil/pseudoblue)!\n\nIn order to explain how the algorithm works we have to first understand how the Bayer indices work (originally I also had an explanation of interleaved gradient noise, but it wasn't really necessary to understand pseudo-blue noise. I left it [as an appendix for the curious though](#AppendixA)).\n\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1895,"value":"$33","pinned":false,"hidden":false,"mode":"ojs"},{"id":1954,"value":"viewof bayerSamples = Inputs.range([0, 8 * 8], {\n value: 1,\n step: 1,\n label: \"Bayer Samples\"\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":1947,"value":"render(\n 8,\n 8,\n (x, y) => toU32Gray(get8x8threshold(x, y) < bayerSamples * 4 ? 0 : 255),\n Math.min(width, 400)\n)","pinned":false,"hidden":false,"mode":"ojs"},{"id":4555,"value":"Let's have look at the actual 4 × 4 and 8 × 8 matrices, and try to follow how the indices step through the grid. Remember, the numbers represent the order in which the pixels should flip: (and also remember that programmers start counting at zero)","pinned":false,"hidden":false,"mode":"md"},{"id":1873,"value":"\\begin{bmatrix}\n0 & 2 \\\\\n3 & 1 \\\\\n\\end{bmatrix}","pinned":false,"hidden":false,"mode":"tex"},{"id":1878,"value":"\\begin{bmatrix}\n 0 & 8 & 2 & 10 \\\\\n12 & 4 & 14 & 6 \\\\\n 3 & 11 & 1 & 9 \\\\\n15 & 7 & 13 & 5 \\\\\n\\end{bmatrix}","pinned":false,"hidden":false,"mode":"tex"},{"id":1883,"value":"\\begin{bmatrix}\n 0 & 32 & 8 & 40 & 2 & 34 & 10 & 42 \\\\\n48 & 16 & 56 & 24 & 50 & 18 & 58 & 26 \\\\\n12 & 44 & 4 & 36 & 14 & 46 & 6 & 38 \\\\\n60 & 28 & 52 & 20 & 62 & 30 & 54 & 22 \\\\\n 3 & 35 & 11 & 43 & 1 & 33 & 9 & 41 \\\\\n51 & 19 & 59 & 27 & 49 & 17 & 57 & 25 \\\\\n15 & 47 & 7 & 39 & 13 & 45 & 5 & 37 \\\\\n63 & 31 & 55 & 23 & 61 & 29 & 53 & 21 \\\\\n\\end{bmatrix}","pinned":false,"hidden":false,"mode":"tex"},{"id":1920,"value":"That might still be a bit tricky to follow, so I'll animate the jumps from index to index in sequences of four:","pinned":false,"hidden":false,"mode":"md"},{"id":1917,"value":"FileAttachment(\"Bayer Matrix 1.gif\").image()","pinned":false,"hidden":false,"mode":"ojs"},{"id":4480,"value":"But now look at how the numbers are divided if you split the 4 × 4 or 8 × 8 tables into four sub-grids:","pinned":false,"hidden":false,"mode":"md"},{"id":1925,"value":"FileAttachment(\"Bayer Matrix 2.gif\").image()","pinned":false,"hidden":false,"mode":"ojs"},{"id":1929,"value":"FileAttachment(\"Bayer Matrix 3.png\").image()","pinned":false,"hidden":false,"mode":"ojs"},{"id":1966,"value":"$34","pinned":false,"hidden":false,"mode":"md"},{"id":2585,"value":"viewof maxBayerRecursionSteps = Inputs.range([0, 4], {\n value: 0,\n step: 1,\n label: \"Recursive Steps\"\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":2579,"value":"$35","pinned":false,"hidden":false,"mode":"ojs"},{"id":4495,"value":"How do we make sure that we are able to divide these numbers evenly into sub-grids all the way down to 1 × 1 sections? Well obviously, if the grid is a square with sides that are a power of two we have nothing to worry about. This is why Bayer indices always are powers of two sized. However, there are tricks to make this \"spread indices evenly into sub-grids\" work for other grid sizes however, and we'll get back to that!\n\nHopefully you can also seen now why this way of spreading out the indices creates this perfectly regular recursive square pattern. Which as mentioned can lead to aliasing artifacts.\n\nWe haven't actually discussed any code for creating a Bayer index yet. That is because getting across the general concept of evenly spreading out indices across a surface was more important. There are various ways to implement an algorithm to create these indices. Here is a top-down recursive version:","pinned":false,"hidden":false,"mode":"md"},{"id":2545,"value":"$36","pinned":true,"hidden":false,"mode":"ojs"},{"id":2524,"value":"$37","pinned":false,"hidden":false,"mode":"md"},{"id":2708,"value":"$38","pinned":true,"hidden":false,"mode":"ojs"},{"id":2845,"value":"Now we can have an arbitrary number of indices and still spread it them out as evenly as possible:","pinned":false,"hidden":false,"mode":"md"},{"id":2770,"value":"viewof maxBayer1DRecursionSteps = Inputs.range([0, 8], {\n value: 2,\n step: 1,\n label: \"Recursive Steps\"\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":2788,"value":"viewof bayerLength1D = Inputs.range([2, Math.min(256, width)], {\n value: 64,\n step: 1,\n label: \"Total Indices\"\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":2832,"value":"*As grayscale*","pinned":false,"hidden":false,"mode":"md"},{"id":2772,"value":"$39","pinned":false,"hidden":false,"mode":"ojs"},{"id":2825,"value":"*As a bar chart showing the value of individual indices*","pinned":false,"hidden":false,"mode":"md"},{"id":2829,"value":"$3a","pinned":false,"hidden":false,"mode":"ojs"},{"id":2849,"value":"As we can see this still looks highly regular, but now we have the freedom to introduce some noise into our system. Before we do that that you might be curious about how we are going turn this back into a 2D image though.\n\nWith Bayer ordering we split the image recursively into quadrants. Here we instead split the image into halves, and always do so along the longest edge. If we repeat this often enough then we are guaranteed to end up with a 1 × 1 subgrid. At that point we have the determined the (x,y) position of a single index. \n\n*(we could also try an alternative method using space-filing curves, and I might do so in a follow-up notebook to see what that looks like. However, then we would end up with the same issue as before: most space-filling curves require sides that are a power of two)*","pinned":false,"hidden":false,"mode":"md"},{"id":2985,"value":"$3b","pinned":true,"hidden":false,"mode":"ojs"},{"id":2989,"value":"{\n const indices = split2D(width, 100);\n return label(\n \"Grayscale\",\n render(width, 100, (x, y) =>\n toU32Gray((indices[x + y * width] * 256) / (width * 100))\n )\n );\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":2992,"value":"{\n const indices = split2D(width, 100);\n return label(\n \"Black and white gradient\",\n render(width, 100, (x, y) =>\n toU32Gray(indices[x + y * width] > x * 100 ? 255 : 0)\n )\n );\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":3008,"value":"We do not quite get a Bayer pattern for two reasons. The first is that our initial dimensions are not a power of two. When our subgrids get close to the final positions the rounding errors become quite severe. For example, dividing a 3 × 3 rectangle will result in one corner getting one index, and the opposite corner getting four. This results in a bias that effectively causes a sort of \"moiré\" pattern. When we change to a grid that has sides that are a power of two we should get a more even result:","pinned":false,"hidden":false,"mode":"md"},{"id":3156,"value":"{\n const { ctx, imgData } = getCtx(2 ** (Math.log2(width) | 0), 128);\n const indices = split2D(imgData.width, imgData.height);\n for (let y = 0; y < imgData.height; y++) {\n for (let x = 0; x < imgData.width; x++) {\n const idx = (x + y * imgData.width) * 4;\n const v = (indices[idx >> 2] * 256) / (imgData.width * imgData.height);\n imgData.data[idx] = v;\n imgData.data[idx + 1] = v;\n imgData.data[idx + 2] = v;\n imgData.data[idx + 3] = 255; // make sure image is not transparent\n }\n }\n ctx.putImageData(imgData, 0, 0);\n\n return html`

Grayscale, power of two sized
${ctx.canvas}

`;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":3161,"value":"{\n const { ctx, imgData } = getCtx(2 ** (Math.log2(width) | 0), 128);\n const indices = split2D(imgData.width, imgData.height);\n for (let y = 0; y < imgData.height; y++) {\n for (let x = 0; x < imgData.width; x++) {\n const idx = (x + y * imgData.width) * 4;\n const v = indices[idx >> 2] > x * imgData.height ? 255 : 0;\n imgData.data[idx] = v;\n imgData.data[idx + 1] = v;\n imgData.data[idx + 2] = v;\n imgData.data[idx + 3] = 255; // make sure image is not transparent\n }\n }\n ctx.putImageData(imgData, 0, 0);\n\n return html`

Black and white gradient, power of two-sized
${ctx.canvas}

`;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":3170,"value":"The other reason is that we don't divide our rectangle into four parts at once and cross from top-left to bottom-right in a recursive pattern. Still, I would say it is a decent start for creating an ordered spread of indices. Now we can turn this starting point into a proper noise function!\n\nFirst we have to address something that some of you might have already noticed: when we split an array of indices into two, alternating between left and right, how do you think all the odd or even numbers are divided?\n\n [0, 1, 2, 3, 4, 5, 6, 7] → [[0, 2, 4, 6], [1, 3, 5, 7]]\n\nLet's visualize that: we take the our `split2D` function and turn it into a \"grayscale\" image like above, but we color all odd numbers blue, and all even numbers red:","pinned":false,"hidden":false,"mode":"md"},{"id":4067,"value":"async function colorBit(bit, noise, visibility) {\n await visibility();\n const h = 100;\n const w = 2 ** (Math.log2(width) | 0);\n const indices = noise(w, h);\n return render(w, h, (x, y) => {\n const v0 = indices[x + y * w];\n const v = ((v0 * 256) / (w * h)) | 0;\n return v0 & (1 << bit) ? toU32Color(v, v, 255) : toU32Color(255, v, v);\n });\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1183,"value":"colorBit(0, split2D, visibility)","pinned":false,"hidden":false,"mode":"ojs"},{"id":3179,"value":"Well now... if we can predict where we can find all odd numbers that easily, then our noise function can't be considered all that *random*, can it? And this is true for all powers of two! If we color multiples (ignoring remainder) of four, eight, sixteen, etc. then we get this:","pinned":false,"hidden":false,"mode":"md"},{"id":3182,"value":"colorBit(1, split2D, visibility)","pinned":false,"hidden":false,"mode":"ojs"},{"id":3188,"value":"colorBit(2, split2D, visibility)","pinned":false,"hidden":false,"mode":"ojs"},{"id":3191,"value":"colorBit(3, split2D, visibility)","pinned":false,"hidden":false,"mode":"ojs"},{"id":3198,"value":"Luckily, the solution for this is quite simple: currently, whenever there is a tie (so when we assign an index to both the left and right half in the same loop), the left half goes first. If we randomize this priority it no longer is guaranteed that all odd or even numbers go to one half. Instead both halves will get an equal number of odd and even numbers on average. The interesting part is that we will still create two monotonically increasing gradients (which is what we want to ensure this more or less even-spread property).","pinned":false,"hidden":false,"mode":"md"},{"id":3016,"value":"$3c","pinned":false,"hidden":false,"mode":"ojs"},{"id":4156,"value":"showDiff(\n split2D,\n split2D_v2,\n \"Diff between split2D and split2D_v2\",\n true\n)","pinned":false,"hidden":false,"mode":"ojs"},{"id":3019,"value":"{\n const { ctx, imgData } = getCtx(width, 100);\n const indices = split2D_v2(imgData.width, imgData.height);\n for (let y = 0; y < imgData.height; y++) {\n for (let x = 0; x < imgData.width; x++) {\n const idx = (x + y * imgData.width) * 4;\n const v = (indices[idx >> 2] * 256) / (imgData.width * imgData.height);\n imgData.data[idx] = v;\n imgData.data[idx + 1] = v;\n imgData.data[idx + 2] = v;\n imgData.data[idx + 3] = 255; // make sure image is not transparent\n }\n }\n ctx.putImageData(imgData, 0, 0);\n\n return ctx.canvas;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":3024,"value":"{\n const { ctx, imgData } = getCtx(width, 100);\n const indices = split2D_v2(imgData.width, imgData.height);\n for (let y = 0; y < imgData.height; y++) {\n for (let x = 0; x < imgData.width; x++) {\n const idx = (x + y * imgData.width) * 4;\n const v = indices[idx >> 2] > x * imgData.height ? 255 : 0;\n imgData.data[idx] = v;\n imgData.data[idx + 1] = v;\n imgData.data[idx + 2] = v;\n imgData.data[idx + 3] = 255; // make sure image is not transparent\n }\n }\n ctx.putImageData(imgData, 0, 0);\n\n return ctx.canvas;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":3203,"value":"colorBit(0, split2D_v2, visibility)","pinned":false,"hidden":false,"mode":"ojs"},{"id":3215,"value":"As a bonus it looks like we have already introduced some noise to our function! Nice! Depending on the resolution of your screen however there might still be visible gridlines. As mentioned, those can appear because we currently always round into the same direction when dividing the available space in half. The resulting bias in the placing of our indices can cause visible gridlines. All we have to do to remove that bias is randomize the rounding direction as well:","pinned":false,"hidden":false,"mode":"md"},{"id":3035,"value":"$3d","pinned":false,"hidden":false,"mode":"ojs"},{"id":4178,"value":"showDiff(\n split2D_v2,\n split2D_v3,\n \"Diff between split2D_v2 and split2D_v3\",\n true\n)","pinned":false,"hidden":false,"mode":"ojs"},{"id":3251,"value":"{\n const { ctx, imgData } = getCtx(width, 100);\n const indices = split2D_v3(imgData.width, imgData.height);\n for (let y = 0; y < imgData.height; y++) {\n for (let x = 0; x < imgData.width; x++) {\n const idx = (x + y * imgData.width) * 4;\n const v = (indices[idx >> 2] * 256) / (imgData.width * imgData.height);\n imgData.data[idx] = v;\n imgData.data[idx + 1] = v;\n imgData.data[idx + 2] = v;\n imgData.data[idx + 3] = 255; // make sure image is not transparent\n }\n }\n ctx.putImageData(imgData, 0, 0);\n\n return ctx.canvas;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":3254,"value":"{\n const { ctx, imgData } = getCtx(width, 100);\n const indices = split2D_v3(imgData.width, imgData.height);\n for (let y = 0; y < imgData.height; y++) {\n for (let x = 0; x < imgData.width; x++) {\n const idx = (x + y * imgData.width) * 4;\n const v = indices[idx >> 2] > x * imgData.height ? 255 : 0;\n imgData.data[idx] = v;\n imgData.data[idx + 1] = v;\n imgData.data[idx + 2] = v;\n imgData.data[idx + 3] = 255; // make sure image is not transparent\n }\n }\n ctx.putImageData(imgData, 0, 0);\n\n return ctx.canvas;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":3260,"value":"{\n const { ctx, imgData } = getCtx(width, 100);\n const indices = split2D_v3(imgData.width, imgData.height);\n for (let y = 0; y < imgData.height; y++) {\n for (let x = 0; x < imgData.width; x++) {\n const idx = (x + y * imgData.width) * 4;\n const v0 = indices[idx >> 2];\n const v = ((v0 * 256) / (imgData.width * imgData.height)) | 0;\n if (v0 & 1) {\n imgData.data[idx] = 255;\n imgData.data[idx + 1] = v;\n imgData.data[idx + 2] = v;\n } else {\n imgData.data[idx] = v;\n imgData.data[idx + 1] = v;\n imgData.data[idx + 2] = 255;\n }\n imgData.data[idx + 3] = 255; // make sure image is not transparent\n }\n }\n ctx.putImageData(imgData, 0, 0);\n\n return ctx.canvas;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":3461,"value":"And there we have our final algorithm! Well, I'm sure more refinements are possible (for example, we could try letting go of axis-aligned rectangles and subdivide the spaces into randomly shaped convex polygons, which might fully get rid of all remaining \"grid artifacts\") but I think this is a good point to stop for now.\n\nTo close off, let's compare what a fade-to-black looks like for random noise, pseudo-blue noise, interleaved gradient dithering and Bayer dithering:","pinned":false,"hidden":false,"mode":"md"},{"id":507,"value":"viewof threshold = Inputs.range([0, 1], {\n value: 0.5,\n step: 1 / 64,\n label: \"Gray value threshold\"\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":543,"value":"render(width, 100, (x, y) =>\n toU32Gray(randomNoiseMap[x + y * width] >= threshold * width * 100 ? 255 : 0)\n)","pinned":false,"hidden":false,"mode":"ojs"},{"id":4738,"value":"render(width, 100, (x, y) =>\n toU32Gray(noiseMaps[0][x + y * width] >= threshold * width * 100 ? 255 : 0)\n)","pinned":false,"hidden":false,"mode":"ojs"},{"id":4736,"value":"render(width, 100, (x, y) =>\n toU32Gray(noiseMaps[1][x + y * width] >= threshold * width * 100 ? 255 : 0)\n)","pinned":false,"hidden":false,"mode":"ojs"},{"id":4732,"value":"render(width, 100, (x, y) =>\n toU32Gray(noiseMaps[2][x + y * width] >= threshold * width * 100 ? 255 : 0)\n)","pinned":false,"hidden":false,"mode":"ojs"},{"id":510,"value":"render(width, 100, (x, y) =>\n toU32Gray(noiseMaps[3][x + y * width] >= threshold * width * 100 ? 255 : 0)\n)","pinned":false,"hidden":false,"mode":"ojs"},{"id":525,"value":"render(width, 100, (x, y) =>\n toU32Gray(interleavedGradientNoise(x, y) >= threshold ? 255 : 0)\n)","pinned":false,"hidden":false,"mode":"ojs"},{"id":3457,"value":"render(width, 100, (x, y) =>\n toU32Gray(get8x8threshold(x, y) / 255 >= threshold ? 255 : 0)\n)","pinned":false,"hidden":false,"mode":"ojs"},{"id":4640,"value":"$3e","pinned":true,"hidden":false,"mode":"ojs"},{"id":1632,"value":"$3f","pinned":false,"hidden":false,"mode":"md"},{"id":1975,"value":"AppendixA = md`---\n\n## [⮨](#returnFromAppendixA) Bonus: Interleaved Gradient Noise`","pinned":false,"hidden":false,"mode":"ojs"},{"id":107,"value":"render(\n (width / 3) | 0,\n 20,\n (x, y) =>\n toU32Gray(interleavedGradientNoise(x, y) < x / ((width / 3) | 0) ? 0 : 255),\n Math.min(width, 800)\n)","pinned":false,"hidden":false,"mode":"ojs"},{"id":2062,"value":"Interleaved Gradient Noise (ILG) is a fairly recent invention. It was created by Jorge Jimenez as a really cheap noise function to use in shaders. The main use-case is reducing visible banding in smooth gradients (this is arguably the main use-case for image dithering in applications these days). The code is *impressively* simple:","pinned":false,"hidden":false,"mode":"md"},{"id":2066,"value":"function interleavedGradientNoise(x, y) {\n // http://www.iryoku.com/next-generation-post-processing-in-call-of-duty-advanced-warfare\n // https://bartwronski.com/2016/10/30/dithering-part-three-real-world-2d-quantization-dithering/\n const v = 52.9829189 * (0.06711056 * x + 0.00583715 * y);\n return v - Math.floor(v);\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":2078,"value":"It is easy to verify that for any given `(x,y)` input this function will return a value from 0 to 1, exclusive. That is what the `v - Math.floor(v)` does.\n\n\nIf we map the `(x,y)` coordinates of each pixel of an image to this function, and multiply the output by 256, we get something like this:","pinned":false,"hidden":false,"mode":"md"},{"id":1106,"value":"render(width >> 1, 100, (x, y) =>\n toU32Gray((interleavedGradientNoise(x, y) * 256) | 0)\n)","pinned":false,"hidden":false,"mode":"ojs"},{"id":2121,"value":"Not as regular as the Bayer pattern, but still evenly spread out. How does it work? Let's break down the function by simplifying it. First, what would happen if we only look at the x-input or y-input, and ignore the other half?","pinned":false,"hidden":false,"mode":"md"},{"id":2097,"value":"function interleavedGradientNoiseX(x, y) {\n const v = 52.9829189 * 0.06711056 * x;\n return v - Math.floor(v);\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":2091,"value":"render(\n width >> 1,\n 100,\n (x, y) => toU32Gray((interleavedGradientNoiseX(x, y) * 256) | 0),\n width\n)","pinned":false,"hidden":false,"mode":"ojs"},{"id":2105,"value":"function interleavedGradientNoiseY(x, y) {\n const v = 52.9829189 * 0.00583715 * y;\n return v - Math.floor(v);\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":2094,"value":"render(\n width >> 1,\n 100,\n (x, y) => toU32Gray((interleavedGradientNoiseY(x, y) * 256) | 0),\n width\n)","pinned":false,"hidden":false,"mode":"ojs"},{"id":2151,"value":"What we can see in both cases is that we have two overlapping frequencies: a relatively high one, and a relatively low one. In the normal case we add these two together to get a fast-changing value, set up in such a way that the remainder of that value produces a pseudo-noise pattern.\n\nWhat happens if we remove the `52.9829189` factor?","pinned":false,"hidden":false,"mode":"md"},{"id":2179,"value":"function interleavedGradientNoiseSmallX(x, y) {\n const v = 0.06711056 * x;\n return v - Math.floor(v);\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":2167,"value":"render(\n width >> 1,\n 100,\n (x, y) => toU32Gray((interleavedGradientNoiseSmallX(x, y) * 256) | 0),\n width\n)","pinned":false,"hidden":false,"mode":"ojs"},{"id":2183,"value":"function interleavedGradientNoiseSmallY(x, y) {\n const v = 0.00583715 * y;\n return v - Math.floor(v);\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":2170,"value":"render(\n width >> 1,\n 100,\n (x, y) => toU32Gray((interleavedGradientNoiseSmallY(x, y) * 256) | 0),\n width\n)","pinned":false,"hidden":false,"mode":"ojs"},{"id":2177,"value":"function interleavedGradientNoiseSmall(x, y) {\n const v = 0.06711056 * x + 0.00583715 * y;\n return v - Math.floor(v);\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":2156,"value":"render(\n width >> 1,\n 100,\n (x, y) => toU32Gray((interleavedGradientNoiseSmall(x, y) * 256) | 0),\n width\n)","pinned":false,"hidden":false,"mode":"ojs"},{"id":2208,"value":"So here we see the base of our noise function: two gradients. One reasonably fast (x), one very slow (y). Adding them together produces a slightly slanted gradient. By multiplying either with a bigger value we turn them into a sort of [beat](https://en.wikipedia.org/wiki/Beat_(acoustics)) of overlapping gradients. The sum of all of it is the final resulting noise function.\n\nWe can also see this if we play with the constants ourselves:","pinned":false,"hidden":false,"mode":"md"},{"id":2237,"value":"function interleavedGradientNoiseCustom(x, y) {\n const v = constantC * (constantX * x + constantY * y);\n return v - Math.floor(v);\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":2239,"value":"viewof constantC = Inputs.range([0, 100], {\n value: 52.9829189,\n label: \"default: 52.9829189\"\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":2249,"value":"viewof constantX = Inputs.range([0, 1], {\n value: 0.06711056,\n label: \"default: 0.06711056\"\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":2250,"value":"viewof constantY = Inputs.range([0, 1], {\n value: 0.00583715,\n label: \"default: 0.00583715\"\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":2257,"value":"render(\n width >> 1,\n 100,\n (x, y) => toU32Gray(interleavedGradientNoiseCustom(x, y) * 256) | 0,\n width\n)","pinned":false,"hidden":false,"mode":"ojs"},{"id":2421,"value":"render(width, 150, (x, y) =>\n toU32Gray(interleavedGradientNoiseCustom(x, y) < x / width ? 0 : 255)\n)","pinned":false,"hidden":false,"mode":"ojs"},{"id":2502,"value":"[⮨](#returnFromAppendixA)\n\n---","pinned":false,"hidden":false,"mode":"md"},{"id":2982,"value":"noise = split2D_v4","pinned":true,"hidden":false,"mode":"ojs"},{"id":535,"value":"noiseMaps = [split2D, split2D_v2, split2D_v3, split2D_v4].map(noise => noise(width, 100));","pinned":false,"hidden":false,"mode":"ojs"},{"id":4266,"value":"noiseMap1D = split1DNoise(width * 100)","pinned":true,"hidden":false,"mode":"ojs"},{"id":4270,"value":"$40","pinned":true,"hidden":false,"mode":"ojs"},{"id":546,"value":"randomNoiseMap = {\n const indices = new Uint32Array(width * 100);\n // initiate + shuffle array\n for (let i = 0; i < indices.length; i++) {\n const j = (Math.random() * (i + 1)) | 0;\n indices[i] = indices[j];\n indices[j] = i;\n }\n return indices;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1600,"value":"get8x8threshold = {\n //prettier-ignore\n const thresholds8x8 = [\n 0, 48, 12, 60, 3, 51, 15, 63,\n 32, 16, 44, 28, 35, 19, 47, 31, \n 8, 56, 4, 52, 11, 59, 7, 55,\n 40, 24, 36, 20, 43, 27, 39, 23,\n 2, 50, 14, 62, 1, 49, 13, 61,\n 34, 18, 46, 30, 33, 17, 45, 29,\n 10, 58, 6, 54, 9, 57, 5, 53,\n 42, 26, 38, 22, 43, 25, 37, 21\n ].map((v) => 2 + 4 * v);\n function get8x8threshold(x, y) {\n return thresholds8x8[(x & 7) + ((y & 7) << 3)];\n }\n get8x8threshold.label = \"8×8\";\n get8x8threshold.width = 8;\n get8x8threshold.height = 8;\n return get8x8threshold;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":4461,"value":"import {render, label, getCtx, toU32Color, toU32Gray} from \"@jobleonard/canvas-tools\"","pinned":true,"hidden":false,"mode":"ojs"},{"id":4417,"value":"import {diff} from \"@jobleonard/diff-tool\"","pinned":true,"hidden":false,"mode":"ojs"},{"id":4163,"value":"showDiff = (left, right, summary, open, opts) =>\n summary\n ? html`

${summary}

${diff(\n left,\n right\n )}

`\n : diff(left, right, opts)","pinned":true,"hidden":false,"mode":"ojs"}],"files":"$W41","stdlib":"1","titling":"auto"}