Tomographic Reconstruction Algorithms: Backprojection, DFR, and FBP
Tomographic Reconstruction Algorithms
Backprojection
Load sinogram: load X.mat
Plot the sinogram: imagesc(sinogram)
Determine the value of the constants numOfAngularProjections and numOfParallelProjections. Each one represents the amount of orientations of the X-ray tube (theta) and the number of X-ray detectors (r), respectively. The latter variable (numOfParallelProjections) defines the original size of the image: size(sinogram)
Generate the array thetas that will contain the angular orientations of the X-ray tube.
Generate the matrix that will contain the backprojected image and initialize it to zeros (BPI).
Calculate the amount of padding required. This variable should have size 1×1 (ceil).
Calculate the size of the image after the padding process.
Loop over each projection: for i = 1:numOfAngularProjections
Convert the attenuation profile for this particular orientation into an image: repmat(sinogram(:, i), 1, padded_dim)
Pad the image so it is square. The dimensions should be pad_dim x pad_dim: padarray(Sin_im, dim_pad, 0, 'both')
Rotate the image clockwise: imrotate(Pad_im, -thetas(i), 'bilinear', 'crop')
Crop the image to its original size. The image after this step should have dimensions numOfParallelProjections x numOfParallelProjections: Rot_im_crop = Rot_im(dim_pad + (1:numOfParallelProjections), dim_pad + (1:numOfParallelProjections));
Sum the rotated and cropped image with the previously obtained images: BPI = BPI + Rot_im_crop
Direct Fourier Reconstruction (DFR)
Load sinogram: load X.mat
Define the original image size (N_image) and define the angular projection array (thetas):
N_image = length(sinogram) or size(sinogram, 1)
thetas = 0:180/size(sinogram, 2):(180-180/size(sinogram, 2))
Define the variable for the oversampling_ratio and set it to 2. This will increase the Nyquist frequency and reduce aliasing: oversampling_ratio = 2
Use the provided functions to reconstruct the image using the Direct Fourier Reconstruction: [DFT_sinogram, omega_s] = compute1DFT(sinogram)
[Cartesian_sinogram, omega_xy] = polar_to_rectA(thetas, omega_s, DFT_sinogram, N_image * oversampling_ratio) (polar_to_rectA — To pass from polar coordinates to cartesian)
[Interpolated_sinogram, axis_xy] = inverse_Fourier_2DA(Cartesian_sinogram, omega_xy)
The various Fourier and inverse Fourier transforms may generate spurious imaginary components in the image (the pixels are represented as complex numbers instead of real numbers). Remove the imaginary part of the pixels: Reconstructed_image = real(Y)
Crop image to its original size (N_image x N_image).
Rotate the image 90º clockwise using the command rot90 (Don’t need padding — 90º / 180º / 270º — They decay to another square): Reconstructed_image_final = rot90(Reconstructed_image_crop, -1);
Filtered Backprojection (FBP)
Load sinogram: load C.mat
Determine the value of the constants numOfAngularProjections and numOfParallelProjections. Each one represents the amount of orientations of the X-ray tube (theta) and the number of X-ray detectors (r), respectively: size(sinogram)
Generate the array thetas that will contain the angular orientations of the X-ray tube.
Generate the matrix that will contain the backprojected image and initialize it to zeros (BPI).
Calculate the amount of padding required. This variable should have size 1×1 (A = ceil).
Calculate the size of the image after pad process (pad_dim).
Calculate the filter coefficients: zeros + calculate + abs(fft(filter_FBP))
Iterate for each orientation: for i = 1:numOfAngularProjections
Filter the attenuation profile that corresponds to the current iteration: Y = conv(sinogram(:, i), filter_FBP, 'same')
Convert the attenuation profile for this particular orientation into an image: Z = repmat(Y, 1, pad_dim)
Pad the image so it is square. The dimensions should be pad_dim x pad_dim: W = padarray(Z, A, 0, 'both')
Rotate the image clockwise: B = imrotate(W, -thetas(i), 'bilinear', 'crop')
Crop the image to its original size. The image after this step should have dimensions numOfParallelProjections x numOfParallelProjections: T = B(A + (1:numOfParallelProjections), A + (1:numOfParallelProjections))
Sum the rotated and cropped image with the previously obtained images: BPI = BPI + T