# A Doubly Constrained TV Algorithm for Image Reconstruction

Authors: Zhiwei Qiao; Gage Redler; Boris Epel; Howard Halpern
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

Purpose. The total variation (TV) minimization algorithm is an effective image reconstruction algorithm capable of accurately reconstructing images from sparse and/or noisy data. The TV model consists of two terms: a data fidelity term and a TV regularization term. Two constrained TV models, data divergence-constrained TV minimization (DDcTV) and TV-constrained data divergence minimization (TVcDM), have been successfully applied to computed tomography (CT) and electron paramagnetic resonance imaging (EPRI). In this work, we propose a new constrained TV model, a doubly constrained TV (dcTV) model, which has the potential to further improve the reconstruction accuracy for the two terms which are both of constraint forms. Methods. We perform an inverse crime study to validate the model and its Chambolle-Pock (CP) solver and characterize the performance of the dcTV-CP algorithm in the context of CT. To demonstrate the superiority of the dcTV model, we compare the convergence rate and the reconstruction accuracy with the DDcTV and TVcDM models via simulated data. Results and Conclusions. The performance-characterizing study shows that the dcTV-CP algorithm is an accurate and convergent algorithm, with the model parameters impacting the reconstruction accuracy and the algorithm parameters impacting the convergence path and rate. The comparison studies show that the dcTV-CP algorithm has a relatively fast convergence rate and can achieve higher reconstruction accuracy from sparse projections or noisy projections relative to the other two single-constrained TV models. The knowledge and insights gained in the work may be utilized in the application of the new model in other imaging modalities including divergence-beam CT, magnetic resonance imaging (MRI), positron emission tomography (PET), and EPRI.


---

## Body

## 1. Introduction

Image reconstruction algorithms are critical components of accurate medical imaging [1]. This includes various tomographic imaging modalities such as computed tomography (CT), magnetic resonance imaging (MRI), positron emission tomography (PET), and electron paramagnetic resonance imaging (EPRI) [2]. There are mainly two reconstruction algorithm frameworks, analytic and iterative algorithms [1]. Before 2006, analytic algorithms were the mainstream approach used in commercial imagers. However, in the wake of increased interest in compressed sensing (CS) [3] and the extensive use of graphics processing units (GPU) [4], iterative algorithms gained popularity.In fact, iterative algorithms are based on a discrete-to-discrete (DD) imaging model, i.e., the reconstructed image and the measured data are both discrete. The DD imaging model is a linear system of equations with the coefficient matrix being the system matrix that represents the forward imaging process. Usually, the linear system is large scale, ill-posed, and underdetermined; therefore, the direct inversion of the linear system is impossible. To achieve accurate reconstruction, the reconstructed image can be formulated as the solution of an optimization model, in which prior information may be included. The CS-based optimization model uses the sparsity prior knowledge since the reconstructed image always has a specific sparse transform [5]. The total variation (TV) minimization algorithm is a classical optimization-based algorithm and has been widely and successfully used in medical imaging [6, 7] due to its capability to accurately reconstruct images from sparse or noisy data [8, 9].The DD imaging model can be formulated as(1)g=Au,where u is a vector of size N, representing the reconstructed image; g is a vector of size M, representing the measured discrete data; and A is a matrix of size M×N, representing the system matrix. For 2D parallel-beam CT, A indicates the 2D Radon transform [10]; for 2D fan-beam CT or 3D cone-beam CT, A indicates the ray transform [11]; for MRI, A indicates the Fourier transform [12], and for 3D EPRI, A indicates the 3D Radon transform [13].The TV model consists of two terms, the data fidelity term and the TV regularization term. Four specific TV models of different constraint forms can be formulated as follows:(2)uo=argminu12g−Au22+wuTV,(3)uo=argminuuTV,s.t.g−Au2≤ϵ,(4)uo=argminug−Au22,s.t.uTV≤t1,(5)uo=argminu0,s.t.uTV≤t1andg−Au2≤ϵ.(2) is the unconstrained TV (ucTV) model; (3) is the data divergence-constrained TV (DDcTV) minimization model; (4) is the TV-constrained data divergence minimization (TVcDM) model; and (5) is the newly proposed doubly constrained TV (dcTV) model.In equation (2), 1/2g−Au22 is the data fidelity term and uTV is the TV regularization term. w is the model parameter that balances the importance of data fidelity and TV regularization. This is a frequently used TV model, especially in CT [14–17]. However, this model has a disadvantage in that the model parameter w does not have any physical meaning and therefore is hard to tune for specific reconstruction tasks. When w is very small, the convergence rate will be very slow [18].In equation (3), ϵ is the data tolerance bound. The DDcTV model has an advantage relative to ucTV in that the model parameter ϵ has a clear physical meaning, as it indicates the noise level of the data and the inconsistence level of the linear system. The DDcTV model was first proposed by Sidky et al. in 2006 [8]. An algorithm for solving the model is the famous adaptive steepest descent projection onto convex sets (ASD-POCS) algorithm [9]. From 2006 to 2015, the ASD-POCS algorithm has been extensively studied and its capability to accurately reconstruct images from sparse-view projections or noisy projections in CBCT [19], offset-detector CBCT [11], and EPRI [13] has been demonstrated.In equation (4), t1 is the TV bound. This model parameter is not as sensitive as that in DDcTV, which is why the TVcDM was later proposed as an improved model in 2014 [20] and has since been the subject of extensive studies. The model has been successful in short scan CBCT [21], C-arm CBCT [22], positron emission tomography (PET) [23], and EPRI [24].In equation (5), there are two constraint terms, the data divergence constraint and the TV constraint. Each constraint defines a convex set. The objective function is a 0 function, which has no substantial optimization meaning, and we just use it to formulate the optimization form. Clearly, the solution of the dcTV model is any point in the intersection set of the two convex sets. In fact, this type of optimization problem is the so-called convex feasibility problem [25–27]. The convex feasibility problem may be inconsistent or consistent. If all the convex sets have intersection set, it is consistent. Otherwise, it is inconsistent. By selecting appropriate model parameters, we can ensure model consistency. Most solving algorithms of the convex feasibility problem are based on projection onto convex sets (POCS). From the initial solution estimate, it is sequentially projected onto each convex set. POCS may always achieve a convergent, useful solution. For the consistent case, the POCS algorithm stops when the moving point enters the intersection set. For the inconsistent case, the POCS algorithm may converge to a point that is close enough to each convex set. Clearly, the solution is not unique. But each solution of the convex feasibility problem in the solution set is considered to be equivalent. So, though the specific solution depends on the initial solution estimate and the solving algorithm, it does not impact the utility of the convex feasibility problem.There are two model parameters,t1 and ϵ, in the dcTV model, which control the TV bound and data tolerance bound. In the DDcTV model, only the data tolerance bound is used; thus, the reconstructed image tends to be overly smooth because of the minimization of TV if the data tolerance bound is comparatively large to correspond to higher noise. In the TVcDM model, only TV bound is used; thus, the reconstructed image tends to be too noisy because of the minimization of the data fidelity term if the projections have comparatively high noise. By comparison, the dcTV model uses both model parameters to control the TV value and the data tolerance, respectively, thus resulting in enhanced capability to appropriately balance these effects and provide more accurate image reconstruction. Using these double constraints is equivalent to using two aspects of quantitative prior knowledge. This is the rationale behind why the proposed dcTV model may be superior to the DDcTV and TVcDM models. We may also use Figure 1 to explain the possible superiority of the dcTV model.Figure 1
Solution relationship schematic diagram of the three constrained TV models. The upper ellipse of sizeϵ indicates the convex set determined by the data fidelity constraint. The lower ellipse of size t1 indicates the convex set determined by the TV constraint. The two red points indicate two of the solutions of the dcTV model, the blue point indicates the solution of the TVcDM model, and the green point indicates the solution of the DDcTV model.We name the upper ellipse asϵ-ellipse that indicates the convex set determined by the data fidelity constraint and name the lower ellipse as the t1-ellipse that indicates the convex set determined by the TV constraint.For the DDcTV model, the solution is the image with the minimal TV value in theϵ-ellipse. Clearly, the green point indicates the solution for it has the minimal TV value. For the TVcDM model, the solution is the image with the minimal data divergence in the t1-ellipse. Clearly, the blue point indicates the solution for it has the minimal data divergence. If the ϵ and t1 are both the optimal value, i.e., ϵ=g−Autruth2 and t1=utruthTV. It is clear that the blue point suffers from too small data divergence and the green point suffers from too small TV value. If we may get the tradeoff of the two points, the solution should be more accurate. The dcTV model has the solution set, i.e., the intersection set of the two ellipses, in which each point is a solution of the model and its data divergence is not so small as that of the TVcDM model and its TV value is not so small as that of the DDcTV model. Clearly, the dcTV model has potential to achieve higher accuracy for it may achieve the balance between noise and smoothness.Selection of an algorithm for solving the above TV models is difficult because of the nonsmooth TV term and the large-scale/ill-posed linear system. However, designing an appropriate algorithm for solving these models plays a critical role in successfully applying these models for image reconstruction. There are currently three main solving algorithms: ASD-POCS, alternating direction method of multipliers (ADMM) [14] (split Bregman [28]), and Chambolle-Pock (CP) algorithm [18, 29–33]. ASD-POCS has been used to solve the DDcTV model; ADMM or split Bregman are usually used to solve the ucTV model; the CP algorithm has been used to solve DDcTV and TVcDM models [18, 20–23].The CP algorithm has certain advantages relative to the ADMM or ASD-POCS algorithms: (1) it may solve all the convex optimization models whether or not they are smooth; (2) each subproblem has a closed form; and (3) all of the algorithm parameters in the CP algorithm may be explicitly determined analytically (i.e., via solving equations) rather than empirically selecting suitable values [31]. So, we select CP algorithm to be the solver of the newly proposed dcTV model.In this work, we perform two studies: (1) characterization of the dcTV-CP algorithm performance and (2) comparison of the three constrained TV algorithms (dcTV-CP, DDcTV-CP, and TVcDM-CP algorithms). In the performance characterization study, we conduct an inverse crime study to validate the algorithm and investigate the impact of the model parameters on reconstruction accuracy and the impact of algorithm parameters on convergence rate. In the comparison study, we compare the convergence rate and the reconstruction accuracy of the three algorithms.There are two related works that should be mentioned here. One is on the first-order convex feasibility algorithms for X-ray CT [34], in which the third model, inequality constrained, TV (ICTV) model, is similar to our dcTV model. However, this work just focused on the acceleration effect of the CP2-ICTV algorithm and the iteration property of the algorithm. The other one is regarding the POCS algorithm for solving the multiconstrained TV models [35]. However, this work just focused on the effectiveness of the POCS algorithm. The main aims of our work are different from these two works. We focus on the performance characterization of the dcTV-CP algorithm and the evaluation of the potential superiority of the new algorithm relative to the other two existing constrained TV algorithms, the DDcTV-CP and TVcDM-CP algorithm. This work may continue to contribute to this type of TV algorithm.In Section2, we describe the methods according to the chain: imaging system modeling, optimization program design, CP algorithm instance derivation, model parameter selection, convergence condition metrics, and reconstruction accuracy metrics. In Section 3, we present the algorithm performance characterization and comparison studies. A brief discussion of this work is given in Section 4.


## 2. Methods

This work focuses on a novel dcTV model and its CP solving algorithm. The terms dcTV-CP and dcTV will be used interchangeably. Without loss of generality, we evaluate the algorithm in the context of 2D parallel CT. In this section, we introduce the imaging system model, optimization program design, reconstruction parameters, and metrics to evaluate algorithm convergence and reconstructed image quality.


### 2.1. Imaging System Modeling

The imaging system model of the 2D parallel CT is a linear system shown in equation (1). The system matrix A represents the forward imaging process, i.e., the 2D Radon transform. Here, the system matrix may also be referred to as a projection matrix, since the projection data g is the multiplication of the projection matrix A and the image u. Therefore, modeling of the imaging system, which is the calculation of the system matrix, is also the modeling of the projection method. There are three common projection methods: pixel-driven, ray-driven [36], and distance-driven [37] projection methods. In this work, we use the novel accurate pixel-driven method proposed by our group [38]. The linear system is large scale, ill-posed and often underdetermined, so it is almost impossible to solve it by calculating the pseudoinverse of A. So we continue to model the imaging system as an optimization model.


### 2.2. dcTV Model: An Optimization Model

To accurately reconstruct an image according to the imaging model, one may design the optimization model by incorporating useful prior information. In this work, we propose a new TV model, the dcTV model shown in equation (5). The data fidelity term and the TV regularization term are both of constraint forms. However, the two corresponding convex functions can have very different magnitudes. This may lead to slow convergence rate. In order to avoid this, two algorithm parameters are incorporated into the model to balance their relative magnitudes, resulting in the optimization program given in the following equation:(6)uo=argminu0,s.t.νuTV≤νt1andλg−Au2≤λϵ.The solution to this model is a point in the intersection set of the two convex sets corresponding to the two constraint terms. Note that the two algorithm parameters do not impact the size of the convex sets and therefore do not impact the solution of the model. However, they may impact the convergence rate of the solving algorithm.


### 2.3. dcTV-CP Algorithm

The CP algorithm framework is a powerful optimization framework. The CP algorithm instances of many optimization models, especially TV models, have been derived and have had their convergence validated. The derivation process is more difficult than that of the ADMM algorithm instances for the calculations of convex conjugate functions and the proximal mapping of convex functions. For that reason, we give the detailed derivation of the dcTV-CP algorithm instance as follows.


#### 2.3.1. CP Algorithm Framework

Here, we summarize the CP algorithm framework briefly, as it will be used for the derivation of the dcTV-CP algorithm instance. The CP algorithm framework solves optimization models of the form given in the following equation:(7)xo=argminxFy+Gx,s.t.y=Kx,where x and y are the vectors in vector spaces X and Y; K is a linear transform, indicating a linear mapping from X to Y; F and G are two convex functions of y and x, respectively. The CP algorithm framework requires F and G to be convex but does not require them to be smooth.The CP algorithm framework is shown in Algorithm1 [18].Algorithm 1: Pseudocode for N steps of the CP algorithm framework.
(1)
L
=
K
SV
;
σ
=
1
/
L
;
τ
=
1
/
L
;
θ
=
1
;
n
=
0
(2)
x
0
=
0
;
x
¯
0
=
0
;
y
0
=
0
(3)
repeat
(4)
yn+1=proxσF∗yn+σKx¯n
(5)
xn+1=proxτGxn−τKTyn+1
(6)
x¯n+1=xn+1+θxn+1−xn
(7)
n=n+1
(8)
until n≥NIn Algorithm1, KSV is the largest singular value of matrix K. F∗ is the convex conjugate of F. proxσF∗ and proxτG are two proximal operations. The superscript T indicates the matrix transpose.The convex conjugate of functionHz may be defined as(8)H∗z=maxczTc−Hc.The proximal mapping may be defined as(9)proxaHz=argmincHc+z−c222a.The key to deriving the CP algorithm instance of a specific optimization model is the calculation of the convex conjugate functionF∗ and the two proximal mapping operations, proxσF∗ and proxτG.


#### 2.3.2. Derivation of the dcTV-CP Algorithm Instance

First, we define two indicator functions as follows:(10)δℓ2Ballax≡0,x2≤a,∞,x2>a,(11)δℓ1Ballax≡0,x1≤a,∞,x1>a.An indicator function is equivalent to a convex set constraint, and thus, a constrained optimization model can be effectively converted into an unconstrained one. The TV term in equation (6) is defined as(12)uTV=Dumag1,where D is a matrix of size 2N×N, indicating the discrete gradient transform, and is of the form,(13)D=D1D2,where, D1 and D2 are both matrices of size N×N, indicating the gradient transform along the x- and y-axes, respectively. Assuming that ux,y, where x∈1,nx,y∈1,ny, indicates a pixel of the 2D image; then, the D1 transform is given by(14)D1ux,y=ux,y−ux−1,y,x∈2,nx,0,x=1.And the D2 transform is given by(15)D2ux,y=ux,y−ux,y−1,y∈2,ny,0,y=1.⋅mag in equation (12) is the magnitude of a 2D vector and may be defined with ℓ1 norm or ℓ2 norm. In this work, we use the ℓ2 norm to define ⋅mag as(16)abmag=ab2=a2+b2.Thus,Du is a vector of size 2N, representing the discrete gradient transform of the image u, whereas Dumag is a vector of size N, representing the isotropic gradient magnitude transform of the image, u. So, the TV norm defined in equation (6) is the isotropic TV norm.According to (10), (11), and (12), the dcTV model may be reformulated into an unconstrained form as follows:(17)uo=argminu0+δℓ1Ballνt1νDumag+δℓ2Ballλϵλg−λAu.To derive the CP algorithm instance, we make the following mechanical associations with equation (7):(18)x=u,p=λAu,q=νDu,K=λAνD,y=pq,Fy=F1p+F2q,F1p=δℓ2Ballλϵλg−p,F2q=δℓ1Ballνt1qmag,Gx=0.According to equations (8) and (9), we get(19)F1∗p=λϵp2+p,λg,(20)proxσF1∗p=maxp−σλg2−σλϵ,0×p−σλgp−σλg2.In equation (19), ⋅,⋅ indicates the inner product of vectors. In equation (20), the right hand side is the shrinkage operation or soft threshold operation. The detailed derivation of equations (19) and (20) are not shown here, but similar derivations may be found in Appendix C of reference [18].According to equations (8) and (9), we get(21)F2∗q=νt1qmag∞,(22)proxσF2∗q=q1I−σProjectOntoℓ1Ballνt1qmag/σqmag.In equation (21), ⋅∞ is the ∞-norm of a vector, which selects the largest element in the vector. In equation (22), ProjectOntoℓ1Ballax is a projection operator which may project a vector x onto the ℓ1 ball of radius a [20, 22]. 1I is a “1” vector in space I=RN.According to equation (9), we get(23)proxτGx=x.Substituting equations (18), (20), (22), and (23) into Algorithm 1, we get Algorithm 2, the CP algorithm instance of the dcTV model.Algorithm 2: Pseudocode for (N) steps of the dcTV-CP algorithm.
INPUT: g, A, D,ϵ, t1; λ, b
(1)
ν
A
=
A
sv
/
D
sv; ν=bνA
(2)
L
=
λ
A
νD
sv; σ=1/L,τ=1/L, θ=1,n=0
(3)
u
0
=
0
;
p
0
=
0
;
q
0
=
0
;
u
¯
0
=
0
(4)
Repeat
(5)
a=pn+σλAu¯n−g
(6)
pn+1=maxa2−σλϵ,0×a/a2
(7)
c=qn+σνDu¯n
(8)
s=ProjectOntoℓ1Ballνt1cmag/σ
(9)
q
n
+
1
=
c
1
I
−
σ
s
/
c
mag
(10)
u
n
+
1
=
u
n
−
τ
λ
A
T
p
n
+
1
−
τ
ν
D
T
q
n
+
1
(11)
u
¯
n
+
1
=
u
n
+
1
+
θ
u
n
+
1
−
u
n
(12)
n
=
n
+
1
(13)
Until n≥N
OUTPUT: The designed solution uN.In Algorithm2, u,u¯∈RN, p,a∈RM, and q,c∈R2N. In line 9, 1I∈RN is a ‘1’ vector in image space I. In line 2, the calculation algorithm of L is shown in Algorithm 8 of reference [18]. The algorithm for ProjectOntoℓ1Ball in line 8 is shown in Algorithm 2 of reference [20].


### 2.4. Reconstruction Parameters

The dcTV-CP algorithm involves many parameters, which may be divided into two types: model parameters and algorithm parameters. Model parameters are those parameters that decide the final solution, whereas algorithm parameters are those parameters that cannot impact the final solution but may impact the convergence rate and path.Model parameters of the dcTV-CP algorithm shown in Algorithm2 include pixel size, projection method, data tolerance bound, and TV bound. In this work, we use a normalized DD model, so the pixel size is 1. The projection method used is the accurate pixel driven method proposed by our group. The data tolerance bound and TV bound are two important model parameters in the dcTV model since they determine the size of the convex sets defined by the data fidelity term and TV regularization term. Thus, they have an important impact on the reconstruction accuracy. We will study the approach to selecting these two model parameters in Section 3.Algorithm parameters consist of two types: the introduced algorithm parameters (e.g.,λandν) and the fixed algorithm parameters in the CP algorithm (e.g., σ,τ,andθ). The three fixed algorithm parameters may be determined by equations in line 2 of Algorithm 2. The two introduced algorithm parameters may balance the magnitude of the two convex functions corresponding to the data fidelity term and TV regularization term. Selection of these parameters impacts the convergence rate of the iteration process. We will study the approach to selecting these parameters in Section 3.


### 2.5. Metrics for Practical Convergence Conditions

#### 2.5.1. For Inverse Crime Study

Due to the fact that, by definition, for inverse crime studies, the truth image exists, the linear system is consistent and the projection data are ideal and sufficient, we make use of the three metrics as follows (equations (24)–(26)) to set up the practical convergence condition:(24)NOEun=RMSEun=un−utruth2N,(25)NDEun=g−Aun2g2,(26)NTVEun=unTV−utruthTVutruthTV.The normalized object errorNOE, describes the distance between the reconstructed image and the truth image. The normalized data error, NDE, describes the distance between the projection data generated from the reconstructed image and the measured projection data. The normalized TV error, NTVE, describes the distance between the TV value of the reconstructed image and the truth image TV value. RMSE stands for the root mean square error.For the inverse crime study, we can simply set the convergence conditions such that all the three metrics must become negligibly small.


#### 2.5.2. For the Nonideal Data Reconstruction

Except for the special case of the inverse crime study, all image reconstructions are nonideal. The nonideal factors include insufficient projection data (i.e., sparse data) and noisy projection data. Even for simulation reconstructions from nonideal data, where the truth image exists, the three metrics shown in equations (24)–(26) cannot be used to formulate practical convergence conditions because we do not know what values to expect for those metrics. Thus, the practical convergence conditions may be set up based on the flatness of each metric plotted versus iteration number, shown in the following equations:(27)dNOEun=un−uref2−un−1−uref2uref2,(28)dNDEun=g−Aun2−g−Aun−12g2,(29)dNTVEun=unTV−un−1TVunTV.The differential normalized object error, dNOE, indicates the flatness of the NOE iteration curve. The differential normalized data error, dNDE, indicates the flatness of the NDE iteration curve. And the differential normalized TV error, dNTVE, indicates the flatness of the NTVE iteration curve.uref means the reference image. For simulation reconstruction, it is the truth image, whereas for real-data reconstruction, it may be the image reconstructed by analytic algorithm or other high-quality reference image.


### 2.6. Metrics for Image Quality Evaluation

There exists many metrics to assess image quality. In this work, we choose to use NOE, i.e. RMSE, of the reconstructed image relative to the reference image, to be the metric:(30)NOEuconv=uconv−uref2N.Here,uconv means the convergent solution.


## 2.1. Imaging System Modeling

The imaging system model of the 2D parallel CT is a linear system shown in equation (1). The system matrix A represents the forward imaging process, i.e., the 2D Radon transform. Here, the system matrix may also be referred to as a projection matrix, since the projection data g is the multiplication of the projection matrix A and the image u. Therefore, modeling of the imaging system, which is the calculation of the system matrix, is also the modeling of the projection method. There are three common projection methods: pixel-driven, ray-driven [36], and distance-driven [37] projection methods. In this work, we use the novel accurate pixel-driven method proposed by our group [38]. The linear system is large scale, ill-posed and often underdetermined, so it is almost impossible to solve it by calculating the pseudoinverse of A. So we continue to model the imaging system as an optimization model.


## 2.2. dcTV Model: An Optimization Model

To accurately reconstruct an image according to the imaging model, one may design the optimization model by incorporating useful prior information. In this work, we propose a new TV model, the dcTV model shown in equation (5). The data fidelity term and the TV regularization term are both of constraint forms. However, the two corresponding convex functions can have very different magnitudes. This may lead to slow convergence rate. In order to avoid this, two algorithm parameters are incorporated into the model to balance their relative magnitudes, resulting in the optimization program given in the following equation:(6)uo=argminu0,s.t.νuTV≤νt1andλg−Au2≤λϵ.The solution to this model is a point in the intersection set of the two convex sets corresponding to the two constraint terms. Note that the two algorithm parameters do not impact the size of the convex sets and therefore do not impact the solution of the model. However, they may impact the convergence rate of the solving algorithm.


## 2.3. dcTV-CP Algorithm

The CP algorithm framework is a powerful optimization framework. The CP algorithm instances of many optimization models, especially TV models, have been derived and have had their convergence validated. The derivation process is more difficult than that of the ADMM algorithm instances for the calculations of convex conjugate functions and the proximal mapping of convex functions. For that reason, we give the detailed derivation of the dcTV-CP algorithm instance as follows.


### 2.3.1. CP Algorithm Framework

Here, we summarize the CP algorithm framework briefly, as it will be used for the derivation of the dcTV-CP algorithm instance. The CP algorithm framework solves optimization models of the form given in the following equation:(7)xo=argminxFy+Gx,s.t.y=Kx,where x and y are the vectors in vector spaces X and Y; K is a linear transform, indicating a linear mapping from X to Y; F and G are two convex functions of y and x, respectively. The CP algorithm framework requires F and G to be convex but does not require them to be smooth.The CP algorithm framework is shown in Algorithm1 [18].Algorithm 1: Pseudocode for N steps of the CP algorithm framework.
(1)
L
=
K
SV
;
σ
=
1
/
L
;
τ
=
1
/
L
;
θ
=
1
;
n
=
0
(2)
x
0
=
0
;
x
¯
0
=
0
;
y
0
=
0
(3)
repeat
(4)
yn+1=proxσF∗yn+σKx¯n
(5)
xn+1=proxτGxn−τKTyn+1
(6)
x¯n+1=xn+1+θxn+1−xn
(7)
n=n+1
(8)
until n≥NIn Algorithm1, KSV is the largest singular value of matrix K. F∗ is the convex conjugate of F. proxσF∗ and proxτG are two proximal operations. The superscript T indicates the matrix transpose.The convex conjugate of functionHz may be defined as(8)H∗z=maxczTc−Hc.The proximal mapping may be defined as(9)proxaHz=argmincHc+z−c222a.The key to deriving the CP algorithm instance of a specific optimization model is the calculation of the convex conjugate functionF∗ and the two proximal mapping operations, proxσF∗ and proxτG.


### 2.3.2. Derivation of the dcTV-CP Algorithm Instance

First, we define two indicator functions as follows:(10)δℓ2Ballax≡0,x2≤a,∞,x2>a,(11)δℓ1Ballax≡0,x1≤a,∞,x1>a.An indicator function is equivalent to a convex set constraint, and thus, a constrained optimization model can be effectively converted into an unconstrained one. The TV term in equation (6) is defined as(12)uTV=Dumag1,where D is a matrix of size 2N×N, indicating the discrete gradient transform, and is of the form,(13)D=D1D2,where, D1 and D2 are both matrices of size N×N, indicating the gradient transform along the x- and y-axes, respectively. Assuming that ux,y, where x∈1,nx,y∈1,ny, indicates a pixel of the 2D image; then, the D1 transform is given by(14)D1ux,y=ux,y−ux−1,y,x∈2,nx,0,x=1.And the D2 transform is given by(15)D2ux,y=ux,y−ux,y−1,y∈2,ny,0,y=1.⋅mag in equation (12) is the magnitude of a 2D vector and may be defined with ℓ1 norm or ℓ2 norm. In this work, we use the ℓ2 norm to define ⋅mag as(16)abmag=ab2=a2+b2.Thus,Du is a vector of size 2N, representing the discrete gradient transform of the image u, whereas Dumag is a vector of size N, representing the isotropic gradient magnitude transform of the image, u. So, the TV norm defined in equation (6) is the isotropic TV norm.According to (10), (11), and (12), the dcTV model may be reformulated into an unconstrained form as follows:(17)uo=argminu0+δℓ1Ballνt1νDumag+δℓ2Ballλϵλg−λAu.To derive the CP algorithm instance, we make the following mechanical associations with equation (7):(18)x=u,p=λAu,q=νDu,K=λAνD,y=pq,Fy=F1p+F2q,F1p=δℓ2Ballλϵλg−p,F2q=δℓ1Ballνt1qmag,Gx=0.According to equations (8) and (9), we get(19)F1∗p=λϵp2+p,λg,(20)proxσF1∗p=maxp−σλg2−σλϵ,0×p−σλgp−σλg2.In equation (19), ⋅,⋅ indicates the inner product of vectors. In equation (20), the right hand side is the shrinkage operation or soft threshold operation. The detailed derivation of equations (19) and (20) are not shown here, but similar derivations may be found in Appendix C of reference [18].According to equations (8) and (9), we get(21)F2∗q=νt1qmag∞,(22)proxσF2∗q=q1I−σProjectOntoℓ1Ballνt1qmag/σqmag.In equation (21), ⋅∞ is the ∞-norm of a vector, which selects the largest element in the vector. In equation (22), ProjectOntoℓ1Ballax is a projection operator which may project a vector x onto the ℓ1 ball of radius a [20, 22]. 1I is a “1” vector in space I=RN.According to equation (9), we get(23)proxτGx=x.Substituting equations (18), (20), (22), and (23) into Algorithm 1, we get Algorithm 2, the CP algorithm instance of the dcTV model.Algorithm 2: Pseudocode for (N) steps of the dcTV-CP algorithm.
INPUT: g, A, D,ϵ, t1; λ, b
(1)
ν
A
=
A
sv
/
D
sv; ν=bνA
(2)
L
=
λ
A
νD
sv; σ=1/L,τ=1/L, θ=1,n=0
(3)
u
0
=
0
;
p
0
=
0
;
q
0
=
0
;
u
¯
0
=
0
(4)
Repeat
(5)
a=pn+σλAu¯n−g
(6)
pn+1=maxa2−σλϵ,0×a/a2
(7)
c=qn+σνDu¯n
(8)
s=ProjectOntoℓ1Ballνt1cmag/σ
(9)
q
n
+
1
=
c
1
I
−
σ
s
/
c
mag
(10)
u
n
+
1
=
u
n
−
τ
λ
A
T
p
n
+
1
−
τ
ν
D
T
q
n
+
1
(11)
u
¯
n
+
1
=
u
n
+
1
+
θ
u
n
+
1
−
u
n
(12)
n
=
n
+
1
(13)
Until n≥N
OUTPUT: The designed solution uN.In Algorithm2, u,u¯∈RN, p,a∈RM, and q,c∈R2N. In line 9, 1I∈RN is a ‘1’ vector in image space I. In line 2, the calculation algorithm of L is shown in Algorithm 8 of reference [18]. The algorithm for ProjectOntoℓ1Ball in line 8 is shown in Algorithm 2 of reference [20].


## 2.3.1. CP Algorithm Framework

Here, we summarize the CP algorithm framework briefly, as it will be used for the derivation of the dcTV-CP algorithm instance. The CP algorithm framework solves optimization models of the form given in the following equation:(7)xo=argminxFy+Gx,s.t.y=Kx,where x and y are the vectors in vector spaces X and Y; K is a linear transform, indicating a linear mapping from X to Y; F and G are two convex functions of y and x, respectively. The CP algorithm framework requires F and G to be convex but does not require them to be smooth.The CP algorithm framework is shown in Algorithm1 [18].Algorithm 1: Pseudocode for N steps of the CP algorithm framework.
(1)
L
=
K
SV
;
σ
=
1
/
L
;
τ
=
1
/
L
;
θ
=
1
;
n
=
0
(2)
x
0
=
0
;
x
¯
0
=
0
;
y
0
=
0
(3)
repeat
(4)
yn+1=proxσF∗yn+σKx¯n
(5)
xn+1=proxτGxn−τKTyn+1
(6)
x¯n+1=xn+1+θxn+1−xn
(7)
n=n+1
(8)
until n≥NIn Algorithm1, KSV is the largest singular value of matrix K. F∗ is the convex conjugate of F. proxσF∗ and proxτG are two proximal operations. The superscript T indicates the matrix transpose.The convex conjugate of functionHz may be defined as(8)H∗z=maxczTc−Hc.The proximal mapping may be defined as(9)proxaHz=argmincHc+z−c222a.The key to deriving the CP algorithm instance of a specific optimization model is the calculation of the convex conjugate functionF∗ and the two proximal mapping operations, proxσF∗ and proxτG.


## 2.3.2. Derivation of the dcTV-CP Algorithm Instance

First, we define two indicator functions as follows:(10)δℓ2Ballax≡0,x2≤a,∞,x2>a,(11)δℓ1Ballax≡0,x1≤a,∞,x1>a.An indicator function is equivalent to a convex set constraint, and thus, a constrained optimization model can be effectively converted into an unconstrained one. The TV term in equation (6) is defined as(12)uTV=Dumag1,where D is a matrix of size 2N×N, indicating the discrete gradient transform, and is of the form,(13)D=D1D2,where, D1 and D2 are both matrices of size N×N, indicating the gradient transform along the x- and y-axes, respectively. Assuming that ux,y, where x∈1,nx,y∈1,ny, indicates a pixel of the 2D image; then, the D1 transform is given by(14)D1ux,y=ux,y−ux−1,y,x∈2,nx,0,x=1.And the D2 transform is given by(15)D2ux,y=ux,y−ux,y−1,y∈2,ny,0,y=1.⋅mag in equation (12) is the magnitude of a 2D vector and may be defined with ℓ1 norm or ℓ2 norm. In this work, we use the ℓ2 norm to define ⋅mag as(16)abmag=ab2=a2+b2.Thus,Du is a vector of size 2N, representing the discrete gradient transform of the image u, whereas Dumag is a vector of size N, representing the isotropic gradient magnitude transform of the image, u. So, the TV norm defined in equation (6) is the isotropic TV norm.According to (10), (11), and (12), the dcTV model may be reformulated into an unconstrained form as follows:(17)uo=argminu0+δℓ1Ballνt1νDumag+δℓ2Ballλϵλg−λAu.To derive the CP algorithm instance, we make the following mechanical associations with equation (7):(18)x=u,p=λAu,q=νDu,K=λAνD,y=pq,Fy=F1p+F2q,F1p=δℓ2Ballλϵλg−p,F2q=δℓ1Ballνt1qmag,Gx=0.According to equations (8) and (9), we get(19)F1∗p=λϵp2+p,λg,(20)proxσF1∗p=maxp−σλg2−σλϵ,0×p−σλgp−σλg2.In equation (19), ⋅,⋅ indicates the inner product of vectors. In equation (20), the right hand side is the shrinkage operation or soft threshold operation. The detailed derivation of equations (19) and (20) are not shown here, but similar derivations may be found in Appendix C of reference [18].According to equations (8) and (9), we get(21)F2∗q=νt1qmag∞,(22)proxσF2∗q=q1I−σProjectOntoℓ1Ballνt1qmag/σqmag.In equation (21), ⋅∞ is the ∞-norm of a vector, which selects the largest element in the vector. In equation (22), ProjectOntoℓ1Ballax is a projection operator which may project a vector x onto the ℓ1 ball of radius a [20, 22]. 1I is a “1” vector in space I=RN.According to equation (9), we get(23)proxτGx=x.Substituting equations (18), (20), (22), and (23) into Algorithm 1, we get Algorithm 2, the CP algorithm instance of the dcTV model.Algorithm 2: Pseudocode for (N) steps of the dcTV-CP algorithm.
INPUT: g, A, D,ϵ, t1; λ, b
(1)
ν
A
=
A
sv
/
D
sv; ν=bνA
(2)
L
=
λ
A
νD
sv; σ=1/L,τ=1/L, θ=1,n=0
(3)
u
0
=
0
;
p
0
=
0
;
q
0
=
0
;
u
¯
0
=
0
(4)
Repeat
(5)
a=pn+σλAu¯n−g
(6)
pn+1=maxa2−σλϵ,0×a/a2
(7)
c=qn+σνDu¯n
(8)
s=ProjectOntoℓ1Ballνt1cmag/σ
(9)
q
n
+
1
=
c
1
I
−
σ
s
/
c
mag
(10)
u
n
+
1
=
u
n
−
τ
λ
A
T
p
n
+
1
−
τ
ν
D
T
q
n
+
1
(11)
u
¯
n
+
1
=
u
n
+
1
+
θ
u
n
+
1
−
u
n
(12)
n
=
n
+
1
(13)
Until n≥N
OUTPUT: The designed solution uN.In Algorithm2, u,u¯∈RN, p,a∈RM, and q,c∈R2N. In line 9, 1I∈RN is a ‘1’ vector in image space I. In line 2, the calculation algorithm of L is shown in Algorithm 8 of reference [18]. The algorithm for ProjectOntoℓ1Ball in line 8 is shown in Algorithm 2 of reference [20].


## 2.4. Reconstruction Parameters

The dcTV-CP algorithm involves many parameters, which may be divided into two types: model parameters and algorithm parameters. Model parameters are those parameters that decide the final solution, whereas algorithm parameters are those parameters that cannot impact the final solution but may impact the convergence rate and path.Model parameters of the dcTV-CP algorithm shown in Algorithm2 include pixel size, projection method, data tolerance bound, and TV bound. In this work, we use a normalized DD model, so the pixel size is 1. The projection method used is the accurate pixel driven method proposed by our group. The data tolerance bound and TV bound are two important model parameters in the dcTV model since they determine the size of the convex sets defined by the data fidelity term and TV regularization term. Thus, they have an important impact on the reconstruction accuracy. We will study the approach to selecting these two model parameters in Section 3.Algorithm parameters consist of two types: the introduced algorithm parameters (e.g.,λandν) and the fixed algorithm parameters in the CP algorithm (e.g., σ,τ,andθ). The three fixed algorithm parameters may be determined by equations in line 2 of Algorithm 2. The two introduced algorithm parameters may balance the magnitude of the two convex functions corresponding to the data fidelity term and TV regularization term. Selection of these parameters impacts the convergence rate of the iteration process. We will study the approach to selecting these parameters in Section 3.


## 2.5. Metrics for Practical Convergence Conditions

### 2.5.1. For Inverse Crime Study

Due to the fact that, by definition, for inverse crime studies, the truth image exists, the linear system is consistent and the projection data are ideal and sufficient, we make use of the three metrics as follows (equations (24)–(26)) to set up the practical convergence condition:(24)NOEun=RMSEun=un−utruth2N,(25)NDEun=g−Aun2g2,(26)NTVEun=unTV−utruthTVutruthTV.The normalized object errorNOE, describes the distance between the reconstructed image and the truth image. The normalized data error, NDE, describes the distance between the projection data generated from the reconstructed image and the measured projection data. The normalized TV error, NTVE, describes the distance between the TV value of the reconstructed image and the truth image TV value. RMSE stands for the root mean square error.For the inverse crime study, we can simply set the convergence conditions such that all the three metrics must become negligibly small.


### 2.5.2. For the Nonideal Data Reconstruction

Except for the special case of the inverse crime study, all image reconstructions are nonideal. The nonideal factors include insufficient projection data (i.e., sparse data) and noisy projection data. Even for simulation reconstructions from nonideal data, where the truth image exists, the three metrics shown in equations (24)–(26) cannot be used to formulate practical convergence conditions because we do not know what values to expect for those metrics. Thus, the practical convergence conditions may be set up based on the flatness of each metric plotted versus iteration number, shown in the following equations:(27)dNOEun=un−uref2−un−1−uref2uref2,(28)dNDEun=g−Aun2−g−Aun−12g2,(29)dNTVEun=unTV−un−1TVunTV.The differential normalized object error, dNOE, indicates the flatness of the NOE iteration curve. The differential normalized data error, dNDE, indicates the flatness of the NDE iteration curve. And the differential normalized TV error, dNTVE, indicates the flatness of the NTVE iteration curve.uref means the reference image. For simulation reconstruction, it is the truth image, whereas for real-data reconstruction, it may be the image reconstructed by analytic algorithm or other high-quality reference image.


## 2.5.1. For Inverse Crime Study

Due to the fact that, by definition, for inverse crime studies, the truth image exists, the linear system is consistent and the projection data are ideal and sufficient, we make use of the three metrics as follows (equations (24)–(26)) to set up the practical convergence condition:(24)NOEun=RMSEun=un−utruth2N,(25)NDEun=g−Aun2g2,(26)NTVEun=unTV−utruthTVutruthTV.The normalized object errorNOE, describes the distance between the reconstructed image and the truth image. The normalized data error, NDE, describes the distance between the projection data generated from the reconstructed image and the measured projection data. The normalized TV error, NTVE, describes the distance between the TV value of the reconstructed image and the truth image TV value. RMSE stands for the root mean square error.For the inverse crime study, we can simply set the convergence conditions such that all the three metrics must become negligibly small.


## 2.5.2. For the Nonideal Data Reconstruction

Except for the special case of the inverse crime study, all image reconstructions are nonideal. The nonideal factors include insufficient projection data (i.e., sparse data) and noisy projection data. Even for simulation reconstructions from nonideal data, where the truth image exists, the three metrics shown in equations (24)–(26) cannot be used to formulate practical convergence conditions because we do not know what values to expect for those metrics. Thus, the practical convergence conditions may be set up based on the flatness of each metric plotted versus iteration number, shown in the following equations:(27)dNOEun=un−uref2−un−1−uref2uref2,(28)dNDEun=g−Aun2−g−Aun−12g2,(29)dNTVEun=unTV−un−1TVunTV.The differential normalized object error, dNOE, indicates the flatness of the NOE iteration curve. The differential normalized data error, dNDE, indicates the flatness of the NDE iteration curve. And the differential normalized TV error, dNTVE, indicates the flatness of the NTVE iteration curve.uref means the reference image. For simulation reconstruction, it is the truth image, whereas for real-data reconstruction, it may be the image reconstructed by analytic algorithm or other high-quality reference image.


## 2.6. Metrics for Image Quality Evaluation

There exists many metrics to assess image quality. In this work, we choose to use NOE, i.e. RMSE, of the reconstructed image relative to the reference image, to be the metric:(30)NOEuconv=uconv−uref2N.Here,uconv means the convergent solution.


## 3. Results

We perform two studies: one to characterize the performance of the dcTV-CP algorithm and one to compare the dcTV-CP, DDcTV-CP, and TVcDM-CP algorithms.


### 3.1. dcTV-CP Algorithm Performance Characterization

This section includes 4 parts: (1) an inverse crime study, (2) determination of the impact of algorithm parameters on convergence rate, (3) determination of the impact of the TV bound on reconstruction accuracy, and (4) determination of the impact of the data tolerance bound on reconstruction accuracy.


#### 3.1.1. Inverse Crime Study

Image reconstruction is a classical inverse problem. Inverse crime studies are a commonly used tool for end-to-end validation of a given approach to inverting an inverse problem. The whole reconstruction chain of the dcTV-CP algorithm includes imaging system modeling, optimization model design, solving algorithm design, and its computer implementation. The inverse crime is considered successful if the reconstructed image is accurate enough such that any residual errors are due to computer float error. For 8-bit grayscale images, there are only 256 possible discrete intensity values to assign to a given pixel. Therefore, a successful inverse crime study in this case will produce a reconstructed image with NOE relative to the truth image being less than or equal to10−4.The phantom used for the inverse crime study is the Shepp–Logan phantom of size256×256. The parallel scanning configuration is used. The projection data are also of size 256×256, which means there are 256 projections evenly distributed in the angle range of 0,πand there are 256 measurements on each projection. The origin of the imaging coordinate system is located at 128,128 of the phantom. The axis range of each projection signal is −128:127. The projections are generated by use of Eq. (1); thus ,the linear system is completely consistent.The system matrix is determined by use of the accurate pixel driven projection method. Since the projection data are ideal, the model parameters are set asϵ=0 and t1=utruthTV. To get fast convergence, the two algorithm parameters are set as λ=1andν=0.1νA.The inverse crime sign, i.e., the convergence condition, is set as(31)NOEun≤10−4,NDEun≤10−4,NTVEun≤10−3.The dcTV-CP algorithm for the study stops at iteration2910. Figures 2(a) and 2(b) show the truth image and the reconstructed image, respectively. It can be seen that the two images are visually identical. Figure 2(c) shows the profiles of the vertical center-line of the images. It shows that the reconstructed profile completely overlaps with the truth profile. At the convergence point, NOE<10−4. Both the qualitative and the quantitative evaluations show that the inverse crime succeeded. This means that the design and implementation of the dcTV-CP algorithm are both correct.Figure 2
dcTV-CP reconstruction images of the Shepp–Logan phantom. (a) The Shepp–Logan phantom image (truth). (b) The reconstructed image. (c) Comparison of the profiles on the vertical center-line of the phantom image and the reconstructed image.
(a)
(b)
(c)The iteration trends of the three metrics for the dcTV-CP algorithm are shown in Figures3(a)–3(c), respectively.Figure 3
Convergence behavior of the dcTV-CP algorithm for reconstructing the Shepp–Logan phantom. (a) The object error NOE. (b) The data error NDE. (c) The relative TV error NTVE, each as a function of the iteration number.
(a)
(b)
(c)From Figure3, one may see that the three metrics decrease monotonically with progressive iterations. Even at the convergence point, the three metrics still maintain the decreasing trend. This shows that the dcTV-CP algorithm is convergent. It is noted that all the iteration curves exhibit a slight oscillatory behavior, which can also be observed in other CP algorithms, for example, DDcTV-CP and TVcDM-CP algorithms. Maybe, this is the iteration law of all the CP algorithms. However, this does not impact the overall downward trend towards convergence.


#### 3.1.2. Impact of Algorithm Parameters on Convergence Rate

The two algorithm parameters in the dcTV-CP algorithm areλandν. To investigate their impact on the convergence rate, we varied λandν and compared convergence rates.The Shepp–Logan phantom was used for this study as well. The imaging conditions and parameter selection are the same as those in the inverse crime study, except that the number of projections is 60 and the iteration number is 1500. We used three different combinations of values forλandν :λ=1andν=0.1νA, λ=1andν=νA, and λ=1andν=10νA.Figure4 shows the error images between the reconstructed images and the Shepp–Logan phantom. Figure 5 shows the iteration trend of NOE for the image reconstructions with the three different parameter combinations. From Figure 4, based on error-image observation after a fixed number of iterations, it can be seen that λ=1andν=0.1νA leads to the fastest convergence, whereas λ=1andν=10νA leads to the slowest convergence. This is further evidenced in Figure 5.Figure 4
Error image between the reconstructed images and the truth image. The error image is the absolute value of the difference between the reconstructed image and the truth image. (a)λ=1andν=0.1νA. (b) λ=1andν=νA. (c) λ=1andν=10νA. The display window is [0, 0.01].
(a)
(b)
(c)Figure 5
Iteration curves of NOE for the dcTV-CP algorithms using different combinations of values forλandν.These results show that the two algorithm parameters have significant impact on the convergence rate. For a specific image reconstruction application, one may run the reconstruction algorithm with a different combination of parameter values and select the most appropriate values that provide rapid convergence for that specific application.


#### 3.1.3. Impact of TV Bound on Reconstruction Accuracy

There are two model parameters, the TV bound and the data tolerance bound, whose selection may impact the reconstruction accuracy.For the TVcDM-CP algorithm and DDcTV algorithm, it is known how the model parameters impact the reconstruction accuracy. For the DDcTV algorithm, the model parameter is the data tolerance bound. In reconstruction from noisy data, a large data tolerance value can lead to an overly smooth image because an image with a smaller TV value can be obtained in larger convex set defined by the data fidelity term. Conversely, a small data tolerance value can lead to a noisy image because the noise information in the projection data will be propagated into the reconstructed image. For the TVcDM-CP algorithm, the TV bound is the model parameter. In noisy data reconstruction, too large of a TV bound leads to a noisy image because the TV noise-reduction effect is suppressed, whereas too small of a TV bound leads to an overly smooth image because low TV value corresponds to a smooth image.In order to isolate the effect of each model parameter for the dcTV-CP algorithm, their impact on reconstruction is studied separately.As before, the Shepp–Logan phantom is of size256×256; the projection data are of size 256×60. Gaussian noise with magnitude of 45 dB is added to the projections, resulting in ϵtruth=24.59. The model parameter ϵ is set to be ϵtruth, and the model parameter t1 is set to be 0.6, 0.8, 1, 1.2, and 1.4 times of utruth  TV, to determine how it may impact reconstruction accuracy. We set λ=1andν=0.1νA to get fast convergence. The practical convergence conditions are given in the following equation:(32)dNOEun≤10−3,dNDEun≤10−3,dNTVEun≤10−3.The reconstructed images and the error images are shown in Figures6 and 7, respectively, and the plot of NOE as function of t1 is shown in Figure 8.Figure 6
Images reconstructed using the dcTV-CP algorithm and varyingt1. (a–e) t1 of values 0.6, 0.8, 1, 1.2 and 1.4 times utruthTV. The display window is [0, 1].
(a)
(b)
(c)
(d)
(e)Figure 7
Reconstructed error images corresponding to the images in Figure6. The display window is [0, 0.05].
(a)
(b)
(c)
(d)
(e)Figure 8
Plot of NOE as function of multiple oft1 relative to the truth TV.From Figure6, it can be seen that the image becomes smoother when the TV bound is decreased, i.e., more noisy when the TV bound is increased. It appears that a multiple of 1 may achieve the best tradeoff. The error images in Figure 7 are displayed with a narrow intensity window to further demonstrate this phenomenon. It can be clearly seen that a multiple of 1 achieves the best tradeoff. Figure 8 shows the NOE of the five reconstructed images. A multiple of 1 corresponds to the lowest reconstruction error. Both the qualitative and quantitative analysis, therefore, show that a multiple of 1 may achieve the most accurate reconstruction.


#### 3.1.4. Impact of Data Tolerance Bound on Reconstruction Accuracy

The impact of the data tolerance bound on reconstruction accuracy is further investigated by repeating the reconstructions in the last section except now;t1 is fixed as the truth TV and ϵ is varied as 0.1, 0.2, 0.3, 0.4, 0.5 0.6, 0.8, 1, 1.2 and 1.4 times of ϵtruth.For briefness, we just show the plot of NOE of the reconstructed images as function ofϵ in Figure 9.Figure 9
Plot of NOE as function of multiple ofϵ relative to the truth ϵ.From Figure9, it can be seen that 0.4 times of the truth data tolerance bound achieve the lowest reconstruction error and, therefore, the most accurate reconstruction. The rationale of this phenomenon is that a large data tolerance bound corresponds to a large convex set defined by the data fidelity term; therefore, a smooth image may be selected in the set. A small data tolerance bound corresponds to a small convex set; therefore, the reconstructed image tends to be noisy due to the existence of noise in the projections.


### 3.2. Comparison of the dcTV-CP, DDcTV-CP, and TVcDM-CP Algorithms

In this section, three comparisons are made between the dcTV-CP, DDcTV-CP, and TVcDM-CP algorithms: comparison of convergence rate, comparison of reconstruction capability from sparse projections, and comparison of reconstruction capability from noisy projections. The DDcTV-CP algorithm is shown in Algorithm3, and the TVcDM-CP algorithm is shown in Algorithm 4.Algorithm 3: Pseudocode for N steps of the DDcTV-CP algorithm instance.
INPUT: g, A, Dϵ; λ, b
(1)
ν
A
=
A
sv
/
D
sv; ν=bνA
(2)
L
=
λ
A
νD
sv; σ=1/L,τ=1/L, θ=1,n=0
(3)
u
0
=
0
;
p
0
=
0
;
q
0
=
0
;
u
¯
0
=
0
(4)
Repeat
(5)
a
=
p
n
+
σ
λ
A
u
¯
n
−
g
(6)
p
n
+
1
=
max
a
2
−
σ
λ
ϵ
,
0
×
a
/
a
2
(7)
c
=
q
n
+
σ
ν
D
u
¯
n
(8)
q
n
+
1
=
c
/
max
1
I
,
c
mag
(9)
u
n
+
1
=
u
n
−
τ
λ
A
T
p
n
+
1
−
τ
ν
D
T
q
n
+
1
(10)
u
¯
n
+
1
=
u
n
+
1
+
θ
u
n
+
1
−
u
n
(11)
n
=
n
+
1
(12)
Until n≥N
OUTPUT: The designed solution uN.Algorithm 4: Pseudocode for N steps of the TVcDM-CP algorithm instance.
INPUT: g, A, D, t1; λ, b
(1)
ν
A
=
A
S
V
/
D
S
V
,
ν
=
b
ν
A
(2)
L
=
A
ν
D
S
V; σ=1/L;τ=1/L;θ=1;n=0
(3)
u
0
=
0
;
u
¯
0
=
0
;
p
0
=
0
;
q
0
=
0
(4)
Repeat
(5)
p
n
+
1
=
p
n
+
σ
A
u
¯
n
−
g
/
1
+
σ
/
λ
(6)
a=qn+σνDu¯n
(7)
s=ProjectOntoℓ1Ballνt1amag/σ
(8)
qn+1=a1I−σs/amag
(9)
u
n
+
1
=
u
n
−
τ
A
T
p
n
+
1
+
ν
D
T
q
n
+
1
(10)
u
¯
n
+
1
=
u
n
+
1
+
θ
u
n
+
1
−
u
n
(11)
n
=
n
+
1
(12)
Until n≥N
OUTPUT: The designed solution uN.


#### 3.2.1. Convergence Rate Comparison

We use the FORBILD phantom to perform this comparison. The imaging conditions and parameters are the same as those used in the inverse crime study except that the number of projections is 60. The convergence condition is NOE≤5×10−4.The convergence curves of NOE are shown in Figure10 for the FORBILD reconstruction. It can be seen that the TVcDM-CP algorithm converges the fastest, the DDcTV-CP algorithm converges the slowest, and the newly proposed dcTV-CP algorithm has a slightly slower convergence rate than the TVcDM-CP algorithm.Figure 10
Plot of NOE as function of iteration number for the three algorithms.


#### 3.2.2. Comparison of Reconstruction Accuracy from Sparse Projections

TV-based algorithms can accurately reconstruct images from sparse-view projections since this type of optimization model embodies the idea of compressed sensing. In this section, we will compare the sparse reconstruction capability of the three different TV algorithms. Each algorithm has a different optimization meaning dictated by their different constraint forms. So they may perform different sparse reconstruction capabilities.To demonstrate the capability of each algorithm to handle sparse projection data, images are reconstructed from 10, 20, 30, 40, and 50 projections. The imaging conditions and parameters are the same as in the inverse crime study except that we use the different projection numbers mentioned above. The FORBILD phantom is used to perform this study.Figures11 and 12 show the reconstructed FORBILD images and the error images, respectively. Figure 13 shows the plot of NOE as function of the projection number. From Figure 11, it can be seen that all the three TV algorithms achieve accurate sparse reconstruction, as it is hard to tell the difference between the images reconstructed from 20, 30, 40, and 50 projections and the truth image for any TV algorithm (see images from the second to the fifth column). From the images reconstructed from 10 projections, we may see that the DDcTV image is too smooth, whereas the TVcDM and dcTV images have clear streak artifacts. From Figure 13, we may see that the order of reconstruction accuracy from high to low is dcTV, TVcDM, and then DDcTV if ≥20 projections are used. From Figure 12, it can be seen that the dcTV error image reconstructed from 40 projections is subtle, whereas the features in the DDcTV error image and TVcDM error image are so pronounced that the structure of the phantom may be distinguished. This analysis shows that the dcTV-CP algorithm can achieve the highest reconstruction accuracy.Figure 11
FORBILD phantom images reconstructed using the DDcTV-CP (the first row), TVcDM-CP (the second row), and dcTV-CP (the third row) algorithms from 10, 20, 30, 40, and 50 projections. The display window is [0, 1].Figure 12
FORBILD phantom reconstructed error images from 40 projections using the DDcTV (the first column), TVcDM (the second column), and dcTV (the third column) algorithms. The display window is [0, 0.01].
(a)
(b)
(c)Figure 13
Plots of NOE as function of the projection number.


#### 3.2.3. Comparison of Reconstruction Accuracy from Noisy Projections

TV algorithms can achieve accurate reconstructions not only from sparse-view projections but also from noisy projections. In fact, TV models have been used in image denoising before its wide application in image reconstruction from 2006. The capability to accurately reconstruct images from noisy projection data is very useful for low-dose CT, in which the dose for each view is reduced, resulting in increased projection noise. The DDcTV and TVcDM models have been deeply investigated in low-dose CT and have shown their capability of suppressing the associated noise. Here, we compare the proposed dcTV-CP algorithm’s capability to reconstruct images from noisy projections with that of the established DDcTV and TVcDM models.The FORBILD phantom is used. The imaging conditions and parameters are the same as those in the last section except that the number of projections is fixed at 60 and the SNR of the noisy projections is varied from 30 dB, 35 dB, 40 dB, 45 dB, to 50 dB.Figures14 and 15 show the reconstructed images and the corresponding error images, respectively, and Figure 16 shows the plot of NOE of the reconstructed images as function of SNR.Figure 14
Reconstructed FORBILD images using the DDcTV-CP (the first row), TVcDM-CP (the second row), and dcTV-CP (the third row) algorithms from projections with SNR of 30 dB, 35 dB, 40 dB, 45 dB, and 50 dB. The display window is [0, 1].Figure 15
FORBILD error images reconstructed from projections with SNR of 30 dB using the DDcTV (the first column), TVcDM (the second column), and dcTV (the third column) algorithms. The display window is [0, 0.01].Figure 16
Plots of NOE as a function of projection SNR.From Figure14, it can be seen that the three TV algorithms all have the capability to accurately reconstruct images from projections with SNR of 35 dB to 50 dB. Subtle differences between the three images reconstructed from projections with SNR of 30 dB can be observed: the DDcTV image is smoother, the TVcDM image is noisier, but the dcTV image appears to balance these two characteristics well (i.e., not as smooth as the DDcTV image and not as noisy as the TVcDM image). This can be clearly seen in Figure 15, where the DDcTV error image exhibits strong features near high-frequency edges in the phantom, the TVcDM error image has a large degree of random differences, but the dcTV error image achieves the balance between these two. This is evidenced quantitatively in Figure 16, where the dcTV-CP algorithm achieves the highest reconstruction accuracy in the context of image reconstruction from noisy projections.


## 3.1. dcTV-CP Algorithm Performance Characterization

This section includes 4 parts: (1) an inverse crime study, (2) determination of the impact of algorithm parameters on convergence rate, (3) determination of the impact of the TV bound on reconstruction accuracy, and (4) determination of the impact of the data tolerance bound on reconstruction accuracy.


### 3.1.1. Inverse Crime Study

Image reconstruction is a classical inverse problem. Inverse crime studies are a commonly used tool for end-to-end validation of a given approach to inverting an inverse problem. The whole reconstruction chain of the dcTV-CP algorithm includes imaging system modeling, optimization model design, solving algorithm design, and its computer implementation. The inverse crime is considered successful if the reconstructed image is accurate enough such that any residual errors are due to computer float error. For 8-bit grayscale images, there are only 256 possible discrete intensity values to assign to a given pixel. Therefore, a successful inverse crime study in this case will produce a reconstructed image with NOE relative to the truth image being less than or equal to10−4.The phantom used for the inverse crime study is the Shepp–Logan phantom of size256×256. The parallel scanning configuration is used. The projection data are also of size 256×256, which means there are 256 projections evenly distributed in the angle range of 0,πand there are 256 measurements on each projection. The origin of the imaging coordinate system is located at 128,128 of the phantom. The axis range of each projection signal is −128:127. The projections are generated by use of Eq. (1); thus ,the linear system is completely consistent.The system matrix is determined by use of the accurate pixel driven projection method. Since the projection data are ideal, the model parameters are set asϵ=0 and t1=utruthTV. To get fast convergence, the two algorithm parameters are set as λ=1andν=0.1νA.The inverse crime sign, i.e., the convergence condition, is set as(31)NOEun≤10−4,NDEun≤10−4,NTVEun≤10−3.The dcTV-CP algorithm for the study stops at iteration2910. Figures 2(a) and 2(b) show the truth image and the reconstructed image, respectively. It can be seen that the two images are visually identical. Figure 2(c) shows the profiles of the vertical center-line of the images. It shows that the reconstructed profile completely overlaps with the truth profile. At the convergence point, NOE<10−4. Both the qualitative and the quantitative evaluations show that the inverse crime succeeded. This means that the design and implementation of the dcTV-CP algorithm are both correct.Figure 2
dcTV-CP reconstruction images of the Shepp–Logan phantom. (a) The Shepp–Logan phantom image (truth). (b) The reconstructed image. (c) Comparison of the profiles on the vertical center-line of the phantom image and the reconstructed image.
(a)
(b)
(c)The iteration trends of the three metrics for the dcTV-CP algorithm are shown in Figures3(a)–3(c), respectively.Figure 3
Convergence behavior of the dcTV-CP algorithm for reconstructing the Shepp–Logan phantom. (a) The object error NOE. (b) The data error NDE. (c) The relative TV error NTVE, each as a function of the iteration number.
(a)
(b)
(c)From Figure3, one may see that the three metrics decrease monotonically with progressive iterations. Even at the convergence point, the three metrics still maintain the decreasing trend. This shows that the dcTV-CP algorithm is convergent. It is noted that all the iteration curves exhibit a slight oscillatory behavior, which can also be observed in other CP algorithms, for example, DDcTV-CP and TVcDM-CP algorithms. Maybe, this is the iteration law of all the CP algorithms. However, this does not impact the overall downward trend towards convergence.


### 3.1.2. Impact of Algorithm Parameters on Convergence Rate

The two algorithm parameters in the dcTV-CP algorithm areλandν. To investigate their impact on the convergence rate, we varied λandν and compared convergence rates.The Shepp–Logan phantom was used for this study as well. The imaging conditions and parameter selection are the same as those in the inverse crime study, except that the number of projections is 60 and the iteration number is 1500. We used three different combinations of values forλandν :λ=1andν=0.1νA, λ=1andν=νA, and λ=1andν=10νA.Figure4 shows the error images between the reconstructed images and the Shepp–Logan phantom. Figure 5 shows the iteration trend of NOE for the image reconstructions with the three different parameter combinations. From Figure 4, based on error-image observation after a fixed number of iterations, it can be seen that λ=1andν=0.1νA leads to the fastest convergence, whereas λ=1andν=10νA leads to the slowest convergence. This is further evidenced in Figure 5.Figure 4
Error image between the reconstructed images and the truth image. The error image is the absolute value of the difference between the reconstructed image and the truth image. (a)λ=1andν=0.1νA. (b) λ=1andν=νA. (c) λ=1andν=10νA. The display window is [0, 0.01].
(a)
(b)
(c)Figure 5
Iteration curves of NOE for the dcTV-CP algorithms using different combinations of values forλandν.These results show that the two algorithm parameters have significant impact on the convergence rate. For a specific image reconstruction application, one may run the reconstruction algorithm with a different combination of parameter values and select the most appropriate values that provide rapid convergence for that specific application.


### 3.1.3. Impact of TV Bound on Reconstruction Accuracy

There are two model parameters, the TV bound and the data tolerance bound, whose selection may impact the reconstruction accuracy.For the TVcDM-CP algorithm and DDcTV algorithm, it is known how the model parameters impact the reconstruction accuracy. For the DDcTV algorithm, the model parameter is the data tolerance bound. In reconstruction from noisy data, a large data tolerance value can lead to an overly smooth image because an image with a smaller TV value can be obtained in larger convex set defined by the data fidelity term. Conversely, a small data tolerance value can lead to a noisy image because the noise information in the projection data will be propagated into the reconstructed image. For the TVcDM-CP algorithm, the TV bound is the model parameter. In noisy data reconstruction, too large of a TV bound leads to a noisy image because the TV noise-reduction effect is suppressed, whereas too small of a TV bound leads to an overly smooth image because low TV value corresponds to a smooth image.In order to isolate the effect of each model parameter for the dcTV-CP algorithm, their impact on reconstruction is studied separately.As before, the Shepp–Logan phantom is of size256×256; the projection data are of size 256×60. Gaussian noise with magnitude of 45 dB is added to the projections, resulting in ϵtruth=24.59. The model parameter ϵ is set to be ϵtruth, and the model parameter t1 is set to be 0.6, 0.8, 1, 1.2, and 1.4 times of utruth  TV, to determine how it may impact reconstruction accuracy. We set λ=1andν=0.1νA to get fast convergence. The practical convergence conditions are given in the following equation:(32)dNOEun≤10−3,dNDEun≤10−3,dNTVEun≤10−3.The reconstructed images and the error images are shown in Figures6 and 7, respectively, and the plot of NOE as function of t1 is shown in Figure 8.Figure 6
Images reconstructed using the dcTV-CP algorithm and varyingt1. (a–e) t1 of values 0.6, 0.8, 1, 1.2 and 1.4 times utruthTV. The display window is [0, 1].
(a)
(b)
(c)
(d)
(e)Figure 7
Reconstructed error images corresponding to the images in Figure6. The display window is [0, 0.05].
(a)
(b)
(c)
(d)
(e)Figure 8
Plot of NOE as function of multiple oft1 relative to the truth TV.From Figure6, it can be seen that the image becomes smoother when the TV bound is decreased, i.e., more noisy when the TV bound is increased. It appears that a multiple of 1 may achieve the best tradeoff. The error images in Figure 7 are displayed with a narrow intensity window to further demonstrate this phenomenon. It can be clearly seen that a multiple of 1 achieves the best tradeoff. Figure 8 shows the NOE of the five reconstructed images. A multiple of 1 corresponds to the lowest reconstruction error. Both the qualitative and quantitative analysis, therefore, show that a multiple of 1 may achieve the most accurate reconstruction.


### 3.1.4. Impact of Data Tolerance Bound on Reconstruction Accuracy

The impact of the data tolerance bound on reconstruction accuracy is further investigated by repeating the reconstructions in the last section except now;t1 is fixed as the truth TV and ϵ is varied as 0.1, 0.2, 0.3, 0.4, 0.5 0.6, 0.8, 1, 1.2 and 1.4 times of ϵtruth.For briefness, we just show the plot of NOE of the reconstructed images as function ofϵ in Figure 9.Figure 9
Plot of NOE as function of multiple ofϵ relative to the truth ϵ.From Figure9, it can be seen that 0.4 times of the truth data tolerance bound achieve the lowest reconstruction error and, therefore, the most accurate reconstruction. The rationale of this phenomenon is that a large data tolerance bound corresponds to a large convex set defined by the data fidelity term; therefore, a smooth image may be selected in the set. A small data tolerance bound corresponds to a small convex set; therefore, the reconstructed image tends to be noisy due to the existence of noise in the projections.


## 3.1.1. Inverse Crime Study

Image reconstruction is a classical inverse problem. Inverse crime studies are a commonly used tool for end-to-end validation of a given approach to inverting an inverse problem. The whole reconstruction chain of the dcTV-CP algorithm includes imaging system modeling, optimization model design, solving algorithm design, and its computer implementation. The inverse crime is considered successful if the reconstructed image is accurate enough such that any residual errors are due to computer float error. For 8-bit grayscale images, there are only 256 possible discrete intensity values to assign to a given pixel. Therefore, a successful inverse crime study in this case will produce a reconstructed image with NOE relative to the truth image being less than or equal to10−4.The phantom used for the inverse crime study is the Shepp–Logan phantom of size256×256. The parallel scanning configuration is used. The projection data are also of size 256×256, which means there are 256 projections evenly distributed in the angle range of 0,πand there are 256 measurements on each projection. The origin of the imaging coordinate system is located at 128,128 of the phantom. The axis range of each projection signal is −128:127. The projections are generated by use of Eq. (1); thus ,the linear system is completely consistent.The system matrix is determined by use of the accurate pixel driven projection method. Since the projection data are ideal, the model parameters are set asϵ=0 and t1=utruthTV. To get fast convergence, the two algorithm parameters are set as λ=1andν=0.1νA.The inverse crime sign, i.e., the convergence condition, is set as(31)NOEun≤10−4,NDEun≤10−4,NTVEun≤10−3.The dcTV-CP algorithm for the study stops at iteration2910. Figures 2(a) and 2(b) show the truth image and the reconstructed image, respectively. It can be seen that the two images are visually identical. Figure 2(c) shows the profiles of the vertical center-line of the images. It shows that the reconstructed profile completely overlaps with the truth profile. At the convergence point, NOE<10−4. Both the qualitative and the quantitative evaluations show that the inverse crime succeeded. This means that the design and implementation of the dcTV-CP algorithm are both correct.Figure 2
dcTV-CP reconstruction images of the Shepp–Logan phantom. (a) The Shepp–Logan phantom image (truth). (b) The reconstructed image. (c) Comparison of the profiles on the vertical center-line of the phantom image and the reconstructed image.
(a)
(b)
(c)The iteration trends of the three metrics for the dcTV-CP algorithm are shown in Figures3(a)–3(c), respectively.Figure 3
Convergence behavior of the dcTV-CP algorithm for reconstructing the Shepp–Logan phantom. (a) The object error NOE. (b) The data error NDE. (c) The relative TV error NTVE, each as a function of the iteration number.
(a)
(b)
(c)From Figure3, one may see that the three metrics decrease monotonically with progressive iterations. Even at the convergence point, the three metrics still maintain the decreasing trend. This shows that the dcTV-CP algorithm is convergent. It is noted that all the iteration curves exhibit a slight oscillatory behavior, which can also be observed in other CP algorithms, for example, DDcTV-CP and TVcDM-CP algorithms. Maybe, this is the iteration law of all the CP algorithms. However, this does not impact the overall downward trend towards convergence.


## 3.1.2. Impact of Algorithm Parameters on Convergence Rate

The two algorithm parameters in the dcTV-CP algorithm areλandν. To investigate their impact on the convergence rate, we varied λandν and compared convergence rates.The Shepp–Logan phantom was used for this study as well. The imaging conditions and parameter selection are the same as those in the inverse crime study, except that the number of projections is 60 and the iteration number is 1500. We used three different combinations of values forλandν :λ=1andν=0.1νA, λ=1andν=νA, and λ=1andν=10νA.Figure4 shows the error images between the reconstructed images and the Shepp–Logan phantom. Figure 5 shows the iteration trend of NOE for the image reconstructions with the three different parameter combinations. From Figure 4, based on error-image observation after a fixed number of iterations, it can be seen that λ=1andν=0.1νA leads to the fastest convergence, whereas λ=1andν=10νA leads to the slowest convergence. This is further evidenced in Figure 5.Figure 4
Error image between the reconstructed images and the truth image. The error image is the absolute value of the difference between the reconstructed image and the truth image. (a)λ=1andν=0.1νA. (b) λ=1andν=νA. (c) λ=1andν=10νA. The display window is [0, 0.01].
(a)
(b)
(c)Figure 5
Iteration curves of NOE for the dcTV-CP algorithms using different combinations of values forλandν.These results show that the two algorithm parameters have significant impact on the convergence rate. For a specific image reconstruction application, one may run the reconstruction algorithm with a different combination of parameter values and select the most appropriate values that provide rapid convergence for that specific application.


## 3.1.3. Impact of TV Bound on Reconstruction Accuracy

There are two model parameters, the TV bound and the data tolerance bound, whose selection may impact the reconstruction accuracy.For the TVcDM-CP algorithm and DDcTV algorithm, it is known how the model parameters impact the reconstruction accuracy. For the DDcTV algorithm, the model parameter is the data tolerance bound. In reconstruction from noisy data, a large data tolerance value can lead to an overly smooth image because an image with a smaller TV value can be obtained in larger convex set defined by the data fidelity term. Conversely, a small data tolerance value can lead to a noisy image because the noise information in the projection data will be propagated into the reconstructed image. For the TVcDM-CP algorithm, the TV bound is the model parameter. In noisy data reconstruction, too large of a TV bound leads to a noisy image because the TV noise-reduction effect is suppressed, whereas too small of a TV bound leads to an overly smooth image because low TV value corresponds to a smooth image.In order to isolate the effect of each model parameter for the dcTV-CP algorithm, their impact on reconstruction is studied separately.As before, the Shepp–Logan phantom is of size256×256; the projection data are of size 256×60. Gaussian noise with magnitude of 45 dB is added to the projections, resulting in ϵtruth=24.59. The model parameter ϵ is set to be ϵtruth, and the model parameter t1 is set to be 0.6, 0.8, 1, 1.2, and 1.4 times of utruth  TV, to determine how it may impact reconstruction accuracy. We set λ=1andν=0.1νA to get fast convergence. The practical convergence conditions are given in the following equation:(32)dNOEun≤10−3,dNDEun≤10−3,dNTVEun≤10−3.The reconstructed images and the error images are shown in Figures6 and 7, respectively, and the plot of NOE as function of t1 is shown in Figure 8.Figure 6
Images reconstructed using the dcTV-CP algorithm and varyingt1. (a–e) t1 of values 0.6, 0.8, 1, 1.2 and 1.4 times utruthTV. The display window is [0, 1].
(a)
(b)
(c)
(d)
(e)Figure 7
Reconstructed error images corresponding to the images in Figure6. The display window is [0, 0.05].
(a)
(b)
(c)
(d)
(e)Figure 8
Plot of NOE as function of multiple oft1 relative to the truth TV.From Figure6, it can be seen that the image becomes smoother when the TV bound is decreased, i.e., more noisy when the TV bound is increased. It appears that a multiple of 1 may achieve the best tradeoff. The error images in Figure 7 are displayed with a narrow intensity window to further demonstrate this phenomenon. It can be clearly seen that a multiple of 1 achieves the best tradeoff. Figure 8 shows the NOE of the five reconstructed images. A multiple of 1 corresponds to the lowest reconstruction error. Both the qualitative and quantitative analysis, therefore, show that a multiple of 1 may achieve the most accurate reconstruction.


## 3.1.4. Impact of Data Tolerance Bound on Reconstruction Accuracy

The impact of the data tolerance bound on reconstruction accuracy is further investigated by repeating the reconstructions in the last section except now;t1 is fixed as the truth TV and ϵ is varied as 0.1, 0.2, 0.3, 0.4, 0.5 0.6, 0.8, 1, 1.2 and 1.4 times of ϵtruth.For briefness, we just show the plot of NOE of the reconstructed images as function ofϵ in Figure 9.Figure 9
Plot of NOE as function of multiple ofϵ relative to the truth ϵ.From Figure9, it can be seen that 0.4 times of the truth data tolerance bound achieve the lowest reconstruction error and, therefore, the most accurate reconstruction. The rationale of this phenomenon is that a large data tolerance bound corresponds to a large convex set defined by the data fidelity term; therefore, a smooth image may be selected in the set. A small data tolerance bound corresponds to a small convex set; therefore, the reconstructed image tends to be noisy due to the existence of noise in the projections.


## 3.2. Comparison of the dcTV-CP, DDcTV-CP, and TVcDM-CP Algorithms

In this section, three comparisons are made between the dcTV-CP, DDcTV-CP, and TVcDM-CP algorithms: comparison of convergence rate, comparison of reconstruction capability from sparse projections, and comparison of reconstruction capability from noisy projections. The DDcTV-CP algorithm is shown in Algorithm3, and the TVcDM-CP algorithm is shown in Algorithm 4.Algorithm 3: Pseudocode for N steps of the DDcTV-CP algorithm instance.
INPUT: g, A, Dϵ; λ, b
(1)
ν
A
=
A
sv
/
D
sv; ν=bνA
(2)
L
=
λ
A
νD
sv; σ=1/L,τ=1/L, θ=1,n=0
(3)
u
0
=
0
;
p
0
=
0
;
q
0
=
0
;
u
¯
0
=
0
(4)
Repeat
(5)
a
=
p
n
+
σ
λ
A
u
¯
n
−
g
(6)
p
n
+
1
=
max
a
2
−
σ
λ
ϵ
,
0
×
a
/
a
2
(7)
c
=
q
n
+
σ
ν
D
u
¯
n
(8)
q
n
+
1
=
c
/
max
1
I
,
c
mag
(9)
u
n
+
1
=
u
n
−
τ
λ
A
T
p
n
+
1
−
τ
ν
D
T
q
n
+
1
(10)
u
¯
n
+
1
=
u
n
+
1
+
θ
u
n
+
1
−
u
n
(11)
n
=
n
+
1
(12)
Until n≥N
OUTPUT: The designed solution uN.Algorithm 4: Pseudocode for N steps of the TVcDM-CP algorithm instance.
INPUT: g, A, D, t1; λ, b
(1)
ν
A
=
A
S
V
/
D
S
V
,
ν
=
b
ν
A
(2)
L
=
A
ν
D
S
V; σ=1/L;τ=1/L;θ=1;n=0
(3)
u
0
=
0
;
u
¯
0
=
0
;
p
0
=
0
;
q
0
=
0
(4)
Repeat
(5)
p
n
+
1
=
p
n
+
σ
A
u
¯
n
−
g
/
1
+
σ
/
λ
(6)
a=qn+σνDu¯n
(7)
s=ProjectOntoℓ1Ballνt1amag/σ
(8)
qn+1=a1I−σs/amag
(9)
u
n
+
1
=
u
n
−
τ
A
T
p
n
+
1
+
ν
D
T
q
n
+
1
(10)
u
¯
n
+
1
=
u
n
+
1
+
θ
u
n
+
1
−
u
n
(11)
n
=
n
+
1
(12)
Until n≥N
OUTPUT: The designed solution uN.


### 3.2.1. Convergence Rate Comparison

We use the FORBILD phantom to perform this comparison. The imaging conditions and parameters are the same as those used in the inverse crime study except that the number of projections is 60. The convergence condition is NOE≤5×10−4.The convergence curves of NOE are shown in Figure10 for the FORBILD reconstruction. It can be seen that the TVcDM-CP algorithm converges the fastest, the DDcTV-CP algorithm converges the slowest, and the newly proposed dcTV-CP algorithm has a slightly slower convergence rate than the TVcDM-CP algorithm.Figure 10
Plot of NOE as function of iteration number for the three algorithms.


### 3.2.2. Comparison of Reconstruction Accuracy from Sparse Projections

TV-based algorithms can accurately reconstruct images from sparse-view projections since this type of optimization model embodies the idea of compressed sensing. In this section, we will compare the sparse reconstruction capability of the three different TV algorithms. Each algorithm has a different optimization meaning dictated by their different constraint forms. So they may perform different sparse reconstruction capabilities.To demonstrate the capability of each algorithm to handle sparse projection data, images are reconstructed from 10, 20, 30, 40, and 50 projections. The imaging conditions and parameters are the same as in the inverse crime study except that we use the different projection numbers mentioned above. The FORBILD phantom is used to perform this study.Figures11 and 12 show the reconstructed FORBILD images and the error images, respectively. Figure 13 shows the plot of NOE as function of the projection number. From Figure 11, it can be seen that all the three TV algorithms achieve accurate sparse reconstruction, as it is hard to tell the difference between the images reconstructed from 20, 30, 40, and 50 projections and the truth image for any TV algorithm (see images from the second to the fifth column). From the images reconstructed from 10 projections, we may see that the DDcTV image is too smooth, whereas the TVcDM and dcTV images have clear streak artifacts. From Figure 13, we may see that the order of reconstruction accuracy from high to low is dcTV, TVcDM, and then DDcTV if ≥20 projections are used. From Figure 12, it can be seen that the dcTV error image reconstructed from 40 projections is subtle, whereas the features in the DDcTV error image and TVcDM error image are so pronounced that the structure of the phantom may be distinguished. This analysis shows that the dcTV-CP algorithm can achieve the highest reconstruction accuracy.Figure 11
FORBILD phantom images reconstructed using the DDcTV-CP (the first row), TVcDM-CP (the second row), and dcTV-CP (the third row) algorithms from 10, 20, 30, 40, and 50 projections. The display window is [0, 1].Figure 12
FORBILD phantom reconstructed error images from 40 projections using the DDcTV (the first column), TVcDM (the second column), and dcTV (the third column) algorithms. The display window is [0, 0.01].
(a)
(b)
(c)Figure 13
Plots of NOE as function of the projection number.


### 3.2.3. Comparison of Reconstruction Accuracy from Noisy Projections

TV algorithms can achieve accurate reconstructions not only from sparse-view projections but also from noisy projections. In fact, TV models have been used in image denoising before its wide application in image reconstruction from 2006. The capability to accurately reconstruct images from noisy projection data is very useful for low-dose CT, in which the dose for each view is reduced, resulting in increased projection noise. The DDcTV and TVcDM models have been deeply investigated in low-dose CT and have shown their capability of suppressing the associated noise. Here, we compare the proposed dcTV-CP algorithm’s capability to reconstruct images from noisy projections with that of the established DDcTV and TVcDM models.The FORBILD phantom is used. The imaging conditions and parameters are the same as those in the last section except that the number of projections is fixed at 60 and the SNR of the noisy projections is varied from 30 dB, 35 dB, 40 dB, 45 dB, to 50 dB.Figures14 and 15 show the reconstructed images and the corresponding error images, respectively, and Figure 16 shows the plot of NOE of the reconstructed images as function of SNR.Figure 14
Reconstructed FORBILD images using the DDcTV-CP (the first row), TVcDM-CP (the second row), and dcTV-CP (the third row) algorithms from projections with SNR of 30 dB, 35 dB, 40 dB, 45 dB, and 50 dB. The display window is [0, 1].Figure 15
FORBILD error images reconstructed from projections with SNR of 30 dB using the DDcTV (the first column), TVcDM (the second column), and dcTV (the third column) algorithms. The display window is [0, 0.01].Figure 16
Plots of NOE as a function of projection SNR.From Figure14, it can be seen that the three TV algorithms all have the capability to accurately reconstruct images from projections with SNR of 35 dB to 50 dB. Subtle differences between the three images reconstructed from projections with SNR of 30 dB can be observed: the DDcTV image is smoother, the TVcDM image is noisier, but the dcTV image appears to balance these two characteristics well (i.e., not as smooth as the DDcTV image and not as noisy as the TVcDM image). This can be clearly seen in Figure 15, where the DDcTV error image exhibits strong features near high-frequency edges in the phantom, the TVcDM error image has a large degree of random differences, but the dcTV error image achieves the balance between these two. This is evidenced quantitatively in Figure 16, where the dcTV-CP algorithm achieves the highest reconstruction accuracy in the context of image reconstruction from noisy projections.


## 3.2.1. Convergence Rate Comparison

We use the FORBILD phantom to perform this comparison. The imaging conditions and parameters are the same as those used in the inverse crime study except that the number of projections is 60. The convergence condition is NOE≤5×10−4.The convergence curves of NOE are shown in Figure10 for the FORBILD reconstruction. It can be seen that the TVcDM-CP algorithm converges the fastest, the DDcTV-CP algorithm converges the slowest, and the newly proposed dcTV-CP algorithm has a slightly slower convergence rate than the TVcDM-CP algorithm.Figure 10
Plot of NOE as function of iteration number for the three algorithms.


## 3.2.2. Comparison of Reconstruction Accuracy from Sparse Projections

TV-based algorithms can accurately reconstruct images from sparse-view projections since this type of optimization model embodies the idea of compressed sensing. In this section, we will compare the sparse reconstruction capability of the three different TV algorithms. Each algorithm has a different optimization meaning dictated by their different constraint forms. So they may perform different sparse reconstruction capabilities.To demonstrate the capability of each algorithm to handle sparse projection data, images are reconstructed from 10, 20, 30, 40, and 50 projections. The imaging conditions and parameters are the same as in the inverse crime study except that we use the different projection numbers mentioned above. The FORBILD phantom is used to perform this study.Figures11 and 12 show the reconstructed FORBILD images and the error images, respectively. Figure 13 shows the plot of NOE as function of the projection number. From Figure 11, it can be seen that all the three TV algorithms achieve accurate sparse reconstruction, as it is hard to tell the difference between the images reconstructed from 20, 30, 40, and 50 projections and the truth image for any TV algorithm (see images from the second to the fifth column). From the images reconstructed from 10 projections, we may see that the DDcTV image is too smooth, whereas the TVcDM and dcTV images have clear streak artifacts. From Figure 13, we may see that the order of reconstruction accuracy from high to low is dcTV, TVcDM, and then DDcTV if ≥20 projections are used. From Figure 12, it can be seen that the dcTV error image reconstructed from 40 projections is subtle, whereas the features in the DDcTV error image and TVcDM error image are so pronounced that the structure of the phantom may be distinguished. This analysis shows that the dcTV-CP algorithm can achieve the highest reconstruction accuracy.Figure 11
FORBILD phantom images reconstructed using the DDcTV-CP (the first row), TVcDM-CP (the second row), and dcTV-CP (the third row) algorithms from 10, 20, 30, 40, and 50 projections. The display window is [0, 1].Figure 12
FORBILD phantom reconstructed error images from 40 projections using the DDcTV (the first column), TVcDM (the second column), and dcTV (the third column) algorithms. The display window is [0, 0.01].
(a)
(b)
(c)Figure 13
Plots of NOE as function of the projection number.


## 3.2.3. Comparison of Reconstruction Accuracy from Noisy Projections

TV algorithms can achieve accurate reconstructions not only from sparse-view projections but also from noisy projections. In fact, TV models have been used in image denoising before its wide application in image reconstruction from 2006. The capability to accurately reconstruct images from noisy projection data is very useful for low-dose CT, in which the dose for each view is reduced, resulting in increased projection noise. The DDcTV and TVcDM models have been deeply investigated in low-dose CT and have shown their capability of suppressing the associated noise. Here, we compare the proposed dcTV-CP algorithm’s capability to reconstruct images from noisy projections with that of the established DDcTV and TVcDM models.The FORBILD phantom is used. The imaging conditions and parameters are the same as those in the last section except that the number of projections is fixed at 60 and the SNR of the noisy projections is varied from 30 dB, 35 dB, 40 dB, 45 dB, to 50 dB.Figures14 and 15 show the reconstructed images and the corresponding error images, respectively, and Figure 16 shows the plot of NOE of the reconstructed images as function of SNR.Figure 14
Reconstructed FORBILD images using the DDcTV-CP (the first row), TVcDM-CP (the second row), and dcTV-CP (the third row) algorithms from projections with SNR of 30 dB, 35 dB, 40 dB, 45 dB, and 50 dB. The display window is [0, 1].Figure 15
FORBILD error images reconstructed from projections with SNR of 30 dB using the DDcTV (the first column), TVcDM (the second column), and dcTV (the third column) algorithms. The display window is [0, 0.01].Figure 16
Plots of NOE as a function of projection SNR.From Figure14, it can be seen that the three TV algorithms all have the capability to accurately reconstruct images from projections with SNR of 35 dB to 50 dB. Subtle differences between the three images reconstructed from projections with SNR of 30 dB can be observed: the DDcTV image is smoother, the TVcDM image is noisier, but the dcTV image appears to balance these two characteristics well (i.e., not as smooth as the DDcTV image and not as noisy as the TVcDM image). This can be clearly seen in Figure 15, where the DDcTV error image exhibits strong features near high-frequency edges in the phantom, the TVcDM error image has a large degree of random differences, but the dcTV error image achieves the balance between these two. This is evidenced quantitatively in Figure 16, where the dcTV-CP algorithm achieves the highest reconstruction accuracy in the context of image reconstruction from noisy projections.


## 4. Discussion and Conclusions

The TV-based reconstruction algorithms have the capability to accurately reconstruct images from sparse-view projections and/or noisy projections in CT. This is very useful for low-dose CT as these cases correspond to two low-dose data acquisition patterns. In MRI and EPRI, the TV algorithm may achieve accurate image reconstruction from more rapid data acquisition sequences.Though the three constrained TV algorithms, DDcTV, TVcDM, and dcTV algorithm, can all achieve accurate reconstructions from sparse data and/or noisy data, their reconstruction characteristics and utility are different.In this work, we propose the novel dcTV-CP algorithm and characterize its performance by analyzing its convergence behavior and exploring how to select the model parameters and algorithm parameters. These studies show that dcTV-CP algorithm is convergent, the two introduced algorithm parameters impact the convergence rate significantly and should be chosen to suit a given application, and both the TV bound and data tolerance bound impact the reconstruction accuracy and therefore should be selected to achieve appropriate tradeoffs in image quality. One may always use the optimal TV bound for the TVcDM model and the optimal data tolerance bound for the DDcTV model as the two optimal model parameters for the dcTV model to achieve the highest or near-highest reconstruction accuracy. According to the optimization meaning explanation of the convex-ellipses schematic diagram of the three TV models, we know that this approach for selecting the two model parameters may guarantee the existence of the intersection set of the two convex sets.We have systematically compared the three TV models using the same solving algorithm, i.e., the CP algorithm. Comparison of the convergence rate demonstrated that the dcTV-CP and TVcDM-CP algorithms have faster convergence rates than the DDcTV-CP algorithm. In the context of sparse reconstruction, studies show that the dcTV-CP algorithm achieves the highest accuracy. In the context of nonideal data reconstruction, the studies show that the dcTV-CP algorithm achieves the highest accuracy.The novel dcTV model proposed in this work tries to find a balance between the issues of the DDcTV and TVcDM models. The dcTV model uses not only the TV prior but also the data fidelity prior quantitatively. Therefore, the TV bound controls the TV value quantitatively and the data tolerance bound controls the data divergence value quantitatively. If the TV bound is the truth TV and the data tolerance bound can optimally control the noise and system inconsistency, then the dcTV model can achieve more accurate reconstructions relative to the DDcTV and TVcDM models. The theoretical analysis and the evaluation experiments consistently show that the dcTV-CP algorithm of appropriately selected model parameters and algorithm parameters may be superior to the existing single-constrained TV algorithm, the DDcTV-CP and TVcDM-CP algorithms.The insight and knowledge gained in this work may also be extended to other optimization models in image reconstruction. Though solving algorithms for multiconstraint models are practically difficult to design, we have shown that the CP algorithm framework is completely capable. Also, we suggest the use of a multiconstraint optimization model in image denoising, image restoration, and other image processing task for its potential to further improve the image quality.



---

*Source: 5382106-2020-01-20.xml*