# A compressed sensing-based iterative algorithm for CT reconstruction and its possible application to phase contrast imaging

- Xueli Li
^{1}and - Shuqian Luo
^{1}Email author

**10**:73

https://doi.org/10.1186/1475-925X-10-73

© Li and Luo; licensee BioMed Central Ltd. 2011

**Received: **9 May 2011

**Accepted: **18 August 2011

**Published: **18 August 2011

## Abstract

### Background

Computed Tomography (CT) is a technology that obtains the tomogram of the observed objects. In real-world applications, especially the biomedical applications, lower radiation dose have been constantly pursued. To shorten scanning time and reduce radiation dose, one can decrease X-ray 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 sensing-based (CS-based) 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 CS-base iterative algorithm, we carried out quantitative evaluation studies in imaging of both software Shepp-Logan phantom and real polystyrene sample. The former is completely absorption based and the later is imaged in phase contrast. The results show that the CS-based 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.

## Keywords

## 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 X-ray 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 real-world 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 X-ray exposure time at each projection view. However, the reduction of exposure time would further lower the signal-to-noise 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 CS-based 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 set-up 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 visualization-based and quantitative-metric-based evaluation levels are compared. Finally, the implication of the results will be further discussed in section 5.

### The experimental set-up for phase contrast imaging

Phase contrast X-ray 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.

^{2}per pixel. The schematic set-up 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 X-ray beams. When the highly parallel and monochromatic X-ray 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 high-angle sides of the rocking curve.

## Materials and methods

*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:

*Ψ*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 Shepp-Logan 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.

*256 × 256*image (Figure 2(a)) is 27521, while the number of non-zero 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 parallel-beam projection-data of image

*f*are modeled by a discrete linear system

*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 parallel-beam scanning from the object. Then substituting Ψx for f, g can be written as

*' =*Φψ 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.

where ε is a small error factor.

*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]

*α*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.

*φ*

_{ 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 X-ray

*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 X-rays are parallel-beam, we can take advantage of the symmetrical characteristic of X-rays. In Figure 6, the X-rays are denoted by

*a, b, c*, and

*d*. Actually, the measurements obtained by the integral along the X-rays

*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 X-rays (

*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 x-axis 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 X-ray

*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.

- (1)initialization of the image
*f*:${f}^{\left(0\right)}=0;$ - (2)iterative process:${f}_{j}^{\left(k\right)}={f}_{j}^{\left(k-1\right)}+\lambda \frac{{g}_{i}-\sum _{n=1}^{N}{\varphi}_{i,n}\cdot {f}_{n}^{\left(k-1\right)}}{\sum _{i=1}^{N}{\varphi}_{i,n}^{2}}\cdot {\varphi}_{i,j};$

*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:

- (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 CS-based iterative algorithm for image reconstruction from under-sampled 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 CS-based iterative algorithm performs on image reconstruction from reduced projection data with the parallel-beam configuration under ideal conditions, and the second set of numerical examples aimed to see how the CS-based 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 under-sampled in view angle. We performed the FBP, traditional algebraic reconstruction technique (ART) and CS-based iterative reconstruction algorithms under the condition that the numbers of views were 60 and 30. The image-quality evaluation for each specimen was performed at two different levels, including 1) visualization-based evaluation, and 2) quantitative-metric-based evaluation. Some of the evaluation concerns make a comparison between the reconstructed and original images.

*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 pixel-to-pixel 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 CS-based iterative algorithm, the middle column using ART, and the right column using FBP. The images in row 1 are reconstructed from 60-view data, and in row 2 the images are reconstructed from 30-view data.

Quantitative assessment RMSE of Shepp-Logan phantom images

RMSE | CS-based Iteration | ART | FBP |
---|---|---|---|

60-view number | 16.82 | 22.61 | 37.92 |

30-view number | 17.26 | 25.07 | 40.76 |

Quantitative assessment UQI of Shepp-Logan phantom images

UQI | CS-based Iteration | ART | FBP |
---|---|---|---|

60-view number | 0.942 | 0.894 | 0.570 |

30-view number | 0.938 | 0.822 | 0.490 |

Quantitative assessment CC of Shepp-Logan phantom images

CC | CS-based Iteration | ART | FBP |
---|---|---|---|

60-view number | 0.947 | 0.900 | 0.719 |

30-view number | 0.945 | 0.891 | 0.646 |

Quantitative assessment RMSE of **Polystyrene** phantom images

RMSE | CS-based Iteration | ART | FBP |
---|---|---|---|

60-view number | 7.47 | 17.90 | 12.16 |

30-view number | 9.52 | 19.44 | 17.37 |

Quantitative assessment UQI of **Polystyrene** phantom images

UQI | CS-based Iteration | ART | FBP |
---|---|---|---|

60-view number | 0.897 | 0.648 | 0.797 |

30-view number | 0.817 | 0.584 | 0.613 |

Quantitative assessment CC of **Polystyrene** phantom images

CC | CS-based Iteration | ART | FBP |
---|---|---|---|

60-view number | 0.900 | 0.668 | 0.803 |

30-view number | 0.831 | 0.601 | 0.635 |

## Discussion and Conclusions

In this article, a CS-base iterative algorithm reconstructing images from substantially reduced projection data was presented. Both the digital Shepp-Logan phantom and the real polystyrene sample experiment results show that the CS-based 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 CS-based 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.

## Declarations

### 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.

## Authors’ Affiliations

## References

- Kak AC, Slaney M:
*Principles of computerized tomographic imaging*. IEEE Press; 1988.Google Scholar - Herman GT:
*Image reconstruction from projections - implementation and applications*. Springer Verlag Press; 1975.Google Scholar - 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.MathSciNetView ArticleGoogle Scholar - Barrett JF, Keat N:
**Artifacts in CT: recognition and avoidance.***Radiographics*2004,**24:**1679–1691. 10.1148/rg.246045065View ArticleGoogle Scholar - Hsieh J:
*Computed Tomography - principles, designs, Artifacts, and Recent Advances*. Bellingham, WA: SPIE Press; 2003.Google Scholar - 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.View ArticleGoogle Scholar - Donoho DL:
**Compressed sensing.***IEEE Trans on Information Theory*2006,**52:**1289–1306.MathSciNetView ArticleGoogle Scholar - Candès EJ:
**Compressive sampling.**In*Proceedings of International Congress of Mathematicians*. Madrid, Spain; 2006:1–20.Google Scholar - Davis TJ, Gao D, Gureyev TE:
**Phase-contrast imaging of weakly absorbing materials using hard X-rays.***Nature*1995,**373:**595–598. 10.1038/373595a0View ArticleGoogle Scholar - Momose A, Fukuda J:
**Phase-contrast radiographs of nonstained rat cerebellar specimen.***Med phys*1995,**22:**375–379. 10.1118/1.597472View ArticleGoogle Scholar - Momose A, Takeda T, Itai Y:
**Phase-contrast X-ray computed-tomography for observing biological specimens and organic materials.***Review of Scientific Instruments*1995,**66:**1434–1436. 10.1063/1.1145931View ArticleGoogle Scholar - Snigirev A, Snigireva I:
**On the possibilities of x-ray phase contrast microimaging by coherent high-energy synchrotron radiation.***Review of Scientific Instruments*1995,**66:**5486–5492. 10.1063/1.1146073View ArticleGoogle Scholar - Wilkins SW, Gao D, Pogany A:
**Phase-contrast imaging using polychromatic hard x-rays.***Nature*1996,**384:**335–338. 10.1038/384335a0View ArticleGoogle Scholar - Yoneyama A, Momose A, Koyama I:
**Large-area phase-contrast X-ray imaging using a two-crystal X-ray interferometer.***J Synchrotron Rad*2002,**9:**277–281. 10.1107/S090904950201350XView ArticleGoogle Scholar - Momose A, Fukuda J:
**Demonstration of phase-contrast X-ray computed-tomography using an X-ray interferometer.***Nuclear Instruments and Methods in Physics Research*1995,**A352:**622–628.View ArticleGoogle Scholar - Dilmanian FA, Zhong Z, Ren B:
**Computed tomography of X-ray index of refraction using the diffraction enhanced imaging method.***Phys Med Biol*2000,**45:**933–946. 10.1088/0031-9155/45/4/309View ArticleGoogle Scholar - Nesterets YI, Gureyev TE, Wilkins SW:
**General reconstruction formulas for analyzer-based computed tomography.***Appl Phys Lett*2006,**89:**264103–264105. 10.1063/1.2423325View ArticleGoogle Scholar - Zachariasen WH:
*Theory of X-Ray Diffraction in Crystals*. New York: John Wiley and Sons Press; 1945.Google Scholar - Bellon PL, Lanzavecchia S:
**A direct Fourier method (DFM) for X-ray tomographic reconstructions and the accurate simulation of sinograms.***Int J Biomed Comput*1995,**38:**55–69. 10.1016/0020-7101(94)01033-WView ArticleGoogle Scholar - Sidky EY, Kao CM, Pan XC:
**Accurate image reconstruction from few-view limited-angle data in divergent-beam CT.***Journal of X-Ray Science and Technology*2006,**14:**119–139.Google Scholar - Candès EJ, Wakin MB:
**An introduction to compressive sampling.***IEEE Signal Processing Magazine*2008,**25:**21–30.View ArticleGoogle Scholar - Sidky EY, Pan XC:
**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 - Yang L, Lu Y, Wang G:
**Compressed sensing inspired image reconstruction from overlapped projections.***International Journal of Biomedical Imaging*2010,**2010:**1–8.View ArticleGoogle Scholar - 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/s11263-008-0168-yView ArticleGoogle Scholar - Herman GT, Meger LB:
**Algebraic reconstruction techniques can be made computationally efficient.***IEEE Transactions on Medical Imaging*1993,**12:**600–609. 10.1109/42.241889View ArticleGoogle Scholar - Bian JG, Siewerdsen JH, Han X, Sidky EY, Prince JL, Pelizzari CA, Pan XC:
**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.