Visualize voxel-density uniformity of N random points in a disk
2 changed files+429−0
README.mdadded+38−0View file
@@ -0,0 +1,38 @@
1+# random-points-in-disk
2+
3+How many random points does it take before a disk *looks* like a disk?
4+
5+[**View the live visualizer**](https://concept-collection.github.io/random-points-in-disk/)
6+
7+Draw *N* points uniformly at random from a disk, bin them into a grid of voxels, and
8+compare the resulting density map against the exact one. Each voxel's count is
9+essentially Poisson with mean μ = points per voxel, so its relative noise is 1/√μ —
10+which depends only on how many points land in a voxel, not on the size of the disk or
11+the resolution per se. The app plots three panels side by side:
12+
13+- **Ideal** — exact density, each voxel shaded by the fraction of it inside the disk
14+- **Sampled** — counts from *N* random points, on the same color scale
15+- **Noise** — the error in units of its expected standard deviation, (count − ideal)/√ideal
16+
17+Sliding *N* and the resolution shows the tradeoff directly: quadrupling the resolution
18+quadruples the number of voxels, so it takes 4× the points to hold the same noise level.
19+The stats row reports the measured coefficient of variation next to the 1/√μ prediction,
20+and how many points a target noise level would require.
21+
22+Statistics are computed only over voxels lying *entirely* inside the disk; boundary
23+voxels have a smaller expected count and would otherwise inflate the measured spread.
24+
25+## Motivation
26+
27+This is the discretization question behind isochromat-based MRI simulation: a voxel's
28+signal is a sum over the isochromats that landed in it, so randomly placed isochromats
29+inject noise of order 1/√μ on top of the physics being modeled.
30+
31+## Running locally
32+
33+No build step — it is a single static `index.html`. Open it directly, or serve the
34+directory:
35+
36+```bash
37+python3 -m http.server 8000
38+```
index.htmladded+391−0View file
@@ -0,0 +1,391 @@
1+<!doctype html>
2+<html lang="en">
3+<head>
4+<meta charset="utf-8" />
5+<meta name="viewport" content="width=device-width, initial-scale=1" />
6+<title>Random points in a disk — how many do you need?</title>
7+<meta name="description" content="Voxel-density heatmaps of N uniform random points in a disk: how many samples it takes before the disk looks smooth at a given resolution." />
8+<link rel="icon" href="data:image/svg+xml,%3Csvg xmlns='http://www.w3.org/2000/svg' viewBox='0 0 32 32'%3E%3Crect width='32' height='32' fill='%230b0e14'/%3E%3Ccircle cx='16' cy='16' r='11' fill='%23dd5136'/%3E%3C/svg%3E" />
9+<style>
10+ :root {
11+ --bg: #0b0e14;
12+ --panel: #141924;
13+ --panel-2: #1b2130;
14+ --line: #2a3243;
15+ --fg: #e7ecf5;
16+ --dim: #97a3b8;
17+ --accent: #f08a2a;
18+ }
19+ * { box-sizing: border-box; }
20+ body {
21+ margin: 0;
22+ background: var(--bg);
23+ color: var(--fg);
24+ font: 15px/1.55 system-ui, -apple-system, "Segoe UI", Roboto, sans-serif;
25+ -webkit-font-smoothing: antialiased;
26+ }
27+ .wrap { max-width: 1140px; margin: 0 auto; padding: 28px 20px 56px; }
28+ h1 { font-size: 1.5rem; margin: 0 0 .4rem; letter-spacing: -0.01em; }
29+ .lede { color: var(--dim); margin: 0 0 1.6rem; max-width: 70ch; }
30+ .lede code { color: var(--fg); background: var(--panel-2); padding: 1px 5px; border-radius: 4px; font-size: .9em; }
31+
32+ .card {
33+ background: var(--panel);
34+ border: 1px solid var(--line);
35+ border-radius: 10px;
36+ }
37+
38+ /* ---- controls ---- */
39+ .controls { padding: 16px 18px; display: grid; gap: 18px; grid-template-columns: 2fr 1fr auto; align-items: end; }
40+ @media (max-width: 720px) { .controls { grid-template-columns: 1fr; } }
41+ label { display: block; font-size: .78rem; text-transform: uppercase; letter-spacing: .06em; color: var(--dim); margin-bottom: 6px; }
42+ .val { color: var(--fg); text-transform: none; letter-spacing: 0; font-size: .95rem; font-variant-numeric: tabular-nums; }
43+ input[type=range] { width: 100%; accent-color: var(--accent); margin: 0; }
44+ select, button {
45+ font: inherit; color: var(--fg); background: var(--panel-2);
46+ border: 1px solid var(--line); border-radius: 7px; padding: 7px 10px;
47+ }
48+ select { width: 100%; }
49+ button { cursor: pointer; white-space: nowrap; }
50+ button:hover { border-color: var(--accent); }
51+
52+ /* ---- panels ---- */
53+ .panels { display: grid; grid-template-columns: repeat(3, 1fr); gap: 14px; margin: 14px 0; }
54+ @media (max-width: 900px) { .panels { grid-template-columns: 1fr; max-width: 380px; } }
55+ .panel { padding: 14px; }
56+ .panel h2 { font-size: .95rem; margin: 0 0 2px; }
57+ .panel p { margin: 0 0 10px; font-size: .82rem; color: var(--dim); min-height: 2.6em; }
58+ canvas {
59+ width: 100%; aspect-ratio: 1; display: block; border-radius: 6px;
60+ image-rendering: pixelated; image-rendering: crisp-edges;
61+ background: #000;
62+ }
63+ .bar { height: 9px; border-radius: 3px; margin-top: 10px; border: 1px solid var(--line); }
64+ .ticks { display: flex; justify-content: space-between; font-size: .72rem; color: var(--dim); font-variant-numeric: tabular-nums; margin-top: 3px; }
65+
66+ /* ---- stats ---- */
67+ .stats { padding: 6px 4px; display: grid; grid-template-columns: repeat(4, 1fr); gap: 1px; background: var(--line); border-radius: 10px; overflow: hidden; padding: 0; }
68+ @media (max-width: 820px) { .stats { grid-template-columns: repeat(2, 1fr); } }
69+ .stat { background: var(--panel); padding: 12px 14px; }
70+ .stat .k { font-size: .74rem; text-transform: uppercase; letter-spacing: .06em; color: var(--dim); }
71+ .stat .v { font-size: 1.15rem; font-variant-numeric: tabular-nums; margin-top: 2px; }
72+ .stat .v small { font-size: .72rem; color: var(--dim); font-variant-numeric: normal; }
73+
74+ .note { color: var(--dim); font-size: .88rem; margin-top: 22px; max-width: 78ch; }
75+ .note b { color: var(--fg); font-weight: 600; }
76+ a { color: var(--accent); }
77+</style>
78+</head>
79+<body>
80+<div class="wrap">
81+
82+ <h1>How many random points does a disk need?</h1>
83+ <p class="lede">
84+ Draw <em>N</em> points uniformly at random from a disk and bin them into a grid of
85+ voxels. Each voxel's count is a Poisson-ish random variable with mean
86+ <code>μ</code> = points per voxel, so its relative noise is
87+ <code>1/√μ</code> — independent of how big the disk is. Sweep <em>N</em> and the
88+ resolution below to see when the sampled disk stops looking speckled and starts
89+ looking like a disk.
90+ </p>
91+
92+ <div class="card controls">
93+ <div>
94+ <label for="n">Points <span class="val" id="nVal"></span></label>
95+ <input type="range" id="n" min="200" max="750" step="1" value="500" />
96+ </div>
97+ <div>
98+ <label for="res">Resolution (voxels across diameter)</label>
99+ <select id="res">
100+ <option value="8">8 × 8</option>
101+ <option value="16">16 × 16</option>
102+ <option value="32">32 × 32</option>
103+ <option value="64" selected>64 × 64</option>
104+ <option value="128">128 × 128</option>
105+ <option value="256">256 × 256</option>
106+ <option value="512">512 × 512</option>
107+ </select>
108+ </div>
109+ <div><button id="resample">New sample</button></div>
110+ </div>
111+
112+ <div class="panels">
113+ <div class="card panel">
114+ <h2>Ideal</h2>
115+ <p>Exact density: each voxel shaded by the fraction of it inside the disk.</p>
116+ <canvas id="cIdeal"></canvas>
117+ <div class="bar" id="barSeq"></div>
118+ <div class="ticks"><span>0</span><span>0.8×</span><span>1.6× ideal</span></div>
119+ </div>
120+ <div class="card panel">
121+ <h2>Sampled</h2>
122+ <p>Counts from <span id="nInline"></span> random points, on the same color scale.</p>
123+ <canvas id="cSampled"></canvas>
124+ <div class="bar" id="barSeq2"></div>
125+ <div class="ticks"><span>0</span><span>0.8×</span><span>1.6× ideal</span></div>
126+ </div>
127+ <div class="card panel">
128+ <h2>Noise</h2>
129+ <p>Error in units of the expected standard deviation, <span style="white-space:nowrap">(count − ideal) / √ideal</span>.</p>
130+ <canvas id="cNoise"></canvas>
131+ <div class="bar" id="barDiv"></div>
132+ <div class="ticks"><span>−3σ</span><span>0</span><span>+3σ</span></div>
133+ </div>
134+ </div>
135+
136+ <div class="stats">
137+ <div class="stat"><div class="k">Points per voxel μ</div><div class="v" id="sMu"></div></div>
138+ <div class="stat"><div class="k">Measured noise</div><div class="v" id="sCv"></div></div>
139+ <div class="stat"><div class="k">Predicted 1/√μ</div><div class="v" id="sPred"></div></div>
140+ <div class="stat"><div class="k">Interior voxels</div><div class="v" id="sK"></div></div>
141+ <div class="stat"><div class="k">Count range (interior)</div><div class="v" id="sRange"></div></div>
142+ <div class="stat"><div class="k">Voxel width</div><div class="v" id="sVox"></div></div>
143+ <div class="stat">
144+ <div class="k">Points for <select id="target" style="width:auto;padding:1px 4px;font-size:.74rem;background:transparent;border:none;color:var(--accent)"><option value="0.10">10%</option><option value="0.05" selected>5%</option><option value="0.02">2%</option><option value="0.01">1%</option></select> noise</div>
145+ <div class="v" id="sNeed"></div>
146+ </div>
147+ <div class="stat"><div class="k">Sample time</div><div class="v" id="sTime"></div></div>
148+ </div>
149+
150+ <p class="note">
151+ <b>Why this matters for MRI simulation.</b> Discretizing a continuous object into
152+ isochromats is exactly this experiment: the signal in a voxel is a sum over the
153+ isochromats that landed in it, so random placement injects noise of order
154+ <code>1/√μ</code> on top of the physics. Statistics are computed only over
155+ voxels lying <em>entirely</em> inside the disk — boundary voxels have a smaller
156+ expected count and would otherwise inflate the measured spread. Note that the noise
157+ panel looks like white noise at every setting: what improves with more points is not
158+ its character but its size relative to the signal.
159+ </p>
160+
161+</div>
162+
163+<script>
164+const TAU = Math.PI * 2;
165+
166+// ---------- colormaps ----------
167+const INFERNO = [
168+ [0,0,4],[20,11,53],[66,10,104],[106,23,110],[147,38,103],
169+ [188,55,84],[221,81,58],[249,142,9],[252,255,164],
170+];
171+const DIVERGING = [
172+ [33,102,172],[103,169,207],[236,238,241],[239,138,98],[178,24,43],
173+];
174+
175+function makeLut(stops) {
176+ const lut = new Uint8Array(256 * 3);
177+ const n = stops.length - 1;
178+ for (let i = 0; i < 256; i++) {
179+ const s = (i / 255) * n;
180+ const j = Math.min(n - 1, Math.floor(s));
181+ const f = s - j;
182+ const a = stops[j], b = stops[j + 1];
183+ lut[i * 3 + 0] = a[0] + (b[0] - a[0]) * f;
184+ lut[i * 3 + 1] = a[1] + (b[1] - a[1]) * f;
185+ lut[i * 3 + 2] = a[2] + (b[2] - a[2]) * f;
186+ }
187+ return lut;
188+}
189+const LUT_SEQ = makeLut(INFERNO);
190+const LUT_DIV = makeLut(DIVERGING);
191+
192+function cssGradient(lut) {
193+ const parts = [];
194+ for (let i = 0; i <= 16; i++) {
195+ const k = Math.round((i / 16) * 255) * 3;
196+ parts.push(`rgb(${lut[k]},${lut[k+1]},${lut[k+2]}) ${(i / 16 * 100).toFixed(1)}%`);
197+ }
198+ return `linear-gradient(to right, ${parts.join(',')})`;
199+}
200+
201+// ---------- geometry: exact-ish area fraction of each voxel inside the unit disk ----------
202+// Domain is [-1,1]^2, so the grid spans the full diameter. Sub-sampled per voxel;
203+// cached because it only depends on the resolution.
204+const fracCache = new Map();
205+function areaFractions(G) {
206+ const hit = fracCache.get(G);
207+ if (hit) return hit;
208+ const S = G >= 256 ? 4 : 8; // sub-samples per axis
209+ const frac = new Float32Array(G * G);
210+ const h = 2 / G, sh = h / S;
211+ for (let iy = 0; iy < G; iy++) {
212+ const y0 = -1 + iy * h;
213+ for (let ix = 0; ix < G; ix++) {
214+ const x0 = -1 + ix * h;
215+ let inside = 0;
216+ for (let sy = 0; sy < S; sy++) {
217+ const y = y0 + (sy + 0.5) * sh;
218+ for (let sx = 0; sx < S; sx++) {
219+ const x = x0 + (sx + 0.5) * sh;
220+ if (x * x + y * y <= 1) inside++;
221+ }
222+ }
223+ frac[iy * G + ix] = inside / (S * S);
224+ }
225+ }
226+ fracCache.set(G, frac);
227+ return frac;
228+}
229+
230+// ---------- sampling ----------
231+function sampleCounts(N, G) {
232+ const counts = new Uint32Array(G * G);
233+ const half = G / 2;
234+ for (let i = 0; i < N; i++) {
235+ const r = Math.sqrt(Math.random()); // sqrt for uniform area density
236+ const t = TAU * Math.random();
237+ let ix = (r * Math.cos(t) + 1) * half | 0;
238+ let iy = (r * Math.sin(t) + 1) * half | 0;
239+ if (ix >= G) ix = G - 1;
240+ if (iy >= G) iy = G - 1;
241+ counts[iy * G + ix]++;
242+ }
243+ return counts;
244+}
245+
246+// ---------- drawing ----------
247+const SCALE_MAX = 1.6; // top of the sequential color scale, in units of the ideal density
248+
249+function paint(canvas, G, valueAt) {
250+ canvas.width = G;
251+ canvas.height = G;
252+ const ctx = canvas.getContext('2d');
253+ const img = ctx.createImageData(G, G);
254+ const d = img.data;
255+ for (let p = 0; p < G * G; p++) {
256+ const rgb = valueAt(p);
257+ d[p * 4 + 0] = rgb[0];
258+ d[p * 4 + 1] = rgb[1];
259+ d[p * 4 + 2] = rgb[2];
260+ d[p * 4 + 3] = 255;
261+ }
262+ ctx.putImageData(img, 0, 0);
263+}
264+
265+const OUTSIDE = [22, 26, 35]; // neutral fill for voxels with no expected signal
266+const rgbSeq = (t) => {
267+ const k = Math.max(0, Math.min(255, Math.round(t * 255))) * 3;
268+ return [LUT_SEQ[k], LUT_SEQ[k + 1], LUT_SEQ[k + 2]];
269+};
270+const rgbDiv = (t) => {
271+ const k = Math.max(0, Math.min(255, Math.round(t * 255))) * 3;
272+ return [LUT_DIV[k], LUT_DIV[k + 1], LUT_DIV[k + 2]];
273+};
274+
275+// ---------- formatting ----------
276+function fmtCount(n) {
277+ if (n >= 1e6) {
278+ const m = n / 1e6;
279+ return (m >= 100 ? m.toFixed(0) : m.toFixed(m >= 10 ? 1 : 2)) + 'M';
280+ }
281+ if (n >= 1e4) return (n / 1e3).toFixed(0) + 'k';
282+ return n.toLocaleString('en-US');
283+}
284+function fmtMu(mu) {
285+ if (mu >= 1000) return fmtCount(Math.round(mu));
286+ if (mu >= 10) return mu.toFixed(0);
287+ if (mu >= 1) return mu.toFixed(2);
288+ return mu.toFixed(3);
289+}
290+function fmtPct(x) {
291+ if (!isFinite(x)) return '—';
292+ if (x >= 1) return (x * 100).toFixed(0) + '%';
293+ if (x >= 0.1) return (x * 100).toFixed(1) + '%';
294+ return (x * 100).toFixed(2) + '%';
295+}
296+
297+// ---------- app ----------
298+const el = (id) => document.getElementById(id);
299+const nSlider = el('n'), resSel = el('res'), targetSel = el('target');
300+
301+document.getElementById('barSeq').style.background = cssGradient(LUT_SEQ);
302+document.getElementById('barSeq2').style.background = cssGradient(LUT_SEQ);
303+document.getElementById('barDiv').style.background = cssGradient(LUT_DIV);
304+
305+const sliderToN = (v) => {
306+ const n = Math.pow(10, v / 100);
307+ // round to 3 significant figures so the readout is stable while dragging
308+ const mag = Math.pow(10, Math.floor(Math.log10(n)) - 2);
309+ return Math.max(1, Math.round(n / mag) * mag);
310+};
311+
312+function showN(N) {
313+ el('nVal').textContent = '= ' + N.toLocaleString('en-US');
314+ el('nInline').textContent = fmtCount(N);
315+}
316+
317+let lastElapsed = 0;
318+
319+function update() {
320+ const N = sliderToN(+nSlider.value);
321+ const G = +resSel.value;
322+
323+ showN(N);
324+
325+ const frac = areaFractions(G);
326+ const t0 = performance.now();
327+ const counts = sampleCounts(N, G);
328+ const elapsed = lastElapsed = performance.now() - t0;
329+
330+ // Expected count in a fully interior voxel: N * (voxel area / disk area).
331+ const mu = N * (4 / (G * G)) / Math.PI;
332+
333+ paint(el('cIdeal'), G, (p) => rgbSeq(frac[p] / SCALE_MAX));
334+ paint(el('cSampled'), G, (p) => rgbSeq(counts[p] / mu / SCALE_MAX));
335+ paint(el('cNoise'), G, (p) => {
336+ const expected = mu * frac[p];
337+ if (expected <= 0) return OUTSIDE;
338+ const z = (counts[p] - expected) / Math.sqrt(expected);
339+ return rgbDiv(z / 6 + 0.5); // ±3σ across the scale
340+ });
341+
342+ // Statistics over voxels entirely inside the disk.
343+ let k = 0, sum = 0, sumSq = 0, lo = Infinity, hi = -Infinity;
344+ for (let p = 0; p < G * G; p++) {
345+ if (frac[p] < 1) continue;
346+ const c = counts[p];
347+ k++; sum += c; sumSq += c * c;
348+ if (c < lo) lo = c;
349+ if (c > hi) hi = c;
350+ }
351+ const mean = k ? sum / k : 0;
352+ const variance = k > 1 ? Math.max(0, sumSq / k - mean * mean) * k / (k - 1) : 0;
353+ const cv = mean > 0 ? Math.sqrt(variance) / mean : Infinity;
354+
355+ el('sMu').innerHTML = fmtMu(mu) + ' <small>expected</small>';
356+ el('sCv').textContent = k ? fmtPct(cv) : '—';
357+ el('sPred').textContent = mu > 0 ? fmtPct(1 / Math.sqrt(mu)) : '—';
358+ el('sK').innerHTML = k.toLocaleString('en-US') + ' <small>of ' + (G * G).toLocaleString('en-US') + '</small>';
359+ el('sRange').textContent = k ? lo.toLocaleString('en-US') + ' – ' + hi.toLocaleString('en-US') : '—';
360+ el('sVox').innerHTML = (2 / G).toPrecision(3).replace(/(\.\d*?)0+$/, '$1').replace(/\.$/, '') + ' <small>disk radii</small>';
361+
362+ const c = +targetSel.value;
363+ const need = Math.PI * G * G / (4 * c * c);
364+ el('sNeed').innerHTML = fmtCount(Math.round(need)) + ' <small>points</small>';
365+
366+ el('sTime').innerHTML = (elapsed < 10 ? elapsed.toFixed(1) : elapsed.toFixed(0)) + ' <small>ms</small>';
367+}
368+
369+// Coalesce rapid input into one update per frame; once a sample gets slow enough to
370+// stutter a drag, wait for the slider to settle instead.
371+let pending = false, timer = 0;
372+function schedule() {
373+ clearTimeout(timer);
374+ if (lastElapsed > 120) {
375+ timer = setTimeout(update, 180);
376+ return;
377+ }
378+ if (pending) return;
379+ pending = true;
380+ requestAnimationFrame(() => { pending = false; update(); });
381+}
382+
383+// The readout is cheap, so keep it live even while a heavy recompute is deferred.
384+nSlider.addEventListener('input', () => { showN(sliderToN(+nSlider.value)); schedule(); });
385+resSel.addEventListener('change', schedule);
386+targetSel.addEventListener('change', schedule);
387+el('resample').addEventListener('click', update);
388+update();
389+</script>
390+</body>
391+</html>