A novel scheme to design the filter for CT reconstruction using FBP algorithm
- Hongli Shi^{1} and
- Shuqian Luo^{1}Email author
https://doi.org/10.1186/1475-925X-12-50
© Shi and Luo; licensee BioMed Central Ltd. 2013
Received: 6 December 2012
Accepted: 26 May 2013
Published: 1 June 2013
Abstract
Background
The Filtered Back-Projection (FBP) algorithm is the most important technique for computerized tomographic (CT) imaging, in which the ramp filter plays a key role. FBP algorithm had been derived using the continuous system model. However, it has to be discretized in practical applications, which necessarily produces distortion in the reconstructed images.
Methods
A novel scheme is proposed to design the filters to substitute the standard ramp filter to improve the reconstruction performance for parallel beam tomography. The design scheme is presented under the discrete image model and discrete projection environment. The designs are achieved by constrained optimization procedures. The designed filter can be regarded as the optimal filter for the corresponding parameters in some ways.
Results
Some filters under given parameters (such as image size and scanning angles) have been designed. The performance evaluation of CT reconstruction shows that the designed filters are better than the ramp filter in term of some general criteria.
Conclusions
The 2-D or 3-D FBP algorithms for fan beam tomography used in most CT systems, are obtained by modifying the FBP algorithm for parallel beam tomography. Therefore, the designed filters can be used for fan beam tomography and have potential applications in practical CT systems.
Keywords
Background
where $\widehat{f}(x,y)$ denotes the reconstructed image; t=x cosθ+y sinθ; |ω| is known as “ramp” filter in the frequency domain. It is well-known that $\widehat{f}(x,y)$ will be identical with f(x,y) almost everywhere according to the properties of FT and IFT.
where N is a positive even integer denotes the number of projection data; L is an even integer that is equal to or larger than the maximum number of the discrete projection data at all directions; ⌊x⌉ denotes the nearest integer of x; θ_{ j },j∈ [ 1,⋯,K], denote the discretized scanning angle, and K is the number of the scanning angles. The discretized ramp filter, $\left|\frac{k}{L}\right|$, is named as the reconstruction filter in this paper.
For the continuous systems, Radon and inverse Radon transforms are solid and perfect in the mathematics principle [8–10]. However, it necessarily produces non-negligible degradation when the projection data are discrete (finite) and Radon and inverse Radon transforms have to be discretized in calculation. Many scholars have studied this problem. In [15], a multilevel back-projection method had been presented to improve the computational speed. The point-spread-function (PSF) convolution techniques had been proposed to reduce blurring. By those approaches the image quality was similar with or superior to that using the standard FBP technique. In [16], the spline interpolation and ramp filtering had been combined to improve the standard FBP algorithm, by which the image quality could also be improved somewhat.
The question can be summarized as how to reduce the degradation caused by the discretizing process. Since the degradation cannot be removed completely, the question can also be simplified as how to design the optimal reconstruction filter for the discrete inverse Radon transform. In this paper, we try to solve this question in a quite different way. First, the discrete image model and discrete projection model are employed in simulating the scanning procedure. The DFT of projection data is regarded as a special 2-D DFT of the original image using non-uniform frequency sampling. Then, the reconstructed image is regarded as a convolution of the original image and a particular kernel. The kernel is constructed by 2-D DFT, IDFT and the reconstruction filter to be designed. If the kernel become a 2-D Dirac delta function, the reconstructed image will be identical with the original image. So, the optimal reconstruction filter can be obtained by an optimization procedure that make the kernel approach the 2-D Dirac delta function as near as possible.
In order to make the idea behind the design scheme clear, a similar question for 1-D signal is proposed at first, and then it is extended to 2-D situations to solve the corresponding question in CT reconstruction.
Methods
The reconstruction filter for 1-D signal
where ⊙ denotes the periodic convolution operator; x= [x(0),⋯,x(N−1)]; $\widehat{\mathbf{x}}=\phantom{\rule{1em}{0ex}}\left[\widehat{x}\right(0),\cdots \phantom{\rule{0.3em}{0ex}},\widehat{x}(N-1\left)\right]$; h= [h(−N+1),⋯,h(0),⋯,h(N−1)], which is referred to as the convolution kernel.
then $\widehat{x}\left(n\right)=x\left(n\right)$. For example, it can be proven that if ${\theta}_{k}=\frac{2k\pi}{N}$ (the standard DFT), it results h(l)=δ(l), and $\widehat{x}\left(n\right)=x\left(n\right)$. It also means the more h(l) is near δ(l), the more $\widehat{x}\left(n\right)$ is near x(n).
Remark 1
From (6), if h(l)≠δ(l), the aliasing distortion will be produced. The more h(l) is far from δ(l), the more aliasing distortion will be produced. F(k) acts as a correction term to make h(l) approach δ(l) as near as possible. However, the residual difference between h(l) and δ(l) not only depends on F(·) but also depends on the original kernel $\frac{1}{N}{\sum}_{k=0}^{N-1}{e}^{\mathrm{i}{\theta}_{k}l}$. Generally, the difference can be reduced and cannot be removed. The farther the original kernel is different from δ(l), the more residual difference between h(l) and δ(l) may remain.
The coefficients of Taylor series expansion of F ( k )( N =16, θ _{ k } =(2 π k / N ) ^{ 0.8) }
a _{0} | 0.001585600279081 | a _{9} | 0.604131093429063 |
a _{1} | -0.030192172552822 | a _{10} | 0.108486066097029 |
a _{2} | 0.213462131921848 | a _{11} | -0.404865317972757 |
a _{3} | -0.686735955366388 | a _{12} | -0.674935753170146 |
a _{4} | 1.122118847097341 | a _{13} | -0.298926140115388 |
a _{5} | -0.387130280716676 | a _{14} | -0.084322132559712 |
a _{6} | -0.834889525296690 | a _{15} | 0.339643128448386 |
a _{7} | -0.005116780492169 | a _{16} | 0.662450316746720 |
a _{8} | 0.655962903880333 | a _{17} | 0.239001712406937 |
a _{18} | -0.529727742063988 |
In this paper, F(k) is referred to as the reconstruction filter. It can improve the reconstruction performance significantly, which is illustrated by the same example above. The new reconstructed signal is plotted in Figure 2 (c), which is quite similar with the original signal.
The example illustrates the reconstruction performance can be improved by an additional digital filter designed properly. In the next section, a similar scheme is used in designing the reconstruction filter for CT reconstruction.
The reconstruction filter for CT reconstruction
where ⊙ denotes the 2-D periodic convolution operator here; H is the 2-D convolution kernel.
where M is the degree of Fourier series; a_{ m } and b_{ m } are the coefficients to be estimated.
Remark 2
From (16), it is obvious that F(k) depends on N (the image size) and θ_{ j } (the scanning angles). The images of same size with same scanning process will have a same optimal F(k), vice versa. Generally, F(k) can reduce the difference between h_{ r }(n,m) and δ(n,m), however, it cannot remove the difference completely.
Remark 3
From (10), the projection and reconstruction model are relatively simple in this paper. There are many factors are not considered, such as the quantization error (the error that when a pixel is not at any projection line and has to be split between the two nearest projection lines) and projection noise. Even though for such a model, the optimization may be rather complicated and difficult. For example, the objective function is very complicated when the image size is large, and/or the number of projection angles is large.
Results
Example 1
The coefficients of Fourier series expansion of F ( k )( N = 256, θ _{ k } = 0 ^{ o , } 3 ^{ o } ,⋯,177 ^{ o } )
a _{0} | 0.122730836679312 | ||
a _{1} | 0.916483009109464 | b _{1} | -0.497982711366623 |
a _{2} | 0.462544582824757 | b _{2} | -0.120545376607102 |
a _{3} | 0.128736533815128 | b _{3} | 0.040355538652453 |
a _{4} | 0.167851665809050 | b _{4} | -0.182518259771096 |
a _{5} | 0.315690517777738 | b _{5} | -0.094362761222605 |
a _{6} | 0.480769918991562 | b _{6} | 0.285895915083407 |
a _{7} | -0.189184291984498 | b _{7} | 0.563061220832784 |
a _{8} | -0.335163337055193 | b _{8} | -0.086997019632140 |
a _{9} | 0.023005162904716 | b _{9} | -0.087383132662148 |
UQI measures the pixel-to-pixel similarity between the reconstructed f^{ i } and reference image f^{0}. It was designed by modeling any image distortion as a combination of three factors: loss of correlation, luminance distortion, and contrast distortion. Its value ranges between 0 and 1. The closer to 1 the UQI value is, the more similar to the reference image the reconstructed image is.
where $p\left({f}_{k}^{i}\right)$ and $p\left({f}_{k}^{0}\right)$ denote the marginal densities of f^{ i } and f^{0}, respectively, which are calculated using the corresponding histograms; the joint density $p({f}_{k}^{i},{f}_{k}^{0})$ is estimated from the joint histogram of f^{ i } and f^{0}; N^{′} denotes the number of bins in the histogram. UQI measures the pixel-to-pixel dependence of the reconstructed image on the reference image, MI measures the histogram correlation between them. MI can be highly sensitive to small differences between visually similar images.
The results of performance evaluation in Example 1
Criteria | MSE | UQI | MI |
---|---|---|---|
The ramp filter | 1.237E-3 | 0.3429 | 0.2758 |
The designed filter | 9.129E-4 | 0.4751 | 0.4563 |
Example 2
The coefficients of Fourier series expansion of F ( k ) ( N = 128, θ _{ k } =4 k ^{ o } , k = 0,⋯,44)
a _{0} | 0.167768278296600 | ||
a _{1} | 0.742071367715883 | b _{1} | -0.050516729362635 |
a _{2} | -0.007222808863187 | b _{2} | -0.257039580419913 |
a _{3} | 0.205121481459339 | b _{3} | -0.107883292543924 |
a _{4} | 0.080764672022586 | b _{4} | 0.021273467611896 |
a _{5} | -0.013240204304737 | b _{5} | -0.106948545272204 |
a _{6} | 0.297619345807562 | b _{6} | -0.114131247626016 |
a _{7} | 0.184808294809679 | b _{7} | 0.324654368672888 |
a _{8} | -0.198785426563876 | b _{8} | 0.153168707602131 |
a _{9} | -0.059263999435747 | b _{9} | -0.059403297538579 |
The coefficients of Fourier series expansion of F ( k ) ( N = 256, θ _{ k } = k ^{ o } , k = 0,⋯,179)
a _{0} | 0.209555620787293 | ||
a _{1} | 0.693546995671764 | b _{1} | -0.201010002024803 |
a _{2} | 0.099842536865856 | b _{2} | -0.101104690137033 |
a _{3} | -0.152784751604807 | b _{3} | -0.180901677852717 |
a _{4} | 0.248221333497748 | b _{4} | -0.353671098582667 |
a _{5} | 0.308013835831055 | b _{5} | -0.032248278459190 |
a _{6} | 0.377740474429379 | b _{6} | 0.104300291184828 |
a _{7} | 0.075385012026711 | b _{7} | 0.465974969229510 |
a _{8} | -0.294404584616438 | b _{8} | 0.114991610901604 |
a _{9} | -0.056064644119820 | b _{9} | -0.088997280339292 |
The results of performance evaluation in Example 2
(Case 1: S i z e=128∗128,θ_{ k }=4k^{ o }, k=0,⋯,44) | |||
---|---|---|---|
Criteria | MSE | UQI | MI |
The ramp filter | 1.399E-3 | 0.6454 | 0.8022 |
The designed filter | 1.229E-3 | 0.6914 | 0.8518 |
(Case 2: S i z e =256∗256, θ _{ k } = k ^{ o } , k =0,⋯,179) | |||
Criteria | MSE | UQI | MI |
The ramp filter | 2.645E-4 | 0.9219 | 0.9289 |
The designed filter | 2.415E-4 | 0.9384 | 0.9301 |
where a is a parameter selected in region [ 0.5,2]. For example, a=1.65 for the Example 1, and a=1 and a=0.6 for the filters in Example 2, respectively.
Generally, the more complicated model F(k) is of, the better properties it has. However, a very much complicated model F(k) means very much heavy burden of calculation, which may cause the optimal parameters even the valid parameters cannot be found. On the other hand, the difference between h_{ r }(n,m) and δ(n,m) can only be reduced and cannot be removed. So a moderate complicated model F(k), such as (7) with S=18 and (15) with M=9, is an appropriate choice.
Conclusion
For the continuous image model and scanning, FBP algorithm is perfect in the mathematics principle, in which the ramp filter |ω| (it is named as the reconstruction filter in this paper) plays an important role. However, it necessarily causes degradation in the reconstructed images when the continuous image model is discretized and the continuous scanning is substituted by the discrete (finite) scanning in the practical calculation. It means the standard discrete version of ramp filter, |k/N|, is not the optimal for the discrete FBP algorithm. In this paper, a novel scheme is proposed to design the new reconstruction filters to substitute the ramp filter. According to analysis, the reconstructed image can be regarded as the 2-D IDFT of 2-D DFT of the original image using non-uniform frequency sampling. It is also equivalent to a 2-D periodic convolution of the original image and a special 2-D kernel (it is named as the reconstruction matrix in this paper). The more the reconstruction matrix is close to a δ-matrix (a matrix whose elements all are zeros except the center element is one), the more the reconstructed image is close to the original image. Therefore, the reconstruction filters are designed to drive the reconstruction matrixes approach δ-matrix as near as possible. The designs are achieved by the constrained optimization procedures. Some simulation examples have been finished. The results demonstrate the filters designed can make the reconstruction matrixes more close to δ-matrix than the ramp filter. The performance evaluation of CT reconstruction also shows the designed filters have outstanding superiority over the ramp filter in the term of three general criteria, such as MSE, UQI and MI.
Declarations
Acknowledgements
This work have been supported by National Natural Science Foundation (NSFC) of China under Grant No. 60972156 and No. 61227802, and Marie Curie International Research Staff Exchange Scheme (IRSES) actions under the 7th Framework Programme of the European Community under Grant No. PIRSES-GA-2009-269124.
The authors greatly appreciate all the anonymous reviewers for their comments, suggestions, questions, even some revisions on language, which helped to improve the quality of this paper.
Authors’ Affiliations
References
- Kak Avinash C, Slaney M: Principles of Computerized Tomographic Imaging. New York: IEEE Press; 1988.Google Scholar
- Hsieh J: Computed Tomography-Principles, Designs, Artifacts, and Recent Advances. (Bellingham: SPIE Press); 2003.Google Scholar
- Lemmens C, Faul D, Nuyts J: Suppression of metal artifacts in ct using a reconstruction procedure that combines map and projection completion. IEEE Trans Med Imag 2009, 28: 250–260.View ArticleGoogle Scholar
- Jingyan Xu, Taguchi K, Tsui B M W: Statistical projection completion in X-ray CT using consistency conditions. IEEE Trans Med Imag 2010, 29: 1528–1540.View ArticleGoogle Scholar
- Thibault JB, Sauer K, Bouman CA, Heieh J: A three-dimensional statistical approach to improved image quality for multi slice helical CT. Med Phys 2007, 34: 4526–4544. 10.1118/1.2789499View ArticleGoogle Scholar
- Sidky EY, Chartrand R, Pan X: Image reconstruction from few views by non-convex optimization. In IEEE Nuclear Science Symposium Conference Record. Honolulu, HI; 2007:3526–3530.Google Scholar
- Sidky EY, Pan X: Image reconstruction in circular cone-beam computed tomography by constrained, total-variation minimization. Phys Med Biol 2008, 53: 4777–4807. 10.1088/0031-9155/53/17/021View ArticleGoogle Scholar
- Faridani A: Introduction to the mathematics of computed tomography. Inside Out: Inverse Problems, MSRI Publications 2003, 47: 1–46.MathSciNetGoogle Scholar
- Natterer F: The Mathematics of Computerized Tomography. Philadelphia: SIAM; 2001.View ArticleGoogle Scholar
- Riedert A, Faridanit A: The semidiscrete filtered backprojection algorithm is optimal for tomographic inversion. SIAM J Numer Anal 2003, 41: 869–892. 10.1137/S0036142902405643MathSciNetView ArticleGoogle Scholar
- Ye Y, Wang Q: Filtered backprojection formula for exact image reconstruction from cone-beam data along a general scanning curve. Med Phys 2005, 32: 42–48. 10.1118/1.1828673MathSciNetView ArticleGoogle Scholar
- Zamyatin AA, Taguchi K, Silver MD: Practical hybrid convolution algorithm for Helical CT reconstruction. IEEE Trans Nuc Sci 2006, 53: 167–174.View ArticleGoogle Scholar
- Borsdorf A, Raupach R, et al.: Wavelet based noise reduction in CT-images using correlation analysis. IEEE Trans Med Imag 2008, 27: 1685–1703.View ArticleGoogle Scholar
- Ivakhnenko VI: A novel quasi-linearization method for CT image reconstruction in scanners with a multi-energy detector system. IEEE Trans Nuc Sci 2010, 57: 870–879.View ArticleGoogle Scholar
- Brandty A, Jordan M, Brodskiy M, Galun M: A fast and accurate multilevel inversion Of the radon transform. SIAM J Appl Math 1999, 60(2):437–462.View ArticleGoogle Scholar
- Horbelt S, Liebling M, Unser M: Filter design for filtered back-projection guided by the interpolation model. Proc. SPIE 4684, Medical Imaging 2002: Image Processing, 806 (May 15,2002) [http://proceedings.spiedigitallibrary.org/proceeding.aspx?articleid=1313610]
- Xu Q, Yu H, et al.: Low-dose X-ray CT reconstruction via dictionary learning. IEEE Trans Med Imag 2012, 31(9):1682–1697.View ArticleGoogle Scholar
- Borsdorf A, Raupach R, Flohr T,etal: Wavelet based noise reduction in CT-images using correlation analysis. IEEE Trans Med Imag 2008, 27(12):1685–1703.View ArticleGoogle Scholar
- Bian J, Siewerdsen JH, Han X,etal.: Evaluation of sparse-view reconstruction from flat-panel-detector cone-beam CT Phys. Med Biol 2010, 55: 6575–6599. 10.1088/0031-9155/55/22/001View ArticleGoogle Scholar
Copyright
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License(http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.