# Combined Energy Minimization for Image Reconstruction from Few Views

Authors: Wei Wei; Xiao-Lin Yang; Bin Zhou; Jun Feng; Pei-Yi Shen
Journal: Mathematical Problems in Engineering
 (2012)

---

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


---

## Body

## 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 thel2 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(2.1)u=(u1,1u1,2…u1,cu2,1u2,2…u2,c…………ur,1ur,2…ur,c).
Here r and c mean the size of image.The projection can be denoted as the following equations:(2.2)M·f=g.
Here f is the vector form rearranged form u. g:=(g1,g2,…,gk)T means the projection data. More exactly, k is the product of the number of views and the number of detector’s pixels. ℳ:=(M1,M2,…,Mk)T is the projection matrix which can be precomputed. Mi (i=1,2,…,k) is the same size as f. 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 (k<r·c). In practice, it is more often that k≪r·c. From the linear algebraic theory, the solution is not unique and the traditional methods cannot be applied. In fact, even if k≥r·c 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:(2.3)argmin∥f∥2subjecttoM·f=g.
Here ||f||2 means the l2 norm of f. 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:
(2.4)argmin∥f∥TV=∥∇f∥1subjecttoM·f=g.
Here ||∇f||1 means the l1 norm of |∇f|=(|∇ft|) which is the rearranged form of matrix |∇u|=(|∇ui,j|). i,j,t satisfy i+(j-1)*r=t. |∇ui,j| can be computed as
(2.5)12(ui+1,j-ui-1,j)2+(ui,j+1-ui,j-1)2.
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.


## 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(2.1)u=(u1,1u1,2…u1,cu2,1u2,2…u2,c…………ur,1ur,2…ur,c).
Here r and c mean the size of image.The projection can be denoted as the following equations:(2.2)M·f=g.
Here f is the vector form rearranged form u. g:=(g1,g2,…,gk)T means the projection data. More exactly, k is the product of the number of views and the number of detector’s pixels. ℳ:=(M1,M2,…,Mk)T is the projection matrix which can be precomputed. Mi (i=1,2,…,k) is the same size as f. 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 (k<r·c). In practice, it is more often that k≪r·c. From the linear algebraic theory, the solution is not unique and the traditional methods cannot be applied. In fact, even if k≥r·c 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:(2.3)argmin∥f∥2subjecttoM·f=g.
Here ||f||2 means the l2 norm of f. 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:
(2.4)argmin∥f∥TV=∥∇f∥1subjecttoM·f=g.
Here ||∇f||1 means the l1 norm of |∇f|=(|∇ft|) which is the rearranged form of matrix |∇u|=(|∇ui,j|). i,j,t satisfy i+(j-1)*r=t. |∇ui,j| can be computed as
(2.5)12(ui+1,j-ui-1,j)2+(ui,j+1-ui,j-1)2.
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 thel2 norm of image gradient is considered as shown in following:
(3.1)argmin∥∇f∥22subjecttoM·f=g.
But this l2 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 thel2 norm and TV directly. The combined energy can be denoted by
(3.2)εCTV(f)=∥∇f∥22+α∥f∥TV.
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 l2 result and TV result.To improve the TV result of image denoising, Chambolle and Lions [27] proposed a combined functional (CL energy)
(3.3)ECL(u)=12∫|∇u|<β‍|∇u|2dx+β∫|∇u|≥β‍|∇u|-β2dx=∫Ω‍F(|∇u|)dx,
where
(3.4)F(s)={12s2,s<β,βs-12β2,s≥β.
Here β is a fixed positive number and it is a threshold of |∇u|. 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 vectorf, we can get
(3.5)εCL(f)=12∑t∈T(f)-‍|∇ft|2+β∑t∈T(f)+‍(|∇ft|-β2).
Here T(f)-:={t∣|∇ft|<β}, T(f)+:={t∣|∇ft|≥β}. It can be found that the l2 energy of image gradient is considered in the noise part T(f)- while the TV energy is computed in features part T(f)+. Then the new model for image reconstruction can be denoted as follows:
(3.6)argminεCL(f)subjecttoM·f=g.With the Lagrange method applied, this constraint optimization problem can be rewritten as an unconstraint optimization problem of following combined Chambolle-Lions (CCL) energy:(3.7)argminεCCL(f)=∥M·f-g∥22+λ·εCL(f).


### 3.3. Conjugate Gradient Descend Algorithm

The gradient ofεCCL(f) can be computed as
(3.8)∇εCCL(f)=2M*(M·f-g)+λ·∇εCL(f),
where
(3.9)∇εCL(f)={-∇·(∇f),|∇f|<β,-β∇·(∇f|∇f|),|∇f|≥β.
In practice, the parameter λ is set to be 10-2.Set a initial valuef0, the conjugate gradient descend algorithm can be given as Algorithm 1. There the time step is denoted by τ.Algorithm 1:Iteration algorithm for CL minimization reconstruction.
% Initialization
maxGrad =10-3; maxIter = 100; maxTau = 10-3; s = 4; k = 0;
τ
0 = 1; g0 = ∇εCCL(f0); Δf0 = −g0;
% Iterations
while (∥gk∥2  > maxGrad and k < maxIter and τ > maxTau){
%Linear search
minE = εCCL(fk); opt = 0; i = −s;
while (i <= s){
τ=2-iτ0;
if (εCCL(fk + τΔfk) < minE){minE = εCCL(fk + τΔfk); opt = τ; }
i=i+1; }
τ0 = opt;
fk+1 = fk+τ0Δfk; gk+1 = ∇εCCL (fk+1); γ = ∥gk+1∥2∥Δfk∥2;
Δfk+1 = −gk+1 + γΔfk;
k = k + 1; }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 l2 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
(3.10)PSNR=10·log(2552/MSE)log10,
where MSE is the mean square error of the gray image.


## 3.1. A Combined Energy of Image

The natural idea is to combine thel2 norm and TV directly. The combined energy can be denoted by
(3.2)εCTV(f)=∥∇f∥22+α∥f∥TV.
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 l2 result and TV result.To improve the TV result of image denoising, Chambolle and Lions [27] proposed a combined functional (CL energy)
(3.3)ECL(u)=12∫|∇u|<β‍|∇u|2dx+β∫|∇u|≥β‍|∇u|-β2dx=∫Ω‍F(|∇u|)dx,
where
(3.4)F(s)={12s2,s<β,βs-12β2,s≥β.
Here β is a fixed positive number and it is a threshold of |∇u|. 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 vectorf, we can get
(3.5)εCL(f)=12∑t∈T(f)-‍|∇ft|2+β∑t∈T(f)+‍(|∇ft|-β2).
Here T(f)-:={t∣|∇ft|<β}, T(f)+:={t∣|∇ft|≥β}. It can be found that the l2 energy of image gradient is considered in the noise part T(f)- while the TV energy is computed in features part T(f)+. Then the new model for image reconstruction can be denoted as follows:
(3.6)argminεCL(f)subjecttoM·f=g.With the Lagrange method applied, this constraint optimization problem can be rewritten as an unconstraint optimization problem of following combined Chambolle-Lions (CCL) energy:(3.7)argminεCCL(f)=∥M·f-g∥22+λ·εCL(f).


## 3.3. Conjugate Gradient Descend Algorithm

The gradient ofεCCL(f) can be computed as
(3.8)∇εCCL(f)=2M*(M·f-g)+λ·∇εCL(f),
where
(3.9)∇εCL(f)={-∇·(∇f),|∇f|<β,-β∇·(∇f|∇f|),|∇f|≥β.
In practice, the parameter λ is set to be 10-2.Set a initial valuef0, the conjugate gradient descend algorithm can be given as Algorithm 1. There the time step is denoted by τ.Algorithm 1:Iteration algorithm for CL minimization reconstruction.
% Initialization
maxGrad =10-3; maxIter = 100; maxTau = 10-3; s = 4; k = 0;
τ
0 = 1; g0 = ∇εCCL(f0); Δf0 = −g0;
% Iterations
while (∥gk∥2  > maxGrad and k < maxIter and τ > maxTau){
%Linear search
minE = εCCL(fk); opt = 0; i = −s;
while (i <= s){
τ=2-iτ0;
if (εCCL(fk + τΔfk) < minE){minE = εCCL(fk + τΔfk); opt = τ; }
i=i+1; }
τ0 = opt;
fk+1 = fk+τ0Δfk; gk+1 = ∇εCCL (fk+1); γ = ∥gk+1∥2∥Δfk∥2;
Δfk+1 = −gk+1 + γΔfk;
k = k + 1; }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 l2 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
(3.10)PSNR=10·log(2552/MSE)log10,
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 Figure1(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.Reconstruction of Shepp-Logan phantom from 72 views.
(a)
Original
(b)
Art result
(c)
TV result
(d)
CL result
(e)
Horizontal gray
(f)
Vertical grayResult analysis of Example4.1.
(a)
Evolution of PSNR
(b)
Evolution ofτ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 Example4.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.Reconstruction of Shepp-Logan phantom from 24 views.
(a)
Original
(b)
Art result
(c)
TV result
(d)
CL result
(e)
Horizontal gray
(f)
Vertical grayResult analysis of Example4.2.
(a)
Evolution of PSNR
(b)
Evolution ofτExample 4.3.
Reconstruction of fruits image from 72 views.
A fruits image with size256×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.Reconstruction of fruits from 72 views.
(a)
Original
(b)
Art result
(c)
TV result
(d)
CL result
(e)
Horizontal gray
(f)
Vertical grayResult analysis of Example4.3.
(a)
Evolution of PSNR
(b)
Evolution ofτExample 4.4.
Reconstruction of fruits image from 30 views.
The same fruits image as Example4.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.Reconstruction of fruits from 30 views.
(a)
Original
(b)
Art result
(c)
TV result
(d)
CL result
(e)
Horizontal gray
(f)
Vertical grayResult analysis of Example4.4.
(a)
Evolution of PSNR
(b)
Evolution ofτExample 4.5.
Reconstruction of a synopsis phantom from 72 views.
A synopsis phantom with size256×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.Reconstruction of a synopsis phantom from 72 views.
(a)
Original
(b)
Art result
(c)
TV result
(d)
CL result
(e)
Horizontal gray
(f)
Vertical grayResult analysis of Example4.5.
(a)
Evolution of PSNR
(b)
Evolution ofτExample 4.6.
Reconstruction of a synopsis phantom from 20 views.
The same synopsis phantom as Example4.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.Reconstruction of a synopsis phantom from 20 views.
(a)
Original
(b)
Art result
(c)
TV result
(d)
CL result
(e)
Horizontal gray
(f)
Vertical grayResult analysis of Example4.6.
(a)
Evolution of PSNR
(b)
Evolution ofτ


## 5. Conclusion

In this paper, a novel model for image reconstruction from few views in parallel-beam data is proposed. First, thel2 energy of the image gradient and the total variation of the image are combined to the CL energy. The l2 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.



---

*Source: 154630-2012-10-31.xml*