Computational and Mathematical Methods in Medicine

Volume 2014 (2014), Article ID 862910, 8 pages

http://dx.doi.org/10.1155/2014/862910

Research Article

A CT Reconstruction Algorithm Based on L1/2 Regularization

Mianyi Chen, Deling Mi, Peng He, Luzhen Deng, and Biao Wei

The Key Lab of Optoelectronic Technology and Systems of the Education Ministry of China, Chongqing University, Chongqing 400044, China

Received 17 February 2014; Accepted 13 March 2014; Published 16 April 2014

Academic Editor: Peng Feng

Copyright © 2014 Mianyi Chen et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

Computed tomography (CT) reconstruction with low radiation dose is a significant research point in current medical CT field. Compressed sensing has shown great potential reconstruct high-quality CT images from few-view or sparse-view data. In this paper, we use the sparser L1/2 regularization operator to replace the traditional L₁ regularization and combine the Split Bregman method to reconstruct CT images, which has good unbiasedness and can accelerate iterative convergence. In the reconstruction experiments with simulation and real projection data, we analyze the quality of reconstructed images using different reconstruction methods in different projection angles and iteration numbers. Compared with algebraic reconstruction technique (ART) and total variance (TV) based approaches, the proposed reconstruction algorithm can not only get better images with higher quality from few-view data but also need less iteration numbers.

1. Introduction

Since computed tomography (CT) technique was born in 1973, CT has been widely applied in medical diagnose, industrial nondestructive detection, and so forth [1]. In medical CT field, CT reconstruction with low radiation dose is a significant research problem, which needs the reconstruction of high-quality CT images from few-view or sparse-view data. Recently, compressed sensing (CS) [2] theory has been applied in CT image reconstruction; it is possible to reconstruct high-quality images from few-view data under the frame of CS. In CS theory, sparse signal is often regularized with regularization; L₀ regularization is the sparest and most ideal regularization norm ( regularization denotes the number of nonzero signal elements). However, L₀ regularization is susceptible to noise interference and it is difficult to solve equations, and regularization is usually used as the regularization operator. Theoretically, if regularization is closer to regularization could get higher-quality CT images in CT reconstruction.

Recently, Xu et al. proposed regularization [3] and a soft thresholding algorithm [4], and Xu and Wang proposed a hybrid soft thresholding algorithm [5]. regularization has some theoretical properties, such as unbiasedness, sparse regularization, and Oracle. It is easier to produce sparser solution compared with the currentregular operator used in CT reconstruction. Theoretically, CT reconstruction method based on regularization will reconstruct higher-quality CT images with few-view data.

On solving regularization problem, Goldstein and Osher proposed Split Bregman method [6] which is derived from Bregman iterative algorithm [7]. Split Bregman method uses an intermediate variable to split regularization and regularization into two equations; regularization equation can be solved by gradient descent method, and regularization equation can be solved by thresholding algorithm. Split Bregman method can accelerate iterative convergence and produce better results. Based on Split Bregman method, Vandeghinste et al. proposed Split Bregman-based sparse-view CT reconstruction [8] and iterative CT reconstruction using shearlet-based regularization [9]. Chu et al. proposed multienergy CT reconstruction based on low rank and sparsity with the Split Bregman method (MLRSS) [10]. Chang et al. proposed a few-view reweighted sparsity hunting (FRESH) method for CT image reconstruction [11].

In this paper, we propose a CT reconstruction algorithm based on regularization and Split Bregman method. In the following section, regularization, Split Bregman method, and the proposed algorithm will be introduced. In the third section, we will use the proposed algorithm to analyze numerical phantom and real projection data. In the last section, we will discuss relevant issues and conclude the paper.

2. Materials and Methods

2.1. L1/2 Regularization

Generally, CT reconstruction algorithm can be divided into analytic reconstruction algorithm and iterative reconstruction algorithm; the current typical analytic reconstruction algorithm is filter back-projection (FBP), and iterative reconstruction algorithm contains algebraic reconstruction technique (ART) [12], simultaneous algebraic reconstruction technique (SART) [13], expectation-maximization (EM) [14], and so forth. Mathematical model [15] of CT image reconstruction can be expressed as

where is the projection matrix, is the projection data, and is the reconstruction image.

For sparse-view data, it is difficult to reconstruct high-quality images using the conventional CT image reconstruction algorithms, especially for analytic reconstruction algorithms which require high completeness of data. Meanwhile, there are also some artifacts in the reconstruction images using the conventional iterative reconstruction algorithms. In 2006, Donoho put forward the compressed sensing (CS) theory [2]; its main idea is that most of signals are sparse in the proper orthogonal transform domain, which means most of signal transformation coefficients are close to zero or equal to zero in orthogonal transformation, such as gradient transformation [16] and shearlet transformation [9]. In CS theory, an image can be reconstructed from a rather limited amount of data as long as an underlying image can be sparsely represented in an appropriate domain and determined from these data. CT sampling signal is a typical sparse signal; the x-ray attenuation coefficients of some regions of tested objects (i.e., human body) are similar or equal. Thus, CT reconstruction approaches based on compressed sensing can reconstruct high-quality CT images from sparse-view data.

CT reconstruction problem can be converted to a constrained optimization problem

where is the regularization function, usually denoted by norm of wavelet, gradient, and so forth. In order to simplify (2), we can use the penalty function method in the optimization method to convert it into an unconstraint optimization problem

where is the weight coefficient, can be denoted by total variation (TV) which is a research hotspot in current CT research field. TV method has been widely used in sparse-view and limited angle CT reconstruction [7].

In compressed sensing theory, norm is the most ideal regularization norm, but it is difficult to solve equations with norm, and regularization is easily interfered by noise in CT reconstruction, so norm is commonly replaced by norm. Theoretically, using a regularization norm which is closer to norm will reconstruct higher-quality CT images. The definitions of , , and regularization norm are

where and denotes the nonzero elements number in matrix .

As shown in Figure 1, regularization norm is closer to L₀ regularization norm than L₁ regularization norm.

Figure 1: , L₁, and regularization norm. Solid line and dotted line represent and regularization norm, respectively; imaginary point line denotes regularization norm.

Xu et al. proposed a regularization fast algorithm [4]; it can be expressed as

where is the regularization coefficient.

If is the optimal solution of (5), then the optimal condition of (5) will be denoted by

where is the gradient of ; multiplying coefficient and adding at the both sides of (6), then we have

The definition of operator is

Then the optimal solution can be represented as

where the operator is

Please see [4] for more proof details.

2.2. Split Bregman Method

In order to solve (3), Goldstein and Osher proposed Split Bregman method [6], using an intermediate variable to split L₁ regularization and regularization into two equations; L₂ regularization equation can be solved by gradient descent method and regularization equation can be solved by thresholding algorithm. Then the unconstrained optimal problem of (3) can be converted into

where is the sparse transform and the common used sparse transform contains gradient, wavelet, shearlet, and so forth.

Using an intermediate variable , (11) can be converted into

where is the coefficient. Then (12) can be converted into two unconstrained optimal Bregman problems; it can be expressed as

Equation (13) can split into two equations:

There are several advantages of Split Bregman method. Firstly, Split Bregman method can accelerate iterative convergence and calculate better results. Secondly, Split Bregman method can be widely used in CT reconstruction; it can not only solve regularization problem but also solve other regularization problems.

2.3. CT Reconstruction Algorithm Based on L1/2 Regularization

According to aforementioned methods, we propose a CT reconstruction algorithm based on regularization, where norm is used as the regularization norm and gradient as the sparse conversion; then (11) can be expressed as

Combine with Split Bregman method to solve (17) as follows.

Step 1. One has

Step 2. One has

Step 3. One has

Step 2 can be solved by the method in Section 2.1 and Step 3 can be solved directly. To solve Step 1, we use gradient descent method

Equation (21) is derivable and derivation of regularization as follows:

To derivate (21), we have

Then

where represents iteration numbers and parameter can be acquired by the following equation:

3. Experimental Study

3.1. Numerical Simulation

In this section, we study the ART algorithm, TV based ART algorithm (ART-TV), and regularization based Split Bregman method (SpBr-) and analyze the reconstructed images. In this paper, we test Shepp-Logan phantom as shown in Figure 2, and the size of phantom image is 256 × 256. We assume that the CT system was viewed as a typical parallel-beam geometry, and the scanning range was from 0° to 180° with a angular increment; projection angles can be indicated as

where is the number of projection angles.

Figure 2: Shepp-Logan phantom.

We will compare the reconstruction results from noise-free and noise data and projection numbers . In the simulation, we add 0.01% Gaussian noise to noise-free projection data, and iteration number for every reconstruction algorithm is 50; the reconstruction results are shown in Figure 3.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 3: The reconstructed images using three different reconstruction algorithms from the noise-free and noise data. (a)–(c) Reconstructed images from noise-free data: (a) reconstructed image using ART method, (b) reconstructed image using ART-TV method, and (c) reconstructed image using SpBr- method; (d)–(f) reconstructed images from noise data: (d) reconstructed image using ART method, (e) reconstructed image using ART-TV method, and (f) reconstructed image using SpBr- method.

From Figures 3 and 4, we can see that the reconstructed images using ART and ART-TV methods contain a lot of noise and artifacts from noise-free and noise data, while the reconstructed images using SpBr- method include less noise and artifacts and have clearer edges.

(a)

(b)

(c)

(d)

Figure 4: The comparison between reconstructed images using three different reconstruction algorithms and original Shepp-Logan image. (a) The profiles of line 128 in reconstructed images using ART and ART-TV methods from noise-free data and original Shepp-Logan image, (b) the profiles of line 128 in reconstructed images using SpBr- methods from noise-free data and original Shepp-Logan image, (c) the profiles of line 128 in reconstructed images using ART and ART-TV methods from noise data and original Shepp-Logan image, and (d) the profiles of line 128 in reconstructed images using SpBr- methods from noise data and original Shepp-Logan image.

In order to evaluate the quality of reconstructed images, we use root mean square errors (RMSE) as the evaluation index. The definition of RMSE is

where and are the reconstructed image and original image, respectively, and the image size is .

As shown in Table 1, the RMSE of reconstructed image using SpBr- method is much smaller than reconstructed images using ART and ART-TV methods, which means SpBr- can reconstruct higher-quality images.

Table 1: The RMSE of reconstructed images using three different algorithms from noise-free and noise data.

Methods
ART
ART-TV
SpBr-L1/2

Noise-free
0.0305
0.0104
0.0044

Noise
0.0388
0.0274
0.0102

As shown in Figure 5, the SpBr- algorithm can reconstruct high-quality images at less iteration numbers, which indicates that SpBr- algorithm can accelerate iterative convergence. From Figure 6, we can see that SpBr- algorithm can reconstruct high-quality images at less projection numbers.

(a)

(b)

Figure 5: The RMSE line of reconstructed images with different reconstruction algorithms at 60 projection angles and different iteration numbers, and the iteration numbers range from 1 to 100. (a) The RMSE of reconstructed images from noise-free projection data and (b) the RMSE of reconstructed images from noise projection data.

(a)

(b)

Figure 6: The RMSE line of reconstructed images with different reconstruction algorithms at 50 iteration numbers and different projection angles; the projection angles range from 1 to 180. (a) The RMSE of reconstructed images from noise-free projection data and (b) the RMSE of reconstructed images from noise projection data.

3.2. Real Data Study

In this section, we reconstruct oral images using three different algorithms with real projection data, where projection numbers are 90 and iteration numbers are 50 while the original projection numbers are 360. And as shown in Figure 7(a), we evaluate three images with reconstructed image using ART with original projection data as a standard image. The reconstructed images and RMSE of three reconstructed images are shown in Figure 7 and Table 2.

Table 2: The RMSE of reconstructed images using three different algorithms with real projection data.

Methods
ART
ART-TV
SpBr-L1/2

RMSE
0.0270
0.0256
0.0236

(a)

(b)

(c)

(d)

Figure 7: The reconstructed images using three algorithms from real projection data; iteration numbers are 50. (a) Reconstructed image using ART with original projection data, (b) the reconstructed image using ART method, (c) the reconstructed image using ART-TV method, and (d) the reconstructed image using SpBr- method.

As shown in Figure 7, the reconstructed images using ART and ART-TV method have more noise and artifacts, while the reconstructed image using SpBr- method has less noise and artifacts. And from Table 2, the RMSE of reconstructed image using SpBr- method is smaller than that of reconstructed images using ART and ART-TV methods, which means SpBr- can reconstruct higher-quality images with clearer details. From the comparison of reconstructed images using different algorithms, we can see that SpBr- can reconstruct high-quality images at less projection angles and less iteration numbers.

4. Discussions and Conclusion

There are several issues worth further discussion in the reconstruction study. Firstly, the thresholding algorithm was not applied to solve the regularization problem. There are two reasons. First, if the thresholding algorithm is applied to solve the regularization problem, the edge and details of reconstructed images will not be clear, while the SpBr- without thresholding algorithm reconstructs images with clear edges. Then, if the threshold value is settled to solve regularization problem, the reconstruction speed will be reduced a lot. Secondly, there are some artifacts in the reconstructed image using SpBr-method from noise data. The reason is that regularization is closer to regularization, and the ability of denoising of regularization is bad; therefore the ability of denoising of regularization is not good enough. Thirdly, in the reconstruction from real projection data, the reconstructed image using ART-TV method lost some details, while details of reconstructed image using SpBr- method are much clearer.

In the further research, we will try to use SpBr- algorithm in interior CT and study the region of interest (ROI) reconstruction, which will reduce radiation dose as much as possible.

In conclusion, we proposed a CT reconstruction algorithm based on regularization; the reconstructed results demonstrate that SpBr- can reconstruct high-quality images from few-view data. In the urgent demand of radiation reduction, SpBr- algorithm will have great potential in clinical application.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

This work was partially supported by the National Natural Science Foundation of China Grants (61201346 and 61171157) and the Fundamental Research Funds for the Central Universities (106112013CDJZR120020). The authors would like to thank real projection data provider Dr. Zhichao Li from the Third Military Medical University, Chongqing, China.

References

G. Wang, H. Yu, and B. de Man, “An outlook on x-ray CT research and development,” Medical Physics, vol. 35, no. 3, pp. 1051–1064, 2008.

D. L. Donoho, “Compressed sensing,” IEEE Transactions on Information Theory, vol. 52, no. 4, pp. 1289–1306, 2006.

Z. B. Xu, H. Zhang, Y. Wang, X. Y. Chang, and Y. Liang, “L1/2 regularization,” Science China Information Sciences, vol. 53, no. 6, pp. 1159–1169, 2010.

Z. B. Xu, X. Y. Chang, F. M. Xu, and H. Zhang, “L1/2regularization: a thresholding representation theory and a fast solver,” IEEE Transactions on Neural Networks and Learning Systems, vol. 23, no. 7, pp. 1013–1027, 2012.

F. M. Xu and S. H. Wang, “A hybrid simulated annealing threshoding algorithm for compressed sensing,” Signal Processing, vol. 93, no. 6, pp. 1577–1585, 2013.

T. Goldstein and S. Osher, “The split Bregman method for L1-regularized problems,” SIAM Journal on Imaging Sciences, vol. 2, no. 2, pp. 323–343, 2009.

L. M. Bregman, “The relaxation method of finding the common point of convex sets and its application to the solution of problems in convex programming,” USSR Computational Mathematics and Mathematical Physics, vol. 7, no. 3, pp. 200–217, 1967.

B. Vandeghinste, B. Goossens, J. de Beenhouwer et al., “Split-Bregman-based sparse-view CT reconstruction,” in Proceedings of the 11th International meeting on Fully Three-Dimensional Image Reconstruction in Radiology and Nuclear Medicine (Fully 3D '11), pp. 431–434, 2011.

B. Vandeghinste, B. Goossens, R. van Holen et al., “Iterative CT reconstruction using shearlet-based regularization,” in Proceedings of the International Society for Optics and Photonics in SPIE Medical Imaging, pp. 83133I–83133I, February 2012.

J. Chu, L. Li, Z. Chen, G. Wang, and H. Gao, “Multi-energy CT reconstruction based on Low Rank and Sparsity with the Split-Bregman Method (MLRSS),” in Proceedings of the Nuclear Science Symposium and Medical Imaging Conference (NSS/MIC '12), pp. 2411–2414, 2012.

M. Chang, L. Li, Z. Chen, Y. Xiao, L. Zhang, and G. Wang, “A few-view reweighted sparsity hunting (FRESH) method for CT image reconstruction,” Journal of X-Ray Science and Technology, vol. 21, no. 2, pp. 161–176, 2013.

R. Gordon, R. Bender, and G. T. Herman, “Algebraic Reconstruction Techniques (ART) for three-dimensional electron microscopy and X-ray photography,” Journal of Theoretical Biology, vol. 29, no. 3, pp. 471–481, 1970.

A. H. Andersen and A. C. Kak, “Simultaneous algebraic reconstruction technique (SART): a superior implementation of the art algorithm,” Ultrasonic Imaging, vol. 6, no. 1, pp. 81–94, 1984.

L. A. Shepp and Y. Vardi, “Maximum likelihood reconstruction for emission tomography,” Medical Imaging, vol. 1, no. 2, pp. 113–122, 1982.

G. Wang and M. Jiang, “Ordered-subset simultaneous algebraic reconstruction techniques (OS-SART),” Journal of X-Ray Science and Technology, vol. 12, no. 3, pp. 169–177, 2004.

E. Y. Sidky, C.-M. Kao, and X. Pan, “Accurate image reconstruction from few-views and limited-angle data in divergent-beam CT,” Journal of X-Ray Science and Technology, vol. 14, no. 2, pp. 119–139, 2006.