Image Reconstruction: Backprojection & Fourier

Posted on Mar 8, 2025 in Architecture

Backprojection

This section demonstrates a basic backprojection algorithm for image reconstruction.

Steps:

• Determine the dimensions of the sinogram:

[numOfParallelProjections, numOfAngularProjections] = size(sinogram);

• Initialize a matrix for the backprojected image (BPI):

BPI = zeros(numOfParallelProjections);
d = numOfParallelProjections * sind(45);
c = d * 2;
dim_pad = ceil((c - numOfParallelProjections) / 2);

• Calculate the padded image dimensions:

padded_dim = numOfParallelProjections + 2 * dim_pad;

• Iterate through each projection angle:

for i = 1:numOfAngularProjections
    % Replicate the attenuation profile to create an image.
    Sin_im = repmat(sinogram(:,i), 1, padded_dim);

    % Pad the image to make it square.
    Pad_im = padarray(Sin_im, dim_pad, 0, 'both');

    % Rotate the image clockwise.
    Rot_im = imrotate(Pad_im, -thetas(i), 'bilinear', 'crop');

    % Crop the image to its original size.
    Rot_im_crop = Rot_im(dim_pad + (1:numOfParallelProjections), dim_pad + (1:numOfParallelProjections));

    % Accumulate the rotated and cropped image.
    BPI = Rot_im_crop + BPI;

    imagesc(BPI, "XData", 0:50:150, "YData", 0:20:180)
    colormap('bone')
    drawnow
end

Direct Fourier Reconstruction

This section implements the Direct Fourier Reconstruction method.

• Define the oversampling ratio:

oversampling_ratio = 2;

• Compute the 1D Discrete Fourier Transform (DFT) of the sinogram:

[DFT_sinogram, omega_s] = compute1DFT(sinogram);

• Convert from polar to Cartesian coordinates:

[Cartesian_sinogram, omega_xy] = polar_to_rectA(thetas, omega_s, DFT_sinogram, numOfAngularProjections * oversampling_ratio);

• Perform the 2D inverse Fourier transform:

[Interpolated_sinogram, axis_xy] = inverse_Fourier_2DA(Cartesian_sinogram, omega_xy);

IM_DFR = Interpolated_sinogram;

• Crop and rotate the reconstructed image:

Reconstructed_image = real(IM_DFR);
Reconstructed_image_crop = Reconstructed_image(round((length(Reconstructed_image) - numOfAngularProjections) / 2) : round((length(Reconstructed_image) - numOfAngularProjections) / 2) + numOfAngularProjections - 1, round((length(Reconstructed_image) - numOfAngularProjections) / 2) : round((length(Reconstructed_image) - numOfAngularProjections) / 2) + numOfAngularProjections - 1);
Reconstructed_image_final = rot90(Reconstructed_image_crop, -1);

imagesc(Reconstructed_image_final, "XData", 0:50:150, "YData", 0:20:180)
colormap('bone')

Filtered Backprojection

This section details the Filtered Backprojection (FBP) algorithm.

• Calculate filter coefficients:

filter_FBP = zeros(1, 2 * numOfParallelProjections + 1);
c = 1;
for k = -(numOfParallelProjections + 1):(numOfParallelProjections + 1)
    if k == 0
        filter_FBP(c) = 1/4;
    elseif mod(k, 2) ~= 0
        filter_FBP(c) = (-1) / (pi^2 * (k)^2);
    elseif mod(k, 2) == 0
        filter_FBP(c) = 0;
    end
    c = c + 1;
end

filter_FBP_spec = abs(fft(filter_FBP));

• Iterate through each projection angle and apply the filter:

for i = 1:numOfAngularProjections
    % Filter the attenuation profile.
    filteredProfile = conv(sinogram(:, i), filter_FBP, 'same');
    Sin_im = repmat(filteredProfile, 1, padded_dim);
    Pad_im = padarray(Sin_im, dim_pad, 0, 'both');

    % Rotate the image clockwise.
    Rot_im = imrotate(Pad_im, -thetas(i), 'bilinear', 'crop');

    % Crop the image to its original size.
    Rot_im_crop = Rot_im(dim_pad + (1:numOfParallelProjections), dim_pad + (1:numOfParallelProjections));

    % Sum the rotated and cropped image.
    BPI = BPI + Rot_im_crop;

    colormap('bone')
    imagesc(BPI, "XData", 0:50:150, "YData", 0:20:180)
end