A Novel Model and ADMM Algorithm for MR Image Reconstruction

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 augmentedLagrangian 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 stepsize 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.


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 samplespace of the spatial Fourier transformation of the object [1]. Therefore, many MR techniques [1][2][3][4][5][6][7][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][10][11][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][16][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][19][20][21][22][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: A estimation of normal vector field from initial or intermediate images; B 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 where is the partially scanned space data, |∇n| is the total variation of the unit normal vector n fl ( 1 , 2 ) for the level curves of image , ‖n‖ = ‖√ 2 1 + 2 2 ‖ 2 , and is a regularization parameter. (ii) The second step, we can obtain the reconstruction image by solving where n * is the solution of (1), 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 finding n and not efficiently utilizing the information of the operator . In this paper, we couple this two steps into one model by introducing a new variable to replace n with n := (cos , sin ) as follows: where ‖∇ ‖ 1 = ∫ Ω |∇ | . Here, ∫ Ω |∇ | replacing ∫ Ω |∇(cos , sin )| = ∫ Ω |∇n| 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 i.e., 1 − cos( − 0 ) is the equivalent infinitesimal quantity of ( − 0 ) 2 /2, where (cos 0 , sin 0 ) = ∇ /|∇ |. 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 Section 3. Some concluding remarks are given in the last section.
Some symbols are introduced in this paper. ∈ R is the objective image consisting of pixels. Φ ⊤ = ( 1 , 2 , . . . , ) ⊤ ∈ C × is usually a proper orthogonal matrix (e.g., Haar wavelet) sparsifying the underlying image u. ∈ C with ≤ is the partially scanned data. ∈ R × is the binary matrix representing the sampling pattern. Considering that the Fourier transform F is in fact × unitary matrix, = F is the partial Fourier transform operator.

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 ) ⊤ = (∇ , ∇ ) ⊤ . Then using another auxiliary variable z transforms Φ ⊤ out of the nondifferentiable norm ‖Φ ⊤ ‖ 1 ; (3) thus can be turned into the following problem: min , , , Mathematical Problems in Engineering The following augmented Lagrangian functional can be used to cope with the constrains in (5): where b = ( 1 , 2 ) ⊤ ∈ C ×2 and ∈ C 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: By applying the alternating minimization method [24], we can divide the minimization problem in (7) into the following four subproblems: Since the updates of the multipliers b and in (7) are not too complicated, in the forthcoming subsections, only how to efficiently solve each of subproblems (8) to (11) will be discussed.
As pointed out in [26], if ⊤ 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 of ( ) = (1/2)‖ − ‖ 2 2 ( , ) = ( ) + ⟨∇ ( ) , − ⟩ which is a linearization of ( ) at point 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 replacing ( ) in (14) with its quadratic approximation ( , ) and considering ∇ ( ) = ⊤ ( − ), the minimization subproblem (14) can be converted to the form By means of the Barzilai-Borwein (BB) step size selection scheme proposed in [25], the iterative step size +1 can be updated by which can be given by Considering ∇ ⊤ ∇ + ∇ ⊤ ∇ = −Δ and ΦΦ ⊤ = , the -subproblem (19) can be solved by its necessary optimality condition, given in the following equation: where 1 = 2 + − 1 Δ, Under periodic boundary condition, since the Laplace operator is block circular and can be diagonalized by Fourier transform F with the property F ⊤ = F −1 , the solution of the equation (22) can be computed by The solution of the -subproblem (15) can be obtained by using a generalized shrinkage formula, as done in [14], and is given by where

-Subproblem.
The next subproblem is -subproblem in (9), to accelerate the efficiency of the solution of , we use to replace 0 in (9), and then the problem would be changed as follows: Considering the optimality condition about , problem (27) can be solved by the following equation: To simplify the solution of , we use sin , cos 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 wherê2 = − 1 FΔF ⊤ , Mathematical Problems in Engineering 5 Solving (29) is also very simple.

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. While ‖ − −1 ‖ > ,

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-1 C 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 = 10 log 10 ( ) .
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 .

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 solving 1 model in [28] (called SBA-1 ), the alternating direction method of multipliers for solving 2 1 model proposed in [27] (called ADMM-2 1 ), and the proposed algorithm in this paper: ADMM-1 C. 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 Figure 1, we show the reconstructed results of ADMM-1 C, SBA-1 , and ADMM-2 1 algorithms from 40% -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 -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.
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-1 C outperforms SBA-1 and ADMM-2 1 in the ability to recover the textures and alleviate the staircase effect.
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

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 Figure 3. We selected the brain MR images as k-space sample data. The first row images of Figure 3 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.

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-1 and ADMM-2 1 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.

Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.