 Research
 Open Access
 Published:
A compressed sensingbased iterative algorithm for CT reconstruction and its possible application to phase contrast imaging
BioMedical Engineering OnLine volume 10, Article number: 73 (2011)
Abstract
Background
Computed Tomography (CT) is a technology that obtains the tomogram of the observed objects. In realworld applications, especially the biomedical applications, lower radiation dose have been constantly pursued. To shorten scanning time and reduce radiation dose, one can decrease Xray exposure time at each projection view or decrease the number of projections. Until quite recently, the traditional filtered back projection (FBP) method has been commonly exploited in CT image reconstruction. Applying the FBP method requires using a large amount of projection data. Especially when the exposure speed is limited by the mechanical characteristic of the imaging facilities, using FBP method may prolong scanning time and cumulate with a high dose of radiation consequently damaging the biological specimens.
Methods
In this paper, we present a compressed sensingbased (CSbased) iterative algorithm for CT reconstruction. The algorithm minimizes the l _{ 1 }  norm of the sparse image as the constraint factor for the iteration procedure. With this method, we can reconstruct images from substantially reduced projection data and reduce the impact of artifacts introduced into the CT reconstructed image by insufficient projection information.
Results
To validate and evaluate the performance of this CSbase iterative algorithm, we carried out quantitative evaluation studies in imaging of both software SheppLogan phantom and real polystyrene sample. The former is completely absorption based and the later is imaged in phase contrast. The results show that the CSbased iterative algorithm can yield images with quality comparable to that obtained with existing FBP and traditional algebraic reconstruction technique (ART) algorithms.
Discussion
Compared with the common reconstruction from 180 projection images, this algorithm completes CT reconstruction from only 60 projection images, cuts the scan time, and maintains the acceptable quality of the reconstructed images.
Background
Computed Tomography (CT), which obtains a series of projection data of objects concerned from several view angles, can get the tomograms of the objects through the technology of image reconstruction. From a purely mathematical standpoint, the solution to the problem of how to reconstruct a function from its projections dated back to the paper by Radon in 1917. The current excitement in tomographic imaging originated with Hounsfield's invention of the Xray computed tomographic scanner for which he received a Nobel Prize in 1972 [1]. The algorithms of image reconstruction from projections can be divided into two classes: the analytical method and the algebraic method [2]. The advantages of the analytical method, such as filtered back projection (FBP) method, are relatively high computational speed and short computational time. When the projection data are densely sampled, images can be reconstructed accurately with analytic methods [3]. Thus, this method is widely used in the commercial CT systems. However, if data containing a reduced number of projections sparsely sampled over an angular range are considered, the analytic algorithms will yield reconstructed images with severe aliasing artifacts such as sharp streaks [4]. The iterative algorithms, on the contrary, can reconstruct images from relatively less projection data. But, it will take much longer time with iterative algorithms versus analytic algorithms.
In realworld applications, especially the biomedical applications, higher temporal resolution and lower radiation dose have been constantly pursued. One can reduce scanning time and radiation dose by decreasing Xray exposure time at each projection view. However, the reduction of exposure time would further lower the signaltonoise ratio (SNR) of the projection images and consequently lower the reconstructed images' quality [5]. The other approach to decrease scanning time and radiation dose is to reduce the number of projections.
Compressed sensing theory, also known as compressive sampling or CS, suggested by Candès, Romberg, Tao and Donoho in 2006 [6, 7], states that one can reconstruct images accurately from a number of samples that are far smaller than the desired resolution of the images [8]. Inspired by the theory's success in signal recovery, we have anticipated that a CSbased algorithm may be used to reconstruct images from substantially reduced projection data. The algorithm minimizes the l _{ 1 } norm of the sparse image as the constraint factor for the iteration procedure. This work focuses on reconstructing images from significantly reduced projection data, shortening scanning time, minimizing radiation dose without reducing image quality, and employing this algorithm in phase contrast imaging experiments.
The paper is organized as follows. In section 2, the experimental setup will be introduced for our polystyrene sample imaging. In section 3, the materials and methods for data scanning and image reconstruction will be described, in section 4, numerical results under different conditions are reported and the reconstructed and reference images at the visualizationbased and quantitativemetricbased evaluation levels are compared. Finally, the implication of the results will be further discussed in section 5.
The experimental setup for phase contrast imaging
Phase contrast Xray imaging [9–13] enables the observation of light samples, such as biological soft tissue, without a contrast agent, because the phase shift cross sections of light elements are much larger than their absorption cross sections [14]. However, in the research of phase contrast imaging, the FBP method has been commonly exploited in CT reconstruction [15–17]. Applying the FBP method requires a large amount of projection data, which can prolong scanning time and cumulate with a high dose of radiation potentially damaging the biological specimens.
Our experiment was performed at the 4W1A beamline of Beijing Synchrotron Radiation Facilities (BSRF). The synchrotron Xray beam was 20 mm in width by 10 mm in height. The Xray beam energy was set at 15keV in this experiment. The distance between the synchrotron radiation source and the sample was approximately 43 m. The Xray CCD was the Photonic Science Xray FDI18 mm camera system with 1300 × 1030 pixels, and 10.9 × 10.9 μm^{2} per pixel. The schematic setup of this analyzer based imaging (ABI) system is shown in Figure 1. It composes of a monochromator crystal, an analyzer crystal, one sample rotation stage, and one image detector. The monochromator crystal is used to produce highly parallel and monochromatic Xray beams. When the highly parallel and monochromatic Xray beams travel through the object, their directions change due to refraction and scattering. According to the Brag diffraction theory, the analyzer crystal only reflects photons coming from a particular angle. Thus, if the analyzer crystal is rotated about an axis perpendicular to the meridian plane, the diffracted intensity will trace out a 'rocking curve' [18]. We obtain the projection data of our polystyrene phantom when the analyzer crystal was set to the half intensity points on the highangle sides of the rocking curve.
Materials and methods
Both software phantom and real sample were used to test our algorithm. The software phantom was a discrete 256 × 256 SheppLogan phantom (see Figure 2(a)). We generated SheepLogan phantom using Matlab command sentence "phantom(256)". We suppose that it is the desired CT image and each pixel value presents an attenuation coefficient. To generate the projection data, we simulated the procedure of an Xray scan, computed line integrals across the image, and eventually, obtained the projection data as shown in Figure 3. More popularly, this projection data is named 'sinogram' [19]. The horizontal and vertical axes represent the detectorbin and viewangle coordinates. In our sample, the number of the detectorbin was 256, and the number of the viewangle had two selections which were 60 and 30. Since minimal levels of noise were introduced to assess the stability of the algorithm, such studies were supposed to give an upper bound to the performance of the CT reconstruction algorithms. The real sample was a polystyrene hexahedron. The length, width, and height of this hexahedron were about 4, 4, and 2 mm respectively. On one 4 mm by 4 mm surface, several lines of holes were punched by a laser gun. We placed the other 4 mm by 4 mm surface on the sample stage. By scanning the sample over 180° with the angular step of 1°, we obtained one hundred and eighty projection images of 1300 × 1030 pixels. Several pieces of projection images in different angles of rotation are shown in Figure 4. Before reconstruction, we moved the background of the projection images. And since the sample was not placed at the absolute center of the sample rotation stage, the projection images should be cut to a suitable size to compose 'sinogram'.
Now, we describe our CSbased iterative algorithm for image reconstruction in parallelbeam CT. A successful application of CS requires that the desired image should have a sparse representation in a known transform domain [7]. Consider an image f, which can be viewed as an N × 1 column vector in R^{N}, whose individual elements f _{ j } , j = 1, 2, ... N are N pixel values of the image. Expand vector f in an orthonormal basis Ψ as follows:
where Ψ is the N by N matrix [ψ _{1}..., ψ _{ N }] with the N × 1 vectors {\left\{{\psi}_{i}\right\}}_{i=1}^{N} as columns and x is also an N × 1 column vector. If all but a few of entries in vector x are zero or almost zero, we will say that f is sparse in the Ψ domain and x is its sparse representation. For example, the SheppLogan phantom in Figure 2(a) and its gradient counterpart in 2(b), we denote the intensity of pixel of a 2D image as f _{ h, w }, where h = 1, 2... H; w = 1, 2... W; H and W are the height and width of the 2D image respectively and W × H = N. If the pixel values are labeled by f _{ h, w }, the gradient modulus is as follows.
We refer to this quantity as the gradient image [20]. The number of nonzero pixels in this 256 × 256 image (Figure 2(a)) is 27521, while the number of nonzero pixels in its gradient image (Figure 2(b)) is only 2182, which is much less than the pixel number of the image. That means clearly, f and x are equivalent representations of the image, with f in the space domain and x in the Ψ domain. In realistic CT imaging, suppose the sampled parallelbeam projectiondata of image f are modeled by a discrete linear system
where vector g has length M with individual measurements referred to as g _{ i } , i = 1, 2, ... M and Φ is the M by N system matrix [21] that yields the discrete set of projection measurements for parallelbeam scanning from the object. Then substituting Ψx for f, g can be written as
where Φ' = Φψ is an M by N matrix. For a sparse image, since M < N in Eq. (4) there are infinitely many \stackrel{\u0303}{x} that satisfy g=\Phi \prime \stackrel{\u0303}{x}. Therefore, the image reconstruction is aimed at finding the vector x in transform domain by solving the linear program [6–8]
Define the l _{ 1 }  norm of the vector x as {\u2225x\u2225}_{{l}_{1}}={\sum}_{i=1}^{N}\left{x}_{i}\right[8]. In this paper, specifically, x represents the gradient image vector of the desired image.
The constraints {\Phi}^{\prime}\stackrel{\u0303}{x}=g can be satisfied under the circumstance that measurements contain no noise. But it is unattainable, so we consider the constraints that
where ε is a small error factor.
To minimize the l _{ 1 }  norm of the gradient image [20, 22, 23], a basic gradient descent method was employed. Usually, the gradient descent method is to reduce the objective function {\chi}^{2}={\u2225\mathsf{\text{g}}\Phi f\u2225}^{2} by iteratively moving the image along the gradient [24]
where α is constant to control the descent speed and \overrightarrow{\Delta} is related to the gradient of the l _{ 1 }  norm of the gradient image which can also be thought of an image. Each pixel value of it is expressed as the following partial derivative [20]
To avoid the condition that the denominator vanishes, a small positive number ε is added in each radical sign.
In the discrete setting, the parallelbeam projectiondata vector \overrightarrow{g} can be written as weighted sums over the pixels traversed by the Xray as
The weight component φ _{ i, j } of the system matrix Φ is computed by the intersection length of the i th ray through the j th pixel. Using a sketch we can understand it clearly. In Figure 5, the image is composed of four pixels f _{0}, f _{1}, f _{2}, and f _{3}; g _{1} is a measurement; and the Xray l _{1} passed pixels f _{0}, f _{1}, and f _{3}; the lengths passed these three pixels are φ _{1, 0}, φ _{1, 1} and φ _{1, 3} respectively. The computation of the weight φ is complex. The immediate computation of each φ will prolong the reconstruction time. Especially with iteration times' increasing, the computation of weights will repeat again and again. A solution is to save the weights in a file in advance and read the weights from this file during the iteration. Since the Xrays are parallelbeam, we can take advantage of the symmetrical characteristic of Xrays. In Figure 6, the Xrays are denoted by a, b, c, and d. Actually, the measurements obtained by the integral along the Xrays a _{ 1 } , b _{ 1 } , c _{ 1 } and d _{ 1 } are detected by the same detector in different rotate angles which are α, (90 α), (90+ α) and (180 α) degrees respectively. The a _{ 1 } , a _{ 2 },..., a _{ m } is a set of parallel beams. These eight Xrays (a _{ 1 } , a _{ m } , b _{ 1 } , b _{ m } , c _{ 1 } , c _{ m } , d _{ 1 } and d _{ m }) in Figure 6 have such symmetrical characteristics that line y = x is the symmetrical axis of a _{ 1 } and b _{ 1 }; the xaxis is the symmetrical axis of b _{ 1 } and c _{ 1 }, also a _{ 1 } and d _{ 1 }; and the parallel beams a _{ 1 } and a _{ m } , b _{ 1 } and b _{ m } , c _{ 1 } and c _{ m } , d _{ 1 } and d _{ m } , are symmetrical to the origin. So, if we have a weight value in Xray a _{ 1 } (the wide segment in line a _{ 1 }), we can gain the other seven weight values (the wide segment in the other seven lines) using the symmetrical characteristics.
In the following, we give the steps of the reconstruction algorithm in the form of a pseudocode and abbreviated notation.

(1)
initialization of the image f:
{f}^{\left(0\right)}=0; 
(2)
iterative process:
{f}_{j}^{\left(k\right)}={f}_{j}^{\left(k1\right)}+\lambda \frac{{g}_{i}\sum _{n=1}^{N}{\varphi}_{i,n}\cdot {f}_{n}^{\left(k1\right)}}{\sum _{i=1}^{N}{\varphi}_{i,n}^{2}}\cdot {\varphi}_{i,j};
where the relaxation parameter λ [25] is a positive real number to adjust the iterative, and k is from 1 to M. When k = M , a complete iteration period was finished. The next iteration will enforce the estimated image to the constraint \left{\Phi}^{\prime}\stackrel{\u0303}{x}g\right<\epsilon by the gradient descent iteration.

(3)
initialization of the gradient descent image:
{\widehat{f}}^{\left(0\right)}={f}^{\left(\mathsf{\text{M}}\right)}; 
(4)
gradient descent iteration:
{\widehat{f}}^{\left(l\right)}={\widehat{f}}^{\left(l1\right)}\alpha \cdot \overrightarrow{\Delta},, and
In this iteration, the end time we selected is 5.

(5)
Initialize the next iterative step:
{f}^{\left(0\right)}={\widehat{f}}^{\left(end\right)},
then we repeat step (2)  (5) until the difference between the current f^{(M)} and the previous f^{(M)} is smaller than the threshold we set or the iteration number is more than 1000.
About the control parameters, we selected λ = 1.0, ε = 0.0001, and α = 0.5 respectively. The threshold value to stop iteration was set as 0.001. These presetting parameters and coefficients only appear to alter the convergence rate.
Results
To demonstrate this CSbased iterative algorithm for image reconstruction from undersampled projection data, we performed two sets of studies: the first set of studies were designed in such a way as to acquire some theoretical understanding of how the CSbased iterative algorithm performs on image reconstruction from reduced projection data with the parallelbeam configuration under ideal conditions, and the second set of numerical examples aimed to see how the CSbased iterative algorithm could be applied to phase contrast CT image reconstruction.
Recall Eq. 3, in the situation that the measurement data g contain no noise and the full scan views data are used, one might expect to reconstruct images accurately. However, in the present studies, the projection data were undersampled in view angle. We performed the FBP, traditional algebraic reconstruction technique (ART) and CSbased iterative reconstruction algorithms under the condition that the numbers of views were 60 and 30. The imagequality evaluation for each specimen was performed at two different levels, including 1) visualizationbased evaluation, and 2) quantitativemetricbased evaluation. Some of the evaluation concerns make a comparison between the reconstructed and original images.
Visual inspection of reconstructions in Figure 7 suggests that under the conditions of fewview (60 and 30view number) projection data, the CSbased iterative algorithm can effectively suppress streak artifacts and noise observed in images obtained with the FBP and traditional ART algorithms, thus yielding images with a higher visual similarity to the SheppLogan phantom image (see Figure 2(a)) than those obtained with other algorithms.
In addition to visualizationbased evaluation, the following three metrics were employed to quantitatively assess the similarity between reconstructed images and the original phantom image: 1) the root mean squared error (RMSE), 2) the universal quality index (UQI) [26], and the correlation coefficient (CC), which are defined as
where vector f _{ r } and f _{ 0 } denote the reconstructed and original images of N pixels, and
The RMSE is widely used for measuring reconstruction accuracy, whereas the UQI and CC can be used for evaluating the pixeltopixel similarity between reconstructed and original images. When assessing the image's quality, we demand the RMSE index to be as small as possible, while expecting the UQI and CC to have the contrary results.
In Figure 7, the images in the left column are the reconstructed images using CSbased iterative algorithm, the middle column using ART, and the right column using FBP. The images in row 1 are reconstructed from 60view data, and in row 2 the images are reconstructed from 30view data.
From the digital SheppLogan phantom and reconstructed images, we computed their RMSEs, UQIs, and CCs and summarized them in Table 1, 2, and 3 respectively. Results of these three metrics suggest that the CSbased iterative algorithm yields images more similar to the original image than the FBP and traditional ART algorithms.
In addition to the simulated experiments with SheppLogan phantom, phase contrast Xray imaging of a real polystyrene phantom was also performed. We collected 180 radiographs at 180 views. From this full data set, we applied the conventional FBP algorithm to yield the reference CT image in the 44th slice (see Figure 8). Then we extracted the 60 and 30view subsets of data evenly distributed over πview to simulate data collected at a reduced number of projection views. The reconstructed images are shown in Figure 9, and the three quantitative metrics RMSEs, UQIs, and CCs are computed and summarized in Table 4, 5, and 6 respectively. The results suggest that the CSbased iterative algorithm can yield the most similar images to the reference image. Different from the results of the reconstructed SheppLogan images, images reconstructed from the traditional ART algorithm seem to have the least similarity to the reference image. The reason might be because the reference image is reconstructed from the FBP algorithm as well.
Discussion and Conclusions
In this article, a CSbase iterative algorithm reconstructing images from substantially reduced projection data was presented. Both the digital SheppLogan phantom and the real polystyrene sample experiment results show that the CSbased iterative algorithm can yield images with quality comparable to that obtained with existing FBP and traditional ART algorithms. However, when the number of gradient descent iterations is increased, the smoothing artifact in the reconstructed images will be more obvious. To improve this situation, reducing the number of gradient descent iterations or the step size is an alternative, but if the gradient descent is too small, the algorithm will reduce to the traditional ART algorithm.
Because the system matrix for sampling is saved in a matrix file, we must create a correlating matrix file for different images in different dimensions at first. The size of such a file, however, relates to the projection data's number and images' dimensions. For this reason, if either of these two elements is rather big, the matrix file will be quite large. Therefore, either compressing the matrix file further or diminishing the region of interest (ROI) in the reconstructed images may be one direction of our further work.
In conclusion, the paper aims to reduce the impact of artifacts introduced into the CT reconstructed image by insufficient projection information. The feasibility of this method which enforces the gradient descent for constraints in traditional iterative algorithms has been demonstrated by both the simulated phantom and the real polystyrene sample experiments. The results show that the CSbased iterative algorithm can yield images with quality comparable to that obtained with existing FBP and traditional ART algorithms. Further research will be performed to improve algorithm efficiency. Moreover, applying this algorithm to a less "sparse" sample such as the real biological soft tissues of small animals and studying how effective the method would be is our future concern.
References
Kak AC, Slaney M: Principles of computerized tomographic imaging. IEEE Press; 1988.
Herman GT: Image reconstruction from projections  implementation and applications. Springer Verlag Press; 1975.
Zou Y, Pan XC, Sidky EY: Theory and algorithms for image reconstruction on chords and within regions of interest. J Opt Soc Am 2005, A22: 2372–2384.
Barrett JF, Keat N: Artifacts in CT: recognition and avoidance. Radiographics 2004, 24: 1679–1691. 10.1148/rg.246045065
Hsieh J: Computed Tomography  principles, designs, Artifacts, and Recent Advances. Bellingham, WA: SPIE Press; 2003.
Candès EJ, Romberg J, Tao T: Robust uncertainty principles: exact signal reconstruction from highly incomplete frequency information. IEEE Trans on Information Theory 2006, 52: 489–509.
Donoho DL: Compressed sensing. IEEE Trans on Information Theory 2006, 52: 1289–1306.
Candès EJ: Compressive sampling. In Proceedings of International Congress of Mathematicians. Madrid, Spain; 2006:1–20.
Davis TJ, Gao D, Gureyev TE: Phasecontrast imaging of weakly absorbing materials using hard Xrays. Nature 1995, 373: 595–598. 10.1038/373595a0
Momose A, Fukuda J: Phasecontrast radiographs of nonstained rat cerebellar specimen. Med phys 1995, 22: 375–379. 10.1118/1.597472
Momose A, Takeda T, Itai Y: Phasecontrast Xray computedtomography for observing biological specimens and organic materials. Review of Scientific Instruments 1995, 66: 1434–1436. 10.1063/1.1145931
Snigirev A, Snigireva I: On the possibilities of xray phase contrast microimaging by coherent highenergy synchrotron radiation. Review of Scientific Instruments 1995, 66: 5486–5492. 10.1063/1.1146073
Wilkins SW, Gao D, Pogany A: Phasecontrast imaging using polychromatic hard xrays. Nature 1996, 384: 335–338. 10.1038/384335a0
Yoneyama A, Momose A, Koyama I: Largearea phasecontrast Xray imaging using a twocrystal Xray interferometer. J Synchrotron Rad 2002, 9: 277–281. 10.1107/S090904950201350X
Momose A, Fukuda J: Demonstration of phasecontrast Xray computedtomography using an Xray interferometer. Nuclear Instruments and Methods in Physics Research 1995, A352: 622–628.
Dilmanian FA, Zhong Z, Ren B: Computed tomography of Xray index of refraction using the diffraction enhanced imaging method. Phys Med Biol 2000, 45: 933–946. 10.1088/00319155/45/4/309
Nesterets YI, Gureyev TE, Wilkins SW: General reconstruction formulas for analyzerbased computed tomography. Appl Phys Lett 2006, 89: 264103–264105. 10.1063/1.2423325
Zachariasen WH: Theory of XRay Diffraction in Crystals. New York: John Wiley and Sons Press; 1945.
Bellon PL, Lanzavecchia S: A direct Fourier method (DFM) for Xray tomographic reconstructions and the accurate simulation of sinograms. Int J Biomed Comput 1995, 38: 55–69. 10.1016/00207101(94)01033W
Sidky EY, Kao CM, Pan XC: Accurate image reconstruction from fewview limitedangle data in divergentbeam CT. Journal of XRay Science and Technology 2006, 14: 119–139.
Candès EJ, Wakin MB: An introduction to compressive sampling. IEEE Signal Processing Magazine 2008, 25: 21–30.
Sidky EY, Pan XC: Image reconstruction in circular conebeam computed tomography by constrained, totalvariation minimization. Phys Med Biol 2008, 53: 4777–4807. 10.1088/00319155/53/17/021
Yang L, Lu Y, Wang G: Compressed sensing inspired image reconstruction from overlapped projections. International Journal of Biomedical Imaging 2010, 2010: 1–8.
Stefan K, Josien PW, Marius S: Adaptive stochastic gradient descent optimization for image registration. International Journal of Computer Vision 2009, 81: 227–239. 10.1007/s112630080168y
Herman GT, Meger LB: Algebraic reconstruction techniques can be made computationally efficient. IEEE Transactions on Medical Imaging 1993, 12: 600–609. 10.1109/42.241889
Bian JG, Siewerdsen JH, Han X, Sidky EY, Prince JL, Pelizzari CA, Pan XC: Evaluation of sparseview reconstruction from flatpaneldetector conebeam CT. Phys Med Biol 2010, 55: 6575–6599. 10.1088/00319155/55/22/001
Acknowledgements
This work was supported by the National Natural Science Foundation of China (Grant Nos. 60532090 and 30770593) and Shanghai Light Source SSRF project (No. 10sr0013). The contents of this paper are solely the responsibility of the authors. The authors would like to thank colleagues of Medical Image Laboratory of Capital Medical University for assistance with data collection and for valuable discussions of studies.
Author information
Authors and Affiliations
Corresponding author
Additional information
Competing interests
The authors declare that they have no competing interests.
Authors' contributions
XL worked on the algorithm design and implementation, and wrote the paper; SL contributed to discussion and suggestions throughout this topic, including the manuscript writing. All authors read and approved the final manuscript.
Authors’ original submitted files for images
Below are the links to the authors’ original submitted files for images.
Rights and permissions
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.
About this article
Cite this article
Li, X., Luo, S. A compressed sensingbased iterative algorithm for CT reconstruction and its possible application to phase contrast imaging. BioMed Eng OnLine 10, 73 (2011). https://doi.org/10.1186/1475925X1073
Received:
Accepted:
Published:
DOI: https://doi.org/10.1186/1475925X1073
Keywords
 Root Mean Square Error
 Projection Data
 Filter Back Projection
 Phase Contrast Imaging
 Algebraic Reconstruction Technique