A Novel Iterative CT Reconstruction Approach Based on FBP Algorithm

The Filtered Back-Projection (FBP) algorithm and its modified versions are the most important techniques for CT (Computerized tomography) reconstruction, however, it may produce aliasing degradation in the reconstructed images due to projection discretization. The general iterative reconstruction (IR) algorithms suffer from their heavy calculation burden and other drawbacks. In this paper, an iterative FBP approach is proposed to reduce the aliasing degradation. In the approach, the image reconstructed by FBP algorithm is treated as the intermediate image and projected along the original projection directions to produce the reprojection data. The difference between the original and reprojection data is filtered by a special digital filter, and then is reconstructed by FBP to produce a correction term. The correction term is added to the intermediate image to update it. This procedure can be performed iteratively to improve the reconstruction performance gradually until certain stopping criterion is satisfied. Some simulations and tests on real data show the proposed approach is better than FBP algorithm or some IR algorithms in term of some general image criteria. The calculation burden is several times that of FBP, which is much less than that of general IR algorithms and acceptable in the most situations. Therefore, the proposed algorithm has the potential applications in practical CT systems.


Introduction
Today the medical images provide very important health information in clinical diagnosis, in which X-ray CT image is one of most important modalities. The quality of CT images heavily depends upon the reconstruction algorithms, especially in the case that only less and less projection data is available in order to reduce the radiation dose. The current algorithms can be roughly divided into two categories: 1) analytical reconstruction algorithms, 2) iterative reconstruction (IR) algorithms. Some of IR algorithms, such as the algorithms based on dictionary learning or compressed sensing(CS) theory, have become the research focuses and have been used in a few special fields, which can solve some reconstruction problems where projection data is far from the requirement of Shannon theorem [1][2][3][4][5][6]. However, these algorithms suffer from heavy calculation burden, poor convergence speed and other drawbacks. For example, the iteration number may be several hundreds of thousands [5], and the average computational time of one iteration may be several seconds [7]. Generally, IR algorithms are based on some hypothetical conditions which are not always satisfied in practice. For example, the algorithm using total variation minimization may cause the details to be weakened or removed. On other hand, the analytical schemes, i.e. FBP algorithm and its modified versions such as FDK (Feldkamp-Davis-Kress) algorithm, are much simpler and faster [8][9][10][11][12][13][14]. They have been used in almost all the fields of straight ray tomography, such as X-ray CT and PET (Positron Emission Tomography) [15][16][17].
For the continuous systems, Radon and inverse Radon transforms using FBP algorithm are in a close form in the mathematics principle [10,11,13]. However, it necessarily produces non-negligible aliasing degradation when the projection data and FBP algorithm have to be discretized. Many approaches had been proposed to deal with this problem. As the description in [18], a multilevel back-projection method had been presented to improve the computational speed while a point-spread-function (PSF) convolution techniques had been proposed to reduce blurring in reconstruction. As a result, the image quality was similar with or superior to that using the classic FBP technique. In [19], the spline interpolation and ramp filtering had been combined to improve the classic FBP algorithm, by which the image quality could also be improved somewhat. In [20], a new filter has been designed to substitute the classic ramp filter to improve reconstruction performance.
In this paper, an iterative CT reconstruction approach based on FBP algorithm is proposed, which is designed to hold both advantages while reduce both disadvantages of analytical and IR algorithms. In the algorithm, the classic FBP is utilized to obtain the initial reconstructed image, which is treated as the intermediate image. It is then projected along the original projection directions to produce the reprojection vectors. The difference between the real and reprojection vector is filtered by a special digital filter to produce the corrected projection term, which is then performed the inverse Radon transform with same parameters to produce a correction term in image domain. By adding the correction term to the previous one, the new intermediate image is obtained. The digital filter is designed to make the new reprojection vectors approach to the real projection vectors as much as possible. This procedure can be performed iteratively until certain stoping criterion is achieved.
Generally, the projection vectors (together with projection angles) are all the information that we have obtained about the image to be reconstructed. If the reprojection data obtained using the proposed approach better resembles the real projection data than the one obtained using the classic FBP, it is reasonable to regard the proposed approach is better than the latter one.
The digital filter is very important for the proposed approach. In order to make the idea behind the design scheme for the filter clear, the reason that produces aliasing degradation by the discretizing process is analyzed at first. Then, a correction scheme is proposed.

Methods
At first of this section, FBP algorithm is introduced in brief, and then the aliasing degradation caused by the discretizing process in FBP algorithm is analyzed. At last, an iterative scheme is proposed to reduce the aliasing degradation.

FBP algorithm
There are two projection modes for the general tomography: parallel beam and fan beam projection. Since a reconstruction problem for the latter can be easily transformed as a problem for the former, only the former is studied in this paper.
The derivation of FBP algorithm for parallel beam tomography is rather simple and straightforward [10,11]. The projection procedure can be instantiated in Fig 1. The Fourier slice theorem shows that the two-dimensional (2-D) Fourier transform (FT) at frequency sample (ω cos θ, ω sin θ), F(ω cos θ, ω sin θ), can be obtained by the one-dimensional (1-D) FT of projection vector at angle θ (at same time, the expression of beam lines in Cartesian coordinates is transformed as the expression in polar coordinates). Therefore, the 1-D inverse FT (IFT) can substitute 2-D IFT, by which the original image can be reconstructed, namely, the inverse Radon transform is achieved. The Jacobian matrix of Cartesian-polar transformation becomes as the ramp filter.
In the calculation, the projection data and reconstructed images have to be discretized. For a discrete projection vector, p θ i (t), t 2 [−L/2, Á Á Á, L/2] (denoted by p i (t) in the next), the discrete FT (DFT) and inverse DFT (IDFT) are employed to approximate FT and IFT, respectively. The FBP algorithm becomes [8] wheref ðn; mÞ denotes the reconstructed image; L is a positive even integer denotes the size of projection vectors; N is an even integer that is equal to or larger than the maximum size of the projection vectors at all directions; j k N j is the "ramp" filter in the frequency domain; bxe denotes the nearest integer of x; θ i , i 2 [1, Á Á Á, K], denote the discretized scanning angle, and K is the number of the scanning angles. Generally, FFT (fast Fourier transform) is utilized to speed up reconstruction.

The aliasing degradation analysis
In this subsection, the implementation of FBP is analyzed in detail to show the cause of aliasing degradation. As shown in Eq (1), the implementation of FBP algorithm can be split into two steps: filtering and interpolating. The filtering can also be expressed in spatial domain as following after IDFT and some simplifications.
When N is large enough such as N ! 128, β N (t) will be substituted by β(t) from now on, where the error caused by substitution will be very small and can be ignored. Since β(t) is symmetrical, the Eq (2) can also be expressed as the circular convolution of the original projection vectors and a kernel h n q i ðtÞ % p i ðtÞ h n ðtÞ ð3Þ where denotes the circular convolution operator, In the second step of FBP algorithm, for the filtered projection q i (t), the intensity of a pixel whose position (n, m) satisfy t − 1 < n cos θ i + m sin θ i < t+1 will be added with the weighted q i (t), which is rather similar with the discrete projection procedure. In that procedure, each pixel's contribution is proportionally split into two nearest projection lines according to the distances between the projection lines and the pixel. The linear weighting coefficient is the common choice in this step, i.e., the intensity of point (n, m) will be increased by the linear interpolation of q i (t), t = −L/2, Á Á Á L/2, which is shown in Fig 2. Finally, the reconstructed image is obtained by accumulating the contributions of all projection angles. Now, suppose the reconstructed image is performed Radon transform along the identical angles again. The reprojection vector along θ i with a distance t from the rotation center,p i ðtÞ, can be expressed asp where S denotes the set of pixels that lies between x cos θ i + y sin θ i = t − 1 and x cos θ i + y sin θ i = t + 1; s(t) indicates the weight of a pixel contributes to the projection linep i ðtÞ. Since each pixel's contribution is proportionally split into the two projection lines that sandwich the pixel according to the distances between the projection lines and the pixel, s(t) is a linear function of t. By substitutingf ðn; mÞ with Eq (1), the reprojection vector can be expressed as following after some simplificationsp where α i (t) denotes the new coefficient, which is also the linear function of t. Consider Eq (3), the Eq (4) can further be expressed as following for briefness.
where Y i ðtÞ ¼ P ðn;mÞ2S P K j¼1 a i ðtÞp i ðtÞ. Remark: The Eq (5) is a summary conclusion about the original projection and reprojection vectors. It can also be comprehended in the following way. Since q i (t) % p i (t) h n (t), the interpolation in the reconstruction and the follow-up reprojection are all linear process when the linear interpolation mode is chosen in the reconstruction process, the reprojection vector can be expressed as the circular convolution of the ramp filter h n and a term that is the linear combination of original projection vectors.

The correction scheme
The original projection vectors (and the corresponding projection angles) are the best knowledge we have got about the unknown object. The Eq (5) shows there is a distinct difference between an original projection vector and the corresponding reprojection one, which indicates the imperfection degree of reconstruction. We seek to take the advantage of the original projection data again in an iteration way to minimize the difference, in other words, improve the reconstructed image quality step by step. The proposed iterative scheme is shown as the flowchart in Fig 3. The idea behind the processing flow is to make the reprojection vectors gradually approach the original projection vectors. First, obtain the initial reconstructed imagef by FBP algorithm. f is then performed Radon transform using the identical setting to obtain the reprojection vectorsp i . Second, calculate the difference between p i andp i , Dp i ¼ p i Àp i , which is filtered to produce a filtered projection residual Δ f p i = Δp i F. At last, an image correction term is obtained by Δ f p i , which is add tof to updatef . The whole flow can be performed iteratively.
The design of the filter F becomes a very key issue. The filter should be designed to meet the following requirement YðtÞ h n ðtÞ þ ð½pðtÞ À YðtÞ h n ðtÞ FðtÞÞ h n ðtÞ % pðtÞ ð6Þ where p(t) denotes an original projection vector; Θ(t) h n (t) denotes the corresponding reprojection vector; p(t) − Θ(t) h n (t) denotes the projection residual, (p(t) − Θ(t) h n (t)) F(t) denotes the filtered version; ((p(t) − Θ(t) h n (t)) F(t)) h n denotes the correction term, which is added to Θ(t) h n (t) to make the updated reprojection vector approach the original projection vector p(t). It requires h n ðtÞ FðtÞ ¼ dðtÞ ð7Þ where δ denotes the discrete Dirac delta function.
In the design, suppose F(t) is symmetrical and its length is identical with that of h n . For example,

Results
In this section, some numerical simulations and tests on real data are shown to demonstrate the performance of the proposed reconstruction algorithm together with the designed filter. All animal experiments and procedures carried out on the animals are approved by the animal welfare committee of Capital Medical University and the approval ID is SCXK-(Army) 2013-X-99. Example 1. This example is employed to show the filter designed can make the reprojection vectors approach the original projection vectors very well. The original image I selected is the head phantom generated by Matlab function phantom with 128 × 128 pixels. It is performed Radon transform using radon with the angle vector θ = [0°,1°, Á Á Á, 179°] to produce p 0 . The reconstructed imageÎ is obtained by iradon (where the classic FBP is utilized) using the linear interpolation mode.Î is performed Radon transform again to get the reprojection vectors p 1 .
Then, calculate Δp f = (p 0 −p 1 ) F, the correction imageR and the corrected imageÎ þR. Iteratively, a same Radon transform is performed to using the corrected image to produce new projection vectors, and so on. This procedure is iterated twice following the flowchart in Fig 3 (it is only a part of the proposed algorithm as the decision process is absent). Since they are matrixes (2d), only one vector of each p 0 , p 1 and p 3 (the reprojection after twice iteration) is shown in Fig 4, whose angle is θ = 47°.
The results show that the difference between the corrected reprojection vector p 3 47 and p 0 47 is much smaller than that of the reprojection vector p 1 47 and p 0 47 . In order to assess the efficiency of the designed filter, MSE (mean square error) is employed as the criterion. Let where p 0 i , p 1 i and p 3 i denote the original, reprojection and corrected projection vectors for the projection angle θ i ; L is the size of a vector, and K denotes the number of projection angles. The results for two examples are shown in Table 1.
The results illustrate the proposed scheme together with the designed filter make the corrected reprojection vectors approach the true projection vectors very well.
Example 2. The original image I is also the phantom image of 128 × 128 pixels. The angle vector is selected θ = [0°,0.3°, Á Á Á, 179.7°], along which I is projected to produce the original projection P. It is reconstructed by the classic FBP and the proposed iterative FBP. The original image and its reconstructed versions by different schemes are shown in Fig 5. As [20], three criteria, MSE (Mean Square Error), UQI (Universal Quality Index) and MI (mutual information), are employed to assess the efficiency of the proposed approach, which are defined as following [21,22].
where I r k and I 0 k denote the pixels of the reconstructed image I r and reference image I°, respectively; M is the total number of pixels. UQI is defined as following where I and σ 2 denote the image mean and variance, respectively; Cov(I r , I°) denote the covariance of the reconstructed image I r and reference image I°. The mean, variance and covariance are defined as the following When the reconstructed image and reference image are interpreted as "stochastic processes",  MI is used for measuring their mutual dependence.
where pðI r k Þ and pðI o k Þ denote the marginal densities of I r and I°, respectively, which are calculated using the corresponding histograms; the joint density pðI r k ; I o k Þ is estimated from the joint histogram of I r and I°; N 0 denotes the number of bins in the histogram.
UQI measures the pixel-to-pixel similarity between the reconstructed I r and reference image I°, while MI measures the histogram correlation between them. The closer to 1 the UQI value is, the more similar the two images are. Similarly, the higher the MI value is, the more similar the two images are. Since it is a simulation example, the original image is known and selected as the reference image.
The results in the form of the three criteria are showed in Table 2, in which the phantom images with different sizes and projection settings are employed. The results illustrate the proposed scheme have better reconstruction performance than the classic FBP algorithm.
Example 3. In order to compare the iterative FBP with the current iterative schemes in the terms of speed and accuracy, some reconstruction problems are implemented by SIRT (Simultaneous Iterative Reconstruction Technique), SART (Simultaneous Algebraic Reconstruction Technique) and MAP-EM (Maximum A Posteriori estimation Expectation Maximization) algorithms [23][24][25]. Unlike the iterative schemes such as the scheme based on minimization of  total variation (TV) regularization or dictionary-learning, these algorithms do not make any artificial assumption on the image to be reconstructed. For the schemes based on TV minimization, there is an assumption that the TV of the underlying image is minimum. For the schemes based on dictionary-learning, the assumption is that the underlying image can be expressed by a linear combination of atoms, the high-dimensional vectors, in a dictionary, the vector set. SIRT, SART and MAP-EM schemes are identical with the proposed scheme in the regard, which depend only on the projection data rather than relying on a prior information or assumption. Thus, they are chosen to compare with the proposed scheme. The head phantom with 512 × 512 pixels is selected as the original image, which is performed Radon transform with the angles θ = [0°,0.5°, Á Á Á, 179.5°] to produce original projection p 0 . The image is reconstructed using the three schemes, and the three criteria, MSE, MI and UQI are employed to evaluated the reconstruction performance. For all the schemes, the image reconstructed using the classic FBP is chosen as the initial image of iteration. The results is showed in Table 3 and Fig 6. The results show that the iterative FBP algorithm has the similar reconstruction performance with SART, SIRT and MAP-EM schemes, and a much smaller operating time. We think the superiority of iterative FBP in the term of speed mainly roots in the fact that it does not employ the linear expression of projection procedure. The expression refers to a very huge matrix (see the Eq (1) or Eq (2) in [25] or in page 277 [8]) and will slow down the calculation, however, it is necessary for SIRT, SART and MAP-EM schemes. When the image size and the number of projection angle is large (as that in this example), the matrix becomes very huge and has to be split into many blocks. Every block has been saved and loaded at least one time, which spends many time.
In order to analyze the robustness of iterative FBP, suppose the original projection is polluted by the white Gaussian noise, i.e., where p 0 denotes the original projection vectors, which is the Radon transform of the head phantom image of 256 × 256 pixels with the projection angles θ = [0°,1°, Á Á Á, 179°]; Z is the white Gaussian noise with same size as p 0 , whose mean is zero and variance is one; σ is the magnitude of noise, which is selected as 0.1max(p 0 ). The image is reconstructed using the three schemes, SART, SIRT and MAP-EM, and Iterative FBP scheme. In all the schemes, the image that is reconstructed using the classic FBP is chosen as the initial image of iteration. The some results are showed in Fig 7. The results show that the iterative FBP algorithm has the poor but acceptable reconstruction performance, while SART, SIRT and MAP-EM have very bad performance even after a rather long running time. SIRT and MAP-EM have worse tendency than SART with the iterative recurrence proceeds. We think the weaknesses mainly roots in the fact that they depend upon Iterative FBP the linear expression of projection procedure, which make the reconstruction performance is very sensitive to the noise in the projection data. Example 4. In this example, the proposed algorithm is used in the practical application. The sample is a mouse liver, and the experiment is performed in the X-ray Imaging and Biomedical Application Beamline station (BL13W1) at the SSRF(Shanghai Synchrotron Radiation Facility, China). All animal experiments and procedures carried out on the animals are approved by the animal welfare committee of Capital Medical University and the approval ID is SCXK-(Army) 2013-X-99. The setup can ensure x-ray beam to be near parallel. The energy is about 21 keV. A CCD camera (the size of pixel is 13μm × 13μm) is used as the detector, comprising 2,588 × 458 pixels. During scanning, the sample is rotated on a turntable around its cylindrical axis by 180°a t step of 0.4411°. The rotation speed is about 0.25°/s and the exposure time was 11 millisecond. Before and after scan, the background images (the sample is absent) and dark field images (the x-ray beam is closed) are recorded for preprocessing. The background images are used for Iterative FBP normalization, the dark field images are used to reduce the noise of various kinds(mainly camera noise). Finally, the logarithm transform is employed to enhance the image contrast, where p b i , p d i and p o i denote i-th column of background image, dark field image and projection image. The images are reconstructed by the classic FBP and the proposed iterative FBP, which are shown in Fig 8. Since there is not a reference image, it is difficult to assess the image quality by some quantitative indexes. The difference between the real projection vectors and the reprojection vectors, i.e., s 1 and s 2 defined in Eq (9), are employed as the criterion. In this example, s 1 = 299.4172, s 2 = 93.4903.

Discussion
As is shown in Fig 2 and Eq (5), the projection and reconstruction model are relatively simple in this paper. There are some factors are not considered, such as the quantization error (for example in radon function in Matlab, every pixel is divided into four sub-pixel and projects each sub-pixel separately) and projection noise. For such a model, the optimization of Eq (7) or Eq (8) will be rather complicated and difficult.
As previously shown in the flow chart (Fig 3), the computation burden induced by the iterative procedure mainly comes from the calculation of the additional reconstruction and filtering procedures. Obviously, the reconstruction part is the dominant portion. Generally, the number of iterative loops is unnecessary to be larger than four, which means the extra calculation burden is only twice or triple that of the classic FBP. To the best of our knowledge, the calculation burden is much less than that of the most IR algorithms, and is acceptable in general.
In the practical applications, since there are many factors that may degrade the quality of reconstructed image, such as the eccentricity of cylindrical axis and measurement noise, the image quality cannot be improved greatly. The eccentricity of cylindrical axis is the most serious factor for the iterative FBP algorithm because the distortion caused by the offset of projection data will be accumulated as the iterative procedure proceeds. The larger the eccentricity of cylindrical axis is, the less times the iterative procedure should be performed. According to our experience, for the experiment setups such as in SSRF, once iterative procedure is an appropriate choice.

Conclusion
FBP algorithm is perfect for the continuous image model and scanning, however, the aliasing degradation cannot be avoided when the image and scanning model are discretized. Generally, there is a significant difference between the original projection data and the reprojection data using the classic FBP algorithm. Since the original projection data are the best information we have, a proper way is to take advantage of the original projection again to make the reprojection data approach the original projection data while improve reconstruction performance.
According to the analysis, a reprojection vector can be approximately expressed as the circular convolution of a certain kernel and a linear combination of the original projection vectors. The difference between the original projection and reprojection vector is served as the motivation to improve the reconstruction performance. It is filtered by a digital filter to compensate the difference caused by FBP, where the filter is designed to remove the filtering effect of the ramp filter. As a result, the sum of the reprojection vector and the filtered difference will approach the original projection vector. In the meantime, the filtered difference is performed inverse Radon transform and the result is added to the reconstruction image of last step to update it. This process can be repeated for several times until certain stopping criterion is satisfied. The results of some numerical simulations and practical applications demonstrate the proposed scheme have better superiority over the classic FBP algorithm in the term of MSE, UQI and MI. The results also show the proposed scheme have much better superiority over some iterative algorithms that depend only on the projection data, such as SIRT, SART and MAP-EM, in the term of reconstruction convergence speed and robustness.