Adaptive clearing passes

/** 1-D squared distance transform (Felzenszwalb & Huttenlocher 2012) of f into d (length n), with buffers v, z. */
function dt1(f, n, d, v, z) {
  let k = 0;
  v[0] = 0; z[0] = -INF; z[1] = INF;
  for (let q = 1; q < n; q++) {
    let s;
    for (;;) { const p = v[k]; s = ((f[q] + q * q) - (f[p] + p * p)) / (2 * q - 2 * p); if (s <= z[k] && k > 0) k--; else break; }
    if (s <= z[k]) { v[0] = q; z[0] = -INF; z[1] = INF; k = 0; continue; }
    k++; v[k] = q; z[k] = s; z[k + 1] = INF;
  }
  k = 0;
  for (let q = 0; q < n; q++) { while (z[k + 1] < q) k++; const p = v[k]; d[q] = (q - p) * (q - p) + f[p]; }
}

/** Euclidean distance (in cells) from every cell to the nearest cell where src is set. */
function edt(src, cols, rows) {
  const n = Math.max(cols, rows), f = new Float64Array(n), d = new Float64Array(n), v = new Int32Array(n), z = new Float64Array(n + 1);
  const g = new Float64Array(cols * rows);
  for (let x = 0; x < cols; x++) {
    for (let y = 0; y < rows; y++) f[y] = src[y * cols + x] ? 0 : INF;
    dt1(f, rows, d, v, z);
    for (let y = 0; y < rows; y++) g[y * cols + x] = d[y];
  }
  for (let y = 0; y < rows; y++) {
    for (let x = 0; x < cols; x++) f[x] = g[y * cols + x];
    dt1(f, cols, d, v, z);
    for (let x = 0; x < cols; x++) g[y * cols + x] = Math.sqrt(d[x]);
  }
  return g;
}

/** The raster of a region: cell centres, inside mask, exact distance to the walls. */
function rasterOf(loops, R, h0, maxCells) {
  const [bx0, by0, bx1, by1] = O.boxOfLoops(loops);
  const W = bx1 - bx0, Hh = by1 - by0;
  let h = h0;
  while ((W / h + 3) * (Hh / h + 3) > maxCells) h *= 1.25;
  const cols = Math.ceil(W / h) + 3, rows = Math.ceil(Hh / h) + 3;
  const x0 = bx0 - h, y0 = by0 - h;
  const inside = new Uint8Array(cols * rows), wall = new Float64Array(cols * rows);
  for (let y = 0; y < rows; y++) for (let x = 0; x < cols; x++) {
    const q = [x0 + x * h, y0 + y * h], i = y * cols + x;
    inside[i] = O.inside(loops, q) ? 1 : 0;
    wall[i] = inside[i] ? O.distToLoops(q, loops) : 0;
  }
  return { h, cols, rows, x0, y0, inside, wall, at: (x, y) => [x0 + x * h, y0 + y * h] };
}

/** Closed iso-contours {F = 0} of a cell-centre field (F > 0 outside the region of interest; the border is > 0). */
function contours(F, cols, rows) {
  const segs = new Map();   // edge key → [other edge key, point] pairs via adjacency
  const pts = new Map(), adj = new Map();
  const ek = (x, y, dir) => (y * cols + x) * 2 + dir;   // dir 0: edge to (x+1, y); 1: edge to (x, y+1)
  const cross = (x, y, dir) => {
    const k = ek(x, y, dir);
    if (!pts.has(k)) {
      const a = F[y * cols + x], b = dir ? F[(y + 1) * cols + x] : F[y * cols + x + 1], t = a / (a - b);
      pts.set(k, dir ? [x, y + t] : [x + t, y]);
    }
    return k;
  };
  const link = (a, b) => { (adj.get(a) || adj.set(a, []).get(a)).push(b); (adj.get(b) || adj.set(b, []).get(b)).push(a); };
  for (let y = 0; y + 1 < rows; y++) for (let x = 0; x + 1 < cols; x++) {
    const v = [F[y * cols + x], F[y * cols + x + 1], F[(y + 1) * cols + x + 1], F[(y + 1) * cols + x]];
    const c = (v[0] <= 0 ? 1 : 0) | (v[1] <= 0 ? 2 : 0) | (v[2] <= 0 ? 4 : 0) | (v[3] <= 0 ? 8 : 0);
    if (c === 0 || c === 15) continue;
    const E = [() => cross(x, y, 0), () => cross(x + 1, y, 1), () => cross(x, y + 1, 0), () => cross(x, y, 1)];   // bottom, right, top, left
    const T = {
      1: [[3, 0]], 2: [[0, 1]], 3: [[3, 1]], 4: [[1, 2]], 6: [[0, 2]], 7: [[3, 2]], 8: [[2, 3]], 9: [[2, 0]], 11: [[2, 1]], 12: [[1, 3]], 13: [[1, 0]], 14: [[0, 3]],
    }[c];
    if (T) { for (const [a, b] of T) link(E[a](), E[b]()); continue; }
    const mid = (v[0] + v[1] + v[2] + v[3]) / 4;
    if (c === 5) { if (mid <= 0) { link(E[3](), E[2]()); link(E[1](), E[0]()); } else { link(E[3](), E[0]()); link(E[1](), E[2]()); } }
    else if (c === 10) { if (mid <= 0) { link(E[0](), E[3]()); link(E[2](), E[1]()); } else { link(E[0](), E[1]()); link(E[2](), E[3]()); } }
  }
  void segs;
  const out = [], seen = new Set();
  for (const start of adj.keys()) {
    if (seen.has(start)) continue;
    const loop = [];
    let prev = null, cur = start;
    for (let guard = 0; guard < 1e7; guard++) {
      seen.add(cur); loop.push(pts.get(cur));
      const nx = (adj.get(cur) || []).find((k) => k !== prev && !seen.has(k)) ?? null;
      if (nx === null) break;
      prev = cur; cur = nx;
    }
    if (loop.length >= 3) out.push(loop);
  }
  return out;
}

// F ≤ 0: the cleared side at ≥ R − w from the material, and allowed
for (let i = 0; i < cols * rows; i++) Fv[i] = allowed[i] ? (R - w) - dM[i] * cell : Math.max(1, (reach + margin) - wall[i]) + cell;
const found = contours(Fv, cols, rows).map((l) => l.map(([x, y]) => [ras.x0 + x * cell, ras.y0 + y * cell]));

// never nearer the walls than R + leave (the raster's interpolation): nudge such points inward
lp = lp.map((p) => {
  const d = O.distToLoops(p, loops);
  if (d >= reach) return p;
  const near = loops.flat().reduce((bq, pr) => { const cq = O.closestOnPrim(p, pr); return !bq || cq.d < bq.d ? cq : bq; }, null);
  const ux = p[0] - near.p[0], uy = p[1] - near.p[1], ul = Math.hypot(ux, uy) || 1;
  return [near.p[0] + ux / ul * (reach + 1e-6), near.p[1] + uy / ul * (reach + 1e-6)];
});

/** Douglas–Peucker on a closed loop. */
function simplify(loop, tol) {
  if (loop.length < 4) return loop;
  const keep = new Uint8Array(loop.length);
  const rec = (a, b) => {
    let worst = -1, wi = -1;
    const A = loop[a], B = loop[b], dx = B[0] - A[0], dy = B[1] - A[1], L = Math.hypot(dx, dy) || 1;
    for (let i = a + 1; i < b; i++) { const d = Math.abs((loop[i][0] - A[0]) * dy - (loop[i][1] - A[1]) * dx) / L; if (d > worst) { worst = d; wi = i; } }
    if (worst > tol) { keep[wi] = 1; rec(a, wi); rec(wi, b); }
  };
  const far = loop.reduce((bi, p, i) => (Math.hypot(p[0] - loop[0][0], p[1] - loop[0][1]) > Math.hypot(loop[bi][0] - loop[0][0], loop[bi][1] - loop[0][1]) ? i : bi), 0);
  keep[0] = keep[far] = 1;
  rec(0, far);
  const tail = [...loop.slice(far), loop[0]];
  const keep2 = new Uint8Array(tail.length);
  const rec2 = (a, b) => {
    let worst = -1, wi = -1;
    const A = tail[a], B = tail[b], dx = B[0] - A[0], dy = B[1] - A[1], L = Math.hypot(dx, dy) || 1;
    for (let i = a + 1; i < b; i++) { const d = Math.abs((tail[i][0] - A[0]) * dy - (tail[i][1] - A[1]) * dx) / L; if (d > worst) { worst = d; wi = i; } }
    if (worst > tol) { keep2[wi] = 1; rec2(a, wi); rec2(wi, b); }
  };
  rec2(0, tail.length - 1);
  const out = loop.filter((_, i) => keep[i] && i <= far);
  for (let i = 1; i < tail.length - 1; i++) if (keep2[i]) out.push(tail[i]);
  return out;
}

/** Mark every cell within r of the polyline (closed when `closed`). */
function sweep(mask, ras, poly, r, closed) {
  const { h, cols, rows, x0, y0 } = ras, rr = r * r;
  const n = poly.length, m = closed ? n : n - 1;
  for (let i = 0; i < m; i++) {
    const a = poly[i], b = poly[(i + 1) % n];
    const xa = Math.max(0, Math.floor((Math.min(a[0], b[0]) - r - x0) / h)), xb = Math.min(cols - 1, Math.ceil((Math.max(a[0], b[0]) + r - x0) / h));
    const ya = Math.max(0, Math.floor((Math.min(a[1], b[1]) - r - y0) / h)), yb = Math.min(rows - 1, Math.ceil((Math.max(a[1], b[1]) + r - y0) / h));
    const dx = b[0] - a[0], dy = b[1] - a[1], L2 = dx * dx + dy * dy;
    for (let y = ya; y <= yb; y++) for (let x = xa; x <= xb; x++) {
      const i2 = y * cols + x;
      if (mask[i2]) continue;
      const px = x0 + x * h - a[0], py = y0 + y * h - a[1];
      const t = L2 > 0 ? Math.max(0, Math.min(1, (px * dx + py * dy) / L2)) : 0;
      const ex = px - t * dx, ey = py - t * dy;
      if (ex * ex + ey * ey <= rr) mask[i2] = 1;
    }
  }
}

/** Re-simulate the passes on a separate, finer raster: the engagement (material area removed per mm of path) and what is left. */
export function engagementCheck(loops, { R, entry, passes, closed = null, links = null, leave = 0, finish = null }, { h = R / 20 } = {}) {
  const ras = rasterOf(loops, R, h, 2e6), { cols, rows, inside, wall } = ras, cell = ras.h;
  const allowed = new Uint8Array(cols * rows);
  for (let i = 0; i < cols * rows; i++) allowed[i] = inside[i] && wall[i] >= R + leave ? 1 : 0;
  const dA = edt(allowed, cols, rows);
  const mat = new Uint8Array(cols * rows);
  for (let i = 0; i < cols * rows; i++) mat[i] = inside[i] ? 1 : 0;
  const clearDisk = (p, r) => {
    let n = 0;
    const xa = Math.max(0, Math.floor((p[0] - r - ras.x0) / cell)), xb = Math.min(cols - 1, Math.ceil((p[0] + r - ras.x0) / cell));
    const ya = Math.max(0, Math.floor((p[1] - r - ras.y0) / cell)), yb = Math.min(rows - 1, Math.ceil((p[1] + r - ras.y0) / cell));
    for (let y = ya; y <= yb; y++) for (let x = xa; x <= xb; x++) { const i = y * cols + x; if (!mat[i]) continue; const dx = ras.x0 + x * cell - p[0], dy = ras.y0 + y * cell - p[1]; if (dx * dx + dy * dy <= r * r) { mat[i] = 0; n++; } }
    return n * cell * cell;
  };
  clearDisk(entry.c, R + entry.r);
  // the cutter's rim, sampled: the arc in material before each step cuts = the engagement angle there
  const NR = 360, rim = Array.from({ length: NR }, (_, k) => [Math.cos(2 * Math.PI * k / NR), Math.sin(2 * Math.PI * k / NR)]);
  const inMat = (q) => { const x = Math.round((q[0] - ras.x0) / cell), y = Math.round((q[1] - ras.y0) / cell); return x >= 0 && y >= 0 && x < cols && y < rows && mat[y * cols + x] === 1; };
  const rr = R - cell;   // just inside the rim, so the cells the step is about to cut are the ones counted
  let angle = 0, sumA = 0, nA = 0, gouge = Infinity, at = null, pos = [entry.c[0] + entry.r, entry.c[1]], inCut = 0, pathLen = 0;
  const step = cell;
  const along = (a, b, pi) => {
    const L = Math.hypot(b[0] - a[0], b[1] - a[1]), n = Math.max(1, Math.ceil(L / step));
    for (let k = 1; k <= n; k++) {
      const p = [a[0] + (b[0] - a[0]) * k / n, a[1] + (b[1] - a[1]) * k / n];
      let m = 0;
      for (const [cx, cy] of rim) if (inMat([p[0] + cx * rr, p[1] + cy * rr])) m++;
      const ang = 360 * m / NR;
      if (ang > angle) { angle = ang; at = { pass: pi, p }; }
      sumA += ang; nA++; if (ang > 0) inCut += L / n; pathLen += L / n;
      clearDisk(p, R);
    }
  };
  for (const [pi, lp] of passes.entries()) {
    const isClosed = closed ? closed[pi] : true;
    if (!links || pi === 0 || links[pi - 1] !== 'up') along(pos, lp[0], pi); else clearDisk(lp[0], R);
    for (let i = 0; i < (isClosed ? lp.length : lp.length - 1); i++) {
      const a = lp[i], b = lp[(i + 1) % lp.length];
      gouge = Math.min(gouge, O.distToLoops(a, loops) - R - leave);
      along(a, b, pi);
    }
    pos = isClosed ? lp[0] : lp[lp.length - 1];
  }
  // the wall pass (exact, at R + leave) the operation ends a level with
  for (const wl of finish || []) { for (const pr of wl) { const n = Math.max(1, Math.ceil(O.lengthOf(pr) / (cell / 2))); for (let k = 0; k <= n; k++) clearDisk(O.pointAt(pr, k / n), R); } }
  let uncut = 0;
  for (let i = 0; i < cols * rows; i++) if (mat[i] && dA[i] * cell <= R - 1.5 * cell) uncut++;
  return { angle, width: R * (1 - Math.cos(Math.min(180, angle) * Math.PI / 180)), meanAngle: nA ? sumA / nA : 0, uncut: uncut * cell * cell, gouge, at, pathLen, inCut };
}