/**
 * Robuste FFT-Periodenschätzung für Fliesen:
 * - Welch-Averaging (überlappende Segmente) -> weniger Rauschen/Verdeckungen
 * - Hann-Fenster + Zero-Padding zur nächsten 2er-Potenz (fft.js-Anforderung)
 * - Subpixel-Peak (parabolische Interpolation)
 * - Phase über Kamm-Summe (Comb-Score)
 * - Fugenbreite via FWHM an Peakpositionen der Projektion
 */
import FFT from 'fft.js';

type FftOk = {
  ok: true;
  periodX: number; periodY: number;
  phaseX: number;  phaseY: number;
  groutPx: number;
  vlines: number[]; hlines: number[];
};
type FftErr = { ok:false; reason:string };
type FftOpts = {
  debug?: boolean;
  weights?: Float32Array | null;
};

function clamp(v:number, a:number, b:number){ return Math.max(a, Math.min(b, v)); }
function nextPow2(n:number){ return n<=1 ? 2 : 1<<Math.ceil(Math.log2(n)); }

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

function realFFTmag(inpReal: Float64Array){
  let n = inpReal.length|0;
  // Defensive: ensure FFT size is a power of two and >= 2
  const isPow2 = (x:number)=> x>1 && (x & (x-1)) === 0;
  if (!isPow2(n)) {
    const n2 = nextPow2(Math.max(2, n));
    const buf = new Float64Array(n2);
    buf.set(inpReal.subarray(0, Math.min(n2, n)));
    inpReal = buf;
    n = n2;
  }
  const fft = new FFT(n);
  const out = fft.createComplexArray() as unknown as Float64Array;
  fft.realTransform(out, inpReal);
  fft.completeSpectrum(out);
  const mags = new Float64Array(n/2 + 1);
  for (let k=0;k<=n/2;k++){
    const re = out[2*k] || 0, im = out[2*k+1] || 0;
    mags[k] = Math.hypot(re, im);
  }
  return mags;
}

/** Welch: mittelt Spektren aus überlappenden Segmenten gleicher Länge. */
function welchAveragedMag(proj: Float32Array, 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);                       // für FFT
  const hop = Math.max(1, Math.floor(segLen * (1-overlap)));
  const numStarts = Math.max(1, Math.floor((N - segLen) / hop) + 1);

  const avg = new Float64Array(segLen/2 + 1);
  const buf = new Float64Array(segLen);

  let used = 0;
  for (let s=0; s<numStarts; s++){
    const start = s*hop;
    if (start+segLen > N) break;
    hannFill(buf, proj, start, segLen);
    const mag = realFFTmag(buf);
    for (let k=0;k<mag.length;k++) avg[k]+=mag[k];
    used++;
  }
  if (used===0){
    // Fallback: einmal über gesamtes Signal
    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 };
}

/** Subpixel: Parabel um (k-1,k,k+1) -> Versatz delta in [-0.5,0.5] */
function subpixelPeak(mags: Float64Array, k:number){
  const a = mags[k-1] ?? 0, b = mags[k] ?? 0, c = mags[k+1] ?? 0;
  const denom = (a - 2*b + c);
  if (denom === 0) return 0;
  const delta = 0.5 * (a - c) / denom;
  return clamp(delta, -0.5, 0.5);
}

function combScore(signal: Float32Array, P:number, off:number){
  const N = signal.length;
  const step = Math.max(2, Math.round(P));
  let s=0, cnt=0;
  for (let k=off; k<N; k+=step){
    const i = Math.round(k);
    const i0 = clamp(i-1,0,N-1), i1=i, i2=clamp(i+1,0,N-1);
    s += Math.max(signal[i0], signal[i1], signal[i2]);
    cnt++;
  }
  return cnt ? s/cnt : 0;
}
function bestPhase(signal: Float32Array, P:number){
  const Pi = Math.max(2, Math.round(P));
  let bestOff=0, bestS=-Infinity;
  for (let off=0; off<Pi; off++){
    const s = combScore(signal, P, off);
    if (s>bestS){ bestS=s; bestOff=off; }
  }
  return bestOff;
}

/** FWHM um Peaks herum -> typische Fugenbreite (Median aus Stichproben) */
function fwhmAt(signal: Float32Array, center: number){
  const N = signal.length;
  const ic = clamp(Math.round(center), 0, N-1);
  const peak = signal[ic];
  if (!(peak>0)) return 1;

  const half = peak/2;
  let L=ic, R=ic;
  while (L>0 && signal[L] > half) L--;
  while (R<N-1 && signal[R] > half) R++;
  return Math.max(1, R - L);
}
function groutFromProjection(signal: Float32Array, P: number, phase: number, samples=8){
  const N = signal.length;
  const step = Math.max(2, Math.round(P));
  const widths:number[] = [];
  for (let n=0; n<samples; n++){
    const x = phase + n*step;
    if (x>=N) break;
    widths.push(fwhmAt(signal, x));
  }
  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 findPeriodWithWelch(proj: Float32Array, 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 as const, reason:'no valid k-range', P:NaN, k:NaN };

  let k0 = kMin, 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 as const, reason:'no peak', P:NaN, k:NaN };

  // Subpixel verfeinern
  const delta = subpixelPeak(mags, clamp(k0, kMin+1, kMax-1));
  const kRef = k0 + delta;
  const P = nFFT / kRef;
  return { ok:true as const, P, k: kRef };
}

/** Gradienten-Projektionen (|gx| für vertikale Linien, |gy| für horizontale). */
function gradientProjections(gray: Float32Array, w:number, h:number, weights?: Float32Array | null){
  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 weight = weights ? weights[i] ?? 1 : 1;
      const gx = (gray[i+1]-gray[i-1]) * weight;
      const gy = (gray[i+w]-gray[i-w]) * weight;
      projX[x]+=Math.abs(gx);
      projY[y]+=Math.abs(gy);
    }
  }
  // leichte Glättung (3-Tap)
  const smooth = (a: Float32Array) => {
    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) };
}

export function analyzeFFTPeriod(
  gray: Float32Array,
  w: number,
  h: number,
  opts?: FftOpts
): FftOk|FftErr {
  try{
    if (!Number.isFinite(w) || !Number.isFinite(h) || w<8 || h<8) {
      return { ok:false, reason:'image too small/invalid' };
    }

    const { projX, projY } = gradientProjections(gray, w, h, opts?.weights ?? null);

    const minP=8, 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 = bestPhase(projX, fx.P);
    const phaseY = bestPhase(projY, fy.P);

    // Fugenbreite aus FWHM der Peakform (Median)
    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:number[]=[]; for (let x=phaseX; x<w; x+=stepX) vlines.push(Math.round(x));
    const hlines:number[]=[]; 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
    };
  }catch(e:any){
    return { ok:false, reason: e?.message || String(e) };
  }
}