Lloyd-Relaxed Stippling

// Stippling via Lloyd's relaxation on a procedural density field.
// Each iteration each point moves toward the density-weighted centroid
// of its (approximate) Voronoi cell. Converges to blue noise whose
// local site density tracks the underlying field.

const N = 900;
const GW = 96, GH = 64;          // density grid
const BW = 12, BH = 8, BCAP = 48; // spatial buckets for nearest-site lookup
const STEP = 0.55;

let W, H, cellW, cellH;
let density, densitySum;
let px, py, qx, qy, qw;
let owner, buckets, bucketCount;
let iter, needReseed, buf, bctx;

function hash(ix, iy, seed) {
  let h = ix * 374761393 + iy * 668265263 + seed * 1013904223;
  h = (h ^ (h >> 13)) * 1274126177;
  h = h ^ (h >> 16);
  return ((h >>> 0) % 10000) / 10000;
}
const smooth = (t) => t * t * (3 - 2 * t);

function vnoise(x, y, s, seed) {
  const sx = x * s, sy = y * s;
  const x0 = Math.floor(sx), y0 = Math.floor(sy);
  const fx = sx - x0, fy = sy - y0;
  const a = hash(x0, y0, seed), b = hash(x0 + 1, y0, seed);
  const c = hash(x0, y0 + 1, seed), d = hash(x0 + 1, y0 + 1, seed);
  const ux = smooth(fx), uy = smooth(fy);
  return (a * (1 - ux) + b * ux) * (1 - uy) + (c * (1 - ux) + d * ux) * uy;
}

function fbm(x, y, seed) {
  let amp = 0.55, freq = 1, sum = 0, norm = 0;
  for (let o = 0; o < 4; o++) {
    sum += amp * vnoise(x, y, freq * 2.2, seed + o * 17);
    norm += amp; amp *= 0.55; freq *= 2.05;
  }
  return sum / norm;
}

function buildDensity(seed) {
  density = new Float32Array(GW * GH);
  let total = 0;
  for (let gy = 0; gy < GH; gy++) {
    for (let gx = 0; gx < GW; gx++) {
      const nx = gx / GW, ny = gy / GH;
      let n = fbm(nx, ny, seed);
      const dx = nx - 0.5, dy = ny - 0.5;
      n = Math.max(0, Math.min(1, (n - 0.35) * 1.7 + Math.exp(-(dx * dx + dy * dy) * 6) * 0.35));
      n = n * n * (3 - 2 * n);
      density[gy * GW + gx] = n;
      total += n;
    }
  }
  densitySum = total;
}

function sampleDensity(x, y) {
  const gx = Math.max(0, Math.min(GW - 1, (x / W) * GW | 0));
  const gy = Math.max(0, Math.min(GH - 1, (y / H) * GH | 0));
  return density[gy * GW + gx];
}

function reseed(seed) {
  buildDensity(seed);
  px = new Float32Array(N); py = new Float32Array(N);
  let placed = 0, guard = 0;
  while (placed < N && guard < N * 200) {
    const x = Math.random() * W, y = Math.random() * H;
    if (Math.random() < sampleDensity(x, y) * 0.95 + 0.05) {
      px[placed] = x; py[placed] = y; placed++;
    }
    guard++;
  }
  for (; placed < N; placed++) { px[placed] = Math.random() * W; py[placed] = Math.random() * H; }
  qx = new Float32Array(N); qy = new Float32Array(N); qw = new Float32Array(N);
  owner = new Int16Array(GW * GH);
  buckets = new Int16Array(BW * BH * BCAP);
  bucketCount = new Int16Array(BW * BH);
  iter = 0;
}

function rebuildBuckets() {
  bucketCount.fill(0);
  for (let i = 0; i < N; i++) {
    const bx = Math.max(0, Math.min(BW - 1, (px[i] / W) * BW | 0));
    const by = Math.max(0, Math.min(BH - 1, (py[i] / H) * BH | 0));
    const b = by * BW + bx, c = bucketCount[b];
    if (c < BCAP) { buckets[b * BCAP + c] = i; bucketCount[b] = c + 1; }
  }
}

function relaxStep() {
  rebuildBuckets();
  qx.fill(0); qy.fill(0); qw.fill(0);
  for (let gy = 0; gy < GH; gy++) {
    const cy = (gy + 0.5) * cellH;
    const by = Math.max(0, Math.min(BH - 1, (cy / H) * BH | 0));
    for (let gx = 0; gx < GW; gx++) {
      const w = density[gy * GW + gx];
      if (w < 0.001) { owner[gy * GW + gx] = -1; continue; }
      const cx = (gx + 0.5) * cellW;
      const bx = Math.max(0, Math.min(BW - 1, (cx / W) * BW | 0));
      let best = Infinity, bi = -1, radius = 1;
      while (bi < 0 && radius <= Math.max(BW, BH)) {
        const bx0 = Math.max(0, bx - radius), bx1 = Math.min(BW - 1, bx + radius);
        const by0 = Math.max(0, by - radius), by1 = Math.min(BH - 1, by + radius);
        for (let yy = by0; yy <= by1; yy++) {
          for (let xx = bx0; xx <= bx1; xx++) {
            const b = yy * BW + xx, cnt = bucketCount[b], base = b * BCAP;
            for (let k = 0; k < cnt; k++) {
              const i = buckets[base + k];
              const dx = cx - px[i], dy = cy - py[i], d = dx * dx + dy * dy;
              if (d < best) { best = d; bi = i; }
            }
          }
        }
        radius++;
      }
      if (bi >= 0) {
        owner[gy * GW + gx] = bi;
        qx[bi] += cx * w; qy[bi] += cy * w; qw[bi] += w;
      } else owner[gy * GW + gx] = -1;
    }
  }
  for (let i = 0; i < N; i++) {
    if (qw[i] > 1e-6) {
      px[i] += (qx[i] / qw[i] - px[i]) * STEP;
      py[i] += (qy[i] / qw[i] - py[i]) * STEP;
    } else {
      px[i] = ((Math.random() * GW | 0) + 0.5) * cellW;
      py[i] = ((Math.random() * GH | 0) + 0.5) * cellH;
    }
  }
  iter++;
}

function init({ ctx, width, height }) {
  W = width; H = height; cellW = W / GW; cellH = H / GH;
  reseed((Math.random() * 1e9) | 0);
  needReseed = false;
  buf = new OffscreenCanvas(GW, GH);
  bctx = buf.getContext('2d');
  ctx.fillStyle = '#f4efe4';
  ctx.fillRect(0, 0, W, H);
}

function tick({ ctx, width, height, input }) {
  if (width !== W || height !== H) { W = width; H = height; cellW = W / GW; cellH = H / GH; }

  if (input.consumeClicks().length > 0) needReseed = true;
  if (needReseed) { reseed((Math.random() * 1e9) | 0); needReseed = false; }

  relaxStep();

  // Paper background.
  ctx.fillStyle = '#f4efe4';
  ctx.fillRect(0, 0, W, H);

  // Subtle density wash beneath the stipples.
  const wash = bctx.createImageData(GW, GH);
  for (let i = 0; i < GW * GH; i++) {
    const d = density[i], v = 244 - d * 70;
    wash.data[i * 4] = v; wash.data[i * 4 + 1] = v - 5;
    wash.data[i * 4 + 2] = v - 16; wash.data[i * 4 + 3] = 90;
  }
  bctx.putImageData(wash, 0, 0);
  ctx.imageSmoothingEnabled = true;
  ctx.drawImage(buf, 0, 0, W, H);

  // Stipples: dot radius tapers with local density.
  ctx.fillStyle = '#161310';
  for (let i = 0; i < N; i++) {
    const r = 0.9 + sampleDensity(px[i], py[i]) * 1.7;
    ctx.beginPath();
    ctx.arc(px[i], py[i], r, 0, Math.PI * 2);
    ctx.fill();
  }

  ctx.fillStyle = 'rgba(22,19,16,0.85)';
  ctx.font = '12px ui-monospace, Menlo, monospace';
  ctx.textBaseline = 'top';
  ctx.fillText(`N = ${N}   iter = ${iter}   click: reseed`, 10, 10);
}