How can I obtain filtered back projection using 1D Fourier Transform

I have an assignment which consist of two part. First part is forward projection of an image and the second part is filtered back projection from the projection data. I managed the first part. In the second part, I can obtain the reconstructed image without filter; however, I cannot apply the filter. The built-in functions such as iradon is forbidden. I am storing the projection data in a matrix MAT, which is MXN. M = length(angles); N = length(beam_number); I have tried the following:

% Triangular Filter

triang_filter = zeros(1,length(tValues));

for ind8 = 1:1:((length(tValues)+1)/2)   %positive slope part of the triangular wave
    
    triang_filter(ind8) = 2 * ((ind8 - 1) / (length(tValues) - 1));
    
end

for ind9 = ((length(tValues)+1)/2):1:(length(tValues))  %negative slope part of the triangular wave
        
    triang_filter(ind9) = 2 * ((length(tValues) - ind9) / (length(tValues) - 1));
end

ft_p_theta_of_t = fft(p_theta_of_t');
ft_p_theta_of_t = ft_p_theta_of_t';

filtered_p_theta_of_t = [];

for ind11 = 1:1:length(thetaValues)     %multiply the fft of projection with filter
    
    filtered_p_theta_of_t(ind11,:) = triang_filter .* ft_p_theta_of_t(ind11,:);

end

inv_filtered_p_theta_of_t = ifft(filtered_p_theta_of_t');
inv_filtered_p_theta_of_t = inv_filtered_p_theta_of_t';

Then, I am using back projection algorithm; however, unfiltered image has less blurring than filtered image.

Thank you in advance :)

Answer 1

[Ross Hanna]

For future reference for people seeing this, two suggestions.

For the triangle filter, either use the triang() function;

https://uk.mathworks.com/help/signal/ref/triang.html

or in one line

filt = (-abs(linspace(-1,1,pixNum))+1);     % Creates a triangle (Ram-Lak) filter

% tmp = linspace(-1,1,pixNum)       % creates a single line from -1 to +1 intercepting x-axis at (0,0) with length "pixNum"
% tmp = abs(tmp)                    % finds the absolute (all positive) of the previous tmp variable (triangle, with point at bottom)
% tmp = -tmp                        % translates filter around y = 0, (triangle with point at top, but in negative domain)
% filt = tmp + 1                    % shifts filter up by 1

this creates the Ram-Lak (triangle filter) in one line as opposed to two loops.

For applying the filter, I'm reconstructing using *.tif images for the x-ray projections, hence the variable name "tf".

*.tif images are imported to matlab, converted to double() and stored in "A" (300 x 370 Matrix)

pixNum = 300;       % number of pixels in reconstruction
Nrot = 370;         % number of rotation steps (360 degrees rotation in 370 steps)

filt = (-abs(linspace(-1,1,pixNum))+1);

for i = 1:Nrot
    
    tf = A(:,rot);      % pick the 1D values to backproject
    tf = fft(tf);       % move into frequency domain
    tf = tf .* filt;    % multiply frequency domain values by the Ram-Lak filter. (multiplication in frequency domain is the same as convolution in spatial domain)
    tf = ifft(tf);      % bring back into spatial domain
    tf = real(tf);      % pick the real part of the signal (this could be debated, something I'm currently investigating)

    Output = backProject(tf);   %Not an actual function, just an example that you then back project the filtered data for each rotation step

end

Hopefully this helps

EDIT - just tried running your example. it throws an error if using an even number for the tValues variable. as it then gives non-integer indices.

Answer 2

[Matt J]

It doesn't look like you are doing any zero-padding.

[Yunus Yilmaz]

I have a 50X50 image and I am trying to obtain 50X50 image reconstructed. How should I process?

[Matt J]

When you convolve A and B (whether directly, or using inverse Fourier transforms), the result is supposed to be of length length(A)+length(B)-1. Your filtered_p_theta_of_t must therefore be large enough to hold the full length of the result.