147 lines
3.5 KiB
JavaScript
147 lines
3.5 KiB
JavaScript
// 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));
|
|
}
|