A list,

Implementation and Application of Filtering Back Projection Reconstruction Algorithm (MATLAB)

  1. Principle of Filtering back projection reconstruction Algorithm Filtering back projection reconstruction algorithm is commonly used in CT imaging reconstruction, the mathematical principle behind it is Fourier transform: one-dimensional Fourier transform of the projection is equivalent to two-dimensional Fourier transform of the original image. (Fourier central section theorem)

CT reconstruction algorithm can be roughly divided into analytical reconstruction algorithm and iterative reconstruction algorithm. With the development of CT technology, reconstruction algorithm has become diversified, each with its own characteristics. In this paper, the most widely used reconstruction algorithm — filtering back projection algorithm (FBP) is used as the basic algorithm of the model. FBP algorithm is a kind of spatial processing technology based on Fourier transform theory. The feature of this algorithm is to convolve the projection under each acquisition projection Angle before backprojection, so as to improve the shape artifact caused by point diffusion function and reconstruct the image with good quality.



The figure above should clearly describe the process of Fourier’s central section theorem: The one-dimensional Fourier transform of the projection is equivalent to the two-dimensional Fourier transform of the original image

The meaning of the Fourier slice theorem is that by performing the Fourier transform on each projection, you can get the two-dimensional Fourier transform from each projection. Thus, the reconstruction of projected images can be solved by the following methods: Collect enough projections at different times (generally 180), solve the one-dimensional Fourier transform of each projection, gather the above slices into the two-dimensional Fourier transform of the image, and then use the Inverse Fourier transform to obtain the reconstructed image. 2. Algorithm process of filtering back projection reconstruction (taking parallel beam as an example) The process of projection reconstruction is that the projection data obtained from the linear array detector is first carried out a one-dimensional Fourier transform, and then convolved with the filter function to obtain the projection data after convolution filtering in all directions; Then, they are back-projected in all directions, that is, they are evenly distributed to each matrix element according to the original path, and the CT value of each matrix element is obtained after overlapping. After proper processing, the algorithm steps of the tomography image of the scanned object are as follows:

  1. Perform a one-dimensional Fourier transform on the original projection
  2. A suitable filter is designed and the original projection P (X_R,φ_i) is obtained at the Angle of φ_i for convolution filtering and the filtered projection is obtained.
  3. The filtered projection is back-projected to obtain the density of the original image in the direction x_r=r cos⁡ ((θ – φ_i)).
  4. Add all the back projections together to get the reconstructed projection.
  5. The selection of filter (filter function) and interpolation function will have two bad influences on experimental results due to the direct use of back projection algorithm:

Inaccurate data reconstruction of the image will produce all kinds of artifacts. The projected data is naturally discrete and can produce large errors if not handled properly. The common filters are R-S filter function and S-L filter function. R-l filter function filtering calculation is simple, avoid a lot of sine and cosine calculation, the obtained sampling sequence is currently divided by sections, and does not significantly reduce the image quality, so the reconstructed image contour is clear, spatial resolution is high.

The most common interpolation methods are nearest neighbor interpolation and double line interpolation. Nearest neighbor interpolation is to replace the missing value in the middle of the discrete point with the projected value at the nearest integer.

Ii. Source code

 
% load a phantom image
im = phantom(256);
% and flatten it to a vector
x = im(:);
 
%% Setting up the geometry
% projection geometry
proj_geom = astra_create_proj_geom('parallel'.1.256, linspace2(0,pi,180));
% object dimensions
vol_geom  = astra_create_vol_geom(256.256);
 
%% Generate projection data
% Create the Spot operator for ASTRA using the GPU.
W = opTomo('cuda', proj_geom, vol_geom);
 
p = W*x;
 
% reshape the vector into a sinogram
sinogram = reshape(p, W.proj_size);
imshow(sinogram, []);
 
 
%% Reconstruction
% We use a least squares solver lsqr from Matlab to solve the 
% equation W*x = p.
% Max number of iterations is 100, convergence tolerance of 1e-6.
y = lsqr(W, p, 1e-6.100);
 
% the output is a vector, so we reshape it into an image
reconstruction = reshape(y, W.vol_size);
 
subplot(1.3.1);
imshow(reconstruction, []);
title('Reconstruction');
 
subplot(1.3.2);
imshow(im, []);
title('Ground truth');
 
% The transpose of the operator corresponds to the backprojection.
backProjection = W'*p;
subplot(1.3.3);
imshow(reshape(backProjection, W.vol_size), []);
title('Backprojection');
Copy the code

3. Operation results

Fourth, note

Version: 2014 a