31:Td6e,function* TorusDiskSample(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.69 * (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(width * rng(), height * 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 + width) % width; // wrap const y = (parent[1] + r*dy + height) % height; // Accept candidates that are inside the allowed extent // and farther than 2 * radius to all existing samples. if (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) { let dx = grid[i] - x; let dy = grid[i + 1] - y; dx -= width * Math.round(dx * (1 / width)); dy -= height * Math.round(dy * (1 / height)); if (dx * dx + dy * dy < radius2) return false; } } return true; } function populate_grid(i, j, x, y) { const jpad = j + 2, ipad = i + 2 const index = 2 * (gridWidth * jpad + ipad); grid[index] = x; grid[index + 1] = y; } function flip(i, w) { return i - w + 2*w*(i < (w >>> 1)); } function sample(x, y) { const s = [x, y]; const w = gridWidth - 4, h = gridHeight - 4; const xg = (x * cellSize) | 0; const yg = (y * cellSize) | 0; populate_grid(xg, yg, x, y); // if we are near an edge of the grid, also save // the sample in the grid padding on the opposite side const xg_near_side = (xg >= w - 2) | (xg < 2); const yg_near_side = (yg >= h - 2) | (yg < 2); if (xg_near_side) { populate_grid(flip(xg, w), yg, x, y); if (yg_near_side) populate_grid(flip(xg, w), flip(yg, h), x, y); } if (yg_near_side) populate_grid(xg, flip(yg, h), x, y); queue.push(s); return s; } }32: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; }33: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; }34:[] 2f:{"title":"Disk Sampling on a Flat Torus","cells":[{"id":0,"value":"# Disk Sampling on a Flat Torus","pinned":false,"hidden":false,"mode":"md"},{"id":1225,"value":"With a few tweaks, we can take the algorithm from [@fil/poisson-distribution-generators](/@fil/poisson-distribution-generators) and make it “wrap” around the edges of a rectangle from top to bottom and side to side. That is, we can sample a [flat torus](https://en.wikipedia.org/wiki/Torus#Flat_torus).","pinned":false,"hidden":false,"mode":"md"},{"id":289,"value":"show(c, visibility)","pinned":false,"hidden":false,"mode":"ojs"},{"id":1454,"value":"viewof replay = {\n await visibility();\n return Inputs.button(\"Randomize and re-animate\");\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":263,"value":"$31","pinned":true,"hidden":false,"mode":"ojs"},{"id":273,"value":"c = replay, Array.from(TorusDiskSample(samplewidth, sampleheight, N))","pinned":true,"hidden":false,"mode":"ojs"},{"id":573,"value":"---\n*Some supporting code*","pinned":false,"hidden":false,"mode":"md"},{"id":364,"value":"$32","pinned":false,"hidden":false,"mode":"ojs"},{"id":1053,"value":"$33","pinned":false,"hidden":false,"mode":"ojs"},{"id":197,"value":"rng = Math.random // import {…} with {rng} to change the random number generator","pinned":true,"hidden":false,"mode":"ojs"},{"id":14,"value":"N = 1000","pinned":true,"hidden":false,"mode":"ojs"},{"id":16,"value":"height = 400","pinned":true,"hidden":false,"mode":"ojs"},{"id":1355,"value":"samplewidth = width/2","pinned":true,"hidden":false,"mode":"ojs"},{"id":1357,"value":"sampleheight = height/2","pinned":true,"hidden":false,"mode":"ojs"},{"id":270,"value":"async function* show(z, visibility) {\n await visibility();\n const context = DOM.context2d(width, height);\n yield context.canvas;\n const blue = \"rgba(20, 150, 240, 0.6)\";\n const red = \"rgba(230, 50, 10, 0.6)\";\n context.fillStyle = blue;\n\n const circle = (x, y, i, j) =>\n context.arc(x + i*samplewidth, y + j*sampleheight, 0.65 * radius, 0, 2 * Math.PI)\n \n let i = 0;\n const radius = Math.sqrt((samplewidth * sampleheight) / N / 2.5);\n for (const p of z) {\n context.fillStyle = blue;\n [[-.5, .5], [.5, -.5], [.5, 1.5], [1.5, .5]].map(([tx, ty]) => {\n context.beginPath(), circle(p[0], p[1], tx, ty), context.fill();\n });\n context.fillStyle = red;\n [[-.5, -.5], [-.5, 1.5], [ .5, .5], [1.5, -.5], [1.5, 1.5]].map(([tx, ty]) => {\n context.beginPath(), circle(p[0], p[1], tx, ty), context.fill();\n });\n \n if (i++ % 10 === 0) yield context.canvas;\n // if (i > N * .6) break; // for making the thumbnail\n }\n}","pinned":true,"hidden":false,"mode":"ojs"}],"files":"$W34","stdlib":"1","titling":"auto"}