Combined Energy Minimization for Image Reconstruction from Few Views

Abstract

Reconstruction from few views is an important problem in medical imaging and applied mathematics. In this paper, a combined energy minimization is proposed for image reconstruction. 𝑙₂ energy of the image gradient is introduced in the lower density region, and it can accelerate the reconstruction speed and improve the results. Total variation of the image is introduced in the higher density region, and the image features can be preserved well. Nonlinear conjugate gradient method is introduced to solve the problem. The efficiency and accuracy of our method are shown in several numerical experiments.

1. Introduction

Computed tomography (CT) is one of the most important advance in diagnostic radiology in recent decades. CT uses multiple X-ray images to build up cross-sectional and 3D pictures of structures inside the human body which enable doctors to view internal organs with unprecedented precision. However, the use of ionizing radiation in CT may induce cancer in the exposed individual after a latent period [1–3]. Cancer induction by ionizing radiation is a probabilistic process. Reduction of radiation dose used in CT will therefore lead to a reduction in the number of induced cancer cases.

Some ways can be used to reduce the radiation dose from CT such as decreasing intensity of X-ray beam, handling scattered radiation, restricting exposure area. Reducing the X-ray exposing time is a simple one. Here we focus on the low-dose X-ray imaging strategy that only a limited number of projection images are taken for the reconstruction, which is called limited-view reconstruction [4–6]. Image reconstruction from few views would enable rapid scanning with a reduced X-ray dose delivered to the patient. But it is possibe to build up blurring artifacts and lose critical spatial resolution in the reconstructed image. So the trade-off needs to be carefully examined for intended use. However, it is beyond the scope of this paper. In a word, we will consider image reconstruction from projection data at few views.

CT reconstruction methods can roughly be categorized as analytic reconstruction methods and iterative reconstruction methods. The analytic reconstruction methods, such as filtered back-projection (FBP) methods [7, 8], require sufficient projection data with low noise level. As the limited-view reconstruction is considered, the analytic methods may induce more noise and produce significant artifacts. The iterative reconstruction methods, such as the algebraic reconstruction technique (ART) [9, 10], require less data than FBP methods and are more robust to the effects of noise, but need much more computation.

Recently, the minimization of the image total variation (TV) has been introduced to divergent-beam CT and an new iterative image reconstruction algorithm was presented [11]. Lustig et al. applied compressed sensing theory to rapid magnetic resonance imaging [12]. Many approaches have been presented based on this [13–16].

In this paper, a novel image reconstruction model is proposed. The image total variation and the 𝑙₂ energy of the image gradient are combined to a new energy functional. Then the functional is introduced to the constrained optimization problem for the reconstruction from 2D parallel-beam data at few views. Our algorithm are performed with various insufficient data problems in fan-beam CT and the numerical results show the efficiency and accuracy of proposed method. The algorithm can be generalized to fan-beam CT and cone-beam CT as well as other tomographic imaging modalities.

2. Reconstruction from Few Views and Total Variation Minimization

There are many approaches about tomographic reconstruction from limited views projection data [17–19]. Algebraic reconstruction technique (ART) [9, 20, 21] and the expectation-maximization (EM) algorithm [22, 23] have been widely used in this field. As the image is discretized on the grids, each projection is regarded as a linear equation of the discrete density distribution. Then a system of simultaneous equations can be obtained and ART tends to solve it via iterative method. ART algorithm can find the image that is consistent with the projection data and the sum-of-squares of the density values is minimized. The EM algorithm applies to positive integral equations, seeking to minimize the Kullback-Liebler distance between the measured data and the projection of the estimated image [11].

However, it is known that the ray does harm to human body and abundant irradiation may lead to cancer [17, 24]. So researchers begin to study the tomographic reconstruction with projection data as little as possible.

2.1. Reconstruction from Few Views

Tomographic reconstruction from few views projection data is an efficient way to reduce the harm caused by ray irradiation, and there are some approaches about it [11, 12].

As the gray image to be reconstructed can be denoted by below 𝑢𝑢=^1,1𝑢_1,2…𝑢_1,𝑐𝑢_2,1𝑢_2,2…𝑢_2,𝑐𝑢…………_𝑟,1𝑢_𝑟,2…𝑢_𝑟,𝑐. (2.1) Here 𝑟 and 𝑐 mean the size of image.

The projection can be denoted as the following equations: 𝑀⋅𝑓=𝑔. (2.2) Here 𝑓 is the vector form rearranged form 𝑢. 𝑔∶=(𝑔₁,𝑔₂,…,𝑔_𝑘)^𝑇 means the projection data. More exactly, 𝑘 is the product of the number of views and the number of detector’s pixels. ℳ∶=(𝑀₁,𝑀₂,…,𝑀_𝑘)^𝑇 is the projection matrix which can be precomputed. 𝑀_𝑖 (𝑖=1,2,…,𝑘) is the same size as 𝑓. The reconstruction procedure equals to solve (2.2).

Unfortunately, this equations are indeterminate if the reconstruction was based on few views. In other words, the number of the equations are less than the number of variables (𝑘<𝑟⋅𝑐). In practice, it is more often that 𝑘≪𝑟⋅𝑐. From the linear algebraic theory, the solution is not unique and the traditional methods cannot be applied. In fact, even if 𝑘≥𝑟⋅𝑐 is satisfied, it is still compromised to deal with the consistency of the projection data and lead to artifacts in the reconstructed image.

The ART can be applied to solve this equation and it means to solve the following problem: argmin𝑓₂subjectto𝑀⋅𝑓=𝑔. (2.3) Here ||𝑓||₂ means the 𝑙₂ norm of 𝑓. Because of the serious insufficiency of the projection data, ART algorithm can hardly provide satisfactory result. The same as to EM algorithm. So some other models should be discussed.

2.2. Total Variation Minimization

The total variation (TV) was first introduced by Rudin et al. [25], and it can be utilized in image processing for images denoising while edges preserved. Candès et al. applied it to image reconstruction with insufficient parallel-beam data [26]. More exactly, they considered the following problem: argmin𝑓_TV=∇𝑓₁subjectto𝑀⋅𝑓=𝑔. (2.4) Here ||∇𝑓||₁ means the 𝑙₁ norm of |∇𝑓|=(|∇𝑓_𝑡|) which is the rearranged form of matrix |∇𝑢|=(|∇𝑢_𝑖,𝑗|). 𝑖,𝑗,𝑡 satisfy 𝑖+(𝑗−1)∗𝑟=𝑡. |∇𝑢_𝑖,𝑗| can be computed as 12𝑢_𝑖+1,𝑗−𝑢_{𝑖−1,𝑗}²+𝑢_𝑖,𝑗+1−𝑢_{𝑖,𝑗−1}². (2.5) Based on these, Sidky et al. developed an iterative image reconstruction algorithm for fan-beam CT in [11].

The TV minimization can efficiently reduce errors and preserve features in the image reconstruction. In next section, we will concentrate on developing a new model to improve the convergence speed and reduce errors based on this TV model.

3. Minimization of a Combined Energy for Image Reconstruction

It is known that the convergence speed will be enhanced when the 𝑙₂ norm of image gradient is considered as shown in following: argmin∇𝑓²₂subjectto𝑀⋅𝑓=𝑔. (3.1) But this 𝑙₂ result can also blur the image features. So some combined energies can be considered.

3.1. A Combined Energy of Image

The natural idea is to combine the 𝑙₂ norm and TV directly. The combined energy can be denoted by 𝜀_CTV𝑓=∇𝑓²₂+𝛼𝑓_TV. (3.2) However, the new result cannot be improved much more than the TV result though the convergence speed may be accelerated in some sense. In fact, the new result is a weighted sum of the 𝑙₂ result and TV result.

To improve the TV result of image denoising, Chambolle and Lions [27] proposed a combined functional (CL energy) 𝐸_CL(𝑢)=12_{|∇𝑢|<𝛽}|∇𝑢|²𝑑𝑥+𝛽_{|∇𝑢|≥𝛽}||−𝛽∇𝑢2𝑑𝑥=_Ω𝐹(|∇𝑢|)𝑑𝑥, (3.3) where 𝐹(𝑠)=12𝑠²1,𝑠<𝛽,𝛽𝑠−2𝛽²,𝑠≥𝛽. (3.4) Here 𝛽 is a fixed positive number and it is a threshold of |∇𝑢|. In some way, it means an approximation of the critical value which can be used to distinguish image features and noise.

3.2. The Reconstruction Model

Based on the discrete form of the CL energy and the rearranged vector 𝑓, we can get 𝜀_CL𝑓=12_{𝑓𝑡∈𝑇−}∇𝑓_𝑡²+𝛽_{𝑓𝑡∈𝑇+}∇𝑓_𝑡−𝛽2. (3.5) Here 𝑇(𝑓)⁻∶={𝑡∣|∇𝑓_𝑡|<𝛽}, 𝑇(𝑓)⁺∶={𝑡∣|∇𝑓_𝑡|≥𝛽}. It can be found that the 𝑙² energy of image gradient is considered in the noise part 𝑇(𝑓)⁻ while the TV energy is computed in features part 𝑇(𝑓)⁺. Then the new model for image reconstruction can be denoted as follows: argmin𝜀_CL𝑓subjectto𝑀⋅𝑓=𝑔. (3.6)

With the Lagrange method applied, this constraint optimization problem can be rewritten as an unconstraint optimization problem of following combined Chambolle-Lions (CCL) energy: argmin𝜀_CCL𝑓=𝑀⋅𝑓−𝑔²₂+𝜆⋅𝜀_CL𝑓. (3.7)

3.3. Conjugate Gradient Descend Algorithm

The gradient of 𝜀_CCL(𝑓) can be computed as ∇𝜀_CCL𝑓=2𝑀^∗𝑀⋅𝑓−𝑔+𝜆⋅∇𝜀_CL𝑓, (3.8) where ∇𝜀_CL𝑓=−∇⋅∇𝑓,∇𝑓<𝛽,−𝛽∇⋅∇𝑓,∇𝑓∇𝑓≥𝛽. (3.9) In practice, the parameter 𝜆 is set to be 10⁻².

Set a initial value 𝑓₀, the conjugate gradient descend algorithm can be given as Algorithm 1. There the time step is denoted by 𝜏.

% Initialization maxGrad = 10−3; maxIter = 100; maxTau = 10−3; 𝑠 = 4; 𝑘 = 0; 𝜏0 = 1; 𝑔0 = ∇𝜀CCL(𝑓0); Δ𝑓0 = −𝑔0; % Iterations while (‖𝑔𝑘‖2 > maxGrad and 𝑘 < maxIter and 𝜏 > maxTau){   %Linear search   min𝐸 = 𝜀CCL(𝑓𝑘); opt = 0; 𝑖 = −𝑠;   while (𝑖 <= s){     𝜏=2−𝑖𝜏0;     if (𝜀CCL(𝑓𝑘 + 𝜏Δ𝑓𝑘) < min𝐸){min𝐸 = 𝜀CCL(𝑓𝑘 + 𝜏Δ𝑓𝑘); opt = 𝜏; }     𝑖=𝑖+1; }   𝜏0 = opt;   𝑓𝑘+1 = 𝑓𝑘+𝜏0Δ𝑓𝑘; 𝑔𝑘+1 = ∇𝜀CCL (𝑓𝑘+1); 𝛾 = ‖𝑔𝑘+1‖2‖Δ𝑓𝑘‖2;   Δ𝑓𝑘+1 = −𝑔𝑘+1 + 𝛾Δ𝑓𝑘;   𝑘 = 𝑘 + 1; }

Algorithm 1 

Iteration algorithm for CL minimization reconstruction.

From the experiments results, we will find some advantages of the proposed model. These are chiefly due to the different optimization problems and it is related with the algorithm. More exactly, the first term in (3.5) can help to enhance the convergence speed and reduce some artifacts in the smooth region. But the efficiency role of this 𝑙² energy of image gradient is depended on the well define of 𝛽. It is related with the characters of the image to be reconstructed. In our experiments, it is set to be 0.01 times the range of phantoms. The advanced researches about this will be approached in our next work. The general metric peak signal to noise ratio (PSNR) is introduced to evaluate the results and it can be computed as PSNR=10⋅log255²/MSE,log10 (3.10) where MSE is the mean square error of the gray image.

4. Numerical Experiments

Example 4.1. Reconstruction of Shepp-Logan phantom from 72 views.
The true image is taken to be the Shepp-Logan image shown in Figure 1(a) discretized on a 256×256 pixel grid. The computational parameters are set as shown in the algorithm, the same to the following experiments. The reconstruction from 72 views is completed after 31 iterations. Figures 1(b) and 1(c) show the ART result and TV result while Figure 1(d) shows our CL result. It can be found that the ART result is enhanced by TV much more. Many artifacts have been removed or slighted. Figures 1(e) and 1(f) show the gray distributions of row 128 and column 128. There are few differences between the reconstructed horizontal gray and the real one. It is similar to the reconstructed vertical gray. The evolutions of PSNR and 𝜏 are shown in Figures 2(a) and 2(b). The PSNR has been improved from 46.4040 (TV result) to 50.5664 (CL result). Though the TV result is almost accurate, our CL result improves it significantly.

(a) Original

(b) Art result

(c) TV result

(d) CL result

(e) Horizontal gray

(f) Vertical gray

Figure 1 

Reconstruction of Shepp-Logan phantom from 72 views.

(a) Evolution of PSNR

(b) Evolution of 𝜏

Figure 2 

Result analysis of Example 4.1.

Example 4.2. Reconstruction of Shepp-Logan phantom from 24 views.
The true image is still taken to be the same Shepp-Logan image as Example 4.1 while the reconstruction views are reduced acutely from 72 to 24. This reconstruction is completed after 100 iterations. Figures 3(b), 3(c), and 3(d) show the ART result, TV result and our CL result. The gray distributions of row 128 and column 128 are shown in Figures 3(e) and 3(f). Figures 4(a) and 4(b) show the evolutions of PSNR and 𝜏. The PSNR has been improved from 21.3451 (TV result) to 34.4123 (CL result). It can be found that there are a lot of serious line artifacts and block artifacts in TV result. Most of these artifacts has been eliminated and the rest has been slighted significantly in our CL result.

(a) Original

(b) Art result

(c) TV result

(d) CL result

(e) Horizontal gray

(f) Vertical gray

Figure 3 

Reconstruction of Shepp-Logan phantom from 24 views.

(a) Evolution of PSNR

(b) Evolution of 𝜏

Figure 4 

Result analysis of Example 4.2.

Example 4.3. Reconstruction of fruits image from 72 views.
A fruits image with size 256×256 is taken to be the true one. The reconstruction views are set to be 72. It costs 25 iterations to finish the reconstruction. Figures 5(b), 5(c), and 5(d) show the ART result, TV result, and our CL result. The gray distributions of row 128 and column 128 are shown in Figures 5(e) and 5(f). Figure 6 shows the evolutions of PSNR and 𝜏. The PSNR has been improved from 31.8559 (TV result) to 32.2934 (CL result). There are few differences between TV result and our CL result though CL result is more satisfactory than TV result in numerical value.

(a) Original

(b) Art result

(c) TV result

(d) CL result

(e) Horizontal gray

(f) Vertical gray

Figure 5 

Reconstruction of fruits from 72 views.

(a) Evolution of PSNR

(b) Evolution of 𝜏

Figure 6 

Result analysis of Example 4.3.

Example 4.4. Reconstruction of fruits image from 30 views.
The same fruits image as Example 4.3 is taken to be the true one while the views are reduced from 72 to 30. This reconstruction is completed after 59 iterations. ART result, TV result, and our CL result are shown in Figures 7(b), 7(c), and 7(d). The gray distributions of row 128 and column 128 are shown in Figures 7(e) and 7(f). The evolutions of PSNR and 𝜏 are shown in Figures 8(a) and 8(b). The PSNR has been improved from 25.9273 (TV result) to 27.3870 (CL result). It can be found that there are many artifacts in TV result and a lot of them have been eliminated or slighted in our CL result.

(a) Original

(b) Art result

(c) TV result

(d) CL result

(e) Horizontal gray

(f) Vertical gray

Figure 7 

Reconstruction of fruits from 30 views.

(a) Evolution of PSNR

(b) Evolution of 𝜏

Figure 8 

Result analysis of Example 4.4.

Example 4.5. Reconstruction of a synopsis phantom from 72 views.
A synopsis phantom with size 256×256 is taken to be the true image. The reconstruction views are set to be 72. After 24 iterations, the reconstruction is finished. Figures 9(b), 9(c), and 9(d) show the ART result, TV result, and our CL result. The gray distributions of row 128 and column 128 are shown in Figures 9(e) and 9(f). Figure 10 shows the evolutions of PSNR and 𝜏. The PSNR has been improved from 48.5852 (TV result) to 49.6270 (CL result). It can be found that CL result is more satisfactory than TV result in numerical value, but there are few differences between TV result and our CL result from vision terms.

(a) Original

(b) Art result

(c) TV result

(d) CL result

(e) Horizontal gray

(f) Vertical gray

Figure 9 

Reconstruction of a synopsis phantom from 72 views.

(a) Evolution of PSNR

(b) Evolution of 𝜏

Figure 10 

Result analysis of Example 4.5.

Example 4.6. Reconstruction of a synopsis phantom from 20 views.
The same synopsis phantom as Example 4.5 is taken to be the true one while the views are reduced from 72 to 20. This reconstruction is completed after 68 iterations. ART result, TV result, and our CL result are shown in Figures 11(c), 11(b), and 11(d). The gray distributions of row 128 and column 128 are shown in Figures 11(e) and 11(f). The evolutions of PSNR and 𝜏 are shown in Figures 12(a) and 12(b). The PSNR has been improved from 31.7194 (TV result) to 36.9629 (CL result). It can be found that there are more artifacts in TV result than in our CL result.

(a) Original

(b) Art result

(c) TV result

(d) CL result

(e) Horizontal gray

(f) Vertical gray

Figure 11 

Reconstruction of a synopsis phantom from 20 views.

(a) Evolution of PSNR

(b) Evolution of 𝜏

Figure 12 

Result analysis of Example 4.6.

5. Conclusion

In this paper, a novel model for image reconstruction from few views in parallel-beam data is proposed. First, the 𝑙₂ energy of the image gradient and the total variation of the image are combined to the CL energy. The 𝑙₂ energy is applied in the lower density region, and it can accelerate the reconstruction speed. The total variation is applied in the higher density region, and it can preserve the image features well. Contributed to the lagrange method and nonlinear conjugate gradient algorithm, the related optimization problem can be solved.

Acknowledgments

The authors would like to thank the anonymous reviewers for their constructive feedback and valuable input. Due thanks are for the supports to their program from the TI, the XILINX, and the Software School of Xidian University. This program is partially supported by NSFC (Grant nos. 61072105, 61007011) and by the Open Projects Program of National Laboratory of Pattern Recognition. The Project is also partially supported by Natural Science Basic Research Plan in Shaanxi Province of China (Program no. 2010JM8005) and Scientific Research Program Funded by Shaanxi Provincial Education Department (Program no. 11JK0504).