Search code examples
matlabimage-processingback-projection

Implementing a filtered backprojection algorithm using the central slice theorem in Matlab

I'm working on a filtered back projection algorithm using the central slice theorem for a homework assignment and while I understand the theory on paper, I've run into an issue implementing it in Matlab. I was provided with a skeleton to follow to do it but there is a step that I think I'm maybe misunderstanding. Here is what I have:

function img = sampleFBP(sino,angs)

% This step is necessary so that frequency information is preserved: we pad
% the sinogram with zeros so that this is ensured.
sino = padarray(sino, floor(size(sino,1)/2), 'both');

% diagDim should be the length of an individual row of the sinogram - don't
% hardcode this!
diagDim = size(sino, 2);

% The 2DFT (2D Fourier transform) of our image will start as a matrix of
% all zeros.
fimg = zeros(diagDim);

% Design your 1-d ramp filter. 
rampFilter_1d  = abs(linspace(-1, 1, diagDim))';

rowIndex = 1;
for nn = angs
    % Each contribution to the image's 2DFT will also begin as all zero.
    imContrib = zeros(diagDim);

    % Get the current row of the sinogram - use rowIndex.
    curRow = sino(rowIndex,:);

    % Take the 1D Fourier transform the current row - be careful, as it's
    % necessary to perform ifftshift and fftshift as Matlab tends to
    % place zero-frequency components of a spectrum at the edges.
    fourierCurRow = fftshift(fft(ifftshift(curRow)));


    % Place the Fourier-transformed sino row and place it at the center of
    % the next image contribution. Add the ramp filter in Fourier domain.
    imContrib(floor(diagDim/2), :) = fourierCurRow;
    imContrib = imContrib * fft(rampFilter_1d);


    % Rotate the current image contribution to be at the correct angle on
    % the 2D Fourier-space image.
    imContrib = imrotate(imContrib, nn, 'crop');

    % Add the current image contribution to the running representation of
    % the image in Fourier space!
    fimg = fimg + imContrib;

    rowIndex = rowIndex + 1;
end

% Finally, just take the inverse 2D Fourier transform of the image! Don't
% forget - you may need an fftshift or ifftshift here.
rcon = fftshift(ifft2(ifftshift(fimg)));

The sinogram I'm inputting is just the output of the radon function on a Shepp-Logan phantom from 0 to 179 degrees. Running the code as it is now gives me a black image. I think I'm missing something in the loop where I add the FTs of rows to the image. From my understanding of the central slice theorem, what I think should be happening is this:

  • Initialize an array the same size as the what the 2DFT will be (i.e., diagDim x diagDim). This is the Fourier space.

  • Take a row of the sinogram which corresponds to the line integral information from a single angle and apply a 1D FT to it

  • According to the Central Slice Theorem, the FT of this line integral is a line through the Fourier domain that passes through the origin at an angle that corresponds to the angle at which the projection was taken. So to emulate that, I take the FT of that line integral and place it in the center row of the diagDim x diagDim matrix I created

  • Next I take the FT of the 1D ramp filter I created and multiply it with the FT of the line integral. Multiplication in the Fourier domain is equivalent to a convolution in the spatial domain so this convolves the line integral with the filter.

  • Now I rotate the entire matrix by the angle the projection was taken at. This should give me a diagDim x diagDim matrix with a single line of information passing through the center at an angle. Matlab increases the size of the matrix when it is rotated but since the sinogram was padded at the beginning, no information is lost and the matrices can still be added

  • If all of these empty matrices with a single line through the center are added up together, it should give me the complete 2D FT of the image. All that needs to be done is take the inverse 2D FT and the original image should be the result.

If the problem I'm running into is something conceptual, I'd be grateful if someone could point out where I messed up. If instead this is a Matlab thing (I'm still kind of new to Matlab), I'd appreciate learning what it is I missed.

Solution

  • The code that you have posted is a pretty good example of filtered backprojection (FBP) and I believe could be useful to people who wanted to learn the basis of FBP. One can use the function iradon(...) in MATLAB (see here) to perform FBP using a variety of filters. In your case of course, the point is to learn the basis of the central slice theorem and so finding a short cut is not the point. I have also learned a lot and refreshed my knowledge through answering to your question!

    Now your code has been perfectly commented and describes the steps that need to be taken. There are a couple of subtle [programming] issues that need to be fixed so that the code works just fine.

    First, your image representation in Fourier domain may end up having a missing array due to floor(diagDim/2) depending on the size of the sinogram. I would change this to round(diagDim/2) to have complete dataset in fimg. Be aware that this may lead to an error for certain sinogram sizes if not handled correctly. I would encourage you to visualize fimg to understand what that missing array is and why it matters.

    Second issue is that your sinogram needs to be transposed to be consistent with your algorithm. Hence an addition of sino = sino'. Again, I do encourage you to try the code without this to see what happens! Note that zero padding must be happened along the views to avoid artifacts due to aliasing. I will demonstrate an example for this in this answer.

    Third and most importantly, imContrib is a temporary holder for an array along fimg. Therefore, it must maintain the same size as fimg, so 

    imContrib = imContrib * fft(rampFilter_1d);

    should be replaced with 

    imContrib(floor(diagDim/2), :) = imContrib(floor(diagDim/2), :)' .* rampFilter_1d;

    Note that the Ramp filter is linear in frequency domain (thanks to @Cris Luengo for correcting this error). Therefore, you should drop the fft in fft(rampFilter_1d) as this filter is applied in the frequency domain (remember fft(x) decomposes the domain of x, such as time, space, etc to its frequency content).

    Now a complete example to show how it works using the modified Shepp-Logan phantom:

    angs = 0:359; % angles of rotation 0, 1, 2... 359
init_img = phantom('Modified Shepp-Logan', 100); % Initial image 2D [100 x 100]
sino = radon(init_img, angs); % Create a sinogram using radon transform

% Here is your function ....

% This step is necessary so that frequency information is preserved: we pad
% the sinogram with zeros so that this is ensured.
sino = padarray(sino, floor(size(sino,1)/2), 'both');

% Rotate the sinogram 90-degree to be compatible with your codes definition of view and radial positions
% dim 1 -> view
% dim 2 -> Radial position
sino = sino'; 

% diagDim should be the length of an individual row of the sinogram - don't
% hardcode this!
diagDim = size(sino, 2);

% The 2DFT (2D Fourier transform) of our image will start as a matrix of
% all zeros.
fimg = zeros(diagDim);

% Design your 1-d ramp filter. 
rampFilter_1d  = abs(linspace(-1, 1, diagDim))';

rowIndex = 1;
for nn = angs

    % fprintf('rowIndex = %g => nn = %g\n', rowIndex, nn);
    % Each contribution to the image's 2DFT will also begin as all zero.
    imContrib = zeros(diagDim);

    % Get the current row of the sinogram - use rowIndex.
    curRow = sino(rowIndex,:);

    % Take the 1D Fourier transform the current row - be careful, as it's
    % necessary to perform ifftshift and fftshift as Matlab tends to
    % place zero-frequency components of a spectrum at the edges.
    fourierCurRow = fftshift(fft(ifftshift(curRow)));


    % Place the Fourier-transformed sino row and place it at the center of
    % the next image contribution. Add the ramp filter in Fourier domain.
    imContrib(round(diagDim/2), :) = fourierCurRow;
    imContrib(round(diagDim/2), :) = imContrib(round(diagDim/2), :)' .* rampFilter_1d; % <-- NOT fft(rampFilter_1d)


    % Rotate the current image contribution to be at the correct angle on
    % the 2D Fourier-space image.
    imContrib = imrotate(imContrib, nn, 'crop');

    % Add the current image contribution to the running representation of
    % the image in Fourier space!
    fimg = fimg + imContrib;

    rowIndex = rowIndex + 1;
end

% Finally, just take the inverse 2D Fourier transform of the image! Don't
% forget - you may need an fftshift or ifftshift here.
rcon = fftshift(ifft2(ifftshift(fimg)));

    Note that your image has complex value. So, I use imshow(abs(rcon),[]) to show the image. A couple of helpful images (food for thought) with the final reconstructed image rcon:

    Reconstruction steps and images

    And here is the same image if you comment out the zero padding step (i.e. comment out sino = padarray(sino, floor(size(sino,1)/2), 'both');):

    Reconstructed image without zero padding sinogram

    Note the different object size in the reconstructed images with and without zero padding. The object shrinks when the sinogram is zero padded since the radial contents are compressed.