Reputation: 75
This question is probably part FFT knowledge and part programming knowledge, but figured I'd post it here to see what you think. I'm trying to implement a ramp filter in JavaScript using Project Nayuki's code and can't quite imitate what I've already done in C++ (FFTW) and Octave/MATLAB. I am zero-padding the initial data array of 672 to 2048 and creating the ramp filter in the spatial domain. Here's images of the data before and after the ramp filter, using Octave's FFT:
And here's the Octave code:
% This script checks my FBP reconstruction code
clear
BaseFolder = "/home/steven/C++/TestSolutio/";
a = textread([BaseFolder "proj.txt"],"%f");
b = textread([BaseFolder "norm.txt"],"%f");
p = zeros(size(a));
for n = 0:499
p((672*n+1):(672*n+672)) = -log(a((672*n+1):(672*n+672)) ./ b);
end
dfan = (2.0*0.0625)/(80.0);
FilterSize = (2*(672-1)) + 1;
Np = 2048;
FilterPadding = (Np-FilterSize-1)/2;
FilterOriginal = zeros(FilterSize, 1);
for f = 1:FilterSize
nf = (-672+1) + f - 1;
if(nf == 0)
FilterOriginal(f) = 1.0 / (8.0 * dfan^2);
else
if(mod(nf,2) == 0) FilterOriginal(f) = 0;
else FilterOriginal(f) = -0.5 / (pi*sin(nf*dfan))^2;
endif
endif
end
RampFilter = zeros(Np, 1);
for f = 1:Np
if(f <= FilterPadding || f > (FilterSize+FilterPadding)) RampFilter(f) = 0;
else RampFilter(f) = FilterOriginal(f-FilterPadding);
endif
end
Filter = abs(fft(RampFilter));
proj_id = 0;
ProjBuffer = zeros(Np,1);
ProjPadding = (Np-672)/2;
for f = 1:Np
if(f <= ProjPadding || f > (672+ProjPadding)) ProjBuffer(f) = 0;
else ProjBuffer(f) = p(672*proj_id+f-ProjPadding);
endif
end
ProjFilter = fft(ProjBuffer);
ProjFilter = ProjFilter .* Filter;
Proj = ifft(ProjFilter);
ProjFinal = Proj((ProjPadding+1):(ProjPadding+672));
plot(1:672, p((672*proj_id+1):(672*proj_id+672)))
axis([1 672 -5 10])
figure
plot(1:Np, Filter)
figure
plot(1:672, ProjFinal)
When I try to do this using JavaScript, it looks as if half the signal is flipped and added to the other half, but I don't know what's really happening:
Here's the JS function:
function filterProj(proj){
// Initialization variables
var padded_size = 2048;
var n_channels = 672;
var d_fan = (2.0*0.0625) / 80.0;
// Create ramp filter
var filter_size = (2*(n_channels-1))+1;
var filter_padding = (padded_size - filter_size - 1)/2;
var ramp_filter = new Array();
var nf;
for(f = 0; f < filter_size; f++){
nf = (-n_channels+1) + f;
if(nf == 0) ramp_filter.push(1.0 / (8.0*Math.pow(d_fan,2.0)));
else {
if(nf % 2 == 0) ramp_filter.push(0.0);
else ramp_filter.push(-0.5 / Math.pow((Math.PI*Math.sin(nf*d_fan)),2.0));
}
}
// Pad filter with zeros & transform
var filter_real = new Array();
var filter_img = new Array();
var filter = new Array();
for(f = 0; f < padded_size; f++){
if(f < filter_padding || f > (filter_size+filter_padding-1)){
filter_real.push(0.0);
}
else {
filter_real.push(ramp_filter[(f-filter_padding)]);
}
filter_img.push(0.0);
}
transform(filter_real, filter_img);
for(f = 0; f < padded_size; f++){
filter_real[f] = Math.abs(filter_real[f]);
}
// For each projection:
// Pad with zeros, take FFT, multiply by filter, take inverse FFT, and remove padding
var proj_padding = (padded_size - n_channels)/2;
for(n = 0; n < 500; n++){
var proj_real = new Array();
var proj_img = new Array();
for(f = 0; f < padded_size; f++){
if(f < proj_padding || f >= (n_channels+proj_padding)){
proj_real.push(0.0);
}
else {
proj_real.push(proj[(n_channels*n + (f-proj_padding))]);
}
proj_img.push(0.0);
}
transform(proj_real, proj_img);
for(f = 0; f < padded_size; f++){
proj_real[f] *= filter_real[f];
}
inverseTransform(proj_real, proj_img);
for(f = 0; f < n_channels; f++){
proj[(n_channels*n+f)] = (d_fan*proj_real[(proj_padding+f)])/padded_size;
}
}
}
Any help/advice would be greatly appreciated!
As an example, here's the ramp filters after FFT, using the same spatial-domain ramp filter as input:
Upvotes: 3
Views: 346
Reputation: 75
Well after some time and investigation I found that my error was not FFT related but a simple complex math error. In C++/Matlab/Octave, the complex data types overload the abs() function to calculate the complex magnitude. However, the Nayuki code splits the input/output data into its real and complex parts, so the complex magnitude needs to be calculated manually. In summary, replacing this line:
filter_real[f] = Math.abs(filter_real[f]);
with this one:
filter_real[f] = Math.sqrt(Math.pow(filter_real[f],2) + Math.pow(filter_img[f],2));
solved all my problems.
Upvotes: 1