Mean Curvature Flow

// Mean curvature flow on a closed planar curve.
// Each point moves inward along the inward normal at speed equal to
// the local curvature: d/dt gamma = kappa * N.
// Gage-Hamilton-Grayson: any embedded simple closed curve flows to a
// round point.

const N_POINTS = 400;
const STEP_H = 0.35;       // per-substep time scale (pixels^2-ish)
const SUBSTEPS = 2;        // explicit Euler substeps per frame
const RESET_AREA_FRAC = 0.0015; // reset when area shrinks below this fraction of W*H
const TRAIL_COUNT = 16;    // how many faded past outlines to draw
const TRAIL_INTERVAL = 6;  // store an outline every N frames

let W, H, cx, cy, R0;
let pts;                   // Float64Array length 2*N_POINTS, [x0,y0,x1,y1,...]
let work;                  // scratch buffer for new positions
let trails;                // ring buffer of Float64Array snapshots
let trailHead;
let trailFilled;
let frameCounter;
let resetCooldown;         // frames before reset-check is allowed
let prevMouseDown;

function rand() {
  return Math.random();
}

function buildNoisyCurve() {
  // Circle radius R0 perturbed by sum of cos modes 3..9 with random phases & amps.
  const modes = [3, 4, 5, 6, 7, 8, 9];
  const phases = modes.map(() => rand() * Math.PI * 2);
  const amps = modes.map(() => (rand() * 0.5 - 0.25) * R0 * 0.18); // up to ~18% of R0 per mode
  // also a small global wobble
  const baseScale = 0.85 + rand() * 0.2;

  for (let i = 0; i < N_POINTS; i++) {
    const theta = (i / N_POINTS) * Math.PI * 2;
    let r = R0 * baseScale;
    for (let m = 0; m < modes.length; m++) {
      r += amps[m] * Math.cos(modes[m] * theta + phases[m]);
    }
    // safety: keep radius positive
    if (r < R0 * 0.15) r = R0 * 0.15;
    pts[i * 2] = cx + r * Math.cos(theta);
    pts[i * 2 + 1] = cy + r * Math.sin(theta);
  }
}

function curveArea() {
  // Shoelace
  let a = 0;
  for (let i = 0; i < N_POINTS; i++) {
    const j = (i + 1) % N_POINTS;
    const x1 = pts[i * 2];
    const y1 = pts[i * 2 + 1];
    const x2 = pts[j * 2];
    const y2 = pts[j * 2 + 1];
    a += x1 * y2 - x2 * y1;
  }
  return Math.abs(a) * 0.5;
}

function curvePerimeter() {
  let s = 0;
  for (let i = 0; i < N_POINTS; i++) {
    const j = (i + 1) % N_POINTS;
    const dx = pts[j * 2] - pts[i * 2];
    const dy = pts[j * 2 + 1] - pts[i * 2 + 1];
    s += Math.sqrt(dx * dx + dy * dy);
  }
  return s;
}

function stepMCF(h) {
  // Discrete Laplacian: Delta p_i ~= p_{i-1} - 2 p_i + p_{i+1}.
  // For an arc-length parameterized curve this equals kappa * N (inward).
  // So move p_i <- p_i + h * Delta p_i  (this is the implicit "inward kappa N" direction).
  for (let i = 0; i < N_POINTS; i++) {
    const im = (i - 1 + N_POINTS) % N_POINTS;
    const ip = (i + 1) % N_POINTS;
    const lx = pts[im * 2] - 2 * pts[i * 2] + pts[ip * 2];
    const ly = pts[im * 2 + 1] - 2 * pts[i * 2 + 1] + pts[ip * 2 + 1];
    work[i * 2] = pts[i * 2] + h * lx;
    work[i * 2 + 1] = pts[i * 2 + 1] + h * ly;
  }
  // swap
  const tmp = pts;
  pts = work;
  work = tmp;
}

function resample() {
  // Resample N_POINTS evenly along arc length to keep numerical Laplacian stable.
  const perim = curvePerimeter();
  if (!isFinite(perim) || perim < 1e-3) return;
  const seg = perim / N_POINTS;

  // build cumulative arc lengths
  const cum = new Float64Array(N_POINTS + 1);
  cum[0] = 0;
  for (let i = 0; i < N_POINTS; i++) {
    const j = (i + 1) % N_POINTS;
    const dx = pts[j * 2] - pts[i * 2];
    const dy = pts[j * 2 + 1] - pts[i * 2 + 1];
    cum[i + 1] = cum[i] + Math.sqrt(dx * dx + dy * dy);
  }

  let k = 0;
  for (let i = 0; i < N_POINTS; i++) {
    const target = i * seg;
    while (k < N_POINTS && cum[k + 1] < target) k++;
    const t = (target - cum[k]) / Math.max(1e-9, cum[k + 1] - cum[k]);
    const a = k % N_POINTS;
    const b = (k + 1) % N_POINTS;
    work[i * 2] = pts[a * 2] + t * (pts[b * 2] - pts[a * 2]);
    work[i * 2 + 1] = pts[a * 2 + 1] + t * (pts[b * 2 + 1] - pts[a * 2 + 1]);
  }
  const tmp = pts;
  pts = work;
  work = tmp;
}

function snapshotTrail() {
  const snap = new Float64Array(N_POINTS * 2);
  snap.set(pts);
  trails[trailHead] = snap;
  trailHead = (trailHead + 1) % TRAIL_COUNT;
  if (trailHead === 0) trailFilled = true;
}

function clearTrails() {
  for (let i = 0; i < TRAIL_COUNT; i++) trails[i] = null;
  trailHead = 0;
  trailFilled = false;
}

function reset() {
  buildNoisyCurve();
  clearTrails();
  frameCounter = 0;
  resetCooldown = 30;
}

function init({ ctx, width, height }) {
  W = width;
  H = height;
  cx = W * 0.5;
  cy = H * 0.5;
  R0 = Math.min(W, H) * 0.36;

  pts = new Float64Array(N_POINTS * 2);
  work = new Float64Array(N_POINTS * 2);
  trails = new Array(TRAIL_COUNT).fill(null);
  trailHead = 0;
  trailFilled = false;
  frameCounter = 0;
  resetCooldown = 30;
  prevMouseDown = false;

  buildNoisyCurve();

  ctx.fillStyle = '#08090c';
  ctx.fillRect(0, 0, W, H);
}

function drawCurve(ctx, arr, strokeStyle, lineWidth) {
  ctx.beginPath();
  ctx.moveTo(arr[0], arr[1]);
  for (let i = 1; i < N_POINTS; i++) {
    ctx.lineTo(arr[i * 2], arr[i * 2 + 1]);
  }
  ctx.closePath();
  ctx.strokeStyle = strokeStyle;
  ctx.lineWidth = lineWidth;
  ctx.stroke();
}

function tick({ ctx, width, height, input }) {
  if (width !== W || height !== H) {
    W = width;
    H = height;
    cx = W * 0.5;
    cy = H * 0.5;
    R0 = Math.min(W, H) * 0.36;
  }

  // Click to reset
  if (input && typeof input.consumeClicks === 'function') {
    if (input.consumeClicks() > 0) {
      reset();
    }
  }

  // Fade background slightly each frame for a soft afterimage
  ctx.globalCompositeOperation = 'source-over';
  ctx.fillStyle = '#08090c';
  ctx.fillRect(0, 0, W, H);

  // Advance flow
  for (let s = 0; s < SUBSTEPS; s++) {
    stepMCF(STEP_H);
  }
  // Resample occasionally to keep arc length uniform
  if (frameCounter % 3 === 0) resample();

  // Periodically store a trail snapshot
  if (frameCounter % TRAIL_INTERVAL === 0) snapshotTrail();

  // Draw trails (oldest -> newest, faint)
  const count = trailFilled ? TRAIL_COUNT : trailHead;
  for (let i = 0; i < count; i++) {
    const idx = (trailHead - 1 - i + TRAIL_COUNT) % TRAIL_COUNT;
    const arr = trails[idx];
    if (!arr) continue;
    const age = i / TRAIL_COUNT; // 0 = freshest trail, 1 = oldest
    const alpha = (1 - age) * 0.18 + 0.02;
    drawCurve(ctx, arr, `rgba(220, 225, 235, ${alpha.toFixed(3)})`, 1);
  }

  // Draw current curve crisp white
  drawCurve(ctx, pts, 'rgba(245, 247, 252, 0.95)', 1.6);

  // Reset condition: area below threshold
  if (resetCooldown > 0) {
    resetCooldown--;
  } else {
    const area = curveArea();
    if (area < RESET_AREA_FRAC * W * H || !isFinite(area)) {
      reset();
    }
  }

  frameCounter++;
}