// Helper: Next power of 2
function nextPow2(n) {
  return 1 << (32 - Math.clz32(n - 1));
}

// Basic Cooley-Tukey FFT (real/imag arrays)
function fft1D(re, im) {
  const N = re.length;
  if (N <= 1) return;

  // Bit-reversal permutation
  const rev = new Uint32Array(N);
  let logN = Math.log2(N);
  for (let i = 0; i < N; i++) {
    let x = i,
      y = 0;
    for (let j = 0; j < logN; j++) {
      y <<= 1;
      y |= x & 1;
      x >>= 1;
    }
    rev[i] = y;
  }
  for (let i = 0; i < N; i++) {
    if (i < rev[i]) {
      [re[i], re[rev[i]]] = [re[rev[i]], re[i]];
      [im[i], im[rev[i]]] = [im[rev[i]], im[i]];
    }
  }

  // FFT
  for (let s = 1; s <= logN; s++) {
    const m = 1 << s;
    const m2 = m >> 1;
    const wAngle = (-2 * Math.PI) / m;
    const cosW = Math.cos(wAngle);
    const sinW = Math.sin(wAngle);
    for (let k = 0; k < N; k += m) {
      let wr = 1,
        wi = 0;
      for (let j = 0; j < m2; j++) {
        const tRe = wr * re[k + j + m2] - wi * im[k + j + m2];
        const tIm = wr * im[k + j + m2] + wi * re[k + j + m2];
        const uRe = re[k + j];
        const uIm = im[k + j];
        re[k + j] = uRe + tRe;
        im[k + j] = uIm + tIm;
        re[k + j + m2] = uRe - tRe;
        im[k + j + m2] = uIm - tIm;
        const tempWr = wr;
        wr = wr * cosW - wi * sinW;
        wi = tempWr * sinW + wi * cosW;
      }
    }
  }
}

// Inverse FFT
function ifft1D(re, im) {
  // Conjugate
  for (let i = 0; i < re.length; i++) im[i] = -im[i];
  fft1D(re, im);
  // Normalize and re-conjugate
  const N = re.length;
  for (let i = 0; i < N; i++) {
    re[i] /= N;
    im[i] = -im[i] / N;
  }
}

// Apply ramp filter in frequency domain
export function applyRampFilter(projection) {
  const N = nextPow2(projection.length);
  const re = new Float32Array(N);
  const im = new Float32Array(N);
  for (let i = 0; i < projection.length; i++) {
    re[i] = projection[i];
  }

  fft1D(re, im);

  for (let i = 0; i < N / 2; i++) {
    const freq = i / N;
    re[i] *= freq;
    im[i] *= freq;
    re[N - i - 1] *= freq;
    im[N - i - 1] *= freq;
  }

  ifft1D(re, im);
  return Array.from(re.slice(0, projection.length));
}

export function applyFilter(projection, type = "ramp") {
  const N = nextPow2(projection.length);
  const re = new Float32Array(N);
  const im = new Float32Array(N);
  for (let i = 0; i < projection.length; i++) {
    re[i] = projection[i];
  }

  fft1D(re, im);

  for (let i = 0; i < N / 2; i++) {
    const normFreq = i / N;
    const ramp = normFreq;

    let filter = 1;

    switch (type) {
      case "ramp":
        filter = ramp;
        break;
      case "shepp-logan":
        const sinc =
          normFreq === 0
            ? 1
            : Math.sin(Math.PI * normFreq) / (Math.PI * normFreq);
        filter = ramp * sinc;
        break;
      case "hann":
        filter = ramp * (0.5 + 0.5 * Math.cos(Math.PI * normFreq));
        break;
      case "hamming":
        filter = ramp * (0.54 + 0.46 * Math.cos(Math.PI * normFreq));
        break;
      case "cosine":
        filter = ramp * Math.cos((Math.PI * normFreq) / 2);
        break;
      case "none":
        filter = 1;
        break;
      default:
        console.warn(`Unknown filter type: ${type}. Defaulting to ramp.`);
        filter = ramp;
    }

    re[i] *= filter;
    im[i] *= filter;
    re[N - i - 1] *= filter;
    im[N - i - 1] *= filter;
  }

  ifft1D(re, im);
  return Array.from(re.slice(0, projection.length));
}