# Weighting Algorithm and Relaxation Strategies of the Landweber Method for Image Reconstruction

Authors: Guanghui Han; Gangrong Qu; Qian Wang
Journal: Mathematical Problems in Engineering
 (2018)

---

## Abstract

The iterative approach is important for image reconstruction with ill-posed problem, especially for limited angle reconstruction. Most of iterative algorithms can be written in the general Landweber scheme. In this context, appropriate relaxation strategies and appropriately chosen weights are critical to yield reconstructed images of high quality. In this paper, based on reducing the condition number of matrixATA, we find one method of weighting matrices for the general Landweber method to improve the reconstructed results. For high resolution images, the approximate iterative matrix is derived. And the new weighting matrices and corresponding relaxation strategies are proposed for the general Landweber method with large dimensional number. Numerical simulations show that the proposed weighting methods are effective in improving the quality of reconstructed image for both complete projection data and limited angle projection data.


---

## Body

## 1. Introduction

Image reconstruction plays a major role in many fields such as X-ray computed tomography (CT). Many imaging processes can be modeled as the following linear systems:(1)Ax=bwhere x=x1,x2,⋯,xNT denotes the reconstructed image, b=b1,b2,⋯,bMT is measurement data, and A=(aij)M×N, where ai,j equals the intersection length of the ith projection ray with the jth pixel xj. Therefore, the reconstructed image can be obtained by solving the system of linear equations. If there is a solution of (1), it is called consistent. Otherwise, it is called inconsistent. Iterative algorithm is important for image reconstruction. In this paper, we use a superscript T to denote the transpose of a vector or matrix and use · to denote the 2-norm of a vector or matrix. We use KN to denote the Hilbert space with the canonical inner product ·,·, where the number field K can be the reals R or the complexes C. The orthogonal component of a subspace E in KN is denoted as E⊥. For a matrix A, the null space and the range of A are denoted by N(A) and R(A), respectively. The orthogonal projection from KN to the null space N(A) is denoted by PA.Because of the rapid increase in computing power of computers, iterative reconstruction algorithms become more and more effective. The general Landweber iterative method is the extensively used scheme which is written as(2)xn+1=xn+λnV-1ATWb-Axn,where λn is the relaxation coefficient and the weighting matrices V and W are symmetric positive definite [1]. The weighted and relaxation strategies are important for the quality of the reconstructed image. Formula (2) becomes some specific iterative methods for the different weighted matrices. If V=IN×N where I denotes the identity matrix and W=1/Mdiag(1/A12,⋯,1/AM2) where Ai denotes the i-th row of A, (2) is the Cimmino’s iterative method [2, 3]. If V-1=diag(M/Nj) and W=diag1/∑jnNjaij2 where Nj(j=1,2,⋯,N) is the number of nonzeros in the jth column of A, (2) becomes the simultaneous version of the diagonally relaxed orthogonal projection (DROP) method [4, 5]. The component averaging (CAV) method is obtained with V=IN×N and W=diag1/∑jnNjaij2 [6]. With V=IN×N and W=IM×M, the Landweber method is as follows:(3)xn+1=xn+λnATb-Axn.Convergence conditions of the method (2) have been established [7–9]. The semiconvergence behavior means that initially the iteration vector approaches a regularized solution, while continuing the iteration often leads to iteration vectors corrupted by noise [10]. Based on the analysis of the semiconvergence behavior, T. Elfving et al. proposed two relaxation strategies [4]. Later, these two strategies were improved by them [11]. If the matrix ATA is positive definite, as a special case of the Richardson iteration, based on minimizing the spectral radius of iterative matrix, the relaxation coefficient of the Landweber iteration (3) is given by(4)λn=2μ1+μminn=1,2,⋯,where μ1 and μmin are the largest and the smallest eigenvalues of ATA [12, 13]. However, if matrix ATA is semipositive definite and singular, it is proved that the iteration x(n) generated by the Landweber scheme (3) with(5)λn=2μ1+μminn=1,2,⋯is divergent because of μmin=0 [8]. In our previous work, based on minimizing the spectral radius of the new iterative matrix, a new relaxation coefficient(6)λn=2μ1+μmn=1,2,⋯,and an accelerating convergence strategy λn∈[1/μ1,2/μ1)(n=1,2,⋯) are proposed for the semipositive definite matrix ATA, where μ1,μm are the largest and smallest positive eigenvalues of ATA, respectively [9].For high resolution images, the number of discretization points of original image is large. The image reconstruction problem is ill-posedness and the coefficient matrix of the reconstructing linear system has a high condition number. Therefore, the motivation of this work is to improve the condition number by weighting applicable matrices. In this paper, we derive an appropriate weighting matrix that can be used to improve the condition number. In addition, for suitable high resolution of reconstructed image, the grid number needs to be increased and the dimensional number of original discrete image becomes large. Given an adequately large number of projections, the discrete original image converges to the continuous original image andATA approaches R∗R, where R is Radon transform and R∗ is the adjoint operator. Here, we do not give the rigorous theoretical proof. Therefore, the smaller eigenvalues of matrix ATA tend to zero. The discrete original image can be represented as a linear combination of the eigenvectors of ATA, and the remainder in this combination tends to zero vector. Then we derive the approximate iterative matrix and propose the new weighting algorithms for Landweber method. Based on minimizing the spectral radius of iterative matrix, we obtain the corresponding relaxation strategies for the proposed methods. The numerical simulations show the advantages of the proposed algorithms in both convergence speed and iteration error compared with Landweber method and Cimmino’s method. The proposed methods are applied to the reconstruction from limited angle projection data and improve the condition number of matrix ATA.The organization of the paper is as follows. Section2 analyzes the relation between the condition number and the spectral radius of iterative matrix. In Section 3, we give a method of weighting matrices to improve the condition number. For the higher resolution images, based on the derived approximate iterative matrix, we propose a new weighting algorithm for the general Landweber and the corresponding relaxation strategies. In Section 4, numerical simulations of image reconstruction are carried out. The advantages of proposed algorithms are verified from two aspects: reconstruction by the complete projection data and reconstruction by the limited angle projection data. Finally, some relevant conclusions are discussed.


## 2. Preliminaries

### 2.1. Preparations

In this section, with the new derived iterative representation, the condition number is analyzed in the iterative matrix. For the sake of argument, some notations are explained as follows. Let{μk}k=1m be the distinct positive eigenvalues of ATA which are ordered such that μ1>⋯>μm>0. Each μk has the algebraic multiplicity pk, respectively, for 1≤k≤m. Let uj,k(1≤j≤pk) be the corresponding orthonormal eigenvectors. It follows that(7)ATAuj,k=μkuj,k.If 0 is an eigenvalue of ATA with an algebraic multiplicity q, let vj be the corresponding orthonormal eigenvectors for 1≤j≤q; then ATAvj=ϑ, where ϑ is a zero vector. With the above notations, the 2-norm of the matrix A satisfies(8)A2=μ1.Let Uk=u1,k,⋯,upk,k,k=1,⋯,m,V = ( v 1 , ⋯ , v q ), and U=(U1,⋯,Um,V), and then we have(9)UTATAU=diagμ1Ip1,⋯,μmIpm,Oq,where Ik represents the k-order unit matrix and Oq denotes the q-order null matrix. The least-squares functional(10)fx=12b-Ax2is used to find the minimal 2-norm solution of (1). The gradient of f(x) is(11)∇fx=-ATb-Ax.All the minimizers of f(x) satisfy the equation ∇f(x)=ϑ, i.e.,(12)ATAx=ATb,which is called the normal equation. Equation (12) is always solvable and has a minimal 2-norm solution x+, which is equal to A†b, where A† is the Moore-Penrose inverse of A [14]. If (1) is inconsistent, which is only caused by the noise, b can be decomposed into two parts b0+δb, where b0∈R(A), δb∈R(A)⊥, and δb is caused by noise. By the conclusion of orthogonal decomposition, there is R(A)⊥=N(AT). Thus, for the inconsistent equation, the treatment is as follows:(13)ATAx=ATb0+δb=ATb0.The minimal 2-norm solution x+ can be written as(14)x+=∑k=1m∑j=1pkcj,kuj,k,where cj,k∈K, for 1≤k≤m, 1≤j≤pk. The initial guess x(0) can be decomposed as(15)x0=∑k=1m∑j=1pkcj,k0uj,k+PAx0=x10+PAx0,where cj,k(0)∈K, for 1≤k≤m, 1≤j≤pk. We use the following convergence results proposed [8, 9].Theorem 1.
Forn=1,⋯, the nth iteration x(n) generated by the Landweber method (3) satisfies(16)xn-x+-PAx0=∏s=0n-1I-λsATAx10-x+=Udiag∏s=0n-11-λsμ1Ip1,⋯,∏s=0n-11-λsμmIpm,Oqr,where r=(r1,1,⋯,r1,p1,⋯,rm,1,,⋯,rm,pm,0,⋯,0)T∈KN with rj,k=cj,k-cj,k(0) for 1≤j≤pk, 1≤k≤m, and p=p1+⋯+pm.Let(17)τn=λnμ1,a1,n=min⁡τn,2-τn,where μ1 is the largest eigenvalue of ATA.Theorem 2.
(a) Forn=1,⋯,(18)xn-x+-PAx0=-∑k=1m∏s=0n-11-λsμk∑j=1pkrj,kuj,k≕rn,(19)xn-x++PAx022=∑k=1m∏s=0n-11-λsμk2∑j=1pkrj,k2.(b) Assume that the relaxation coefficients are chosen such that 0<τn<2 for n≥0. Then the iteration x(n) generated by the Landweber method (3) converges to a solution of the normal equation (12) for any b∈KM and any initial guess x(0)∈KN if(20)∑n=0∞a1,n=+∞.Under the conditions, the limit of the iteration x(n) is equal to x++PA(x(0)).For the Landweber method (3), by (16), we define the error vector ϵ(n)(21)ϵn≔xn-x+-PAx0,n≥0,and the iterative matrix G(λn)(22)Gλn=Udiag1-λnμ1Ip1,⋯,1-λnμmIpm,OqUT.Then(23)ϵn=Gλnϵn-1,n=1,2,⋯,and(24)ϵn=∏s=1nGλsϵ0,where ϵ(0)=Ur. The following two relaxation strategies are proposed [9].Proposition 3 (see [9]).
For an initial guessx(0)∈RN, based on minimizing the spectral radius of the iterative matrix G(λn), the relaxation strategy for the Landweber method (3) is(25)arg⁡minλnρGλn=2μm+μ1,forn=1,2,⋯.Lemma 4 (see [9]).
For anyλ(1)∈(0,1/μ1], its symmetric point about λ=1/μ1 is λ(2)=2/μ1-λ(1), where λ(2)∈[1/μ1,2/μ1) and(26)ρGλ1>ρGλ2.For the Landweber method, Proposition3 gives the relaxation strategy and Lemma 4 obtains an accelerating convergence strategy; i.e., λn∈[1/μ1,2/μ1) in the case of only calculating μ1.


### 2.2. Relation between the Condition Number and Spectral Radius

Definition 5.
Thecondition number of the matrix ATA is defined as(27)KATA=μ1μm,where μ1, μm are the largest and smallest positive eigenvalues of matrix ATA, respectively.In the following, we analyze the relation between the condition number and spectral radius. Letτn=λnμ1; then the equivalent form of (3) can be written as(28)xn+1=xn+τn1μ1ATb-Axn.Based on minimizing the spectral radius of the newly derived iterative matrix G(λn) in (22), if μ1 and μm are all known, then the optimal relaxation strategy is(29)τn=21+1/KATA.If only μ1 is known, then the second relaxation strategy is(30)τn∈1,2.For clarity, we useρ(G(λ),A) to represent the spectral radius of the iterative matrix G(λn) in (22) being relative to matrix A. By (22) and Lemma 4, let 1/μ1≤λ<2/μ1; then τ=λμ1∈[1,2) and(31)ρGλ,A=max⁡λμ1-1,1-λμm=max⁡τ-1,1-τμmμ1,for1≤τ<2.If τ-1=1-τ(μm/μ1), then τ=2/(1+μm/μ1). Hence,(32)ρGτ,A=ρGλ,A=1-τ1KATA,for1≤τ≤21+1/KATA;τ-1,for21+1/KATA<τ<2.According to (32), if matrices A1, A2 satisfy K(A1TA1)<K(A2TA2), then(33)1KA1TA1>1KA2TA2,21+1/KA1TA1<21+1/KA2TA2.Therefore, in this case, we get(34)ρGτ,A1=1-τ1KA1TA1,for1≤τ≤21+1/KA1TA1;τ-1,for21+1/KA1TA1<τ≤21+1/KA2TA2;τ-1,for21+1/KA2TA2<τ<2.By (32), (34),(35)ρGτ,A1<ρGτ,A2,for1≤τ≤21+1/KA1TA1;(36)ρGτ,A1=τ-1<1-τ1KA2TA2=ρGτ,A2,for21+1/KA1TA1<τ≤21+1/KA2TA2;(37)ρGτ,A1=τ-1=ρGτ,A2,for21+1/KA2TA2<τ<2.IfA of (28) is equal to A1 and A2, respectively, two iterative algorithms are as follows:(38)xn+1=xn+τn1μ11A1Tb-A1xn,and(39)xn+1=xn+τn1μ12A2Tb-A2xn,where μ1(1), μ1(2) are the largest eigenvalues of A1TA1, A2TA2, respectively. The iterative matrices corresponding to the two algorithms (38), (39) are denoted as G1(τn), G2(τn). By (35), (36), and (37), we have the following theorem.Proposition 6.
Assume thatK(A1TA1)<K(A2TA2).
(1) If the optimal relaxation strategy (29) τn=2/(1+1/K(A1TA1)) is selected, then(40)ρG1τn,A1<ρG2τn,A2,forn=1,2,⋯.(2) For the relaxation strategy (30) τn∈[1,2), if τn∈[1,2/(1+1/K(A2TA2))), then(41)ρG1τn,A1<ρG2τn,A2,forn=1,2,⋯.If τn∈[2/(1+1/K(A2TA2)),2), then(42)ρG1τn,A1=ρG2τn,A2,forn=1,2,⋯.Proof.
Forτn=2/(1+1/K(A1TA1))∈[1,2/(1+1/K(A2TA2))], by (35),(43)ρG1τn,A1<ρG2τn,A2,forn=1,2,⋯.For τn∈[1,2), by (36), (37),(44)ρG1τn,A1<ρG2τn,A2,forτn∈1,21+1/KA2TA2,ρG1τn,A1=ρG2τn,A2,forτn∈21+1/KA2TA2,2,n=1,2,⋯.The theorem is established.


## 2.1. Preparations

In this section, with the new derived iterative representation, the condition number is analyzed in the iterative matrix. For the sake of argument, some notations are explained as follows. Let{μk}k=1m be the distinct positive eigenvalues of ATA which are ordered such that μ1>⋯>μm>0. Each μk has the algebraic multiplicity pk, respectively, for 1≤k≤m. Let uj,k(1≤j≤pk) be the corresponding orthonormal eigenvectors. It follows that(7)ATAuj,k=μkuj,k.If 0 is an eigenvalue of ATA with an algebraic multiplicity q, let vj be the corresponding orthonormal eigenvectors for 1≤j≤q; then ATAvj=ϑ, where ϑ is a zero vector. With the above notations, the 2-norm of the matrix A satisfies(8)A2=μ1.Let Uk=u1,k,⋯,upk,k,k=1,⋯,m,V = ( v 1 , ⋯ , v q ), and U=(U1,⋯,Um,V), and then we have(9)UTATAU=diagμ1Ip1,⋯,μmIpm,Oq,where Ik represents the k-order unit matrix and Oq denotes the q-order null matrix. The least-squares functional(10)fx=12b-Ax2is used to find the minimal 2-norm solution of (1). The gradient of f(x) is(11)∇fx=-ATb-Ax.All the minimizers of f(x) satisfy the equation ∇f(x)=ϑ, i.e.,(12)ATAx=ATb,which is called the normal equation. Equation (12) is always solvable and has a minimal 2-norm solution x+, which is equal to A†b, where A† is the Moore-Penrose inverse of A [14]. If (1) is inconsistent, which is only caused by the noise, b can be decomposed into two parts b0+δb, where b0∈R(A), δb∈R(A)⊥, and δb is caused by noise. By the conclusion of orthogonal decomposition, there is R(A)⊥=N(AT). Thus, for the inconsistent equation, the treatment is as follows:(13)ATAx=ATb0+δb=ATb0.The minimal 2-norm solution x+ can be written as(14)x+=∑k=1m∑j=1pkcj,kuj,k,where cj,k∈K, for 1≤k≤m, 1≤j≤pk. The initial guess x(0) can be decomposed as(15)x0=∑k=1m∑j=1pkcj,k0uj,k+PAx0=x10+PAx0,where cj,k(0)∈K, for 1≤k≤m, 1≤j≤pk. We use the following convergence results proposed [8, 9].Theorem 1.
Forn=1,⋯, the nth iteration x(n) generated by the Landweber method (3) satisfies(16)xn-x+-PAx0=∏s=0n-1I-λsATAx10-x+=Udiag∏s=0n-11-λsμ1Ip1,⋯,∏s=0n-11-λsμmIpm,Oqr,where r=(r1,1,⋯,r1,p1,⋯,rm,1,,⋯,rm,pm,0,⋯,0)T∈KN with rj,k=cj,k-cj,k(0) for 1≤j≤pk, 1≤k≤m, and p=p1+⋯+pm.Let(17)τn=λnμ1,a1,n=min⁡τn,2-τn,where μ1 is the largest eigenvalue of ATA.Theorem 2.
(a) Forn=1,⋯,(18)xn-x+-PAx0=-∑k=1m∏s=0n-11-λsμk∑j=1pkrj,kuj,k≕rn,(19)xn-x++PAx022=∑k=1m∏s=0n-11-λsμk2∑j=1pkrj,k2.(b) Assume that the relaxation coefficients are chosen such that 0<τn<2 for n≥0. Then the iteration x(n) generated by the Landweber method (3) converges to a solution of the normal equation (12) for any b∈KM and any initial guess x(0)∈KN if(20)∑n=0∞a1,n=+∞.Under the conditions, the limit of the iteration x(n) is equal to x++PA(x(0)).For the Landweber method (3), by (16), we define the error vector ϵ(n)(21)ϵn≔xn-x+-PAx0,n≥0,and the iterative matrix G(λn)(22)Gλn=Udiag1-λnμ1Ip1,⋯,1-λnμmIpm,OqUT.Then(23)ϵn=Gλnϵn-1,n=1,2,⋯,and(24)ϵn=∏s=1nGλsϵ0,where ϵ(0)=Ur. The following two relaxation strategies are proposed [9].Proposition 3 (see [9]).
For an initial guessx(0)∈RN, based on minimizing the spectral radius of the iterative matrix G(λn), the relaxation strategy for the Landweber method (3) is(25)arg⁡minλnρGλn=2μm+μ1,forn=1,2,⋯.Lemma 4 (see [9]).
For anyλ(1)∈(0,1/μ1], its symmetric point about λ=1/μ1 is λ(2)=2/μ1-λ(1), where λ(2)∈[1/μ1,2/μ1) and(26)ρGλ1>ρGλ2.For the Landweber method, Proposition3 gives the relaxation strategy and Lemma 4 obtains an accelerating convergence strategy; i.e., λn∈[1/μ1,2/μ1) in the case of only calculating μ1.


## 2.2. Relation between the Condition Number and Spectral Radius

Definition 5.
Thecondition number of the matrix ATA is defined as(27)KATA=μ1μm,where μ1, μm are the largest and smallest positive eigenvalues of matrix ATA, respectively.In the following, we analyze the relation between the condition number and spectral radius. Letτn=λnμ1; then the equivalent form of (3) can be written as(28)xn+1=xn+τn1μ1ATb-Axn.Based on minimizing the spectral radius of the newly derived iterative matrix G(λn) in (22), if μ1 and μm are all known, then the optimal relaxation strategy is(29)τn=21+1/KATA.If only μ1 is known, then the second relaxation strategy is(30)τn∈1,2.For clarity, we useρ(G(λ),A) to represent the spectral radius of the iterative matrix G(λn) in (22) being relative to matrix A. By (22) and Lemma 4, let 1/μ1≤λ<2/μ1; then τ=λμ1∈[1,2) and(31)ρGλ,A=max⁡λμ1-1,1-λμm=max⁡τ-1,1-τμmμ1,for1≤τ<2.If τ-1=1-τ(μm/μ1), then τ=2/(1+μm/μ1). Hence,(32)ρGτ,A=ρGλ,A=1-τ1KATA,for1≤τ≤21+1/KATA;τ-1,for21+1/KATA<τ<2.According to (32), if matrices A1, A2 satisfy K(A1TA1)<K(A2TA2), then(33)1KA1TA1>1KA2TA2,21+1/KA1TA1<21+1/KA2TA2.Therefore, in this case, we get(34)ρGτ,A1=1-τ1KA1TA1,for1≤τ≤21+1/KA1TA1;τ-1,for21+1/KA1TA1<τ≤21+1/KA2TA2;τ-1,for21+1/KA2TA2<τ<2.By (32), (34),(35)ρGτ,A1<ρGτ,A2,for1≤τ≤21+1/KA1TA1;(36)ρGτ,A1=τ-1<1-τ1KA2TA2=ρGτ,A2,for21+1/KA1TA1<τ≤21+1/KA2TA2;(37)ρGτ,A1=τ-1=ρGτ,A2,for21+1/KA2TA2<τ<2.IfA of (28) is equal to A1 and A2, respectively, two iterative algorithms are as follows:(38)xn+1=xn+τn1μ11A1Tb-A1xn,and(39)xn+1=xn+τn1μ12A2Tb-A2xn,where μ1(1), μ1(2) are the largest eigenvalues of A1TA1, A2TA2, respectively. The iterative matrices corresponding to the two algorithms (38), (39) are denoted as G1(τn), G2(τn). By (35), (36), and (37), we have the following theorem.Proposition 6.
Assume thatK(A1TA1)<K(A2TA2).
(1) If the optimal relaxation strategy (29) τn=2/(1+1/K(A1TA1)) is selected, then(40)ρG1τn,A1<ρG2τn,A2,forn=1,2,⋯.(2) For the relaxation strategy (30) τn∈[1,2), if τn∈[1,2/(1+1/K(A2TA2))), then(41)ρG1τn,A1<ρG2τn,A2,forn=1,2,⋯.If τn∈[2/(1+1/K(A2TA2)),2), then(42)ρG1τn,A1=ρG2τn,A2,forn=1,2,⋯.Proof.
Forτn=2/(1+1/K(A1TA1))∈[1,2/(1+1/K(A2TA2))], by (35),(43)ρG1τn,A1<ρG2τn,A2,forn=1,2,⋯.For τn∈[1,2), by (36), (37),(44)ρG1τn,A1<ρG2τn,A2,forτn∈1,21+1/KA2TA2,ρG1τn,A1=ρG2τn,A2,forτn∈21+1/KA2TA2,2,n=1,2,⋯.The theorem is established.


## 3. New Weighting Algorithms

### 3.1. Reduction of Condition Number

According to Proposition6, we consider making K(ATA) become smaller by the weighting matrix V, W in (2). Since μ1>μ2>⋯>μm>0 is decreasing, by (9), using U as the transformation matrix of similarity diagonalization, the positive eigenvalues of the diagonalization matrix constitute an increasing sequence. With the same U as in (9), if ω>μ1/2, we obtain(45)UT2ωI-ATAU=diag2ω-μ1Ip1,⋯,2ω-μmIpm,Oqwhere 2ω-μ1<2ω-μ2<⋯<2ω-μm. Therefore we propose that the weighted matrix V-1 in (2) is taken as (2ωI-ATA).(46)xn+1=xn+λn2ωI-ATAATb-Axn,where ω>μ1/2 and the matrix (2ωI-ATA) is symmetric positive definite. By (16), (22), and Theorem 1, the iterative matrix is(47)Gλn=Udiag⁡1-λn2ω-μ1μ1Ip1,⋯,1-λn2ω-μmμmIpm,OqUT.Proposition 7.
The condition number of the matrix(2ωI-ATA)ATA is denoted as K((2ωI-ATA)ATA), a function of ω. Then
(1) it is monotonically decreasing forμ1/2<ω<(μ1+μm)/2 and monotonically increasing for (μ1+μm)/2≤ω≤μ1;
(2)K((2ωI-ATA)ATA)<K(ATA) for μ1/2<ω≤μ1.Proof.
Using the sameU as in (9), we obtain(48)UT2ωI-ATAATAU=Diag2ω-μ1μ1Ip1,⋯,2ω-μmμmIpm,Oq,where 2ω-μ1μ1,⋯,(2ω-μm)μm are the positive eigenvalues of (2ωI-ATA)ATA and the algebraic multiplicity of every positive eigenvalue is denoted as pk(1≤k≤m), respectively. For μm≤t≤μ1, ω>μ1/2, let(49)lt,ω=2ω-tt=-t-ω2+ω2.It is easy to get(50)arg⁡maxμm≤t≤μ1⁡lt,ω=ω,maxμm≤t≤μ1⁡lt,ω=ω2.For μ1/2<ω≤(μm+μ1)/2 and μm≤t≤μ1, since t=μ1 is the farthest point from t=ω, we obtain(51)arg⁡minμm≤t≤μ1⁡lt,ω=μ1,minμm≤t≤μ1⁡lt,ω=-μ1-ω2+ω2.We need to find the minimum point of the fraction(52)maxμm≤t≤μ1⁡lt,ωminμm≤t≤μ1⁡lt,ω=ω2-μ1-ω2+ω2or the maximum point of its reciprocal(53)-μ1-ω2+ω2ω2=-μ1ω-12+1.For μ1/2<ω≤(μm+μ1)/2<μ1, -(μ1/ω-1)2+1is a monotonically increasing function of ω,(54)arg⁡maxμ1/2<ω≤μm+μ1/2⁡-μ1-ω2+ω2ω2=arg⁡maxμ1/2<ω≤μm+μ1/2⁡-μ1ω-12+1=μm+μ12,and(55)arg⁡minμ1/2<ω≤μm+μ1/2ω2-μ1-ω2+ω2=μm+μ12,(56)minμ1/2<ω≤μm+μ1/2⁡ω2-μ1-ω2+ω2=μm+μ1/22μ1μm.Therefore,(57)Kμ1+μmI-ATAATA=μ1+μm-μk0μk0μ1μm,where(58)μk0=arg⁡min1≤k≤mμ1+μm2-μk.
Similarly, for(μm+μ1)/2≤ω<μ1, t=μmis the farthest point from t=ω,(59)arg⁡maxμm≤t≤μ1⁡lt,ω=ω,maxμm≤t≤μ1⁡lt,ω=ω2,and(60)arg⁡minμm≤t≤μ1⁡lt,ω=μm,minμm≤t≤μ1⁡lt,ω=-μm-ω2+ω2.Therefore, we need to find the minimum point of the fraction(61)maxμm≤t≤μ1⁡lt,ωminμm≤t≤μ1⁡lt,ω=ω2-μm-ω2+ω2or the maximum point of its reciprocal(62)-μm-ω2+ω2ω2=-μmω-12+1.For μm<(μm+μ1)/2≤ω<μ1, -(μm/ω-1)2+1is a monotonically decreasing function of ω,(63)arg⁡maxμm+μ1/2≤ω<μ1-μm-ω2+ω2ω2=arg⁡maxμm+μ1/2≤ω<μ1-μmω-12+1=μm+μ12.and(64)arg⁡minμm+μ1/2≤ω<μ1ω2-μm-ω2+ω2=μm+μ12,(65)minμm+μ1/2≤ω<μ1⁡ω2-μm-ω2+ω2=μm+μ1/22μ1μm.Therefore,(66)Kμ1+μmI-ATAATA=μ1+μm-μk0μk0μ1μm,where(67)μk0=arg⁡min1≤k≤mμ1+μm2-μk.
Forω≥μ1, the function l(t,ω) is monotonically increasing for μm≤t≤μ1, so we get(68)arg⁡maxμm≤t≤μ1lt,ω=μ1,maxμm≤t≤μ1⁡lt,ω=2ωμ1-μ12,and(69)arg⁡minμm≤t≤μ1lt,ω=μm,minμm≤t≤μ1⁡lt,ω=2ωμm-μm2.Since(70)arg⁡minω≥μ1⁡2ωμ1-μ122ωμm-μm2=μ1,we obtain(71)K2μ1I-ATAATA=2μ1μ1-μ122μ1μm-μm2=μ122μ1-μmμm.By (57), (66), and (71), we can readily compute that(72)Kμ1+μmI-ATAATA=μ1+μm-μk0μk0μ1μm<μ12μ1μm=μ1μm=KATA,and(73)K2μ1I-ATAATA=μ12μ1-μm·μ1μm<μ1μm=KATA.Since(74)Kμ1+μmI-ATAATA-K2μ1I-ATAATA<0,we obtain(75)Kμ1+μmI-ATAATA<K2μ1I-ATAATA.This completes the proof.


### 3.2. Weighting and Relaxation Strategies

#### 3.2.1. Analysis for Dimensional Number Increases

In this subsection, we analyze the relations between the original continuous image and discrete image, the Radon transform, and the coefficient matrix of reconstructing algebraic system. Based on these relations, the weighted methods are proposed. The original image is a continuous functionf generally, assumed to be l2-integrable. Algebraic iterative reconstruction is discretizing the original continuous image and Radon transform R to be the discrete image and the linear system, respectively, and then finding an iterative solution of the linear system. To satisfy the suitable high reconstructed resolution for applications, the discretized grid is refined and the grid number is increasing, so the dimensional numbers of the iteration solution and coefficient matrix A become larger. If the projected data are many enough, the discrete original image converges to the continuous original image and ATA approaches R∗R as the grid number increases in a certain form. Since the imaging object is in a bounded domain, the image function is compactly supported. For the Radon transform R in bounded domain, R has singular value decomposition and R∗R is the positive definite compact operator [10]. By the theory of spectrum of the positive definite compact operator, as the dimension of the matrix ATA becomes very large, the smaller positive eigenvalues of ATA tend to be 0 [10]. Since Radon transform R has unique solution, the linear system (1) has unique solution as the projection data is appropriately large [10]. Due to the space limit of the paper, we will not prove these. Then the original discrete image x is(76)x=∑k=1m∑j=1pkdj,kuj,kand(77)x2=∑k=1m∑j=1pkdj,k2,where dj,k∈K, for 1≤k≤m, 1≤j≤pk, and the eigenvectors uj,k,(j=1,…,pk;k=1,…,m) constitute the orthonormal basis of Hilbert space KN. As the grid number increases, the value of m increases and the original discrete image x tends to be the original image f. Assuming that ε is an allowed iterative error. The remainder in (77) is denoted as(78)ϵm0=∑k=m0+1m∑j=1pkdj,k2,where m0 is unknown. If the reconstructed error vector ϵ(n)<ε, there exists m0 such that ϵm0<ε/2. So ϵm0 can be negligible. By Theorem 2,(79)xn-x+-PAx0≈-∑k=1m0∏s=0n-11-λsμk∑j=1pkrj,kuj,k,(80)xn-x++PAx022≈∑k=1m0∏s=0n-11-λsμk2∑j=1pkrj,k2.By (47), the approximate iterative matrix is(81)Gλn≈Udiag⁡1-λn2ω-μ1μ1Ip1,⋯,1-λn2ω-μm0μm0Ipm0,Opm0+1,⋯,Opm,OqUT=Udiag⁡1-λnaμ1-μ1μ1Ip1,⋯,1-λnaμ1-μm0μm0Ipm0,Opm0+1,⋯,Opm,OqUT.


#### 3.2.2. Weighting and Relaxation Strategies

Based on analysis in Section3.2.1 and minimizing the condition number of the approximate iterative matrix (81), by Proposition 7, the parameter ω in (46) is taken as(82)ω=μ1+μm02=μ11+μm0/μ12=aμ12,a=1+μm0μ1,and the corresponding iterative method is(83)NWL-a:xn+1=xn+λnaμ1I-ATAATb-Axn,which is marked as NWL-a in this paper. Since μm0 is unknown, the value of μm0/μ1 is taken in a suitable interval 0.002≤μm0/μ1≤0.5. Then 1.002≤a≤1.5.Now we study the relaxation strategy after determining the weighted matrix. For the iterative method (83), by (25), (47), and Proposition 3,(84)arg⁡minλn⁡ρGλn=2aμ1-μkaμka+aμ1-μmμm>2aμ1-aμ1/2aμ1/2+aμ1-μ1μ1=2a2μ12/4+aμ1-μ1μ1,where(85)μka=arg⁡min1≤k≤m⁡aμ12-μk.We hope that the spectral radius of the iterative matrix is smaller. By the inequality (41) in Proposition 6, the relaxation coefficient is(86)λna=2a2μ12/4+aμ1-μ1μ1.Then λn(a) is a monotonically decreasing function with respect to a and satisfies 2/μ12≤λn(a)≤8/μ12 for 1.002≤a≤1.5.To compare the reconstructed results of different weighting algorithms, by Proposition7, the parameter ω in (46) is taken as ω=μ1 and ω=(μ1+μm)/2, and the corresponding methods are marked as NWL-1 and NWL-2, respectively. By (25) in Proposition 3, for NWL-1, its relaxation coefficient is(87)arg⁡minλn⁡ρG=2μ12+2μ1-μmμm,for NWL-2,(88)arg⁡minλn⁡ρG=2μ1+μm-μk0μk0+μ1μm,by (58),(89)μk0=arg⁡min1≤k≤m⁡μ1+μm2-μk.


## 3.1. Reduction of Condition Number

According to Proposition6, we consider making K(ATA) become smaller by the weighting matrix V, W in (2). Since μ1>μ2>⋯>μm>0 is decreasing, by (9), using U as the transformation matrix of similarity diagonalization, the positive eigenvalues of the diagonalization matrix constitute an increasing sequence. With the same U as in (9), if ω>μ1/2, we obtain(45)UT2ωI-ATAU=diag2ω-μ1Ip1,⋯,2ω-μmIpm,Oqwhere 2ω-μ1<2ω-μ2<⋯<2ω-μm. Therefore we propose that the weighted matrix V-1 in (2) is taken as (2ωI-ATA).(46)xn+1=xn+λn2ωI-ATAATb-Axn,where ω>μ1/2 and the matrix (2ωI-ATA) is symmetric positive definite. By (16), (22), and Theorem 1, the iterative matrix is(47)Gλn=Udiag⁡1-λn2ω-μ1μ1Ip1,⋯,1-λn2ω-μmμmIpm,OqUT.Proposition 7.
The condition number of the matrix(2ωI-ATA)ATA is denoted as K((2ωI-ATA)ATA), a function of ω. Then
(1) it is monotonically decreasing forμ1/2<ω<(μ1+μm)/2 and monotonically increasing for (μ1+μm)/2≤ω≤μ1;
(2)K((2ωI-ATA)ATA)<K(ATA) for μ1/2<ω≤μ1.Proof.
Using the sameU as in (9), we obtain(48)UT2ωI-ATAATAU=Diag2ω-μ1μ1Ip1,⋯,2ω-μmμmIpm,Oq,where 2ω-μ1μ1,⋯,(2ω-μm)μm are the positive eigenvalues of (2ωI-ATA)ATA and the algebraic multiplicity of every positive eigenvalue is denoted as pk(1≤k≤m), respectively. For μm≤t≤μ1, ω>μ1/2, let(49)lt,ω=2ω-tt=-t-ω2+ω2.It is easy to get(50)arg⁡maxμm≤t≤μ1⁡lt,ω=ω,maxμm≤t≤μ1⁡lt,ω=ω2.For μ1/2<ω≤(μm+μ1)/2 and μm≤t≤μ1, since t=μ1 is the farthest point from t=ω, we obtain(51)arg⁡minμm≤t≤μ1⁡lt,ω=μ1,minμm≤t≤μ1⁡lt,ω=-μ1-ω2+ω2.We need to find the minimum point of the fraction(52)maxμm≤t≤μ1⁡lt,ωminμm≤t≤μ1⁡lt,ω=ω2-μ1-ω2+ω2or the maximum point of its reciprocal(53)-μ1-ω2+ω2ω2=-μ1ω-12+1.For μ1/2<ω≤(μm+μ1)/2<μ1, -(μ1/ω-1)2+1is a monotonically increasing function of ω,(54)arg⁡maxμ1/2<ω≤μm+μ1/2⁡-μ1-ω2+ω2ω2=arg⁡maxμ1/2<ω≤μm+μ1/2⁡-μ1ω-12+1=μm+μ12,and(55)arg⁡minμ1/2<ω≤μm+μ1/2ω2-μ1-ω2+ω2=μm+μ12,(56)minμ1/2<ω≤μm+μ1/2⁡ω2-μ1-ω2+ω2=μm+μ1/22μ1μm.Therefore,(57)Kμ1+μmI-ATAATA=μ1+μm-μk0μk0μ1μm,where(58)μk0=arg⁡min1≤k≤mμ1+μm2-μk.
Similarly, for(μm+μ1)/2≤ω<μ1, t=μmis the farthest point from t=ω,(59)arg⁡maxμm≤t≤μ1⁡lt,ω=ω,maxμm≤t≤μ1⁡lt,ω=ω2,and(60)arg⁡minμm≤t≤μ1⁡lt,ω=μm,minμm≤t≤μ1⁡lt,ω=-μm-ω2+ω2.Therefore, we need to find the minimum point of the fraction(61)maxμm≤t≤μ1⁡lt,ωminμm≤t≤μ1⁡lt,ω=ω2-μm-ω2+ω2or the maximum point of its reciprocal(62)-μm-ω2+ω2ω2=-μmω-12+1.For μm<(μm+μ1)/2≤ω<μ1, -(μm/ω-1)2+1is a monotonically decreasing function of ω,(63)arg⁡maxμm+μ1/2≤ω<μ1-μm-ω2+ω2ω2=arg⁡maxμm+μ1/2≤ω<μ1-μmω-12+1=μm+μ12.and(64)arg⁡minμm+μ1/2≤ω<μ1ω2-μm-ω2+ω2=μm+μ12,(65)minμm+μ1/2≤ω<μ1⁡ω2-μm-ω2+ω2=μm+μ1/22μ1μm.Therefore,(66)Kμ1+μmI-ATAATA=μ1+μm-μk0μk0μ1μm,where(67)μk0=arg⁡min1≤k≤mμ1+μm2-μk.
Forω≥μ1, the function l(t,ω) is monotonically increasing for μm≤t≤μ1, so we get(68)arg⁡maxμm≤t≤μ1lt,ω=μ1,maxμm≤t≤μ1⁡lt,ω=2ωμ1-μ12,and(69)arg⁡minμm≤t≤μ1lt,ω=μm,minμm≤t≤μ1⁡lt,ω=2ωμm-μm2.Since(70)arg⁡minω≥μ1⁡2ωμ1-μ122ωμm-μm2=μ1,we obtain(71)K2μ1I-ATAATA=2μ1μ1-μ122μ1μm-μm2=μ122μ1-μmμm.By (57), (66), and (71), we can readily compute that(72)Kμ1+μmI-ATAATA=μ1+μm-μk0μk0μ1μm<μ12μ1μm=μ1μm=KATA,and(73)K2μ1I-ATAATA=μ12μ1-μm·μ1μm<μ1μm=KATA.Since(74)Kμ1+μmI-ATAATA-K2μ1I-ATAATA<0,we obtain(75)Kμ1+μmI-ATAATA<K2μ1I-ATAATA.This completes the proof.


## 3.2. Weighting and Relaxation Strategies

### 3.2.1. Analysis for Dimensional Number Increases

In this subsection, we analyze the relations between the original continuous image and discrete image, the Radon transform, and the coefficient matrix of reconstructing algebraic system. Based on these relations, the weighted methods are proposed. The original image is a continuous functionf generally, assumed to be l2-integrable. Algebraic iterative reconstruction is discretizing the original continuous image and Radon transform R to be the discrete image and the linear system, respectively, and then finding an iterative solution of the linear system. To satisfy the suitable high reconstructed resolution for applications, the discretized grid is refined and the grid number is increasing, so the dimensional numbers of the iteration solution and coefficient matrix A become larger. If the projected data are many enough, the discrete original image converges to the continuous original image and ATA approaches R∗R as the grid number increases in a certain form. Since the imaging object is in a bounded domain, the image function is compactly supported. For the Radon transform R in bounded domain, R has singular value decomposition and R∗R is the positive definite compact operator [10]. By the theory of spectrum of the positive definite compact operator, as the dimension of the matrix ATA becomes very large, the smaller positive eigenvalues of ATA tend to be 0 [10]. Since Radon transform R has unique solution, the linear system (1) has unique solution as the projection data is appropriately large [10]. Due to the space limit of the paper, we will not prove these. Then the original discrete image x is(76)x=∑k=1m∑j=1pkdj,kuj,kand(77)x2=∑k=1m∑j=1pkdj,k2,where dj,k∈K, for 1≤k≤m, 1≤j≤pk, and the eigenvectors uj,k,(j=1,…,pk;k=1,…,m) constitute the orthonormal basis of Hilbert space KN. As the grid number increases, the value of m increases and the original discrete image x tends to be the original image f. Assuming that ε is an allowed iterative error. The remainder in (77) is denoted as(78)ϵm0=∑k=m0+1m∑j=1pkdj,k2,where m0 is unknown. If the reconstructed error vector ϵ(n)<ε, there exists m0 such that ϵm0<ε/2. So ϵm0 can be negligible. By Theorem 2,(79)xn-x+-PAx0≈-∑k=1m0∏s=0n-11-λsμk∑j=1pkrj,kuj,k,(80)xn-x++PAx022≈∑k=1m0∏s=0n-11-λsμk2∑j=1pkrj,k2.By (47), the approximate iterative matrix is(81)Gλn≈Udiag⁡1-λn2ω-μ1μ1Ip1,⋯,1-λn2ω-μm0μm0Ipm0,Opm0+1,⋯,Opm,OqUT=Udiag⁡1-λnaμ1-μ1μ1Ip1,⋯,1-λnaμ1-μm0μm0Ipm0,Opm0+1,⋯,Opm,OqUT.


### 3.2.2. Weighting and Relaxation Strategies

Based on analysis in Section3.2.1 and minimizing the condition number of the approximate iterative matrix (81), by Proposition 7, the parameter ω in (46) is taken as(82)ω=μ1+μm02=μ11+μm0/μ12=aμ12,a=1+μm0μ1,and the corresponding iterative method is(83)NWL-a:xn+1=xn+λnaμ1I-ATAATb-Axn,which is marked as NWL-a in this paper. Since μm0 is unknown, the value of μm0/μ1 is taken in a suitable interval 0.002≤μm0/μ1≤0.5. Then 1.002≤a≤1.5.Now we study the relaxation strategy after determining the weighted matrix. For the iterative method (83), by (25), (47), and Proposition 3,(84)arg⁡minλn⁡ρGλn=2aμ1-μkaμka+aμ1-μmμm>2aμ1-aμ1/2aμ1/2+aμ1-μ1μ1=2a2μ12/4+aμ1-μ1μ1,where(85)μka=arg⁡min1≤k≤m⁡aμ12-μk.We hope that the spectral radius of the iterative matrix is smaller. By the inequality (41) in Proposition 6, the relaxation coefficient is(86)λna=2a2μ12/4+aμ1-μ1μ1.Then λn(a) is a monotonically decreasing function with respect to a and satisfies 2/μ12≤λn(a)≤8/μ12 for 1.002≤a≤1.5.To compare the reconstructed results of different weighting algorithms, by Proposition7, the parameter ω in (46) is taken as ω=μ1 and ω=(μ1+μm)/2, and the corresponding methods are marked as NWL-1 and NWL-2, respectively. By (25) in Proposition 3, for NWL-1, its relaxation coefficient is(87)arg⁡minλn⁡ρG=2μ12+2μ1-μmμm,for NWL-2,(88)arg⁡minλn⁡ρG=2μ1+μm-μk0μk0+μ1μm,by (58),(89)μk0=arg⁡min1≤k≤m⁡μ1+μm2-μk.


## 3.2.1. Analysis for Dimensional Number Increases

In this subsection, we analyze the relations between the original continuous image and discrete image, the Radon transform, and the coefficient matrix of reconstructing algebraic system. Based on these relations, the weighted methods are proposed. The original image is a continuous functionf generally, assumed to be l2-integrable. Algebraic iterative reconstruction is discretizing the original continuous image and Radon transform R to be the discrete image and the linear system, respectively, and then finding an iterative solution of the linear system. To satisfy the suitable high reconstructed resolution for applications, the discretized grid is refined and the grid number is increasing, so the dimensional numbers of the iteration solution and coefficient matrix A become larger. If the projected data are many enough, the discrete original image converges to the continuous original image and ATA approaches R∗R as the grid number increases in a certain form. Since the imaging object is in a bounded domain, the image function is compactly supported. For the Radon transform R in bounded domain, R has singular value decomposition and R∗R is the positive definite compact operator [10]. By the theory of spectrum of the positive definite compact operator, as the dimension of the matrix ATA becomes very large, the smaller positive eigenvalues of ATA tend to be 0 [10]. Since Radon transform R has unique solution, the linear system (1) has unique solution as the projection data is appropriately large [10]. Due to the space limit of the paper, we will not prove these. Then the original discrete image x is(76)x=∑k=1m∑j=1pkdj,kuj,kand(77)x2=∑k=1m∑j=1pkdj,k2,where dj,k∈K, for 1≤k≤m, 1≤j≤pk, and the eigenvectors uj,k,(j=1,…,pk;k=1,…,m) constitute the orthonormal basis of Hilbert space KN. As the grid number increases, the value of m increases and the original discrete image x tends to be the original image f. Assuming that ε is an allowed iterative error. The remainder in (77) is denoted as(78)ϵm0=∑k=m0+1m∑j=1pkdj,k2,where m0 is unknown. If the reconstructed error vector ϵ(n)<ε, there exists m0 such that ϵm0<ε/2. So ϵm0 can be negligible. By Theorem 2,(79)xn-x+-PAx0≈-∑k=1m0∏s=0n-11-λsμk∑j=1pkrj,kuj,k,(80)xn-x++PAx022≈∑k=1m0∏s=0n-11-λsμk2∑j=1pkrj,k2.By (47), the approximate iterative matrix is(81)Gλn≈Udiag⁡1-λn2ω-μ1μ1Ip1,⋯,1-λn2ω-μm0μm0Ipm0,Opm0+1,⋯,Opm,OqUT=Udiag⁡1-λnaμ1-μ1μ1Ip1,⋯,1-λnaμ1-μm0μm0Ipm0,Opm0+1,⋯,Opm,OqUT.


## 3.2.2. Weighting and Relaxation Strategies

Based on analysis in Section3.2.1 and minimizing the condition number of the approximate iterative matrix (81), by Proposition 7, the parameter ω in (46) is taken as(82)ω=μ1+μm02=μ11+μm0/μ12=aμ12,a=1+μm0μ1,and the corresponding iterative method is(83)NWL-a:xn+1=xn+λnaμ1I-ATAATb-Axn,which is marked as NWL-a in this paper. Since μm0 is unknown, the value of μm0/μ1 is taken in a suitable interval 0.002≤μm0/μ1≤0.5. Then 1.002≤a≤1.5.Now we study the relaxation strategy after determining the weighted matrix. For the iterative method (83), by (25), (47), and Proposition 3,(84)arg⁡minλn⁡ρGλn=2aμ1-μkaμka+aμ1-μmμm>2aμ1-aμ1/2aμ1/2+aμ1-μ1μ1=2a2μ12/4+aμ1-μ1μ1,where(85)μka=arg⁡min1≤k≤m⁡aμ12-μk.We hope that the spectral radius of the iterative matrix is smaller. By the inequality (41) in Proposition 6, the relaxation coefficient is(86)λna=2a2μ12/4+aμ1-μ1μ1.Then λn(a) is a monotonically decreasing function with respect to a and satisfies 2/μ12≤λn(a)≤8/μ12 for 1.002≤a≤1.5.To compare the reconstructed results of different weighting algorithms, by Proposition7, the parameter ω in (46) is taken as ω=μ1 and ω=(μ1+μm)/2, and the corresponding methods are marked as NWL-1 and NWL-2, respectively. By (25) in Proposition 3, for NWL-1, its relaxation coefficient is(87)arg⁡minλn⁡ρG=2μ12+2μ1-μmμm,for NWL-2,(88)arg⁡minλn⁡ρG=2μ1+μm-μk0μk0+μ1μm,by (58),(89)μk0=arg⁡min1≤k≤m⁡μ1+μm2-μk.


## 4. Numerical Simulation

In order to illustrate the performance of the proposed algorithms, we simulate the reconstructions of parallel scanning model from noise-free, noise, and limited angle projections. The MATLAB software is used in numerical simulation. All the algorithms are implemented in a PC with 1.70 GHz processor and 4 GB of memory. As the original image, the Shepp-Logan phantom is obtained by the MATLAB function “phantom” and its pixel is64×64. The original image is labeled in a 1D sequential order x¯=(x¯1,x¯2,…,x¯642)T, where x¯ is the exact solution. The initial vector x(0) is selected as zero vector. The projection b=(b1,b2,…,bM)T is obtained by the MATLAB function “radon”. The element aij in A is the segment length of the ith projection ray within the jth pixel xj, which is calculated in terms of the method proposed in [15]. In the process of computing, all iterative reconstruction algorithms are computed by the form(90)xn+1=Qxn+q.For the sake of comparison of different algorithms, four quantitative evaluation indices are used to evaluate the quality of reconstructed images. The relative error is defined as follows:(91)Relative  error≔xn-x¯x¯,and the smaller its value is, the faster the convergence rate is. The other three kinds of image distance indices, d, r, and e, are calculated by the method proposed in [16]. Let tu,v and ru,v denote the element of the row u and the column v, respectively, in the original image T and reconstructed image R. Here T and R are obtained by rearranging the exact solution x¯ and iterative solution x(n), respectively. The formulae of image distance are as follows:(92)d=∑u=1N∑v=1Ntu,v-ru,v2∑u=1N∑v=1Ntu,v-t¯21/2,(93)r=∑u=1N∑v=1Ntu,v-ru,v∑u=1N∑v=1Ntu,v,(94)e=max1≤i≤N/2,1≤j≤N/2⁡Ti,j-Ri,j,where N=64×64,[N/2] denotes the largest integer less than N/2 and(95)t¯=∑i=1N∑j=1Nti,jN,(96)Ti,j=14t2i,2j+t2i+1,2j+t2i,2j+1+t2i+1,2j+1,(97)Ri,j=14r2i,2j+r2i+1,2j+r2i,2j+1+r2i+1,2j+1.The criterion d describes the overall nearness between the original and reconstructed image, the criterion r is the mean relative error of reconstruction, and the criterion e is the maximum difference between the original and reconstructed image. The smaller the values of d,r,e, the better the reconstructed results.Based on the relaxation strategies proposed in Section3, for Landweber, Cimmino, NWL-1, and NWL-2 algorithms, the corresponding relaxation coefficients are obtained by the conclusions (25), (87), and (88), respectively. For the algorithms NWL-a, the values of parameter a and λn satisfy a∈[1.002,1.5], λn∈[2/μ12,8/μ12].The principle of determining the number of iterations is that the relative error does not significantly decrease as the number of iterations increases. According to the results of repeated experiments, the number of iterations is chosen as 200 and 1000, respectively, for the complete projection data and the limited angle projection data.


### 4.1. Complete Projection Data

We use180 projections with 95 rays per projection for the reconstruction of complete projection data. The calculated coefficient matrix A has the dimensions 17100×4096. For the algorithm NWL-a, the corresponding values of a and λn are as follows.The largest eigenvalue ofATA is μ1.Fora=1.002, λn=8/μ12; for a=1.005, λn=7.5/μ12; for a=1.01, λn=7/μ12; and for a=1.2, λn=5/μ12.We compare the performance of different algorithms in terms of reconstruction error and image distances. Table1 shows the computation time for completing 200 iterations by formula (90). It can be seen that the time spent by other algorithms is similar except that the Cimmino algorithm takes a little longer time.Table 1
Computation time for completing 200 iterations.
Algorithms Landweber NWL-a Cimmino DROP CAV Time (s) 2.6839 2.3738 14.4995 2.6840 3.6551


#### 4.1.1. Noise-Free Projection Data

The comparison of relative errors from the proposed algorithms is shown in Figure1, from which we observe that the relative errors obtained from the NWL-a method are less than those obtained from the Cimmino method, the NWL-1 and NWL-2 methods, and the Landweber method. Figure 2 shows the reconstructed phantom images from the noise-free projections. It can be seen that the reconstructed images obtained from the NWL-a method are clearer than those obtained from the Cimmino method, the NWL-1 and NWL-2 methods, and the Landweber method. The image distances for the noise-free projection data are given in Table 2, from which we see that the image distances from NWL-a algorithm are smaller than those from Cimmino and Landweber algorithms. Figure 3 plots the curves of relative error and it shows that the relative errors from the NWL-a method are less than those from the other methods. Figure 4 shows the comparisons of the 32nd column data in the original and the reconstructed images. It is observed that the main features of the reconstructed images obtained from the NWL-a method with a=1.01 are closer to those of the original image. The image distances generated from different algorithms are given in Table 3, from which we see that the image distances from NWL-a algorithm are smaller than those from DROP and CAV algorithms.Table 2
Image distances from different methods for the noise-free projection data.
Algorithms Landweber NWL-a=1.005 NWL-a=1.01 NWL-a=1.2 Cimmino d 0.4873 0.0039 0.0042 0.0052 0.0351 r 0.8743 0.0028 0.0032 0.0046 0.0386 e 0.1460 0.0033 0.0035 0.0039 0.0064Table 3
Comparison of image distances with DROP and CAV algorithms.
Algorithms NWL-a=1.005 NWL-a=1.01 NWL-a=1.2 DROP CAV d 0.0039 0.0042 0.0052 0.3512 0.3489 r 0.0028 0.0032 0.0046 0.5865 0.5893 e 0.0033 0.0035 0.0039 0.0906 0.0862Figure 1
Relative error histories generated from seven iterative algorithms.Figure 2
(a) Original image. (b)-(h) Images reconstructed by seven iterative algorithms from the noise-free projections.
(a) (b) (c) (d) (e) (f) (g) (h)Figure 3
Relative error histories generated from seven iterative algorithms.Figure 4
Comparisons of the 32nd column in the reconstructed and original images. The red line is reconstruction and blue line represents the original image.


#### 4.1.2. Noisy Projection Data

In this section, we test the performance of different algorithms for the projection data with noise. For this purpose, the white Gauss noise with signal-to-noise ratio (SNR) of 10 is added to the projectionsb. The reconstructed results are shown in Figures 5, 6, and 7. The image distances for the noisy projection data are given in Table 4, from which we see that the image distances from NWL-a algorithm are smaller than those from Cimmino and Landweber algorithms. The conclusion is the same as above. The quality of the reconstructed images obtained from the NWL-a method is better than the quality of those obtained from the other methods.Table 4
Image distances from different methods for the noisy projection data.
Algorithms Landweber NWL-a=1.005 NWL-a=1.01 NWL-a=1.2 Cimmino d 0.4954 0.1156 0.1187 0.1199 0.1296 r 0.8745 0.1471 0.1506 0.1517 0.1614 e 0.1617 0.0334 0.0369 0.0394 0.0488Figure 5
(a) Original image. (b)-(h) Images reconstructed by seven iterative algorithms from the projections with noise.
(a) (b) (c) (d) (e) (f) (g) (h)Figure 6
Relative error histories generated from seven iterative algorithms.Figure 7
Comparisons of the 32nd column in the reconstructed and original images. The red line is reconstruction and blue line represents the original image.


### 4.2. Limited Angle Reconstruction

In the theory of continuous case, the problem of limited angle reconstruction has been discussed [17–20]. Assume that f∈L2(R2) is compactly supported and the support of f is supp(f)⊂B={x∈R2:|x|≤c}. In two-dimensional space, S+1={ω=(cos⁡θ,sin⁡θ):0≤θ≤π}. Let Θ0={ω=(cos⁡θ,sin⁡θ):0≤θ≤α<π}. The Radon transform of f is defined by(98)Rfω,p=∫ω·x=pfxds.Let Ω=R×S+1 and Ω0=R×Θ0. If Rf(ω,p) is only known on p∈R and a subset Θ0⊂S+1, the data is called angle-limited. By projection slice theorem,(99)F2fσω=F1Rfσω,and(100)F2fξ=∫Rne-2iπξ·xfxdx.Let PB and PΩ0 be the characteristic functions of B and Ω0, H(ξ)=PΩ0F(ξ). Then(101)Ω0FBf=Gwhere Bf=PBf and Ω0FBf=PΩ0FBf. We let A=Ω0FB. It follows from (101) that the adjoint operator A∗ of A is given by(102)A∗Gx=BF-1Ω0Gx=PBx∫ξ∈Ω0Gξe-i2πξ·xdξ.Multiplying A∗ on both sides of (101), we obtain(103)A∗Af=A∗G.The integral equation (101) can be uniquely solved since A∗A is self-adjoint and positive definite [18]. Therefore the eigenvalues σk of the eigenfunctions of A∗A satisfy 1>σ1≥σ2≥⋯. It turns out that the eigenvalues σk are split into two rather distinct subsets: some are “close to 1” and the rest are “close to 0” [18]. When the range of known limited angles become small, more eigenvalues are “close to 0”, which leads to ill-posed problem becoming serious and degrading the quality of reconstructed image. Here we do not study the theoretical relation between the eigenvalue distributions of derived equation (103) in continuous case and those of discrete algebraic linear systems for limited angle reconstruction problem; we just use our methods to improve the condition number of discrete algebraic linear systems, so as to improve the results of limited angle reconstruction.


#### 4.2.1. Limited Angle 0°−119°

To test the effectiveness of the new algorithms in this aspect, image reconstruction using limited angle projection data is simulated in this section. The projection angle is set from 0°−119° and we use 120 projections with 95 rays per projection. The calculated coefficient matrixA has the dimensions 11400×4096. The number of iterations is 1000. For the algorithms NWL-a, the corresponding values of a and λn are as follows. The largest eigenvalue of ATA is μ1. For a=1.005, λn=8/μ12; for a=1.01, λn=7/μ12; and for a=1.2, λn=5/μ12. Completing 1000 iterations, the Cimmino method and the proposed methods take about 15 seconds.From Figure8, we see that the reconstruction images from the NWL-a method are clearer than those obtained from the other methods. Although the results obtained from the Cimmino method are good, it is time-consuming to calculate the iteration matrix Q of Cimmino method before the iteration begins. Figure 9 shows that the NWL-a method makes the relative error smaller. The comparison of the 32nd column data is shown in Figure 10, from which we observe that, for incomplete projection data, the reconstructed image has a certain deviation from the original image. However the reconstructed images obtained from the NWL-a method are the closest to the original image.Figure 8
(a) Original image. (b)-(h) Images reconstructed by seven iterative algorithms for the limited projections 0°−119°.
(a) (b) (c) (d) (e) (f) (g) (h)Figure 9
Relative error histories generated from seven iterative algorithms for the limited angle 0°−119°.Figure 10
Comparisons of the 32nd column in the reconstructed and original images. The red line is reconstruction and blue line represents the original image.


#### 4.2.2. Four Kinds of Limited Angle Projection Data

In order to further test the performance of new algorithms on limited angle reconstruction, comparing with the Cimmino method, we select four kinds of limited angle scans which cover0∘-109∘, 0∘-99∘, 0∘-89∘, and 0∘-79∘, respectively. The number of iterations is 1000. Corresponding to limited angle 0∘-109∘, 0∘-99∘, 0∘-89∘, and 0∘-79∘, the calculated coefficient matrices A have the dimensions, respectively, 10450×4096, 9500×4096, 8550×4096, and 7600×4096. For the NWL-a algorithm, the corresponding values of a and λn are as follows. For a=1.005, λn=8/μ1a=1.0052; for a=1.01, λn=7/μ1a=1.012; and for a=1.2, λn=5/μ1a=1.22, where μ1a denotes the largest eigenvalue corresponding to the parameter a.Figure11 shows the reconstructed results by the methods FBP (Filter Backprojection), NWL-1, NWL-2, and Cimmino. Figure 12 shows the reconstructed results from the NWL-a method with a=1.005,a=1.01,a=1.2. From these figures, it is observed that the quality of the reconstructed images becomes worse with the narrowing of the limited angle range. Corresponding to Figure 12, the image distances from the NWL-a method are represented by Tables 5 and 6. It can be seen that the image distances d, r, and e obtained by NWL-a method with a = 1.005, 1.01, 1.2 are similar for each kind of limited angle projection data. The values of d, r, and e become larger with the narrowing of the limited angle range, which is consistent with the fact. However the reconstructions of the main features in the original image are clear when the range of the limited angle is at least 90∘. In contrast, the reconstructed images obtained via the NWL-a method are more clear. The reconstructed images are not clear when the range of the limited angle is less than 90∘. Figures 13 and 15 show the comparison of relative errors from the six iterative methods. It is observed that the convergence speed of the NWL-a method is the fastest. From Figures 14 and 16, we see that the reconstructed image obtained from the NWL-a method is closer to the original image. The numerical simulations of limited angle reconstruction show that the NWL-a method has improved the reconstruction results.Table 5
Image distances from the NWL-a method for the limited angle projection data0∘-109∘ and 0∘-99∘.
Limited angle 0 ∘ - 10 9 ∘ 0 ∘ - 9 9 ∘ NWL-a a=1.005 a=1.01 a=1.2 a=1.005 a=1.01 a=1.2 d 0.2112 0.2172 0.2328 0.2269 0.2338 0.2521 r 0.2267 0.2317 0.2443 0.2386 0.2441 0.2584 e 0.1904 0.1979 0.2165 0.1963 0.2033 0.2208Table 6
Image distances from the NWL-a method for the limited angle projection data0∘-89∘ and 0∘-79∘.
Limited angle 0 ∘ - 8 9 ∘ 0 ∘ - 7 9 ∘ NWL-a a=1.005 a=1.01 a=1.2 a=1.005 a=1.01 a=1.2 d 0.3328 0.3393 0.3560 0.4487 0.4549 0.4711 r 0.3288 0.3327 0.3426 0.3830 0.3864 0.3950 e 0.2900 0.2971 0.3149 0.4971 0.5049 0.5244Figure 11
The first to fourth rows are the results reconstructed by FBP, NWL-1, NWL-2, and Cimmino, respectively, corresponding to four kinds of limited angle projection data.Figure 12
The first to third rows are the results reconstructed by NWL-a witha=1.005, a=1.01, and a=1.2, respectively, corresponding to four kinds of limited angle projection data.Figure 13
Relative error histories generated from seven iterative algorithms for the limited angle 0°−89°.Figure 14
Comparisons of the 32nd column in the reconstructed and original images. The red line is reconstruction.Figure 15
Relative error histories generated from seven iterative algorithms for the limited angle 0°−79°.Figure 16
Comparisons of the 32nd column in the reconstructed and original images. The red line is reconstruction.The simulation results show that the proposed weighting method NWL-a (a∈[1.002,1.5]) has better convergence results.


## 4.1. Complete Projection Data

We use180 projections with 95 rays per projection for the reconstruction of complete projection data. The calculated coefficient matrix A has the dimensions 17100×4096. For the algorithm NWL-a, the corresponding values of a and λn are as follows.The largest eigenvalue ofATA is μ1.Fora=1.002, λn=8/μ12; for a=1.005, λn=7.5/μ12; for a=1.01, λn=7/μ12; and for a=1.2, λn=5/μ12.We compare the performance of different algorithms in terms of reconstruction error and image distances. Table1 shows the computation time for completing 200 iterations by formula (90). It can be seen that the time spent by other algorithms is similar except that the Cimmino algorithm takes a little longer time.Table 1
Computation time for completing 200 iterations.
Algorithms Landweber NWL-a Cimmino DROP CAV Time (s) 2.6839 2.3738 14.4995 2.6840 3.6551


### 4.1.1. Noise-Free Projection Data

The comparison of relative errors from the proposed algorithms is shown in Figure1, from which we observe that the relative errors obtained from the NWL-a method are less than those obtained from the Cimmino method, the NWL-1 and NWL-2 methods, and the Landweber method. Figure 2 shows the reconstructed phantom images from the noise-free projections. It can be seen that the reconstructed images obtained from the NWL-a method are clearer than those obtained from the Cimmino method, the NWL-1 and NWL-2 methods, and the Landweber method. The image distances for the noise-free projection data are given in Table 2, from which we see that the image distances from NWL-a algorithm are smaller than those from Cimmino and Landweber algorithms. Figure 3 plots the curves of relative error and it shows that the relative errors from the NWL-a method are less than those from the other methods. Figure 4 shows the comparisons of the 32nd column data in the original and the reconstructed images. It is observed that the main features of the reconstructed images obtained from the NWL-a method with a=1.01 are closer to those of the original image. The image distances generated from different algorithms are given in Table 3, from which we see that the image distances from NWL-a algorithm are smaller than those from DROP and CAV algorithms.Table 2
Image distances from different methods for the noise-free projection data.
Algorithms Landweber NWL-a=1.005 NWL-a=1.01 NWL-a=1.2 Cimmino d 0.4873 0.0039 0.0042 0.0052 0.0351 r 0.8743 0.0028 0.0032 0.0046 0.0386 e 0.1460 0.0033 0.0035 0.0039 0.0064Table 3
Comparison of image distances with DROP and CAV algorithms.
Algorithms NWL-a=1.005 NWL-a=1.01 NWL-a=1.2 DROP CAV d 0.0039 0.0042 0.0052 0.3512 0.3489 r 0.0028 0.0032 0.0046 0.5865 0.5893 e 0.0033 0.0035 0.0039 0.0906 0.0862Figure 1
Relative error histories generated from seven iterative algorithms.Figure 2
(a) Original image. (b)-(h) Images reconstructed by seven iterative algorithms from the noise-free projections.
(a) (b) (c) (d) (e) (f) (g) (h)Figure 3
Relative error histories generated from seven iterative algorithms.Figure 4
Comparisons of the 32nd column in the reconstructed and original images. The red line is reconstruction and blue line represents the original image.


### 4.1.2. Noisy Projection Data

In this section, we test the performance of different algorithms for the projection data with noise. For this purpose, the white Gauss noise with signal-to-noise ratio (SNR) of 10 is added to the projectionsb. The reconstructed results are shown in Figures 5, 6, and 7. The image distances for the noisy projection data are given in Table 4, from which we see that the image distances from NWL-a algorithm are smaller than those from Cimmino and Landweber algorithms. The conclusion is the same as above. The quality of the reconstructed images obtained from the NWL-a method is better than the quality of those obtained from the other methods.Table 4
Image distances from different methods for the noisy projection data.
Algorithms Landweber NWL-a=1.005 NWL-a=1.01 NWL-a=1.2 Cimmino d 0.4954 0.1156 0.1187 0.1199 0.1296 r 0.8745 0.1471 0.1506 0.1517 0.1614 e 0.1617 0.0334 0.0369 0.0394 0.0488Figure 5
(a) Original image. (b)-(h) Images reconstructed by seven iterative algorithms from the projections with noise.
(a) (b) (c) (d) (e) (f) (g) (h)Figure 6
Relative error histories generated from seven iterative algorithms.Figure 7
Comparisons of the 32nd column in the reconstructed and original images. The red line is reconstruction and blue line represents the original image.


## 4.1.1. Noise-Free Projection Data

The comparison of relative errors from the proposed algorithms is shown in Figure1, from which we observe that the relative errors obtained from the NWL-a method are less than those obtained from the Cimmino method, the NWL-1 and NWL-2 methods, and the Landweber method. Figure 2 shows the reconstructed phantom images from the noise-free projections. It can be seen that the reconstructed images obtained from the NWL-a method are clearer than those obtained from the Cimmino method, the NWL-1 and NWL-2 methods, and the Landweber method. The image distances for the noise-free projection data are given in Table 2, from which we see that the image distances from NWL-a algorithm are smaller than those from Cimmino and Landweber algorithms. Figure 3 plots the curves of relative error and it shows that the relative errors from the NWL-a method are less than those from the other methods. Figure 4 shows the comparisons of the 32nd column data in the original and the reconstructed images. It is observed that the main features of the reconstructed images obtained from the NWL-a method with a=1.01 are closer to those of the original image. The image distances generated from different algorithms are given in Table 3, from which we see that the image distances from NWL-a algorithm are smaller than those from DROP and CAV algorithms.Table 2
Image distances from different methods for the noise-free projection data.
Algorithms Landweber NWL-a=1.005 NWL-a=1.01 NWL-a=1.2 Cimmino d 0.4873 0.0039 0.0042 0.0052 0.0351 r 0.8743 0.0028 0.0032 0.0046 0.0386 e 0.1460 0.0033 0.0035 0.0039 0.0064Table 3
Comparison of image distances with DROP and CAV algorithms.
Algorithms NWL-a=1.005 NWL-a=1.01 NWL-a=1.2 DROP CAV d 0.0039 0.0042 0.0052 0.3512 0.3489 r 0.0028 0.0032 0.0046 0.5865 0.5893 e 0.0033 0.0035 0.0039 0.0906 0.0862Figure 1
Relative error histories generated from seven iterative algorithms.Figure 2
(a) Original image. (b)-(h) Images reconstructed by seven iterative algorithms from the noise-free projections.
(a) (b) (c) (d) (e) (f) (g) (h)Figure 3
Relative error histories generated from seven iterative algorithms.Figure 4
Comparisons of the 32nd column in the reconstructed and original images. The red line is reconstruction and blue line represents the original image.


## 4.1.2. Noisy Projection Data

In this section, we test the performance of different algorithms for the projection data with noise. For this purpose, the white Gauss noise with signal-to-noise ratio (SNR) of 10 is added to the projectionsb. The reconstructed results are shown in Figures 5, 6, and 7. The image distances for the noisy projection data are given in Table 4, from which we see that the image distances from NWL-a algorithm are smaller than those from Cimmino and Landweber algorithms. The conclusion is the same as above. The quality of the reconstructed images obtained from the NWL-a method is better than the quality of those obtained from the other methods.Table 4
Image distances from different methods for the noisy projection data.
Algorithms Landweber NWL-a=1.005 NWL-a=1.01 NWL-a=1.2 Cimmino d 0.4954 0.1156 0.1187 0.1199 0.1296 r 0.8745 0.1471 0.1506 0.1517 0.1614 e 0.1617 0.0334 0.0369 0.0394 0.0488Figure 5
(a) Original image. (b)-(h) Images reconstructed by seven iterative algorithms from the projections with noise.
(a) (b) (c) (d) (e) (f) (g) (h)Figure 6
Relative error histories generated from seven iterative algorithms.Figure 7
Comparisons of the 32nd column in the reconstructed and original images. The red line is reconstruction and blue line represents the original image.


## 4.2. Limited Angle Reconstruction

In the theory of continuous case, the problem of limited angle reconstruction has been discussed [17–20]. Assume that f∈L2(R2) is compactly supported and the support of f is supp(f)⊂B={x∈R2:|x|≤c}. In two-dimensional space, S+1={ω=(cos⁡θ,sin⁡θ):0≤θ≤π}. Let Θ0={ω=(cos⁡θ,sin⁡θ):0≤θ≤α<π}. The Radon transform of f is defined by(98)Rfω,p=∫ω·x=pfxds.Let Ω=R×S+1 and Ω0=R×Θ0. If Rf(ω,p) is only known on p∈R and a subset Θ0⊂S+1, the data is called angle-limited. By projection slice theorem,(99)F2fσω=F1Rfσω,and(100)F2fξ=∫Rne-2iπξ·xfxdx.Let PB and PΩ0 be the characteristic functions of B and Ω0, H(ξ)=PΩ0F(ξ). Then(101)Ω0FBf=Gwhere Bf=PBf and Ω0FBf=PΩ0FBf. We let A=Ω0FB. It follows from (101) that the adjoint operator A∗ of A is given by(102)A∗Gx=BF-1Ω0Gx=PBx∫ξ∈Ω0Gξe-i2πξ·xdξ.Multiplying A∗ on both sides of (101), we obtain(103)A∗Af=A∗G.The integral equation (101) can be uniquely solved since A∗A is self-adjoint and positive definite [18]. Therefore the eigenvalues σk of the eigenfunctions of A∗A satisfy 1>σ1≥σ2≥⋯. It turns out that the eigenvalues σk are split into two rather distinct subsets: some are “close to 1” and the rest are “close to 0” [18]. When the range of known limited angles become small, more eigenvalues are “close to 0”, which leads to ill-posed problem becoming serious and degrading the quality of reconstructed image. Here we do not study the theoretical relation between the eigenvalue distributions of derived equation (103) in continuous case and those of discrete algebraic linear systems for limited angle reconstruction problem; we just use our methods to improve the condition number of discrete algebraic linear systems, so as to improve the results of limited angle reconstruction.


### 4.2.1. Limited Angle 0°−119°

To test the effectiveness of the new algorithms in this aspect, image reconstruction using limited angle projection data is simulated in this section. The projection angle is set from 0°−119° and we use 120 projections with 95 rays per projection. The calculated coefficient matrixA has the dimensions 11400×4096. The number of iterations is 1000. For the algorithms NWL-a, the corresponding values of a and λn are as follows. The largest eigenvalue of ATA is μ1. For a=1.005, λn=8/μ12; for a=1.01, λn=7/μ12; and for a=1.2, λn=5/μ12. Completing 1000 iterations, the Cimmino method and the proposed methods take about 15 seconds.From Figure8, we see that the reconstruction images from the NWL-a method are clearer than those obtained from the other methods. Although the results obtained from the Cimmino method are good, it is time-consuming to calculate the iteration matrix Q of Cimmino method before the iteration begins. Figure 9 shows that the NWL-a method makes the relative error smaller. The comparison of the 32nd column data is shown in Figure 10, from which we observe that, for incomplete projection data, the reconstructed image has a certain deviation from the original image. However the reconstructed images obtained from the NWL-a method are the closest to the original image.Figure 8
(a) Original image. (b)-(h) Images reconstructed by seven iterative algorithms for the limited projections 0°−119°.
(a) (b) (c) (d) (e) (f) (g) (h)Figure 9
Relative error histories generated from seven iterative algorithms for the limited angle 0°−119°.Figure 10
Comparisons of the 32nd column in the reconstructed and original images. The red line is reconstruction and blue line represents the original image.


### 4.2.2. Four Kinds of Limited Angle Projection Data

In order to further test the performance of new algorithms on limited angle reconstruction, comparing with the Cimmino method, we select four kinds of limited angle scans which cover0∘-109∘, 0∘-99∘, 0∘-89∘, and 0∘-79∘, respectively. The number of iterations is 1000. Corresponding to limited angle 0∘-109∘, 0∘-99∘, 0∘-89∘, and 0∘-79∘, the calculated coefficient matrices A have the dimensions, respectively, 10450×4096, 9500×4096, 8550×4096, and 7600×4096. For the NWL-a algorithm, the corresponding values of a and λn are as follows. For a=1.005, λn=8/μ1a=1.0052; for a=1.01, λn=7/μ1a=1.012; and for a=1.2, λn=5/μ1a=1.22, where μ1a denotes the largest eigenvalue corresponding to the parameter a.Figure11 shows the reconstructed results by the methods FBP (Filter Backprojection), NWL-1, NWL-2, and Cimmino. Figure 12 shows the reconstructed results from the NWL-a method with a=1.005,a=1.01,a=1.2. From these figures, it is observed that the quality of the reconstructed images becomes worse with the narrowing of the limited angle range. Corresponding to Figure 12, the image distances from the NWL-a method are represented by Tables 5 and 6. It can be seen that the image distances d, r, and e obtained by NWL-a method with a = 1.005, 1.01, 1.2 are similar for each kind of limited angle projection data. The values of d, r, and e become larger with the narrowing of the limited angle range, which is consistent with the fact. However the reconstructions of the main features in the original image are clear when the range of the limited angle is at least 90∘. In contrast, the reconstructed images obtained via the NWL-a method are more clear. The reconstructed images are not clear when the range of the limited angle is less than 90∘. Figures 13 and 15 show the comparison of relative errors from the six iterative methods. It is observed that the convergence speed of the NWL-a method is the fastest. From Figures 14 and 16, we see that the reconstructed image obtained from the NWL-a method is closer to the original image. The numerical simulations of limited angle reconstruction show that the NWL-a method has improved the reconstruction results.Table 5
Image distances from the NWL-a method for the limited angle projection data0∘-109∘ and 0∘-99∘.
Limited angle 0 ∘ - 10 9 ∘ 0 ∘ - 9 9 ∘ NWL-a a=1.005 a=1.01 a=1.2 a=1.005 a=1.01 a=1.2 d 0.2112 0.2172 0.2328 0.2269 0.2338 0.2521 r 0.2267 0.2317 0.2443 0.2386 0.2441 0.2584 e 0.1904 0.1979 0.2165 0.1963 0.2033 0.2208Table 6
Image distances from the NWL-a method for the limited angle projection data0∘-89∘ and 0∘-79∘.
Limited angle 0 ∘ - 8 9 ∘ 0 ∘ - 7 9 ∘ NWL-a a=1.005 a=1.01 a=1.2 a=1.005 a=1.01 a=1.2 d 0.3328 0.3393 0.3560 0.4487 0.4549 0.4711 r 0.3288 0.3327 0.3426 0.3830 0.3864 0.3950 e 0.2900 0.2971 0.3149 0.4971 0.5049 0.5244Figure 11
The first to fourth rows are the results reconstructed by FBP, NWL-1, NWL-2, and Cimmino, respectively, corresponding to four kinds of limited angle projection data.Figure 12
The first to third rows are the results reconstructed by NWL-a witha=1.005, a=1.01, and a=1.2, respectively, corresponding to four kinds of limited angle projection data.Figure 13
Relative error histories generated from seven iterative algorithms for the limited angle 0°−89°.Figure 14
Comparisons of the 32nd column in the reconstructed and original images. The red line is reconstruction.Figure 15
Relative error histories generated from seven iterative algorithms for the limited angle 0°−79°.Figure 16
Comparisons of the 32nd column in the reconstructed and original images. The red line is reconstruction.The simulation results show that the proposed weighting method NWL-a (a∈[1.002,1.5]) has better convergence results.


## 4.2.1. Limited Angle 0°−119°

To test the effectiveness of the new algorithms in this aspect, image reconstruction using limited angle projection data is simulated in this section. The projection angle is set from 0°−119° and we use 120 projections with 95 rays per projection. The calculated coefficient matrixA has the dimensions 11400×4096. The number of iterations is 1000. For the algorithms NWL-a, the corresponding values of a and λn are as follows. The largest eigenvalue of ATA is μ1. For a=1.005, λn=8/μ12; for a=1.01, λn=7/μ12; and for a=1.2, λn=5/μ12. Completing 1000 iterations, the Cimmino method and the proposed methods take about 15 seconds.From Figure8, we see that the reconstruction images from the NWL-a method are clearer than those obtained from the other methods. Although the results obtained from the Cimmino method are good, it is time-consuming to calculate the iteration matrix Q of Cimmino method before the iteration begins. Figure 9 shows that the NWL-a method makes the relative error smaller. The comparison of the 32nd column data is shown in Figure 10, from which we observe that, for incomplete projection data, the reconstructed image has a certain deviation from the original image. However the reconstructed images obtained from the NWL-a method are the closest to the original image.Figure 8
(a) Original image. (b)-(h) Images reconstructed by seven iterative algorithms for the limited projections 0°−119°.
(a) (b) (c) (d) (e) (f) (g) (h)Figure 9
Relative error histories generated from seven iterative algorithms for the limited angle 0°−119°.Figure 10
Comparisons of the 32nd column in the reconstructed and original images. The red line is reconstruction and blue line represents the original image.


## 4.2.2. Four Kinds of Limited Angle Projection Data

In order to further test the performance of new algorithms on limited angle reconstruction, comparing with the Cimmino method, we select four kinds of limited angle scans which cover0∘-109∘, 0∘-99∘, 0∘-89∘, and 0∘-79∘, respectively. The number of iterations is 1000. Corresponding to limited angle 0∘-109∘, 0∘-99∘, 0∘-89∘, and 0∘-79∘, the calculated coefficient matrices A have the dimensions, respectively, 10450×4096, 9500×4096, 8550×4096, and 7600×4096. For the NWL-a algorithm, the corresponding values of a and λn are as follows. For a=1.005, λn=8/μ1a=1.0052; for a=1.01, λn=7/μ1a=1.012; and for a=1.2, λn=5/μ1a=1.22, where μ1a denotes the largest eigenvalue corresponding to the parameter a.Figure11 shows the reconstructed results by the methods FBP (Filter Backprojection), NWL-1, NWL-2, and Cimmino. Figure 12 shows the reconstructed results from the NWL-a method with a=1.005,a=1.01,a=1.2. From these figures, it is observed that the quality of the reconstructed images becomes worse with the narrowing of the limited angle range. Corresponding to Figure 12, the image distances from the NWL-a method are represented by Tables 5 and 6. It can be seen that the image distances d, r, and e obtained by NWL-a method with a = 1.005, 1.01, 1.2 are similar for each kind of limited angle projection data. The values of d, r, and e become larger with the narrowing of the limited angle range, which is consistent with the fact. However the reconstructions of the main features in the original image are clear when the range of the limited angle is at least 90∘. In contrast, the reconstructed images obtained via the NWL-a method are more clear. The reconstructed images are not clear when the range of the limited angle is less than 90∘. Figures 13 and 15 show the comparison of relative errors from the six iterative methods. It is observed that the convergence speed of the NWL-a method is the fastest. From Figures 14 and 16, we see that the reconstructed image obtained from the NWL-a method is closer to the original image. The numerical simulations of limited angle reconstruction show that the NWL-a method has improved the reconstruction results.Table 5
Image distances from the NWL-a method for the limited angle projection data0∘-109∘ and 0∘-99∘.
Limited angle 0 ∘ - 10 9 ∘ 0 ∘ - 9 9 ∘ NWL-a a=1.005 a=1.01 a=1.2 a=1.005 a=1.01 a=1.2 d 0.2112 0.2172 0.2328 0.2269 0.2338 0.2521 r 0.2267 0.2317 0.2443 0.2386 0.2441 0.2584 e 0.1904 0.1979 0.2165 0.1963 0.2033 0.2208Table 6
Image distances from the NWL-a method for the limited angle projection data0∘-89∘ and 0∘-79∘.
Limited angle 0 ∘ - 8 9 ∘ 0 ∘ - 7 9 ∘ NWL-a a=1.005 a=1.01 a=1.2 a=1.005 a=1.01 a=1.2 d 0.3328 0.3393 0.3560 0.4487 0.4549 0.4711 r 0.3288 0.3327 0.3426 0.3830 0.3864 0.3950 e 0.2900 0.2971 0.3149 0.4971 0.5049 0.5244Figure 11
The first to fourth rows are the results reconstructed by FBP, NWL-1, NWL-2, and Cimmino, respectively, corresponding to four kinds of limited angle projection data.Figure 12
The first to third rows are the results reconstructed by NWL-a witha=1.005, a=1.01, and a=1.2, respectively, corresponding to four kinds of limited angle projection data.Figure 13
Relative error histories generated from seven iterative algorithms for the limited angle 0°−89°.Figure 14
Comparisons of the 32nd column in the reconstructed and original images. The red line is reconstruction.Figure 15
Relative error histories generated from seven iterative algorithms for the limited angle 0°−79°.Figure 16
Comparisons of the 32nd column in the reconstructed and original images. The red line is reconstruction.The simulation results show that the proposed weighting method NWL-a (a∈[1.002,1.5]) has better convergence results.


## 5. Conclusion

Based on the analysis of the relation between the condition number and the spectral radius of the iterative matrix, the new weighted matrices have been derived to improve the ill conditioned number. For the higher resolution of reconstructed image, the approximate iterative matrix has been derived for the weighted Landweber method with large dimensional number. According to these, the new weighting algorithm NWL-a and corresponding relaxation strategies for the weighted Landweber method are proposed. Numerical simulations verify the feasibility and the better performance of NWL-a method for complete projection data, especially for limited angle projection data. Fora∈[1.002,1.5], the reconstruction results obtained by NWL-a are similar. Therefore, when using the NWL-a method, a number can be selected from the interval as the value of the parameter a. Since we have no imaging laboratory here, the proposed algorithm cannot be applied to actual experiments, which will be a part of our future work.



---

*Source: 5674647-2018-07-18.xml*