import { orientedProjectionsWithOccluder } from './occluder_mask';

type RefineRes = {
  ok: boolean; reason?: string;
  periodX: number; periodY: number;
  phaseX: number;  phaseY: number;
  groutPx: number;
  vlines: number[]; hlines: number[];
};

function clamp(v:number, lo:number, hi:number){ return Math.max(lo, Math.min(hi, v)); }
function percentile(arr:number[], p:number): number {
  if (!arr.length) return 0;
  const a = arr.slice().sort((x,y)=>x-y);
  const idx = clamp((p/100)*(a.length-1), 0, a.length-1);
  const lo = Math.floor(idx), hi = Math.ceil(idx);
  if (lo===hi) return a[lo];
  const t = idx-lo; return a[lo]*(1-t)+a[hi]*t;
}

function combScore(signal: Float32Array, P: number, off: number){
  // Kamm über ganze Länge, robust: lokal +-1 mit max (kleiner Subpixel-Effekt)
  const N = signal.length;
  let s = 0;
  const step = Math.max(2, Math.round(P));
  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);
    const v = Math.max(signal[i0], signal[i1], signal[i2]);
    s += v;
  }
  return s / Math.max(1, Math.floor(N/step));
}

function bestPeriodAndPhase(signal: Float32Array, approxP: number, searchFrac=0.2){
  if (!(approxP>2)) return {ok:false, P:NaN, off:NaN, reason:'approx period invalid'};
  const minP = Math.max(2, approxP*(1-searchFrac));
  const maxP = approxP*(1+searchFrac);
  let best = {score:-Infinity, P:approxP, off:0};

  // Schrittweite: fein genug, aber schnell
  const stepP = Math.max(0.25, approxP/80);
  for (let P=minP; P<=maxP; P+=stepP){
    const P_i = Math.round(P);
    for (let off=0; off<P_i; off++){
      const s = combScore(signal, P, off);
      if (s>best.score) best = {score:s, P, off};
    }
  }
  return { ok:true, P:best.P, off:best.off };
}

function localPeak(signal: Float32Array, x: number, win=2){
  // suche Maximum in +-win
  const N = signal.length;
  let bestIdx = clamp(Math.round(x), 0, N-1);
  let best = -Infinity;
  const a = clamp(Math.round(x-win), 0, N-1);
  const b = clamp(Math.round(x+win), 0, N-1);
  for (let i=a; i<=b; i++){
    if (signal[i] > best){ best = signal[i]; bestIdx = i; }
  }
  return bestIdx;
}

function robustLineSet(signal: Float32Array, P:number, off:number){
  // Zähne lokalisieren + RANSAC-ähnliches Re-Fit (Regression über akzeptierte Zähne)
  const N = signal.length;
  const teethIdx:number[] = [];
  const teethScore:number[] = [];
  const step = Math.max(2, Math.round(P));
  for (let k=off; k<N; k+=step){
    const idx = localPeak(signal, k, 2);
    teethIdx.push(idx);
    teethScore.push(signal[idx]);
  }
  if (!teethIdx.length) return { idx:[], P, off };

  // Outlier-Rejection: verwerfe schwache Zähne
  const med = percentile(teethScore, 50);
  const keep:number[] = [];
  for (let i=0;i<teethIdx.length;i++){
    if (teethScore[i] >= 0.4*med) keep.push(i);
  }
  if (keep.length<2) return { idx: teethIdx, P, off };

  // Regressiere Position ~ n*P + off
  const xs:number[] = [], ys:number[] = [];
  for (let j=0;j<keep.length;j++){
    const i = keep[j];
    xs.push(j); ys.push(teethIdx[i]);
  }
  // y = a*x + b
  const n = xs.length;
  let sumX=0,sumY=0,sumXY=0,sumXX=0;
  for (let i=0;i<n;i++){
    sumX+=xs[i]; sumY+=ys[i]; sumXY+=xs[i]*ys[i]; sumXX+=xs[i]*xs[i];
  }
  const a = (n*sumXY - sumX*sumY) / Math.max(1e-6, (n*sumXX - sumX*sumX));
  const b = (sumY - a*sumX) / n;
  const Pnew = Math.max(2, a);
  const offNew = clamp(Math.round(b)%Math.round(Pnew), 0, Math.round(Pnew)-1);

  // Rekonstruiere Linien mit neuem P/off
  const idxNew:number[] = [];
  for (let y=offNew; y<N; y+=Math.round(Pnew)){
    idxNew.push(y);
  }
  return { idx: idxNew, P: Pnew, off: offNew };
}

function estimateGroutWidth1D(signal: Float32Array, lines: number[]): number {
  if (lines.length<2) return NaN;
  const N = signal.length;
  const widths:number[] = [];
  for (const id of lines){
    const i = clamp(Math.round(id), 1, N-2);
    const peak = signal[i];
    const half = peak*0.5;
    let L=i, R=i;
    while (L>0 && signal[L]>half) L--;
    while (R<N-1 && signal[R]>half) R++;
    const w = R-L;
    if (w>0 && w<=16) widths.push(w);
  }
  if (!widths.length) return NaN;
  widths.sort((a,b)=>a-b);
  return widths[Math.floor(widths.length/2)];
}

export function refineLatticeWithOccluder(
  gray: Float32Array, w:number, h:number,
  approx: { periodX:number; periodY:number; phaseX:number; phaseY:number },
  opts?: { searchFrac?: number }
): RefineRes {
  if (!(w>8 && h>8)) return {ok:false, reason:'image too small', periodX:NaN,periodY:NaN,phaseX:NaN,phaseY:NaN,groutPx:NaN, vlines:[], hlines:[]};

  // Oriented projections (Occluder gedämpft)
  const { projX, projY } = orientedProjectionsWithOccluder(gray, w, h, { orientGate: 0.55 });

  // Suche um die Approx.-Periode herum
  const fx = bestPeriodAndPhase(projX, approx.periodX, opts?.searchFrac ?? 0.22);
  const fy = bestPeriodAndPhase(projY, approx.periodY, opts?.searchFrac ?? 0.22);
  if (!fx.ok || !fy.ok){
    return { ok:false, reason:'refine failed', periodX:NaN,periodY:NaN,phaseX:NaN,phaseY:NaN,groutPx:NaN, vlines:[], hlines:[] };
  }

  // RANSAC-ähnliches Re-Fit
  const rx = robustLineSet(projX, fx.P, fx.off);
  const ry = robustLineSet(projY, fy.P, fy.off);

  const groutX = estimateGroutWidth1D(projX, rx.idx);
  const groutY = estimateGroutWidth1D(projY, ry.idx);
  let grout = Number.isFinite(groutX) && Number.isFinite(groutY)
    ? Math.min(groutX as number, groutY as number)
    : (Number.isFinite(groutX) ? (groutX as number) :
       (Number.isFinite(groutY) ? (groutY as number) : NaN));
  if (!Number.isFinite(grout)) grout = Math.max(1, Math.round(0.02*Math.min(rx.P, ry.P)));

  return {
    ok: true,
    periodX: rx.P,
    periodY: ry.P,
    phaseX: rx.off,
    phaseY: ry.off,
    groutPx: grout,
    vlines: rx.idx,
    hlines: ry.idx
  };
}