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.
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 sample -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 where is the partially scanned space data, is the total variation of the unit normal vector for the level curves of image , , and is a regularization parameter.(ii)The second step, we can obtain the reconstruction image by solving where 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 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 as follows:where . Here, replacing mainly uses the fact that moving the and 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., is the equivalent infinitesimal quantity of , where . 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. is the objective image consisting of pixels. is usually a proper orthogonal matrix (e.g., Haar wavelet) sparsifying the underlying image u. with is the partially scanned data. is the binary matrix representing the sampling pattern. Considering that the Fourier transform is in fact unitary matrix, 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 . Then using another auxiliary variable z transforms out of the nondifferentiable norm ; (3) thus can be turned into the following problem:
The following augmented Lagrangian functional can be used to cope with the constrains in (5):where and are the Lagrangian multipliers and and 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 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.
2.1. -Subproblem
The first problem is the u-subproblem in (8), which can be solved by the split Bregman iterative scheme. By using variable to replace in (8) subject to and applying the augmented Lagrangian technique, the split Bregman iteration for (8) can be described as follows:where the variable is the Lagrangian multiplier and is positive penalty parameter.
By applying the alternating minimization scheme to the minimization problem (12) and combining a simplified calculation of , problems (12)-(13) are converted to the following form: where .
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 which is a linearization of at point plus a proximity term penalized by parameter . 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 can be updated bywhich can be given by
Considering and , the -subproblem (19) can be solved by its necessary optimality condition, given in the following equation:where ,
Under periodic boundary condition, since the Laplace operator is block circular and can be diagonalized by Fourier transform with the property , the solution of the equation (22) can be computed bywhere . Solving (24) is very simple because is diagonal.
The solution of the -subproblem (15) can be obtained by using a generalized shrinkage formula, as done in [14], and is given bywhere
2.2. -Subproblem
The next subproblem is -subproblem in (9), to accelerate the efficiency of the solution of , we use to replace 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 to replace in (28), respectively, so as to reduce the complexity of calculating . Combining the diagonalization technique of the Fourier transform, (28) is converted to the formwhere ,
Solving (29) is also very simple.
2.3. Shrink-Subproblem
The remaining subproblems would be discussed in this subsection, including (10) about -subproblem and (11) about -subproblem. By a generalized shrinkage formula, as done in [14], we can obtain the solution of the -subproblem (10) as follows:where The z-subproblem (11) can be solved easily by the shrink operator, defined by , in the following:
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-C). Initialize: , , ,;
While 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-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 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 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 , and the parameters , were chosen by and .
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 solving model in [28] (called SBA-), the alternating direction method of multipliers for solving model proposed in [27] (called ADMM-), and the proposed algorithm in this paper: ADMM-C. In the first experiment, the original image is a human abdomen MR image of size 256256, 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 in this experiment.
In Figure 1, we show the reconstructed results of ADMM-C, SBA-, and ADMM- 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 to 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 198198 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-C outperforms SBA- and ADMM- 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 198198 to test the convergence speed of three algorithms under 27% to 36% k-space sampling rate. We use the mean square errors of . 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- and ADMM-.
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 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.
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- and ADMM- 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.
Conflicts of Interest
The authors declare that they have no conflicts of interest.