# A Novel Model and ADMM Algorithm for MR Image Reconstruction

Authors: Bo Zhou; Yu-Fei Yang; Wei-Si Xie
Journal: Mathematical Problems in Engineering
 (2018)

---

## Abstract

Motivated by the ideas from the LOT model and its deformations, we propose a coupling model for the MR image reconstruction and apply the split Bregman iterative method on the proposed model by utilizing the augmented Lagrangian technique. The related minimization problem is then divided into four subproblems by means of the alternating minimization method. And on this basis, by combining the Barzilai-Borwein step size selection scheme, generalized shrinkage formulas, and the shrink operator, we propose an ADMM type algorithm to solve the proposed model. Several numerical examples are implemented; the experimental results demonstrate the feasibility and effectiveness of the proposed model and algorithm.


---

## Body

## 1. Introduction

Magnetic resonance (MR) images are obtained by placing an object in a strong magnetic field and then turning on and off a radio frequency electromagnetic field. Different body parts produce different signals which are detected by a receiver. So the obtained images via this technique are noninvasive and nonionizing and then we can get excellent visualization of both anatomical structure and physiological function. However, this imaging modality is relatively slow because the data in MR imaging are acquired to sequentially samplek-space of the spatial Fourier transformation of the object [1]. Therefore, many MR techniques [1–8] have been developed to reduce the number of data needed for the accurate reconstruction. Due to the inherent instability of the MR image reconstruction, it is necessary to apply regularization in order to compute meaningful reconstructions.In general, the existing regularization schemes can be fell into two categories: the adaptively learned dictionary and the predefined transformation [3, 4, 9–12]. The first type of regularization method exploits the nonlocal similarity and sparse representation. Most of the existing models adopt a two-step iterative scheme, where the sparse representation approximations are found with the dictionary fixed and then the dictionary is subsequently optimized based on the current sparse representation coefficients [4, 9, 12]. For the second category, it is usually related to total variation and wavelet transform. For example, Ma et al. in [11] employed the total variation (TV) penalty [13, 14] and the wavelet transform for the MR imaging problem. This is due to the fact that the total variation term is a powerful technique when the sought solution is required to have sharp edges. However, this term often leads to the stair casing effect and also results in patchy, cartoon-like images which appear unnatural. In order to overcome the drawbacks of the TV penalty, some high-order models such as the classic total generalized variation (TGV) have been proposed in [3, 10, 11, 15–17]. The TGV in nature is equivalent to TV in terms of edge preservation and noise removal, which can also be true of imaging situations where the assumption of what the TV penalty is based on is not effective. It is more precise in describing intensity variation in smooth regions and thus reduces oil painting artifacts while still being able to preserve sharp edges as the TV-norm did. Based on these facts, Guo et al. in [3] recently proposed an outstanding method that combines shearlet transformation and the TGV term, which is able to recover both the texture and the smoothly varied intensities while the other methods such as the shearlet and the TGV only models either return a cartoon image or lose the textures.In order to eliminate the undesired patchy artifacts of the above-mentioned TV-based methods, another high-order derivative model, called two-step reconstruction model, has attracted great attention in the field of the image restoration [18–23]. This kind of model is aimed at eliminating the patchy artifacts by introducing normal direction of isophote lines and then fitting the obtained normal vector field to get the restoration image. In detail, this method involves two procedures in its implementation: ① estimation of normal vector field from initial or intermediate images; ② image reconstruction with data constraints in the normal field, which is derived from the estimated tangent field. The introduction of the normal vector field is corresponding to retain continuous and smooth isophote lines. The second fitting step indicates the more consistency on the isophote lines for continuity and smoothness. Thus, visually pleasant images with smooth regions and continuous boundaries were obtained, where the staircase and patchy artifacts caused by the oversmoothing along the normal directions in the above TV-based approaches were efficiently mitigated. Motivated by the advantages of the above two-step models, this model can be extended to the MR image reconstruction problem:(i)
The first step, we can get the regularized normal vector field by solving(1)minn=1∫Ω∇ndx+α2n-∇f∇f2,wheref is the partially scanned k space data, ∇n is the total variation of the unit normal vector n≔(n1,n2) for the level curves of image u, n=n12+n222, and α is a regularization parameter.(ii)
The second step, we can obtain the reconstruction image by solving(2)minu∫Ω∇u-n∗·∇udx+λΦTu1+β2Au-f22,wheren∗ is the solution of (1), A is involved in the sampling matrix, β is a regularization parameter, and Φ is the wavelet transformation.However, the above model has a few drawbacks such as the numerical difficulty for findingn and not efficiently utilizing the information of the operator A. In this paper, we couple this two steps into one model by introducing a new variable θ to replacen with n:=(cos⁡θ,sin⁡θ) as follows:(3)minθ,u∫Ω∇θdx+λ1∇u1+λΦTu1+α2θ-θ02+β2Au-f22+γ2∇u-cos⁡θ,sin⁡θ·∇u2,where ∇u1=∫Ω‍∇udx. Here, ∫Ω‍∇θdx replacing ∫Ω‍∇cos⁡θ,sin⁡θdx=∫Ω‍∇ndx mainly uses the fact that moving the cos and sin transformations does not change when using the total variation regularization to describe the image structures. In addition, we here also use the fact that (4)cos⁡θ-cos⁡θ02=1-cos⁡θ-θ0~θ-θ022,i.e., 1-cos⁡(θ-θ0) is the equivalent infinitesimal quantity of (θ-θ0)2/2, where (cos⁡θ0,sin⁡θ0)=∇f/∇f. The proposed model (3) is a nonsmoothing optimization problem; here we employ the alternating direction method of multipliers (ADMM) [24] to solve it and then use the Barzilai-Borwein step size selection scheme [25] to improve the efficiency of the numerical computation. The experimental results show that the proposed algorithm is an efficient and fast algorithm and can get better MR reconstructed image.The organization of this paper is as follows. In the next section, based on the augmented Lagrangian scheme, we first apply the split Bregman iterative method to solve the proposed coupling model. The related minimization problem is then divided into four subproblems by means of the alternating minimization method. Furthermore, by combining the Barzilai-Borwein step size selection scheme, generalized shrinkage formulas, and the shrink operator, we propose an ADMM type algorithm for the MR image reconstruction. In order to verify the feasibility and effectiveness of the proposed model and algorithm, we implement some numerical experiments by the proposed algorithm in Section3. Some concluding remarks are given in the last section.Some symbols are introduced in this paper.u∈RN is the objective image consisting of N pixels. Φ⊤=(ϕ1,ϕ2,…,ϕn)⊤∈CN×N is usually a proper orthogonal matrix (e.g., Haar wavelet) sparsifying the underlying image u. f∈CM with M≤N is the partially scanned data. R∈RM×N is the binary matrix representing the sampling pattern. Considering that the Fourier transform F is in fact N×N unitary matrix, A=RF is the partial Fourier transform operator.


## 2. Numerical Method

The nonsmoothing term in the model (3) leads to some numerical difficulties. One of efficient schemes is to introduce some auxiliary variables to transform the original problem into several easily solvable subproblems. Here we employ the split Bregman iterative scheme to solve the above model. The first step is of using an auxiliary variable η to transform ∇θ in (3) with η=(η1,η2)⊤=(∇xθ,∇yθ)⊤. Then using another auxiliary variable z transforms Φ⊤u out of the nondifferentiable norm Φ⊤u1; (3) thus can be turned into the following problem:(5)minu,θ,η,zη1+λ1∇u1+λz1+α2θ-θ022+β2Au-f22+γ2∇u-cos⁡θ,sin⁡θ·∇u22.s.t.η=∇θ,z=Φ⊤u.The following augmented Lagrangian functional can be used to cope with the constrains in (5):(6)Lu,θ,η,z;b,c=η1+ρ12η1-∇xθ-b122+ρ12η2-∇yθ-b222+λz1+ρ22z-Φ⊤u-c22+α2θ-θ022+λ1∇u1+β2Au-f22+γ2∇u-cos⁡θ,sin⁡θ·∇u22,where b=(b1,b2)⊤∈CN×2N and c∈CN are the Lagrangian multipliers and ρ1 and ρ2 are positive penalty parameters. The split Bregman iteration applied to minimize L in (6) then can be written as the following form:(7)uk+1,θk+1,ηk+1,zk+1=argminu,θ,η,,zLu,θ,η,z;bk,ck,b1k+1=b1k+κ1∇xθk+1-η1k+1,b2k+1=b2k+κ1∇yθk+1-η2k+1,ck+1=ck+κ2Φ⊤uk+1-zk+1.By applying the alternating minimization method [24], we can divide the minimization problem in (7) into the following four subproblems:(8)uk+1=argminu2λ1∇u1+λρ2zk-Φ⊤u-ck22+βAu-f22+γ∇u-cos⁡θk,sin⁡θk·∇u22,(9)θk+1=argminθρ1ηk-∇θ-bk2+αθ-θ022+γ∇uk+1-cos⁡θ,sin⁡θ·∇uk+122,(10)η1k+1,η2k+1=argminη1,η2⁡η1+ρ12η-∇θk+1-bk2,(11)zk+1=argminz⁡z1+ρ22z-Φ⊤uk+1-ck22.Since the updates of the multipliersb and c in (7) are not too complicated, in the forthcoming subsections, only how to efficiently solve each of subproblems (8) to (11) will be discussed.


### 2.1.u-Subproblem

The first problem is the u-subproblem in (8), which can be solved by the split Bregman iterative scheme. By using variable ω=(ω1,ω2)⊤ to replace ∇u in (8) subject to ω=∇u and applying the augmented Lagrangian technique, the split Bregman iteration for (8) can be described as follows:(12)uk+1,ωk+1=argminu,ω2λ1ω1+λρ2zk-Φ⊤u-ck22+βAu-f22+γω-cos⁡θk,sin⁡θk·ω22+γ1ω-∇u-dk2,(13)dk+1=dk+κ3∇uk+1-ωk+1,where the variable d=(d1,d2)⊤ is the Lagrangian multiplier and γ1 is positive penalty parameter.By applying the alternating minimization scheme to the minimization problem (12) and combining a simplified calculation of ωk+1, problems (12)-(13) are converted to the following form: (14)uk+1=argminuλρ2zk-Φ⊤u-ck22+βAu-f22+γ1ωk-∇u-dk2,(15)ω1k+1,ω2k+1=argminω1,ω2gk+1ω1+γ1ω1-∇xuk+1-d1k22+γ1ω2-∇yuk+1-d2k22,(16)d1k+1=d1k+κ3∇xuk+1-ω1k+1,(17)d2k+1=d2k+κ3∇yuk+1-ω2k+1,where gk+1=2λ1+γ1-(cos⁡θk,sin⁡θk)·∇uk+1/∇uk+12·∇uk+1.As pointed out in [26], if A⊤A cannot be diagonalized, then utilizing the fast Fourier transform to directly solve (14) would lead to losing efficiency of fast algorithm. To overcome this drawback, as done in [27], we consider quadratic approximation at point uk of H(u)=1/2Au-f22(18)Qδu,uk=Huk+∇Huk,u-uk+δ2u-uk22,which is a linearization of H(u) at point uk plus a proximity term penalized by parameter δ>0. In addition, to speed up convergence, we also adopt the Barzilai-Borwein step size selection scheme [25] to update the parameter δ, as done in [26, 27].By replacingH(u) in (14) with its quadratic approximation Qδ(u,uk) and considering ∇H(u)=A⊤(Au-f), the minimization subproblem (14) can be converted to the form(19)uk+1=argminuλρ2zk-Φ⊤u-c22+βQδu,uk+γ1ωk-∇u-dk2=argminuλρ2zk-Φ⊤u-c22+βδku-uk+δk-1A⊤Auk-f22+γ1ωk-∇u-dk2.By means of the Barzilai-Borwein (BB) step size selection scheme proposed in [25], the iterative step size δk+1 can be updated by(20)δk+1=argminδ⁡∇Huk+1-∇Huk-δuk+1-uk22,which can be given by(21)δk+1=∇Huk+1-∇Huk,uk+1-ukuk+1-uk22=Auk+1-uk22uk+1-uk22.Considering∇x⊤∇x+∇y⊤∇y=-Δ and ΦΦ⊤=I, the u-subproblem (19) can be solved by its necessary optimality condition, given in the following equation:(22)L1kuk+1=rhs1k,where L1k=λρ2I+βδkI-γ1Δ,(23)rhs1k=γ1∇x⊤w1k-d1k+γ1∇y⊤w2k-d2k+λρ2Φzk-ck+βδkuk-βA⊤Auk-f.Under periodic boundary condition, since the Laplace operator is block circular and can be diagonalized by Fourier transformF with the property F⊤=F-1, the solution of the equation (22) can be computed by(24)F⊤L^1kFuk+1=rhs1k,where L^1k=λρ2I+βδkI-γ1FΔF⊤. Solving (24) is very simple because L^1k is diagonal.The solution of theω-subproblem (15) can be obtained by using a generalized shrinkage formula, as done in [14], and is given by(25)ω1k+1=max⁡sk+1gk+1-2gk+1γ1,0·∇xuk+1+d1ksk+1,ω2k+1=max⁡sk+1gk+1-2gk+1γ1,0·∇yuk+1+d2ksk+1,where(26)sk+1=∇xuk+1+d1k2+∇yuk+1+d2k2.


### 2.2.θ-Subproblem

The next subproblem isθ-subproblem in (9), to accelerate the efficiency of the solution of θ, we use θk to replace θ0 in (9), and then the problem would be changed as follows: (27)θk+1=argminθρ1η1k-∇xθ-b1k22+ρ1η2k-∇yθ-b2k22+αθ-θk22+γ∇uk+1-cos⁡θ,sin⁡θ·∇uk+122.Considering the optimality condition aboutθ, problem (27) can be solved by the following equation: (28)αI-ρ1Δθ=γ-sin⁡θ,cos⁡θ·∇uk+1∇uk+1-cos⁡θ,sin⁡θ·∇uk+1+αθk+ρ1∇x⊤η1k-b1k+ρ1∇y⊤η2k-b2k.To simplify the solution ofθ, we use sin⁡θk,cos⁡θk to replace sin⁡θ,cos⁡θ in (28), respectively, so as to reduce the complexity of calculating θ. Combining the diagonalization technique of the Fourier transform, (28) is converted to the form(29)F⊤L^2kFθk+1=rhs2k,where L^2k=αI-ρ1FΔF⊤, (30)rhs2k=γ-sin⁡θk,cos⁡θk·∇uk+1∇uk+1-cos⁡θk,sin⁡θk·∇uk+1+αθk+ρ1∇x⊤η1k-b1k+ρ1∇y⊤η2k-b2k.Solving (29) is also very simple.


### 2.3.(η,z) Shrink-Subproblem

The remaining subproblems would be discussed in this subsection, including (10) about η-subproblem and (11) about z-subproblem. By a generalized shrinkage formula, as done in [14], we can obtain the solution of the η-subproblem (10) as follows:(31)η1k+1=max⁡tk+1-1ρ1,0·∇xθk+1+bxktk+1,η2k+1=max⁡tk+1-1ρ1,0·∇yθk+1+byktk+1,where (32)tk+1=∇xθk+1+b1k2+∇yθk+1+b2k2.The z-subproblem (11) can be solved easily by the shrink operator, defined by shrink(a,b)=a/amax(a-1/b,0), in the following:(33)zk+1=shrinkΦ⊤uk+1+ck,1ρ2.


### 2.4. The Proposed Algorithm

Based on the above discussion, we summarize an ADMM type algorithm for the MR image reconstruction coping model (3) in the following Algorithm 1.Algorithm 1 (ADMM-TVL1C).
Initialize:u0:=F-1f, δ0:=1, b10:=b20:=c0:=d10:=d20:=0;pij0:=z0:=0,κ1=κ2=κ3=1;
Whileuk-uk-1>tol(34)L^1k=λρ2I+βδkI-γ1FΔF⊤,rhs1k=γ1∇x⊤w1k-d1k+γ1∇y⊤w2k-d2k+λρ2Φzk-ck+βδkuk-βA⊤Auk-f,uk+1=F-1L^1k-1Frhs1k,gk+1=2λ1+γ1-cos⁡θk,sin⁡θk·∇uk+1∇uk+12·∇uk+1,sk+1=∇xuk+1+d1k2+∇yuk+1+d2k2,ω1k+1=max⁡sk+1gk+1-2gk+1γ1,0·∇xuk+1+d1ksk+1,ω2k+1=max⁡sk+1gk+1-2gk+1γ1,0·∇yuk+1+d2ksk+1,d1k+1=d1k+κ3∇xuk+1-ω1k+1,d2k+1=d2k+κ3∇yuk+1-ω2k+1,L^2k=αI-ρ1FΔF⊤,rhs2k=γ-sin⁡θk,cos⁡θk·∇uk+1∇uk+1-cos⁡θk,sin⁡θk·∇uk+1+αθk+ρ1∇x⊤η1k-b1k+ρ1∇y⊤η2k-b2k,θk+1=F-1L^2k-1Frhs2k,tk+1=∇xθk+1+b1k2+∇yθk+1+b2k2,η1k+1=max⁡tk+1-1ρ1,0·∇xθk+1+bxktk+1,η2k+1=max⁡tk+1-1ρ1,0·∇yθk+1+byktk+1,zk+1=shrinkΦ⊤uk+1+ck,gk+1ρ2,b1k+1=b1k+κ1∇xθk+1-η1k+1,b2k+1=b2k+κ1∇yθk+1-η2k+1,ck+1=ck+κ2Φ⊤uk+1-zk+1,δk+1=Auk+1-uk22uk+1-uk22,
end


## 2.1.u-Subproblem

The first problem is the u-subproblem in (8), which can be solved by the split Bregman iterative scheme. By using variable ω=(ω1,ω2)⊤ to replace ∇u in (8) subject to ω=∇u and applying the augmented Lagrangian technique, the split Bregman iteration for (8) can be described as follows:(12)uk+1,ωk+1=argminu,ω2λ1ω1+λρ2zk-Φ⊤u-ck22+βAu-f22+γω-cos⁡θk,sin⁡θk·ω22+γ1ω-∇u-dk2,(13)dk+1=dk+κ3∇uk+1-ωk+1,where the variable d=(d1,d2)⊤ is the Lagrangian multiplier and γ1 is positive penalty parameter.By applying the alternating minimization scheme to the minimization problem (12) and combining a simplified calculation of ωk+1, problems (12)-(13) are converted to the following form: (14)uk+1=argminuλρ2zk-Φ⊤u-ck22+βAu-f22+γ1ωk-∇u-dk2,(15)ω1k+1,ω2k+1=argminω1,ω2gk+1ω1+γ1ω1-∇xuk+1-d1k22+γ1ω2-∇yuk+1-d2k22,(16)d1k+1=d1k+κ3∇xuk+1-ω1k+1,(17)d2k+1=d2k+κ3∇yuk+1-ω2k+1,where gk+1=2λ1+γ1-(cos⁡θk,sin⁡θk)·∇uk+1/∇uk+12·∇uk+1.As pointed out in [26], if A⊤A cannot be diagonalized, then utilizing the fast Fourier transform to directly solve (14) would lead to losing efficiency of fast algorithm. To overcome this drawback, as done in [27], we consider quadratic approximation at point uk of H(u)=1/2Au-f22(18)Qδu,uk=Huk+∇Huk,u-uk+δ2u-uk22,which is a linearization of H(u) at point uk plus a proximity term penalized by parameter δ>0. In addition, to speed up convergence, we also adopt the Barzilai-Borwein step size selection scheme [25] to update the parameter δ, as done in [26, 27].By replacingH(u) in (14) with its quadratic approximation Qδ(u,uk) and considering ∇H(u)=A⊤(Au-f), the minimization subproblem (14) can be converted to the form(19)uk+1=argminuλρ2zk-Φ⊤u-c22+βQδu,uk+γ1ωk-∇u-dk2=argminuλρ2zk-Φ⊤u-c22+βδku-uk+δk-1A⊤Auk-f22+γ1ωk-∇u-dk2.By means of the Barzilai-Borwein (BB) step size selection scheme proposed in [25], the iterative step size δk+1 can be updated by(20)δk+1=argminδ⁡∇Huk+1-∇Huk-δuk+1-uk22,which can be given by(21)δk+1=∇Huk+1-∇Huk,uk+1-ukuk+1-uk22=Auk+1-uk22uk+1-uk22.Considering∇x⊤∇x+∇y⊤∇y=-Δ and ΦΦ⊤=I, the u-subproblem (19) can be solved by its necessary optimality condition, given in the following equation:(22)L1kuk+1=rhs1k,where L1k=λρ2I+βδkI-γ1Δ,(23)rhs1k=γ1∇x⊤w1k-d1k+γ1∇y⊤w2k-d2k+λρ2Φzk-ck+βδkuk-βA⊤Auk-f.Under periodic boundary condition, since the Laplace operator is block circular and can be diagonalized by Fourier transformF with the property F⊤=F-1, the solution of the equation (22) can be computed by(24)F⊤L^1kFuk+1=rhs1k,where L^1k=λρ2I+βδkI-γ1FΔF⊤. Solving (24) is very simple because L^1k is diagonal.The solution of theω-subproblem (15) can be obtained by using a generalized shrinkage formula, as done in [14], and is given by(25)ω1k+1=max⁡sk+1gk+1-2gk+1γ1,0·∇xuk+1+d1ksk+1,ω2k+1=max⁡sk+1gk+1-2gk+1γ1,0·∇yuk+1+d2ksk+1,where(26)sk+1=∇xuk+1+d1k2+∇yuk+1+d2k2.


## 2.2.θ-Subproblem

The next subproblem isθ-subproblem in (9), to accelerate the efficiency of the solution of θ, we use θk to replace θ0 in (9), and then the problem would be changed as follows: (27)θk+1=argminθρ1η1k-∇xθ-b1k22+ρ1η2k-∇yθ-b2k22+αθ-θk22+γ∇uk+1-cos⁡θ,sin⁡θ·∇uk+122.Considering the optimality condition aboutθ, problem (27) can be solved by the following equation: (28)αI-ρ1Δθ=γ-sin⁡θ,cos⁡θ·∇uk+1∇uk+1-cos⁡θ,sin⁡θ·∇uk+1+αθk+ρ1∇x⊤η1k-b1k+ρ1∇y⊤η2k-b2k.To simplify the solution ofθ, we use sin⁡θk,cos⁡θk to replace sin⁡θ,cos⁡θ in (28), respectively, so as to reduce the complexity of calculating θ. Combining the diagonalization technique of the Fourier transform, (28) is converted to the form(29)F⊤L^2kFθk+1=rhs2k,where L^2k=αI-ρ1FΔF⊤, (30)rhs2k=γ-sin⁡θk,cos⁡θk·∇uk+1∇uk+1-cos⁡θk,sin⁡θk·∇uk+1+αθk+ρ1∇x⊤η1k-b1k+ρ1∇y⊤η2k-b2k.Solving (29) is also very simple.


## 2.3.(η,z) Shrink-Subproblem

The remaining subproblems would be discussed in this subsection, including (10) about η-subproblem and (11) about z-subproblem. By a generalized shrinkage formula, as done in [14], we can obtain the solution of the η-subproblem (10) as follows:(31)η1k+1=max⁡tk+1-1ρ1,0·∇xθk+1+bxktk+1,η2k+1=max⁡tk+1-1ρ1,0·∇yθk+1+byktk+1,where (32)tk+1=∇xθk+1+b1k2+∇yθk+1+b2k2.The z-subproblem (11) can be solved easily by the shrink operator, defined by shrink(a,b)=a/amax(a-1/b,0), in the following:(33)zk+1=shrinkΦ⊤uk+1+ck,1ρ2.


## 2.4. The Proposed Algorithm

Based on the above discussion, we summarize an ADMM type algorithm for the MR image reconstruction coping model (3) in the following Algorithm 1.Algorithm 1 (ADMM-TVL1C).
Initialize:u0:=F-1f, δ0:=1, b10:=b20:=c0:=d10:=d20:=0;pij0:=z0:=0,κ1=κ2=κ3=1;
Whileuk-uk-1>tol(34)L^1k=λρ2I+βδkI-γ1FΔF⊤,rhs1k=γ1∇x⊤w1k-d1k+γ1∇y⊤w2k-d2k+λρ2Φzk-ck+βδkuk-βA⊤Auk-f,uk+1=F-1L^1k-1Frhs1k,gk+1=2λ1+γ1-cos⁡θk,sin⁡θk·∇uk+1∇uk+12·∇uk+1,sk+1=∇xuk+1+d1k2+∇yuk+1+d2k2,ω1k+1=max⁡sk+1gk+1-2gk+1γ1,0·∇xuk+1+d1ksk+1,ω2k+1=max⁡sk+1gk+1-2gk+1γ1,0·∇yuk+1+d2ksk+1,d1k+1=d1k+κ3∇xuk+1-ω1k+1,d2k+1=d2k+κ3∇yuk+1-ω2k+1,L^2k=αI-ρ1FΔF⊤,rhs2k=γ-sin⁡θk,cos⁡θk·∇uk+1∇uk+1-cos⁡θk,sin⁡θk·∇uk+1+αθk+ρ1∇x⊤η1k-b1k+ρ1∇y⊤η2k-b2k,θk+1=F-1L^2k-1Frhs2k,tk+1=∇xθk+1+b1k2+∇yθk+1+b2k2,η1k+1=max⁡tk+1-1ρ1,0·∇xθk+1+bxktk+1,η2k+1=max⁡tk+1-1ρ1,0·∇yθk+1+byktk+1,zk+1=shrinkΦ⊤uk+1+ck,gk+1ρ2,b1k+1=b1k+κ1∇xθk+1-η1k+1,b2k+1=b2k+κ1∇yθk+1-η2k+1,ck+1=ck+κ2Φ⊤uk+1-zk+1,δk+1=Auk+1-uk22uk+1-uk22,
end


## 3. Numerical Implementation

In this section, we will give some numerical experiments, composed of two kinds of numerical comparisons. One comparison is between the proposed algorithm ADMM-TVL1C and two existing algorithms by using same images including abdominal and brain images; numerical results demonstrated that the proposed algorithm can reconstruct MR images higher effectively and faster than the other two algorithms. Another comparison is between different sampling rates by using a brain image; results illustrated that the proposed algorithm can get higher quality image with low staircase effect.In the following experiments, the ratio of signal to noise (SNR) is defined by(35)SNR=10log10uture22u-uture22.The partially scanned data were generated by nonuniformly randomly selecting some of the k-space samples from the entire k-space data of the tested MR images. We chose ρ1=10-3,ρ2=10-4,γ1=10-3 for all the numerical experiments. The other two parameters α,β related to the sparsity of the underlying of image and the noise level were selected by α=10-3,β=0.9, and the parameters λ, λ1 were chosen by λ1=10-3 and λ=0.5∗10-3.


### 3.1. Comparisons with Other Reconstruction Algorithms

In this subsection, by using the abdominal MR image for the test, we compare the reconstruction effectiveness of three algorithms: the split Bregman algorithm for solvingTVL1 model in [28] (called SBA-TVL1), the alternating direction method of multipliers for solving TV2L1 model proposed in [27] (called ADMM-TV2L1), and the proposed algorithm in this paper: ADMM-TVL1C. In the first experiment, the original image is a human abdomen MR image of size 256×256, we simulate incomplete Fourier data with various sampling rates. In order to demonstrate the performance of the proposed algorithm, the image is corrupted with the Gaussian white noise at level σ=10 in this experiment.In Figure1, we show the reconstructed results of ADMM-TVL1C, SBA-TVL1, and ADMM-TV2L1 algorithms from 40%k-space data for the abdomen MR image. For better visual comparison, we zoom in on the central area to observe the texture features of these reconstructed images. Moreover, the SNRs of the abdomen MR image under different k-space sampling rates from 32% to 50% are given in Table 1, which show that the proposed algorithm is superior to the other two reconstruction algorithms.Table 1
SNRs via different algorithms and sample rates for abdomen image.
Sample rate 32% 34% 36% 38% 40% SBA-TVL1 25.87 25.90 25.91 26.03 26.01 ADMM-TV2L1 37.23 37.69 38.08 38.35 38.74 ADMM-TVL1C 37.92 39.04 39.41 40.08 40.37  Sample rate 42% 44% 46% 48% 50%  SBA-TVL1 26.06 26.14 26.16 26.13 26.28 ADMM-TV2L1 39.23 39.81 40.29 40.46 40.90 ADMM-TVL1C 41.49 41.89 42.37 42.82 43.00Figure 1
Abdomen MR image reconstruction via different algorithms from 40% sample rate.In the second experiment, we use a human brain MR image of size 198×198 to test the three algorithms under different sampling rates, as done in normal MR image reconstruction [29]. The SNRs of the human brain MR image under different k-space sampling rates from 32% to 50% are given in Table 2; the reconstructed images from 40% k-space data and the images zoomed in on their central areas are shown in Figure 2. The results demonstrate that the quality of the reconstructed images is improved as the sampling rate increases and that ADMM-TVL1C outperforms SBA-TVL1 and ADMM-TV2L1 in the ability to recover the textures and alleviate the staircase effect.Table 2
SNRs via different algorithms and sample rates for brain image.
Sample rate 32% 34% 36% 38% 40% SBA-TVL1 28.07 28.11 28.18 28.22 28.30 ADMM-TV2L1 40.41 40.91 41.38 41.44 42.06 ADMM-TVL1C 41.02 42.52 43.70 44.82 46.36  Sample rate 42% 44% 46% 48% 50%  SBA-TVL1 28.42 28.44 28.35 28.57 28.48 ADMM-TV2L1 42.61 42.87 43.11 43.61 43.95 ADMM-TVL1C 46.72 47.51 47.73 48.15 48.55Figure 2
Brain MR image reconstruction via different algorithms from 40% sample rate.In the third experiment, we also use the human brain MR image of size 198×198 to test the convergence speed of three algorithms under 27% to 36% k-space sampling rate. We use the mean square errors of 5∗10-3. The CPU times about three algorithms are shown in Table 3. The results demonstrate that the proposed algorithm is efficient and has fast convergence speed than the SBA-TVL1 and ADMM-TV2L1.Table 3
CPU times via different algorithms and sample rates for brain image.
Sample rate 27% 28% 29% 30% 31% SBA-TVL1 1.80 1.78 1.83 1.57 1.54 ADMM-TV2L1 4.12 4.11 3.88 3.78 3.98 ADMM-TVL1C 0.31 0.13 0.12 0.09 0.11  Sample rate 32% 33% 34% 35% 36%  SBA-TVL1 1.58 1.49 1.45 1.38 1.48 ADMM-TV2L1 3.68 3.80 3.79 3.59 3.54 ADMM-TVL1C 0.11 0.09 0.09 0.10 0.11


### 3.2. Comparisons under Different Sampling Rates

In this subsection, we make a comparison to illustrate the proposal algorithm’s performance under different sampling rates from 26%-42% and list the reconstructed images in Figure3. We selected the brain MR images as k-space sample data.Figure 3
ADMM-TVL1C reconstruction for brain MR image from 26% to 42% sample rates.The first row images of Figure3 demonstrate five reconstructed brain MR images, and the second row images demonstrate five zoomed-in details on the central area of the reconstructed images. The experimental results show that the proposed algorithm is stable and can get better MR reconstructed image with the increase of sampling rate.


## 3.1. Comparisons with Other Reconstruction Algorithms

In this subsection, by using the abdominal MR image for the test, we compare the reconstruction effectiveness of three algorithms: the split Bregman algorithm for solvingTVL1 model in [28] (called SBA-TVL1), the alternating direction method of multipliers for solving TV2L1 model proposed in [27] (called ADMM-TV2L1), and the proposed algorithm in this paper: ADMM-TVL1C. In the first experiment, the original image is a human abdomen MR image of size 256×256, we simulate incomplete Fourier data with various sampling rates. In order to demonstrate the performance of the proposed algorithm, the image is corrupted with the Gaussian white noise at level σ=10 in this experiment.In Figure1, we show the reconstructed results of ADMM-TVL1C, SBA-TVL1, and ADMM-TV2L1 algorithms from 40%k-space data for the abdomen MR image. For better visual comparison, we zoom in on the central area to observe the texture features of these reconstructed images. Moreover, the SNRs of the abdomen MR image under different k-space sampling rates from 32% to 50% are given in Table 1, which show that the proposed algorithm is superior to the other two reconstruction algorithms.Table 1
SNRs via different algorithms and sample rates for abdomen image.
Sample rate 32% 34% 36% 38% 40% SBA-TVL1 25.87 25.90 25.91 26.03 26.01 ADMM-TV2L1 37.23 37.69 38.08 38.35 38.74 ADMM-TVL1C 37.92 39.04 39.41 40.08 40.37  Sample rate 42% 44% 46% 48% 50%  SBA-TVL1 26.06 26.14 26.16 26.13 26.28 ADMM-TV2L1 39.23 39.81 40.29 40.46 40.90 ADMM-TVL1C 41.49 41.89 42.37 42.82 43.00Figure 1
Abdomen MR image reconstruction via different algorithms from 40% sample rate.In the second experiment, we use a human brain MR image of size 198×198 to test the three algorithms under different sampling rates, as done in normal MR image reconstruction [29]. The SNRs of the human brain MR image under different k-space sampling rates from 32% to 50% are given in Table 2; the reconstructed images from 40% k-space data and the images zoomed in on their central areas are shown in Figure 2. The results demonstrate that the quality of the reconstructed images is improved as the sampling rate increases and that ADMM-TVL1C outperforms SBA-TVL1 and ADMM-TV2L1 in the ability to recover the textures and alleviate the staircase effect.Table 2
SNRs via different algorithms and sample rates for brain image.
Sample rate 32% 34% 36% 38% 40% SBA-TVL1 28.07 28.11 28.18 28.22 28.30 ADMM-TV2L1 40.41 40.91 41.38 41.44 42.06 ADMM-TVL1C 41.02 42.52 43.70 44.82 46.36  Sample rate 42% 44% 46% 48% 50%  SBA-TVL1 28.42 28.44 28.35 28.57 28.48 ADMM-TV2L1 42.61 42.87 43.11 43.61 43.95 ADMM-TVL1C 46.72 47.51 47.73 48.15 48.55Figure 2
Brain MR image reconstruction via different algorithms from 40% sample rate.In the third experiment, we also use the human brain MR image of size 198×198 to test the convergence speed of three algorithms under 27% to 36% k-space sampling rate. We use the mean square errors of 5∗10-3. The CPU times about three algorithms are shown in Table 3. The results demonstrate that the proposed algorithm is efficient and has fast convergence speed than the SBA-TVL1 and ADMM-TV2L1.Table 3
CPU times via different algorithms and sample rates for brain image.
Sample rate 27% 28% 29% 30% 31% SBA-TVL1 1.80 1.78 1.83 1.57 1.54 ADMM-TV2L1 4.12 4.11 3.88 3.78 3.98 ADMM-TVL1C 0.31 0.13 0.12 0.09 0.11  Sample rate 32% 33% 34% 35% 36%  SBA-TVL1 1.58 1.49 1.45 1.38 1.48 ADMM-TV2L1 3.68 3.80 3.79 3.59 3.54 ADMM-TVL1C 0.11 0.09 0.09 0.10 0.11


## 3.2. Comparisons under Different Sampling Rates

In this subsection, we make a comparison to illustrate the proposal algorithm’s performance under different sampling rates from 26%-42% and list the reconstructed images in Figure3. We selected the brain MR images as k-space sample data.Figure 3
ADMM-TVL1C reconstruction for brain MR image from 26% to 42% sample rates.The first row images of Figure3 demonstrate five reconstructed brain MR images, and the second row images demonstrate five zoomed-in details on the central area of the reconstructed images. The experimental results show that the proposed algorithm is stable and can get better MR reconstructed image with the increase of sampling rate.


## 4. Conclusion

In this paper, based on the LOT model and its deformations, we proposed a coupling model for MR image reconstruction and proposed an ADMM type algorithm to solve the proposed model. The feasibility and effectiveness of the proposed model and algorithm are illustrated from two aspects. First, to test the effectiveness of the proposed algorithm on MR image reconstruction, we compare it with SBA-TVL1 and ADMM-TV2L1 algorithms under different sampling rates. The numerical results showed that the proposed algorithm is effective on the reconstructed MR images. The reconstructed images under 40% sampling rate were also listed and showed that the proposed algorithm preserved diverse texture better than the other two algorithms. Secondly, we make a comparison under different sampling rates from 26% to 42% and listed the reconstructed images; results demonstrated that the proposed algorithm had stable performance and can get better MR reconstructed image with the increase of sampling rate. The experimental results show that proposed algorithm is efficient and has faster convergence speed. However, we do not obtain the convergence of the proposed algorithm; it will be one of our future research topics to give the convergence analysis.



---

*Source: 5490458-2018-09-05.xml*