A CT Reconstruction Algorithm Based on L1/2 Regularization.

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 L1 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 L(0 ≤ p ≤ 1) regularization; L0 regularization is the sparest and most ideal regularization norm (L0 regularization denotes the number of nonzero signal elements). However, L0 regularization is susceptible to noise interference and it is difficult to solve equations, and L1 regularization is usually used as the regularization operator. Theoretically, if regularization is closer to L0 regularization could get higher-quality CT images in CT reconstruction.

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

On solving L1 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 L1 regularization and L2  regularization into two equations; L2  regularization equation can be solved by gradient descent method, and L1 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 L1/2 regularization and Split Bregman method. In the following section, L1/2 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 A = (a ) is the projection matrix, b = (b 1, b 2,…, b ) ∈ R is the projection data, and u = (u 1, u 2,…, u ) ∈ R 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 E(u) is the regularization function, usually denoted by L1 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, E(u) 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, L0 norm is the most ideal regularization norm, but it is difficult to solve equations with L0 norm, and L0 regularization is easily interfered by noise in CT reconstruction, so L0 norm is commonly replaced by L1 norm. Theoretically, using a regularization norm which is closer to L0 norm will reconstruct higher-quality CT images. The definitions of L0, L1, and L1/2 regularization norm are where U = (u 1, u 2,…, u ) and a denotes the nonzero elements number in matrix U.

As shown in Figure 1, L1/2 regularization norm is closer to L0 regularization norm than L1 regularization norm.

Figure 1

L0, L1, and L1/2 regularization norm. Solid line and dotted line represent L0 and L1 regularization norm, respectively; imaginary point line denotes L1/2 regularization norm.

Xu et al. proposed a L1/2 regularization fast algorithm [4]; it can be expressed as where γ is the regularization coefficient.

If u is the optimal solution of (5), then the optimal condition of (5) will be denoted by where ∇(||u||1/2 1/2) is the gradient of ||u||1/2 1/2; multiplying coefficient μ and adding u at the both sides of (6), then we have

The definition of operator is

Then the optimal solution u 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 L1 regularization and L2 regularization into two equations; L2 regularization equation can be solved by gradient descent method and L1 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 d = Φ(u), (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 L1 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 L1/2 regularization, where L1/2 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 L2 regularization as follows:

To derivate (21), we have

Then where k 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 L1/2 regularization based Split Bregman method (SpBr-L1/2) 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 N view 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 N view = 60. 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.

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-L1/2 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-L1/2 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-L1/2 method include less noise and artifacts and have clearer edges.

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-L1/2 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-L1/2 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 u and u* are the reconstructed image and original image, respectively, and the image size is M × N.

As shown in Table 1, the RMSE of reconstructed image using SpBr-L1/2 method is much smaller than reconstructed images using ART and ART-TV methods, which means SpBr-L1/2 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-L1/2 algorithm can reconstruct high-quality images at less iteration numbers, which indicates that SpBr-L1/2 algorithm can accelerate iterative convergence. From Figure 6, we can see that SpBr-L1/2 algorithm can reconstruct high-quality images at less projection numbers.

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.

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.

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-L1/2 method.

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

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-L1/2 method has less noise and artifacts. And from Table 2, the RMSE of reconstructed image using SpBr-L1/2 method is smaller than that of reconstructed images using ART and ART-TV methods, which means SpBr-L1/2 can reconstruct higher-quality images with clearer details. From the comparison of reconstructed images using different algorithms, we can see that SpBr-L1/2 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 L1/2 regularization problem. There are two reasons. First, if the thresholding algorithm is applied to solve the L1 regularization problem, the edge and details of reconstructed images will not be clear, while the SpBr-L1/2 without thresholding algorithm reconstructs images with clear edges. Then, if the threshold value is settled to solve L1/2 regularization problem, the reconstruction speed will be reduced a lot. Secondly, there are some artifacts in the reconstructed image using SpBr-L1/2  method from noise data. The reason is that L1/2 regularization is closer to L0 regularization, and the ability of denoising of L0 regularization is bad; therefore the ability of denoising of L1/2 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-L1/2 method are much clearer.

In the further research, we will try to use SpBr-L1/2 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 L1/2 regularization; the reconstructed results demonstrate that SpBr-L1/2 can reconstruct high-quality images from few-view data. In the urgent demand of radiation reduction, SpBr-L1/2 algorithm will have great potential in clinical application.