#!/usr/bin/env node
const fs = require('fs');
const path = require('path');
const sharp = require('sharp');
const FFT = require('fft.js');

const DEFAULT_MAX_SIDE = 1200;

function clamp(v, lo, hi) {
  return Math.max(lo, Math.min(hi, v));
}

function median(values) {
  const arr = values.filter(Number.isFinite).slice().sort((a, b) => a - b);
  if (!arr.length) return 0;
  const mid = Math.floor(arr.length / 2);
  return arr.length % 2 ? arr[mid] : 0.5 * (arr[mid - 1] + arr[mid]);
}

async function loadImage(filePath, maxSide = DEFAULT_MAX_SIDE) {
  const input = sharp(filePath).removeAlpha();
  const meta = await input.metadata();
  const srcW = meta.width || 0;
  const srcH = meta.height || 0;
  if (!(srcW > 0 && srcH > 0)) {
    throw new Error(`invalid image dimensions for ${filePath}`);
  }
  let work = input;
  let W = srcW;
  let H = srcH;
  if (Math.max(W, H) > maxSide) {
    const scale = maxSide / Math.max(W, H);
    W = Math.max(1, Math.round(W * scale));
    H = Math.max(1, Math.round(H * scale));
    work = work.resize(W, H, { fit: 'inside' });
  }
  const { data, info } = await work.raw().ensureAlpha().toBuffer({ resolveWithObject: true });
  const rgba = new Uint8Array(data);
  const gray = toGrayFloat32(rgba, info.width, info.height);
  const mag = sobelMag(gray, info.width, info.height);
  return {
    rgba,
    gray,
    mag,
    workWidth: info.width,
    workHeight: info.height,
    srcWidth: srcW,
    srcHeight: srcH,
    scaleX: srcW / info.width,
    scaleY: srcH / info.height,
  };
}

function toGrayFloat32(rgba, w, h) {
  const out = new Float32Array(w * h);
  for (let i = 0, j = 0; i < rgba.length; i += 4, j += 1) {
    out[j] = 0.299 * rgba[i] + 0.587 * rgba[i + 1] + 0.114 * rgba[i + 2];
  }
  return out;
}

function sobelMag(gray, w, h) {
  const out = new Float32Array(w * h);
  const gxk = [-1, 0, 1, -2, 0, 2, -1, 0, 1];
  const gyk = [-1, -2, -1, 0, 0, 0, 1, 2, 1];
  for (let y = 1; y < h - 1; y++) {
    for (let x = 1; x < w - 1; x++) {
      let gx = 0;
      let gy = 0;
      let k = 0;
      for (let dy = -1; dy <= 1; dy++) {
        for (let dx = -1; dx <= 1; dx++, k += 1) {
          const v = gray[(y + dy) * w + (x + dx)];
          gx += gxk[k] * v;
          gy += gyk[k] * v;
        }
      }
      out[y * w + x] = Math.hypot(gx, gy);
    }
  }
  return out;
}

function projectColumns(mag, w, h) {
  const out = new Float32Array(w);
  for (let x = 0; x < w; x++) {
    let acc = 0;
    for (let y = 0; y < h; y++) acc += mag[y * w + x];
    out[x] = acc;
  }
  return out;
}

function projectRows(mag, w, h) {
  const out = new Float32Array(h);
  for (let y = 0; y < h; y++) {
    let acc = 0;
    for (let x = 0; x < w; x++) acc += mag[y * w + x];
    out[y] = acc;
  }
  return out;
}

function autoCorr1D(sig) {
  const n = sig.length;
  const out = new Float32Array(n);
  let mean = 0;
  for (let i = 0; i < n; i++) mean += sig[i];
  mean /= n;
  for (let i = 0; i < n; i++) sig[i] = sig[i] - mean;
  for (let lag = 1; lag < n; lag++) {
    let s = 0;
    for (let i = 0; i + lag < n; i++) s += sig[i] * sig[i + lag];
    out[lag] = s;
  }
  out[0] = 0;
  return out;
}

function bestPeriod(ac, minP, maxP) {
  minP = Math.max(3, minP | 0);
  maxP = Math.max(minP + 1, maxP | 0);
  let best = 0;
  let bestVal = -Infinity;
  for (let p = minP; p <= maxP && p < ac.length; p++) {
    const v = ac[p];
    if (v > bestVal) {
      bestVal = v;
      best = p;
    }
  }
  return best;
}

function phaseForPeriod(proj, period) {
  if (period <= 0) return 0;
  let bestOff = 0;
  let bestScore = -Infinity;
  for (let off = 0; off < period; off++) {
    let score = 0;
    for (let i = off; i < proj.length; i += period) score += proj[i];
    if (score > bestScore) {
      bestScore = score;
      bestOff = off;
    }
  }
  return bestOff;
}

function analyzeProjectionPipeline(prep) {
  const { mag, workWidth: w, workHeight: h, scaleX, scaleY } = prep;
  const projX = projectColumns(mag, w, h);
  const projY = projectRows(mag, w, h);
  const acX = autoCorr1D(projX.slice());
  const acY = autoCorr1D(projY.slice());
  const perX = bestPeriod(acX, Math.max(12, Math.floor(w / 100)), Math.max(24, Math.floor(w / 3)));
  const perY = bestPeriod(acY, Math.max(12, Math.floor(h / 100)), Math.max(24, Math.floor(h / 3)));
  const phaseX = phaseForPeriod(projX, perX);
  const phaseY = phaseForPeriod(projY, perY);
  const grout = Math.max(1, Math.round(0.025 * ((perX + perY) / 2)));
  const vlines = [];
  for (let x = phaseX; x < w; x += perX || w) vlines.push(Math.round(x));
  const hlines = [];
  for (let y = phaseY; y < h; y += perY || h) hlines.push(Math.round(y));
  return {
    tile_w_px: Math.round(perX * scaleX),
    tile_h_px: Math.round(perY * scaleY),
    grout_px: Math.round(grout * Math.sqrt(scaleX * scaleY)),
    grid: {
      vlines: vlines.map((x) => Math.round(x * scaleX)),
      hlines: hlines.map((y) => Math.round(y * scaleY)),
    },
    meta: {
      perX,
      perY,
      phaseX,
      phaseY,
    },
  };
}

function gradientFromGray(gray, w, h) {
  const gx = new Float32Array(w * h);
  const gy = new Float32Array(w * h);
  for (let y = 1; y < h - 1; y++) {
    for (let x = 1; x < w - 1; x++) {
      const i = y * w + x;
      const p00 = gray[(y - 1) * w + (x - 1)];
      const p01 = gray[(y - 1) * w + x];
      const p02 = gray[(y - 1) * w + (x + 1)];
      const p10 = gray[y * w + (x - 1)];
      const p12 = gray[y * w + (x + 1)];
      const p20 = gray[(y + 1) * w + (x - 1)];
      const p21 = gray[(y + 1) * w + x];
      const p22 = gray[(y + 1) * w + (x + 1)];
      gx[i] = -p00 - 2 * p10 - p20 + p02 + 2 * p12 + p22;
      gy[i] = p00 + 2 * p01 + p02 - p20 - 2 * p21 - p22;
    }
  }
  return { gx, gy };
}

function confidenceByLines(vlines, hlines, gx, gy, w, h) {
  const col = new Float64Array(w);
  const row = new Float64Array(h);
  for (let y = 1; y < h - 1; y++) {
    let base = y * w;
    for (let x = 1; x < w - 1; x++) {
      const i = base + x;
      col[x] += Math.abs(gx[i]);
      row[y] += Math.abs(gy[i]);
    }
  }
  const medCol = median(Array.from(col));
  const medRow = median(Array.from(row));
  function pick(idxs, arr) {
    if (!idxs || !idxs.length) return 0;
    const samples = [];
    for (const k of idxs) {
      const x = clamp(Math.round(k), 0, arr.length - 1);
      const a0 = arr[x - 1] ?? arr[x];
      const a1 = arr[x];
      const a2 = arr[x + 1] ?? arr[x];
      samples.push((a0 + 2 * a1 + a2) / 4);
    }
    return median(samples);
  }
  const peakCol = pick(vlines, col);
  const peakRow = pick(hlines, row);
  const sV = peakCol > 0 && medCol > 0 ? peakCol / medCol : 0;
  const sH = peakRow > 0 && medRow > 0 ? peakRow / medRow : 0;
  const map = (s) => clamp((s - 1) / 3, 0, 1);
  return clamp(0.5 * (map(sV) + map(sH)), 0, 1);
}

function hannFill(dst, src, start, len) {
  for (let i = 0; i < len; i++) {
    const s = src[start + i] ?? 0;
    const w = 0.5 * (1 - Math.cos((2 * Math.PI * i) / Math.max(1, len - 1)));
    dst[i] = s * w;
  }
  for (let i = len; i < dst.length; i++) dst[i] = 0;
}

function nextPow2(n) {
  return n <= 1 ? 2 : 1 << Math.ceil(Math.log2(n));
}

function realFFTmag(inp) {
  let n = inp.length | 0;
  const isPow2 = (x) => x > 1 && (x & (x - 1)) === 0;
  if (!isPow2(n)) {
    const n2 = nextPow2(Math.max(2, n));
    const buf = new Float64Array(n2);
    buf.set(inp.subarray(0, Math.min(n2, n)));
    inp = buf;
    n = n2;
  }
  const fft = new FFT(n);
  const out = fft.createComplexArray();
  fft.realTransform(out, inp);
  fft.completeSpectrum(out);
  const mags = new Float64Array(n / 2 + 1);
  for (let k = 0; k <= n / 2; k++) {
    const re = out[2 * k] || 0;
    const im = out[2 * k + 1] || 0;
    mags[k] = Math.hypot(re, im);
  }
  return mags;
}

function welchAveragedMag(proj, targetSegs = 4, overlap = 0.5) {
  const N = proj.length;
  const segLenRaw = Math.max(16, Math.floor(N / Math.max(1, targetSegs)));
  const segLen = nextPow2(segLenRaw);
  const hop = Math.max(1, Math.floor(segLen * (1 - overlap)));
  const avg = new Float64Array(segLen / 2 + 1);
  const buf = new Float64Array(segLen);
  let used = 0;
  for (let s = 0; s * hop + segLen <= N; s++) {
    const start = s * hop;
    hannFill(buf, proj, start, segLen);
    const mag = realFFTmag(buf);
    for (let k = 0; k < mag.length; k++) avg[k] += mag[k];
    used += 1;
  }
  if (!used) {
    const nFFT = nextPow2(N);
    const b2 = new Float64Array(nFFT);
    hannFill(b2, proj, 0, N);
    return { mags: realFFTmag(b2), nFFT };
  }
  for (let k = 0; k < avg.length; k++) avg[k] /= used;
  return { mags: avg, nFFT: segLen };
}

function subpixelPeak(mags, k) {
  const a = mags[k - 1] ?? 0;
  const b = mags[k] ?? 0;
  const c = mags[k + 1] ?? 0;
  const denom = a - 2 * b + c;
  if (denom === 0) return 0;
  return clamp(0.5 * (a - c) / denom, -0.5, 0.5);
}

function findPeriodWithWelch(proj, minP = 8, maxP = 1024) {
  const { mags, nFFT } = welchAveragedMag(proj, 4, 0.5);
  const kMin = Math.max(1, Math.floor(nFFT / (maxP || 1024)));
  const kMax = Math.min(Math.floor(nFFT / 2) - 1, Math.ceil(nFFT / Math.max(minP, 2)));
  if (!(kMax > kMin)) {
    return { ok: false, reason: 'no valid k-range', P: NaN };
  }
  let k0 = kMin;
  let v0 = -Infinity;
  for (let k = kMin; k <= kMax; k++) {
    const v = mags[k];
    if (v > v0) {
      v0 = v;
      k0 = k;
    }
  }
  if (!(v0 > 0)) return { ok: false, reason: 'no peak', P: NaN };
  const delta = subpixelPeak(mags, clamp(k0, kMin + 1, kMax - 1));
  const kRef = k0 + delta;
  return { ok: true, P: nFFT / kRef };
}

function gradientProjections(gray, w, h) {
  const projX = new Float32Array(w);
  const projY = new Float32Array(h);
  for (let y = 1; y < h - 1; y++) {
    const o = y * w;
    for (let x = 1; x < w - 1; x++) {
      const i = o + x;
      const gx = gray[i + 1] - gray[i - 1];
      const gy = gray[i + w] - gray[i - w];
      projX[x] += Math.abs(gx);
      projY[y] += Math.abs(gy);
    }
  }
  function smooth(a) {
    const N = a.length;
    if (N < 3) return a.slice();
    const out = new Float32Array(N);
    out[0] = (a[0] * 2 + a[1]) / 3;
    for (let i = 1; i < N - 1; i++) out[i] = (a[i - 1] + a[i] * 2 + a[i + 1]) / 4;
    out[N - 1] = (a[N - 2] + a[N - 1] * 2) / 3;
    return out;
  }
  return { projX: smooth(projX), projY: smooth(projY) };
}

function bestPhaseComb(signal, P) {
  const Pi = Math.max(2, Math.round(P));
  let bestOff = 0;
  let bestS = -Infinity;
  for (let off = 0; off < Pi; off++) {
    let s = 0;
    for (let k = off; k < signal.length; k += Pi) {
      const i = Math.round(k);
      const i0 = clamp(i - 1, 0, signal.length - 1);
      const i2 = clamp(i + 1, 0, signal.length - 1);
      s += Math.max(signal[i0], signal[i], signal[i2]);
    }
    if (s > bestS) {
      bestS = s;
      bestOff = off;
    }
  }
  return bestOff;
}

function groutFromProjection(signal, P, phase, samples = 8) {
  const step = Math.max(2, Math.round(P));
  const widths = [];
  for (let n = 0; n < samples; n++) {
    const x = phase + n * step;
    if (x >= signal.length) break;
    const idx = clamp(Math.round(x), 0, signal.length - 1);
    const peak = signal[idx];
    if (!(peak > 0)) continue;
    const half = peak / 2;
    let L = idx;
    let R = idx;
    while (L > 0 && signal[L] > half) L -= 1;
    while (R < signal.length - 1 && signal[R] > half) R += 1;
    widths.push(R - L);
  }
  if (!widths.length) return Math.max(1, Math.round(0.04 * P));
  widths.sort((a, b) => a - b);
  return clamp(Math.round(widths[Math.floor(widths.length / 2)]), 1, Math.round(P / 2));
}

function analyzeFFTPeriod(gray, w, h) {
  if (!(w > 8 && h > 8)) return { ok: false, reason: 'image too small' };
  const { projX, projY } = gradientProjections(gray, w, h);
  const minP = 8;
  const maxP = Math.max(32, Math.floor(Math.max(w, h) / 2));
  const fx = findPeriodWithWelch(projX, minP, maxP);
  const fy = findPeriodWithWelch(projY, minP, maxP);
  if (!fx.ok || !fy.ok) {
    return { ok: false, reason: `welch fail x:${fx.ok ? 'ok' : fx.reason} y:${fy.ok ? 'ok' : fy.reason}` };
  }
  const phaseX = bestPhaseComb(projX, fx.P);
  const phaseY = bestPhaseComb(projY, fy.P);
  const groutX = groutFromProjection(projX, fx.P, phaseX);
  const groutY = groutFromProjection(projY, fy.P, phaseY);
  const groutPx = Math.max(1, Math.round((groutX + groutY) / 2));
  const stepX = Math.max(2, Math.round(fx.P));
  const stepY = Math.max(2, Math.round(fy.P));
  const vlines = [];
  for (let x = phaseX; x < w; x += stepX) vlines.push(Math.round(x));
  const hlines = [];
  for (let y = phaseY; y < h; y += stepY) hlines.push(Math.round(y));
  return {
    ok: true,
    periodX: fx.P,
    periodY: fy.P,
    phaseX,
    phaseY,
    groutPx,
    vlines,
    hlines,
  };
}

function sobelFull(gray, w, h) {
  const gx = new Float32Array(w * h);
  const gy = new Float32Array(w * h);
  const mag = new Float32Array(w * h);
  for (let y = 1; y < h - 1; y++) {
    for (let x = 1; x < w - 1; x++) {
      const i = y * w + x;
      const p00 = gray[(y - 1) * w + (x - 1)];
      const p01 = gray[(y - 1) * w + x];
      const p02 = gray[(y - 1) * w + (x + 1)];
      const p10 = gray[y * w + (x - 1)];
      const p12 = gray[y * w + (x + 1)];
      const p20 = gray[(y + 1) * w + (x - 1)];
      const p21 = gray[(y + 1) * w + x];
      const p22 = gray[(y + 1) * w + (x + 1)];
      const sx = -p00 - 2 * p10 - p20 + p02 + 2 * p12 + p22;
      const sy = p00 + 2 * p01 + p02 - p20 - 2 * p21 - p22;
      gx[i] = sx;
      gy[i] = sy;
      mag[i] = Math.hypot(sx, sy);
    }
  }
  return { gx, gy, mag };
}

function percentile(array, q) {
  const vals = [];
  const step = Math.max(1, Math.floor(array.length / 50000));
  for (let i = 0; i < array.length; i += step) {
    const v = array[i];
    if (Number.isFinite(v)) vals.push(v);
  }
  if (!vals.length) return 0;
  vals.sort((a, b) => a - b);
  const p = clamp(q, 0, 1) * (vals.length - 1);
  const lo = Math.floor(p);
  const hi = Math.ceil(p);
  if (lo === hi) return vals[lo];
  const t = p - lo;
  return vals[lo] * (1 - t) + vals[hi] * t;
}

function smooth1d(arr, win = 5) {
  if (win <= 1) return Float64Array.from(arr);
  const n = arr.length;
  const out = new Float64Array(n);
  const r = Math.floor(win / 2);
  for (let i = 0; i < n; i++) {
    let s = 0;
    let c = 0;
    for (let k = -r; k <= r; k++) {
      const j = i + k;
      if (j < 0 || j >= n) continue;
      s += arr[j];
      c += 1;
    }
    out[i] = c ? s / c : arr[i];
  }
  return out;
}

function pickTheta(points, w, h, degStart, degEnd, degStep) {
  const maxR = Math.hypot(w, h);
  const thetas = [];
  for (let d = degStart; d <= degEnd; d += degStep) thetas.push((d * Math.PI) / 180);
  const acc = new Float64Array(thetas.length);
  for (const [x, y, wt] of points) {
    for (let t = 0; t < thetas.length; t++) {
      const th = thetas[t];
      const rho = x * Math.cos(th) + y * Math.sin(th);
      const bin = Math.round(rho + maxR);
      if (bin >= 0 && bin <= Math.ceil(2 * maxR)) acc[t] += wt;
    }
  }
  const accS = smooth1d(acc, 5);
  let bestI = 0;
  let bestV = -Infinity;
  for (let i = 0; i < accS.length; i++) {
    if (accS[i] > bestV) {
      bestV = accS[i];
      bestI = i;
    }
  }
  return thetas[bestI] ?? ((degStart + degEnd) * Math.PI) / 360;
}

function rhoHistogram(points, theta, w, h) {
  const maxR = Math.hypot(w, h);
  const bins = new Float64Array(Math.ceil(2 * maxR) + 1);
  for (const [x, y, wt] of points) {
    const rho = x * Math.cos(theta) + y * Math.sin(theta);
    const b = Math.round(rho + maxR);
    if (b >= 0 && b < bins.length) bins[b] += wt;
  }
  return { bins: smooth1d(bins, 7), offset: maxR };
}

function findPeaks(hist, minDist = 8, minRel = 0.1) {
  const maxVal = hist.length ? Math.max(...hist) : 0;
  const thr = maxVal * minRel;
  const peaks = [];
  for (let i = 1; i < hist.length - 1; i++) {
    if (hist[i] > thr && hist[i] > hist[i - 1] && hist[i] > hist[i + 1]) peaks.push(i);
  }
  peaks.sort((a, b) => hist[b] - hist[a]);
  const keep = [];
  for (const p of peaks) {
    if (keep.every((k) => Math.abs(k - p) >= minDist)) keep.push(p);
  }
  keep.sort((a, b) => a - b);
  return keep;
}

function robustPeriodFromPeaks(peaks) {
  if (peaks.length < 2) return NaN;
  const diffs = [];
  for (let i = 1; i < peaks.length; i++) {
    const d = peaks[i] - peaks[i - 1];
    if (d > 2) diffs.push(d);
  }
  if (!diffs.length) return NaN;
  const bins = new Map();
  const maxD = Math.max(...diffs);
  const minD = Math.min(...diffs);
  for (const d of diffs) bins.set(d, (bins.get(d) || 0) + 1);
  let bestD = minD;
  let bestC = -1;
  for (let d = minD; d <= maxD; d++) {
    const c = bins.get(d) || 0;
    if (c > bestC) {
      bestC = c;
      bestD = d;
    }
  }
  const near = diffs.filter((d) => Math.abs(d - bestD) <= 2);
  return median(near.length ? near : [bestD]);
}

function phaseFromPeaks(peaks, period) {
  if (!(period > 0)) return 0;
  const r = peaks.map((p) => ((p % period) + period) % period);
  return Math.round(median(r) || 0);
}

function projectionsFromGrad(gx, gy, w, h) {
  const projX = new Float32Array(w);
  const projY = new Float32Array(h);
  for (let y = 1; y < h - 1; y++) {
    const row = y * w;
    for (let x = 1; x < w - 1; x++) {
      const i = row + x;
      projX[x] += Math.abs(gx[i]);
      projY[y] += Math.abs(gy[i]);
    }
  }
  const smooth = (a) => {
    const out = new Float32Array(a.length);
    for (let i = 0; i < a.length; i++) {
      const a0 = a[i - 1] ?? a[i];
      const a1 = a[i];
      const a2 = a[i + 1] ?? a[i];
      out[i] = (a0 + 2 * a1 + a2) / 4;
    }
    return out;
  };
  return { projX: smooth(projX), projY: smooth(projY) };
}

function groutFromProjectionSmooth(signal, period, phase, samples = 8) {
  const step = Math.max(2, Math.round(period));
  const widths = [];
  for (let n = 0; n < samples; n++) {
    const x = phase + n * step;
    if (x >= signal.length) break;
    const idx = clamp(Math.round(x), 0, signal.length - 1);
    const peak = signal[idx];
    if (!(peak > 0)) continue;
    const half = peak / 2;
    let L = idx;
    let R = idx;
    while (L > 0 && signal[L] > half) L -= 1;
    while (R < signal.length - 1 && signal[R] > half) R += 1;
    widths.push(R - L);
  }
  if (!widths.length) return Math.max(1, Math.round(0.04 * period));
  widths.sort((a, b) => a - b);
  return clamp(Math.round(widths[Math.floor(widths.length / 2)]), 1, Math.round(period / 2));
}

function analyzeHoughGrid(gray, w, h) {
  if (!(w > 16 && h > 16)) return { ok: false, reason: 'image too small/invalid' };
  const { gx, gy, mag } = sobelFull(gray, w, h);
  const thr = percentile(mag, 0.85);
  const pts = [];
  const stride = Math.max(1, Math.floor(Math.max(w, h) / 800));
  for (let y = 1; y < h - 1; y += stride) {
    for (let x = 1; x < w - 1; x += stride) {
      const i = y * w + x;
      const m = mag[i];
      if (m > thr) pts.push([x, y, m]);
    }
  }
  if (pts.length < 200) return { ok: false, reason: 'too few edges' };
  const thetaV = pickTheta(pts, w, h, -15, 15, 1);
  const thetaH = pickTheta(pts, w, h, 75, 105, 1);
  const { bins: rhoV, offset: offV } = rhoHistogram(pts, thetaV, w, h);
  const { bins: rhoH, offset: offH } = rhoHistogram(pts, thetaH, w, h);
  const peaksV = findPeaks(rhoV, 8, 0.15);
  const peaksH = findPeaks(rhoH, 8, 0.15);
  if (peaksV.length < 2 || peaksH.length < 2) return { ok: false, reason: 'no grid families' };
  const periodX = robustPeriodFromPeaks(peaksV);
  const periodY = robustPeriodFromPeaks(peaksH);
  if (!(periodX > 2) || !(periodY > 2)) return { ok: false, reason: 'no robust period' };
  const phaseX = phaseFromPeaks(peaksV, periodX);
  const phaseY = phaseFromPeaks(peaksH, periodY);
  const vlines = [];
  for (let x = phaseX; x < w; x += Math.round(periodX)) vlines.push(Math.round(x));
  const hlines = [];
  for (let y = phaseY; y < h; y += Math.round(periodY)) hlines.push(Math.round(y));
  if (!vlines.length) {
    for (const p of peaksV) {
      const x = Math.round(((p - offV)) / Math.cos(thetaV));
      if (x >= 0 && x < w) vlines.push(x);
    }
    vlines.sort((a, b) => a - b);
  }
  if (!hlines.length) {
    for (const p of peaksH) {
      const y = Math.round(((p - offH)) / Math.sin(thetaH));
      if (y >= 0 && y < h) hlines.push(y);
    }
    hlines.sort((a, b) => a - b);
  }
  const { projX, projY } = projectionsFromGrad(gx, gy, w, h);
  const groutX = groutFromProjectionSmooth(projX, periodX, phaseX);
  const groutY = groutFromProjectionSmooth(projY, periodY, phaseY);
  const groutPx = Math.max(1, Math.round((groutX + groutY) / 2));
  return {
    ok: true,
    periodX,
    periodY,
    phaseX,
    phaseY,
    groutPx,
    thetaV,
    thetaH,
    vlines,
    hlines,
  };
}

function runPeriodicAnalysis(gray, w, h) {
  const grads = gradientFromGray(gray, w, h);
  const fft = analyzeFFTPeriod(gray, w, h);
  let best = null;
  if (fft.ok) {
    const grout = fft.groutPx || 1;
    const tile_w_px = Math.max(1, Math.round(fft.periodX - grout));
    const tile_h_px = Math.max(1, Math.round(fft.periodY - grout));
    const conf = confidenceByLines(fft.vlines, fft.hlines, grads.gx, grads.gy, w, h);
    best = {
      ok: true,
      tile_w_px,
      tile_h_px,
      grout_px: Math.max(1, Math.round(grout)),
      vlines: fft.vlines,
      hlines: fft.hlines,
      source: 'cv@fft-welch',
      periodX: fft.periodX,
      periodY: fft.periodY,
      phaseX: fft.phaseX,
      phaseY: fft.phaseY,
      confidence: conf,
      thetaV: 0,
      thetaH: Math.PI / 2,
    };
  }
  let hbest = null;
  if (!best || best.confidence < 0.35) {
    const hg = analyzeHoughGrid(gray, w, h);
    if (hg.ok) {
      const grout = Math.max(1, Math.round(hg.groutPx));
      const tile_w_px = Math.max(1, Math.round(hg.periodX - grout));
      const tile_h_px = Math.max(1, Math.round(hg.periodY - grout));
      const conf = confidenceByLines(hg.vlines, hg.hlines, grads.gx, grads.gy, w, h);
      hbest = {
        ok: true,
        tile_w_px,
        tile_h_px,
        grout_px: grout,
        vlines: hg.vlines,
        hlines: hg.hlines,
        source: 'cv@hough-grid',
        periodX: hg.periodX,
        periodY: hg.periodY,
        phaseX: hg.phaseX,
        phaseY: hg.phaseY,
        confidence: conf,
        thetaV: hg.thetaV,
        thetaH: hg.thetaH,
      };
    }
  }
  if (best && hbest) return best.confidence >= hbest.confidence ? best : hbest;
  if (best) return best;
  if (hbest) return hbest;
  return { ok: false, error: 'fft-fail; hough-fail' };
}

function refineGrid1D(lines, extent) {
  const xs = [...new Set(lines.filter((v) => Number.isFinite(v) && v >= 0 && v <= extent))].sort((a, b) => a - b);
  if (xs.length < 2) {
    return { lines: xs, T: 0, a: xs[0] || 0, conf: 0 };
  }
  const diffs = [];
  for (let i = 0; i < xs.length; i++) {
    for (let j = i + 1; j < Math.min(xs.length, i + 10); j++) diffs.push(xs[j] - xs[i]);
  }
  const fallbackT = xs.length > 1 ? xs[1] - xs[0] : 0;
  let T = histPeak(diffs, 1, 6, Math.max(12, Math.floor(extent / 2)));
  if (!(T > 0)) T = median(diffs) || fallbackT;
  const a = phaseFromModulo(xs, T);
  const grid = generateGrid(a, T, extent);
  const tol = Math.max(2, Math.round(0.1 * T));
  let inliers = 0;
  for (const x of xs) {
    if (grid.some((g) => Math.abs(g - x) <= tol)) inliers += 1;
  }
  return { lines: grid, T, a, conf: xs.length ? inliers / xs.length : 0 };
}

function histPeak(values, bin = 1, minT = 6, maxT = 1024) {
  const hist = new Map();
  for (const v of values) {
    if (!Number.isFinite(v) || v < minT || v > maxT) continue;
    const bucket = Math.round(v / bin) * bin;
    hist.set(bucket, (hist.get(bucket) || 0) + 1);
  }
  let best = minT;
  let count = -1;
  for (const [k, v] of hist) {
    if (v > count) {
      best = k;
      count = v;
    }
  }
  return best;
}

function phaseFromModulo(xs, T) {
  if (T <= 0) return 0;
  const remainders = xs.map((x) => ((x % T) + T) % T);
  return median(remainders);
}

function generateGrid(a, T, extent) {
  const result = [];
  if (T <= 0 || extent <= 0) return result;
  let k = Math.floor((0 - a) / T) - 1;
  for (let i = 0; i < 10000; i++) {
    const x = a + k * T;
    if (x > extent) break;
    if (x >= 0) result.push(Math.round(x));
    k += 1;
    if (result.length > 1) {
      const last = result[result.length - 1];
      const prev = result[result.length - 2];
      if (last - prev > 2 * extent) break;
    }
  }
  return Array.from(new Set(result)).sort((x, y) => x - y);
}

function refineGridRansac(w, h, vlines, hlines) {
  const vertical = refineGrid1D(vlines, w);
  const horizontal = refineGrid1D(hlines, h);
  return {
    vlines: vertical.lines,
    hlines: horizontal.lines,
    meta: {
      v: { T: vertical.T, a: vertical.a, conf: vertical.conf },
      h: { T: horizontal.T, a: horizontal.a, conf: horizontal.conf },
    },
  };
}

async function benchmarkImage(filePath, opts) {
  const prep = await loadImage(filePath, opts.maxSide);
  const current = analyzeProjectionPipeline(prep);
  const periodic = runPeriodicAnalysis(prep.gray, prep.workWidth, prep.workHeight);
  let periodicScaled = null;
  if (periodic.ok) {
    periodicScaled = {
      tile_w_px: Math.round(periodic.tile_w_px * prep.scaleX),
      tile_h_px: Math.round(periodic.tile_h_px * prep.scaleY),
      grout_px: Math.round(periodic.grout_px * Math.sqrt(prep.scaleX * prep.scaleY)),
      periodX: periodic.periodX * prep.scaleX,
      periodY: periodic.periodY * prep.scaleY,
      confidence: periodic.confidence,
      source: periodic.source,
      grid: {
        vlines: periodic.vlines.map((x) => Math.round(x * prep.scaleX)),
        hlines: periodic.hlines.map((y) => Math.round(y * prep.scaleY)),
      },
    };
  }
  let refinedCurrent = current.grid;
  let refinedPeriodic = periodicScaled ? periodicScaled.grid : null;
  if (opts.useRansac) {
    const rc = refineGridRansac(prep.srcWidth, prep.srcHeight, current.grid.vlines, current.grid.hlines);
    refinedCurrent = { vlines: rc.vlines, hlines: rc.hlines, meta: rc.meta };
    if (periodicScaled) {
      const rp = refineGridRansac(prep.srcWidth, prep.srcHeight, periodicScaled.grid.vlines, periodicScaled.grid.hlines);
      refinedPeriodic = { vlines: rp.vlines, hlines: rp.hlines, meta: rp.meta };
    }
  }
  return {
    image: path.basename(filePath),
    path: filePath,
    size: { width: prep.srcWidth, height: prep.srcHeight },
    current: {
      tile_w_px: current.tile_w_px,
      tile_h_px: current.tile_h_px,
      grout_px: current.grout_px,
      periodX: current.meta.perX,
      periodY: current.meta.perY,
      grid: refinedCurrent,
    },
    periodic: periodicScaled
      ? {
          ...periodicScaled,
          grid: refinedPeriodic,
        }
      : { ok: false },
  };
}

function parseArgs(argv) {
  const args = { folder: '../testimg', maxSide: DEFAULT_MAX_SIDE, useRansac: true };
  for (let i = 2; i < argv.length; i++) {
    const arg = argv[i];
    if (arg === '--folder' && argv[i + 1]) {
      args.folder = argv[++i];
    } else if (arg === '--max-side' && argv[i + 1]) {
      args.maxSide = Number(argv[++i]) || DEFAULT_MAX_SIDE;
    } else if (arg === '--no-ransac') {
      args.useRansac = false;
    } else if (arg === '--help' || arg === '-h') {
      args.help = true;
    }
  }
  return args;
}

async function main() {
  const args = parseArgs(process.argv);
  if (args.help) {
    console.log('Usage: node scripts/benchmark-analyzers.js [--folder path] [--max-side N] [--no-ransac]');
    process.exit(0);
  }
  const folder = path.resolve(process.cwd(), args.folder);
  const files = fs.readdirSync(folder)
    .filter((f) => /\.(png|jpg|jpeg)$/i.test(f))
    .map((f) => path.join(folder, f));
  if (!files.length) {
    console.error(`No images found in ${folder}`);
    process.exit(1);
  }
  const results = [];
  for (const file of files) {
    try {
      const res = await benchmarkImage(file, args);
      results.push(res);
    } catch (err) {
      console.error(`Failed on ${file}:`, err.message);
    }
  }
  console.log(JSON.stringify({ params: args, results }, null, 2));
}

main().catch((err) => {
  console.error(err);
  process.exit(1);
});