/**
 * Hough-Fallback für Gittererkennung (zwei orthogonale Linienfamilien).
 * - Kanten: Sobel |gx|,|gy| -> Magnitude
 * - Adaptive Schwelle (Perzentil) + optionales Subsampling
 * - θ-Scan in Fenstern um 0° (vertikal) und 90° (horizontal)
 * - Für bestes θ*: ρ-Histogramm, Glättung, Peak-Detektion -> Perioden/Phase
 * - Fugenbreite via FWHM aus Achs-Projektion (wie im FFT-Analyzer)
 *
 * Keine externen Dependencies.
 */

type Ok = {
  ok: true;
  periodX: number; periodY: number;
  phaseX: number;  phaseY: number;
  groutPx: number;
  thetaV: number;  thetaH: number;
  vlines: number[]; hlines: number[];
};
type Err = { ok:false; reason:string };

const toRad = (deg:number)=>deg*Math.PI/180;
const clamp = (v:number,a:number,b:number)=>Math.max(a,Math.min(b,v));
const median = (a:number[])=>{
  if (a.length===0) return NaN;
  const b=[...a].sort((x,y)=>x-y);
  const i=Math.floor(b.length/2);
  return b.length%2?b[i]:(b[i-1]+b[i])/2;
};

function sobel(gray: Float32Array, w:number, h:number, weights?: Float32Array | null){
  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)], p01=gray[(y-1)*w+x], p02=gray[(y-1)*w+(x+1)];
      const p10=gray[y*w+(x-1)],     p12=gray[y*w+(x+1)];
      const p20=gray[(y+1)*w+(x-1)], p21=gray[(y+1)*w+x], p22=gray[(y+1)*w+(x+1)];
      const weight = weights ? weights[i] ?? 1 : 1;
      const sx = (-p00 -2*p10 -p20 + p02 +2*p12 + p22) * weight;
      const sy = ( p00 +2*p01 +p02 - p20 -2*p21 - p22) * weight;
      gx[i]=sx; gy[i]=sy;
      mag[i]=Math.hypot(sx,sy);
    }
  }
  return {gx,gy,mag};
}

function percentile(arr: Float32Array, q:number){
  const n=arr.length;
  const step = Math.max(1, Math.floor(n/50000)); // subsample bei sehr großen Bildern
  const tmp:number[]=[];
  for (let i=0;i<n;i+=step){ const v=arr[i]; if (Number.isFinite(v)) tmp.push(v); }
  tmp.sort((a,b)=>a-b);
  if (tmp.length===0) return 0;
  const p = clamp(q,0,1)*(tmp.length-1);
  const lo=Math.floor(p), hi=Math.ceil(p);
  if (lo===hi) return tmp[lo];
  const t=p-lo;
  return tmp[lo]*(1-t)+tmp[hi]*t;
}

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

function pickThetaFromWindow(points: Array<[number,number,number]>, w:number, h:number, degStart:number, degEnd:number, degStep=1){
  const maxR = Math.hypot(w,h);
  const thetas:number[]=[];
  for (let d=degStart; d<=degEnd; d+=degStep) thetas.push(toRad(d));
  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)); // 1px bins, implizit
      if (bin>=0 && bin<=Math.ceil(2*maxR)) acc[t]+=wt;
    }
  }
  // leichte Glättung über θ
  const accS = smooth1d(acc, 5);
  let bestI=0, bestV=-Infinity;
  for (let i=0;i<accS.length;i++){ if (accS[i]>bestV){ bestV=accS[i]; bestI=i; } }
  const bestTheta = thetas[bestI] ?? toRad((degStart+degEnd)/2);
  return bestTheta;
}

function rhoHistogram(points: Array<[number,number,number]>, theta:number, w:number, h:number){
  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: Float64Array, minDist=8, minRel=0.1){
  // einfache Nachbarschaftsmaxima
  const maxVal = hist.length? Math.max(...hist as any):0;
  const thr = maxVal*minRel;
  const peaks:number[]=[];
  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);
  }
  // Non-maximum suppression (minDist)
  const keep:number[]=[];
  peaks.sort((a,b)=>hist[b]-hist[a]);
  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:number[]){
  if (peaks.length<2) return NaN;
  const diffs:number[]=[];
  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;
  // Modus-Annäherung per 1px-Binning + Median in der Hauptklasse
  const maxD = Math.max(...diffs), minD = Math.min(...diffs);
  const bins = new Array(maxD+1).fill(0);
  diffs.forEach(d=>bins[d]++);
  let bestD=0, bestC=-1;
  for (let d=minD; d<=maxD; d++){ if (bins[d]>bestC){ bestC=bins[d]; bestD=d; } }
  // nähere Werte um bestD sammeln
  const near = diffs.filter(d=>Math.abs(d-bestD)<=2);
  return median(near.length?near:[bestD]);
}

function phaseFromPeaks(peaks:number[], period:number){
  if (!Number.isFinite(period) || period<=0) return 0;
  const r = peaks.map(p=>((p%period)+period)%period);
  return Math.round(median(r) || 0);
}

function projectionsFromGrad(gx: Float32Array, gy: Float32Array, w:number, h:number){
  const projX = new Float32Array(w); // Summe |gx| -> vertikale Kanten
  const projY = new Float32Array(h); // Summe |gy| -> horizontale Kanten
  for (let y=1;y<h-1;y++){
    let 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]);
    }
  }
  // etwas glätten
  const s = (a: Float32Array)=> {
    const out = new Float32Array(a.length);
    for (let i=0;i<a.length;i++){
      const a0=a[i-1]??a[i], a1=a[i], a2=a[i+1]??a[i];
      out[i]=(a0+2*a1+a2)/4;
    }
    return out;
  };
  return { projX: s(projX), projY: s(projY) };
}

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, period:number, phase:number, samples=8){
  const N=signal.length;
  const step=Math.max(2, Math.round(period));
  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*period));
  return clamp(Math.round(median(widths)),1,Math.round(period/2));
}

type HoughOpts = {
  debug?: boolean;
  weights?: Float32Array | null;
};

export function analyzeHoughGrid(gray: Float32Array, w:number, h:number, opts?:HoughOpts): Ok|Err {
  if (!Number.isFinite(w)||!Number.isFinite(h)||w<16||h<16) return {ok:false,reason:'image too small/invalid'};
  const {gx,gy,mag}=sobel(gray,w,h, opts?.weights ?? null);

  // adaptive Schwelle
  const thr = percentile(mag, 0.85);
  const pts:Array<[number,number,number]>=[];
  const stride = Math.max(1, Math.floor(Math.max(w,h)/800)); // subsample für Speed
  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'};

  // beste Orientierungen in zwei Fenstern
  const thetaV = pickThetaFromWindow(pts, w,h, -15, +15, 1);   // nahe 0°  => vertikal
  const thetaH = pickThetaFromWindow(pts, w,h,  75, 105, 1);   // nahe 90° => horizontal

  // ρ-Histogramme
  const {bins:rhoV, offset:offV} = rhoHistogram(pts, thetaV, w,h);
  const {bins:rhoH, offset:offH} = rhoHistogram(pts, thetaH, w,h);

  // Peaks
  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 (!Number.isFinite(periodX)||periodX<=2 || !Number.isFinite(periodY)||periodY<=2) {
    return {ok:false,reason:'no robust period'};
  }

  const phaseX = phaseFromPeaks(peaksV, periodX);
  const phaseY = phaseFromPeaks(peaksH, periodY);

  // Linienlisten als pixelbasierte "x" bzw. "y" Position (bei y=0 bzw. x=0)
  // Für θ≈0 gilt: x ≈ ρ / cosθ  (bei y=0), ρ = bin - off
  const vlines:number[]=[];
  for (let x=phaseX; x<w; x+=Math.round(periodX)) vlines.push(Math.round(x));
  if (!vlines.length){
    for (let 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);
  }

  const hlines:number[]=[];
  for (let y=phaseY; y<h; y+=Math.round(periodY)) hlines.push(Math.round(y));
  if (!hlines.length){
    for (let 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);
  }

  // Fugenbreite aus Gradienten-Projektionen
  const {projX, projY} = projectionsFromGrad(gx, gy, w,h);
  const groutX = groutFromProjection(projX, periodX, phaseX);
  const groutY = groutFromProjection(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
  };
}