2b:null 2c:["$","title",null,{"children":"Poisson disk sampling functions | Fil | Observable"}] 31:T828,A brief history of blue noise 2 dimensions point distributions (a.k.a. Poisson disk sampling): - dark ages: Monte Carlo, aka “dart-throwing”; very slow - 1991: [Don Mitchell](http://mentallandscape.com/) publishes “[Spectrally optimal sampling for distribution ray tracing](http://mentallandscape.com/Papers_siggraph91.pdf),” introducing a sequential sampling algorithm which has a time complexity of ${tex`\mathcal{O}(N^2)`} - 2007: [Robert Bridson](https://www.cs.ubc.ca/~rbridson/) publishes “[Fast Poisson disk sampling in arbitrary dimensions](https://dl.acm.org/doi/10.1145/1278780.1278807)”, which improves on Mitchell’s algorithm—now linear time. - 2012: Mike Bostock implements [Mitchell’s](https://observablehq.com/@mbostock/best-candidate-circles) algorithm (it was the time of d3v2 and bl.ocks!). - 2014: Jason Davies implements [Bridson’s](https://www.jasondavies.com/poisson-disc/) as a replacement of Mitchell’s. (**A**) See Mike Bostock’s 2014 [animation](https://bl.ocks.org/mbostock/dbb02448b0f93e4c82c3). - 2019: Martin Roberts [tweaks Bridson’s for speed](http://extremelearning.com.au/an-improved-version-of-bridsons-algorithm-n-for-poisson-disc-sampling/). (**B**) - 2021: Jacob Rus [tweaks Roberts’ for quality](https://observablehq.com/@jrus/bridson-fork/2); adds another speed improvement in 2022. (**C**) This type of distribution very useful to generate random data which look organic. Applying a [Delaunay](https://observablehq.com/collection/@d3/d3-delaunay) triangulation and some graph magic ([Gabriel](https://observablehq.com/@rpgove/constructing-a-gabriel-graph) or [Urquhart](https://observablehq.com/@mbostock/urquhart-graph)) results in nice graphs; the associated Voronoi tessellation is harmonious. (Note that Jacob has also published a [spherical version](https://observablehq.com/@jrus/spheredisksample).) This notebook compares the last three implementations for speed and quality. Let’s call them A, B and C. We've selected the radius in each of them to get a good approximation of the number of points requested.32:T9f3,// https://observablehq.com/@kemper/voronoi-polygon-smoothing // https://observablehq.com/@mbostock/poisson-disc-distribution // https://www.jasondavies.com/poisson-disc/ function* A(width, height, N) { const radius = Math.sqrt((0.63 * (width * height)) / N); const points = [], max = (3 * (width * height)) / (radius * radius), sampleSites = sampler(width, height, radius); for (var i = 0; i < max; ++i) { let s = sampleSites(); if (s) yield s; } function sampler(width, height, radius) { var k = 40, // maximum number of samples before rejection radius2 = radius * radius, R = 3 * radius2, cellSize = radius * Math.SQRT1_2, gridWidth = Math.ceil(width / cellSize), gridHeight = Math.ceil(height / cellSize), grid = new Array(gridWidth * gridHeight), queue = [], queueSize = 0, sampleSize = 0; return function () { if (!sampleSize) return sample(rng() * width, rng() * height); // Pick a random existing sample and remove it from the queue. while (queueSize) { var i = (rng() * queueSize) | 0, s = queue[i]; // Make a new candidate between [radius, 2 * radius] from the existing sample. for (var j = 0; j < k; ++j) { var a = 2 * Math.PI * rng(), r = Math.sqrt(rng() * R + radius2), x = s[0] + r * Math.cos(a), y = s[1] + r * Math.sin(a); // Reject candidates that are outside the allowed extent, // or closer than 2 * radius to any existing sample. if (0 <= x && x < width && 0 <= y && y < height && far(x, y)) return sample(x, y); } queue[i] = queue[--queueSize]; queue.length = queueSize; } }; function far(x, y) { var i = (x / cellSize) | 0, j = (y / cellSize) | 0, i0 = Math.max(i - 2, 0), j0 = Math.max(j - 2, 0), i1 = Math.min(i + 3, gridWidth), j1 = Math.min(j + 3, gridHeight); for (j = j0; j < j1; ++j) { var o = j * gridWidth; for (i = i0; i < i1; ++i) { if ((s = grid[o + i])) { var s, dx = s[0] - x, dy = s[1] - y; if (dx * dx + dy * dy < radius2) return false; } } } return true; } function sample(x, y) { var s = [x, y]; queue.push(s); grid[gridWidth * ((y / cellSize) | 0) + ((x / cellSize) | 0)] = s; ++sampleSize; ++queueSize; return s; } } }33:T84a,function* B(width, height, N) { const k = 4; // maximum number of samples before rejection const radius2 = (0.755 * (width * height)) / N; const radius = Math.sqrt(radius2); const cellSize = radius * Math.SQRT1_2; const gridWidth = Math.ceil(width / cellSize); const gridHeight = Math.ceil(height / cellSize); const grid = new Array(gridWidth * gridHeight); const queue = []; // Pick the first sample near the center yield sample( 0.1 * width * (4.5 + rng()), 0.1 * height * (4.5 + rng()), null); // Pick a random existing sample from the queue. pick: while (queue.length) { const i = (rng() * queue.length) | 0; const parent = queue[i]; const seed = rng(); const epsilon = 0.0000001; // Make a new candidate. for (let j = 0; j < k; ++j) { const a = 2 * Math.PI * (seed + (1.0 * j) / k); const r = radius + epsilon; const x = parent[0] + r * Math.cos(a); const y = parent[1] + r * Math.sin(a); // Accept candidates that are inside the allowed extent // and farther than 2 * radius to all existing samples. if (0 <= x && x < width && 0 <= y && y < height && far(x, y)) { yield sample(x, y); continue pick; } } // If none of k candidates were accepted, remove it from the queue. const r = queue.pop(); if (i < queue.length) queue[i] = r; } function far(x, y) { const i = (x / cellSize) | 0; const j = (y / cellSize) | 0; const i0 = Math.max(i - 2, 0); const j0 = Math.max(j - 2, 0); const i1 = Math.min(i + 3, gridWidth); const j1 = Math.min(j + 3, gridHeight); for (let j = j0; j < j1; ++j) { const o = j * gridWidth; for (let i = i0; i < i1; ++i) { const s = grid[o + i]; if (s) { const dx = s[0] - x; const dy = s[1] - y; if (dx * dx + dy * dy < radius2) return false; } } } return true; } function sample(x, y) { const s = (grid[gridWidth * ((y / cellSize) | 0) + ((x / cellSize) | 0)] = [ x, y ]); queue.push(s); return s; } }34:Tad7,function* C(width, height, N) { const k = 11; // maximum number of samples before rejection (use 13 or more for higher quality) const m = 4; // a number mutually prime to k const radius2 = (0.7 * (width * height)) / N; const radius = Math.sqrt(radius2); const cellSize = 1 / (radius * Math.SQRT1_2); const gridWidth = Math.ceil(width * cellSize) + 4; // pad grid by 2 on each side const gridHeight = Math.ceil(height * cellSize) + 4; const grid = new Float64Array(2 * gridWidth * gridHeight).fill(Infinity); const queue = []; const rotx = Math.cos((2 * Math.PI * m) / k); const roty = Math.sin((2 * Math.PI * m) / k); // Pick the first sample near the center yield sample(0.1 * width * (4.5 + rng()), 0.1 * height * (4.5 + rng()), null); // Pick a random existing sample from the queue. pick: while (queue.length) { const i = (rng() * queue.length) | 0; const parent = queue[i]; const t = tanpi_2(2 * rng() - 1); const q = 1 / (1 + t * t); const epsilon = 1e-6; let dw, dx, dy; dx = q ? (1 - t * t) * q : -1; // [dx, dy] = random unit vector dy = q ? 2 * t * q : 0; // Make a new candidate. for (let j = 0; j < k; ++j) { dw = dx * rotx - dy * roty; // temporary name for dx dy = dx * roty + dy * rotx; dx = dw; const rand0 = rng(); const rand1 = rng(); const r = (radius + epsilon) * (1 + 0.65 * rand0 * rand1); const x = parent[0] + r * dx; const y = parent[1] + r * dy; // Accept candidates that are inside the allowed extent // and farther than 2 * radius to all existing samples. if ((0 <= x) & (x < width) & (0 <= y) & (y < height) & far(x, y)) { yield sample(x, y); continue pick; } } // If none of k candidates were accepted, remove the point from the queue. const r = queue.pop(); if (i < queue.length) queue[i] = r; } // fast approximation of tan(πx/2) function tanpi_2(a) { let b = 1 - a * a; return a * (-0.0187108 * b + 0.31583526 + 1.27365776 / b); } function far(x, y) { const j0 = (y * cellSize) | 0; const i0 = (x * cellSize) | 0; for (let j = j0; j < j0 + 5; ++j) { const index0 = 2 * (j * gridWidth + i0); for (let i = index0; i < index0 + 10; i += 2) { const dx = grid[i] - x; const dy = grid[i + 1] - y; if (dx * dx + dy * dy < radius2) return false; } } return true; } function sample(x, y) { const s = [x, y]; // offset grid cell by 2 in each direction to account for border const i = (x * cellSize + 2) | 0; const j = (y * cellSize + 2) | 0; const index = 2 * (gridWidth * j + i); grid[index] = x; grid[index + 1] = y; queue.push(s); return s; } }35:T656,function showVoronoi(M, width, height, type = "sides") { const delaunay = d3.Delaunay.from(M(width, height, 2000)); const voronoi = delaunay.voronoi([0, 0, width, height]); const context = DOM.context2d(width, height); const areaColor = d3.scaleSequentialLog(d3.interpolateRdBu).domain([0.5, 2]); const shapeColor = d3.scaleSequential(d3.interpolateRdBu).domain([2.5, 9.5]); const {points, halfedges, triangles} = delaunay; for (let i = 0, n = halfedges.length; i < n; ++i) { const j = halfedges[i]; if (j < i) continue; const ti = triangles[i]; const tj = triangles[j]; context.moveTo(points[ti * 2], points[ti * 2 + 1]); context.lineTo(points[tj * 2], points[tj * 2 + 1]); } const N = delaunay.points.length / 2; const I = d3 .range(N) .filter( (i) => !( delaunay.points[2 * i] < 10 || delaunay.points[2 * i] > width - 10 || delaunay.points[2 * i + 1] < 10 || delaunay.points[2 * i + 1] > height - 10 ) ); const areas = Array.from(d3.range(N), (i) => Math.abs(d3.polygonArea(voronoi.cellPolygon(i))) ); const meanArea = d3.mean(I, (i) => areas[i]); for (const i of I) { context.beginPath(); voronoi.renderCell(i, context); if (type === "area") { context.fillStyle = areaColor(areas[i] / meanArea); } if (type === "shape") { const sides = voronoi.cellPolygon(i).length - 1; context.fillStyle = shapeColor(sides); } context.fill(); } voronoi.render(context); context.lineWidth = 0.6; context.stroke(); return context.canvas; }36:T5d2,function showDelaunay(M, width, height, type = "edges") { const delaunay = d3.Delaunay.from(M(width, height, 2000)); const context = DOM.context2d(width, height); const {points, halfedges, triangles} = delaunay; const startPt = (i) => { let ti = triangles[i]; return points.slice(ti * 2, ti * 2 + 2); } const endPt = (i) => { let tj = triangles[halfedges[i]]; return points.slice(tj * 2, tj * 2 + 2); } const edgeLength = (i) => { const [xi, yi] = startPt(i); const [xj, yj] = endPt(i); return Math.hypot(xj - xi, yj - yi); } const N = halfedges.length; const I = d3 .range(N) .filter(i => halfedges[i] > i) .filter(i => { const [xi, yi] = startPt(i); const [xj, yj] = endPt(i); return !( ((xi < 10) && (xj < 10)) || ((xi > width - 10) && (xj > width - 10)) || ((yi < 10) && (yj < 10)) || ((yi > height - 10) && (yj > height - 10))); }); const edgeLengths = Array.from(I, i => edgeLength(i)); const medianLength = Math.sqrt(d3.mean(edgeLengths.map(x => x*x))); const edgeColor = d3.scaleSequentialLog(d3.interpolateRdBu).domain([1/2, 2]); context.lineWidth = 1.5; for (const i of I) { const [xi, yi] = startPt(i); const [xj, yj] = endPt(i); context.beginPath(); context.moveTo(xi, yi); context.lineTo(xj, yj); context.strokeStyle = edgeColor(edgeLength(i) / medianLength); context.stroke() } return context.canvas; }37:[] 2f:{"id":"783b84bd84437222","slug":"poisson-distribution-generators","title":"Poisson disk sampling functions","thumbnail":"dc456e7b781a5bdbbb2f14dac6be9f0c69a560a5061f048b06a87c0e5bb79bae","access_level":"public","create_time":"$D2021-11-02T23:03:47.489Z","update_time":"$D2025-04-20T07:34:41.557Z","view_time":null,"first_public_time":"$D2022-05-29T20:32:54.938Z","fork":false,"stars":43,"starred":false,"owner":{"type":"team","id":"da2a32154638c457","login":"fil","name":"Fil","avatar_url":"https://avatars-next.observableusercontent.com/d5b1777907117a2bdf2a07eb98600d7a9beb21639fde42b9753e95c88c7731a6","tier":"starter_2024"},"version":1318,"publish_version":null,"role":null,"authored":false,"authors":[{"type":"individual","id":"45a379fcfcb14253","login":"fil","name":"Fil","avatar_url":"https://avatars-next.observableusercontent.com/d5b1777907117a2bdf2a07eb98600d7a9beb21639fde42b9753e95c88c7731a6","tier":null},{"type":"individual","id":"ddd7235babc79ce0","login":"jrus","name":"Jacob Rus","avatar_url":"https://avatars-next.observableusercontent.com/c33dbeed76527c50e64208b0e005e0b6e54a3bef2a14ea6a823f1c7de4867662","tier":null}],"members":[{"type":"individual","id":"45a379fcfcb14253","login":"fil","name":"Fil","avatar_url":"https://avatars-next.observableusercontent.com/d5b1777907117a2bdf2a07eb98600d7a9beb21639fde42b9753e95c88c7731a6","tier":null,"role":"editor"},{"type":"individual","id":"ddd7235babc79ce0","login":"jrus","name":"Jacob Rus","avatar_url":"https://avatars-next.observableusercontent.com/c33dbeed76527c50e64208b0e005e0b6e54a3bef2a14ea6a823f1c7de4867662","tier":null,"role":"editor"}],"comments":[],"custom_thumbnail":true,"trashed":false,"subscribed":null,"forks":7,"collection":{"id":"7f51a09dd1333537","slug":"algorithms","title":"Algorithms","owner":{"id":"da2a32154638c457","login":"fil","name":"Fil"}},"type":"notebook","body":{"title":"Poisson disk sampling functions","cells":[{"id":0,"value":"# Poisson disk sampling functions","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined"},{"id":167,"value":"$31","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined"},{"id":234,"value":"----\n#### *A:* Jason Davies implements Bridson’s.","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined"},{"id":24,"value":"show(a, visibility)","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1225,"value":"Bridson’s algorithm keeps a queue of “live” points. At each step, it picks one of the points *p* from the queue, chooses up to *k* (here *k* = 40) new candidate points of a distance between *r* and 2*r* from *p*, and checks whether each candidate is at least distance *r* away from every previous point. If one of the candidates meets this check, it is added to the list of points and a new step begins. If none of the candidates meets this test, then *p* is removed from the queue. See [@mbostock/poisson-disc-distribution](/@mbostock/poisson-disc-distribution) to watch the algorithm animate one point at a time.\n\nTo speed candidate checks, the algorithm maintains a “grid” of square bins of size chosen so that no more than 1 point can ever lie in the same bin. To check a candidate we do not need to test every previous point, but only those lying within the nearest 5x5 block of grid bins.","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined"},{"id":19,"value":"$32","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":201,"value":"a = Array.from(A(width, height, N))","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":237,"value":"---\n#### *B*: Martin Roberts tweaks Davies’s.","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined"},{"id":226,"value":"show(b, visibility)","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1224,"value":"Martin Roberts realized that Bridson’s algorithm always requires testing at least *k* points per point in the final set, and tried to reduce the value of *k*. Instead of picking *k* = 40 candidates at random inside a ring between radius *r* and 2*r*, he instead picks *k* = 4 candidates equally spaced around a circle of radius *r*, with random starting angle.\n\nThis requires testing on order 10 times fewer points than Bridson’s algorithm, but has the problem that many points end up in chains each of distance exactly *r* apart, with wide gaps between chains. This is a noticeably different point distribution than dart throwing would produce.","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined"},{"id":214,"value":"$33","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":223,"value":"b = Array.from(B(width, height, N))","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":239,"value":"---\n#### *C*: Jacob Rus revises Roberts’s.","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined"},{"id":289,"value":"show(c, visibility)","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1223,"value":"We (Jacob and Fil) discussed the deficits of these two algorithms in 2021, and tried to combine their strengths, balancing performance with distribution quality. Here we choose up to *k* = 11 candidates per starting point, with candidates at equally spaced angles, but with randomized radii in a suitable distribution. Parameters were tweaked by eye until the result looked good (feel free to fiddle with it more yourself).\n\nJacob sped up performance further in 2022 (especially in Firefox) by switching away from using an array of 2-element arrays to store the “grid” in favor of a single large array of floating point coordinates, and padding the grid on the sides to eliminate the need for special-case code near edges.","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined"},{"id":263,"value":"$34","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":273,"value":"c = Array.from(C(width, height, N))","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":805,"value":"---\n#### For comparison: *Dart throwing*\n\nThe code below shows how dart throwing works: we would need to set radius to the maximum value that still allows to reach N points in a reasonable time. But if the radius is too high or unlucky, it gets stuck in a situation where it can’t find any more new points before it has yielded enough of them. Set the radius lower than that, and the quality of the generator sinks.","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined"},{"id":802,"value":"/*\n function* D(width, height, N) {\n let radius = Math.sqrt((width * height) / N) * 0.84;\n const quad = d3.quadtree();\n let steps = 0;\n do {\n const point = [rng() * width, rng() * height];\n if (quad.find(...point, radius) === undefined) {\n quad.add(point);\n yield point;\n }\n } while (steps++ < 1e7);\n}\n*/","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":295,"value":"---\nOne needs a sharp eye to see the difference between these distributions, but we can compute a few metrics:","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined"},{"id":299,"value":"### Speed\n\nWe run each generator 5 times, and plot the sum of time spent *(smaller is faster)*:","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined"},{"id":336,"value":"Plot.barX(times, {\n x: \"time\",\n y: \"generator\",\n stroke: \"white\",\n fill: \"generator\",\n sort: \"time\",\n reverse: true\n}).plot({\n x: { nice: true, label: 'time (ms) →'},\n width\n})","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1208,"value":"*Note: the results depend on the browser; on Firefox, C is the fastest.*","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined"},{"id":339,"value":"viewof replayTimes = {\n await visibility();\n return Inputs.button(\"replay\");\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":304,"value":"times = {\n replayTimes;\n const k = 5; // number of times to repeat each generator\n const timer = (f) => {\n const p = performance.now();\n f();\n return performance.now() - p;\n };\n return Promise.all(\n d3.shuffle(d3.range(k).flatMap(i => [A, B, C])).map(\n (f) =>\n new Promise(async (resolve) => {\n await Promises.delay(500 * rng());\n resolve({\n generator: f.name,\n time: timer(() => Array.from(f(width, height, N)))\n });\n })\n )\n );\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":357,"value":"---\n### Area\n\nEach point in the distribution commands an area (the surface of its voronoi cell); we can check the distribution of these values. Colors are White: mean; Red: small; blue: large.","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined"},{"id":361,"value":"{\n await visibility();\n return html`\nA ${showVoronoi(A, 290, 300, \"area\")}\nB ${showVoronoi(B, 290, 300, \"area\")}\nC ${showVoronoi(C, 290, 300, \"area\")}\n`;\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":449,"value":"The differences in areas between the three contenders is hard to judge by eye. But the quality of C appears to be better, with less deviation in the cells’ areas—in particular fewer large cells:","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined"},{"id":457,"value":"{\n const data = [A, B, C].flatMap((f) =>\n areas(f).map((area) => ({ generator: f.name, area }))\n );\n return Plot.plot({\n facet: { data, x: \"generator\" },\n marks: [\n Plot.rectY(\n data,\n Plot.binX(\n { y: \"count\" },\n { x: \"area\", fill: \"generator\", curve: \"basis\", thresholds: 40 }\n )\n ),\n Plot.ruleY([0])\n ],\n width,\n height: 300\n });\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":607,"value":"Plot.barX(deviations, {\n x: \"deviation\",\n y: \"generator\",\n stroke: \"white\",\n fill: \"generator\",\n reverse: true\n}).plot({\n x: { nice: true }\n})","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":705,"value":"Plot.plot({\n marks: [\n Plot.barX(outliers, {\n x: \"outliers\",\n y: \"generator\",\n fill: \"generator\",\n stroke: \"white\",\n strokeWidth: 2,\n order: \"type\"\n }),\n Plot.textX(\n outliers,\n Plot.stackX({\n x: \"outliers\",\n y: \"generator\",\n fill: \"white\",\n text: \"type\",\n order: \"type\"\n })\n )\n ],\n x: { nice: true }\n})","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1069,"value":"viewof replayCharts = {\n await visibility();\n return Inputs.button(\"replay\")\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":684,"value":"*Comment: C is winning here, with fewer outliers. Furthermore its outliers tend to be small polygons (more clustered points), but they rarely form large gaps—which tend to cripple B.*","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined"},{"id":459,"value":"function areas(M) {\n replayCharts;\n const delaunay = d3.Delaunay.from(M(width, height, 10000, 500));\n const voronoi = delaunay.voronoi([0, 0, width, height]);\n const N = delaunay.points.length / 2;\n const I = d3\n .range(N)\n .filter(\n (i) =>\n !(\n delaunay.points[2 * i] < 10 ||\n delaunay.points[2 * i] > width - 10 ||\n delaunay.points[2 * i + 1] < 10 ||\n delaunay.points[2 * i + 1] > height - 10\n )\n );\n return Array.from(d3.range(N), (i) =>\n Math.abs(d3.polygonArea(voronoi.cellPolygon(i)))\n );\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":614,"value":"deviations = {\n replayCharts;\n return [A, B, C].map((M) => ({\n deviation: d3.deviation(areas(M)),\n generator: M.name\n }));\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":699,"value":"outliers = [A, B, C].flatMap((M) => {\n const areasM = areas(M);\n const mean = d3.mean(areasM);\n return [\n {\n outliers: areasM.filter((d) => d < mean * .7).length,\n type: \"small\",\n generator: M.name\n },\n {\n outliers: areasM.filter((d) => d > mean * 1.35).length,\n type: \"large\",\n generator: M.name\n }\n]})","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1064,"value":"* * *\n\n### Distance\n\nIf we alternately draw the Delaunay triangulation, we can color the edges by distance between adjacent points. Colors are White: root mean square; Red: short; Blue: long.","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined"},{"id":1063,"value":"{\n await visibility();\n return html`\nA ${showDelaunay(A, 290, 300, \"edges\")}\nB ${showDelaunay(B, 290, 300, \"edges\")}\nC ${showDelaunay(C, 290, 300, \"edges\")}\n`;\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1062,"value":"What really distinguishes B is that it has an extreme number of barely long-enough edges. But C also has fewer unusually long edges than B. In tweaking the candidate selection parameters/algorithm for C, there is a trade-off between a smooth and tight distribution for Delaunay edge lengths and a tight distribution for Voronoi cell areas.","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined"},{"id":1061,"value":"{\n const data = [A, B, C].flatMap((f) =>\n lengths(f).map((length) => ({ generator: f.name, length }))\n );\n return Plot.plot({\n facet: { data, x: \"generator\" },\n marks: [\n Plot.frame(),\n Plot.rectY(\n data,\n Plot.binX(\n { y: \"count\" },\n { x: \"length\", fill: \"generator\", interval: 0.25 }\n )\n ),\n Plot.ruleY([0]),\n Plot.text(\n data,\n Plot.selectMaxY(\n // selectMaxY ensures we only ever have 0 or 1 label per facet\n Plot.binX(\n {\n y: (data) => (data.length > 3500 ? 3200 : NaN),\n text: (data) => `← ${d3.format(\",\")(data.length)} edges`\n },\n { x: \"length\", thresholds: 40, textAnchor: \"start\", dx: 6 }\n )\n )\n )\n ],\n width,\n height: 300,\n marginLeft: 45,\n y: { domain: [0, 3500] },\n x: { domain: [4, 14] }\n });\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1060,"value":"Plot.barX(length_deviations, {\n x: \"deviation\",\n y: \"generator\",\n stroke: \"white\",\n fill: \"generator\",\n reverse: true\n}).plot({\n x: { nice: true }\n})","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1059,"value":"length_deviations = {\n replayCharts;\n return [A, B, C].map((M) => ({\n deviation: d3.deviation(lengths(M)),\n generator: M.name\n }));\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1058,"value":"function lengths(M) {\n replayCharts;\n const delaunay = d3.Delaunay.from(M(width, height, 10000, 500));\n const { points, halfedges, triangles } = delaunay;\n\n const startPt = (i) => {\n let ti = triangles[i];\n return points.slice(ti * 2, ti * 2 + 2);\n };\n const endPt = (i) => {\n let tj = triangles[halfedges[i]];\n return points.slice(tj * 2, tj * 2 + 2);\n };\n const edgeLength = (i) => {\n const [xi, yi] = startPt(i);\n const [xj, yj] = endPt(i);\n return Math.hypot(xj - xi, yj - yi);\n };\n\n const N = halfedges.length;\n const I = d3\n .range(N)\n .filter((i) => halfedges[i] > i)\n .filter((i) => {\n const [xi, yi] = startPt(i);\n const [xj, yj] = endPt(i);\n return !(\n (xi < 10 && xj < 10) ||\n (xi > width - 10 && xj > width - 10) ||\n (yi < 10 && yj < 10) ||\n (yi > height - 10 && yj > height - 10)\n );\n });\n\n return Array.from(I, (i) => edgeLength(i));\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":634,"value":"---\n### Shape\n\nWe also want to compare the distribution of voronoi cells for their number of edges: how many cells are hexagons (preferred), pentagons, heptagons (acceptable), or other shapes…?","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined"},{"id":1056,"value":"{\n await visibility();\n return html`\nA ${showVoronoi(A, 290, 300, \"shape\")}\nB ${showVoronoi(B, 290, 300, \"shape\")}\nC ${showVoronoi(C, 290, 300, \"shape\")}\n`;\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":651,"value":"Plot.plot({\n facet: { data: ngons, y: \"generator\" },\n marks: [\n Plot.barX(\n ngons,\n Plot.groupZ({ x: \"proportion-facet\" }, { fill: \"sides\", order: \"z\" })\n ),\n Plot.textX(\n ngons,\n Plot.stackX(\n Plot.groupZ(\n {\n x: \"proportion-facet\",\n text: ([d]) => (d.sides > 3 && d.sides < 9 ? d.sides : \"\")\n },\n { z: \"sides\", order: \"z\" }\n )\n )\n )\n ],\n x: { nice: true, percent: true }\n})","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":690,"value":"*Comment: B has more weird shapes than A. C is winning on hexagons.*","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined"},{"id":638,"value":"function sides(M) {\n const delaunay = d3.Delaunay.from(M(width, height, 10000));\n const voronoi = delaunay.voronoi([0, 0, width, height]);\n const N = delaunay.points.length / 2;\n const I = d3\n .range(N)\n .filter(\n (i) =>\n !(\n delaunay.points[2 * i] < 10 ||\n delaunay.points[2 * i] > width - 10 ||\n delaunay.points[2 * i + 1] < 10 ||\n delaunay.points[2 * i + 1] > height - 10\n )\n );\n return Array.from(d3.range(N), (i) => voronoi.cellPolygon(i).length - 1);\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":644,"value":"ngons = [A, B, C].flatMap((M) =>\n sides(M).map((sides) => ({\n sides,\n generator: M.name\n }))\n)","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":744,"value":"---\n\n_How to decide?_\n\nThe level of order _vs_ entropy counts: in fact, the regular hexagonal grid is\nthe _best_ if we only care about optimizing the above metrics. Speed counts,\ntoo. Maybe we need some other kinds of metrics. Given its higher quality and\ngood speed, _C_ has become the default “poisson” method in our\n[2D point distributions](https://observablehq.com/@fil/2d-point-distributions)\ntool.\n","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined"},{"id":573,"value":"---\n*Some supporting code*","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined"},{"id":364,"value":"$35","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1053,"value":"$36","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":197,"value":"rng = Math.random // import {…} with {rng} to change the random number generator","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":14,"value":"N = 4000","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":16,"value":"height = 400","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":270,"value":"async function* show(z, visibility) {\n await visibility();\n const context = DOM.context2d(width, height);\n yield context.canvas;\n context.fillStyle = \"rgba(20, 150, 240, 0.6)\";\n\n let i = 0;\n const radius = Math.sqrt((width * height) / N / 2.5);\n for (const p of z) {\n context.beginPath();\n context.arc(p[0], p[1], 0.65 * radius, 0, 2 * Math.PI);\n context.fill();\n\n if (i++ % 40 === 0) yield context.canvas;\n }\n}","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"}],"files":"$W37","stdlib":"1","titling":"auto"}} 30:"$2f:body"