Abstract
It is challenging and inspiring for us to achieve high spatiotemporal resolutions in dynamic cardiac magnetic resonance imaging (MRI). In this paper, we introduce two novel models and algorithms to reconstruct dynamic cardiac MRI data from under-sampled space data. In contrast to classical low-rank and sparse model, we use rank-one and transformed sparse model to exploit the correlations in the dataset. In addition, we propose projected alternative direction method (PADM) and alternative hard thresholding method (AHTM) to solve our proposed models. Numerical experiments of cardiac perfusion and cardiac cine MRI data demonstrate improvement in performance.
1. Introduction
Dynamic magnetic resonance imaging (MRI) reconstructs a temporal series of images to resolve the motion or the variation of the imaged object. It is often used in cardiac, perfusion, functional, and gastrointestinal imaging. The dynamic contrast variations in the myocardium are typically imaged 40–60 s in cardiac perfusion imaging; however, most patients cannot maintain a breath-hold for such long durations. The data in MRI, samples in space of the spatial Fourier transform of the object, are acquired sequentially in time. Therefore, achieving high spatiotemporal resolutions is challenging and inspiring in dynamic MRI due to the hardware limitations and the risk of peripheral nerve stimulation.
Compressed sensing (CS) [1, 2] has been successfully applied to accelerate MRI, for example, [3–7] and references therein. Lustig et al. [3] systematically exploited the sparsity which is implicit in MRI and developed a framework for using CS in MRI. Liang [4] utilized low-rank matrix completion to dynamic MRI by considering each temporal frame as a column of a space-time matrix. Lingala et al. [5] and Zhao et al. [6] combined CS and low-rank matrix completion which finds a solution that is both low-rank and sparse. Gao et al. [8] proposed a different method which decomposes dynamic MRI matrix as a superposition of a low-rank component and a sparse component. It is worth noting that the low-rank component can model the temporally correlated background and the sparse component can model the dynamic information that lies on the top of the background. Otazo et al. [7] presented low-rank and sparse reconstructions for dynamic MRI using joint multicoil reconstruction for Cartesian and non-Cartesian -space sampling.
Most of the above work is concerned about convex relaxation; then it comes to an interesting question whether we can deal with dynamic MRI using -minimization. Many authors have made great effort; see [9–13], for example. Blumensath and Davies [9] were the first to propose the iterative hard thresholding (IHT) to solve a type of -regularized problems and showed that the IHT method converges to a local minimizer. Beck and Eldar [10] introduced and analyzed three kinds of necessary optimality conditions for -constrained problems. Bahmani et al. [14] studied -constrained optimization in cases where nonlinear models are involved or the cost function is not quadratic. Recently, Lu and Zhang [11] presented a penalty decomposition method for solving a more general class of -constrained and -regularized minimization problems. Lu [12] considered -regularized convex cone problems and then gave an IHT method and its convergence. More recently, Pan et al. [13] provided the first- and second-order necessary and sufficient optimality conditions for -constrained optimization.
Motivated by the recent advances in the use of CS in MRI and -minimization, we propose rank-one and transformed sparse models and algorithms to significantly accelerate dynamic MRI. In contrast to classical low-rank and transformed sparse decomposition, this paper has the following merits:(i)We introduce the use of rank-one matrix in dynamic cardiac MRI data which avoids singular value decomposition (SVD).(ii)We assume the dynamic component has a sparse representation in some known sparsifying transformation, not itself.(iii)We develop new -constrained model and -regularized model, rather than their convex relaxation.(iv)We design two ADM-based algorithms, projected alternative direction method (PADM) and alternative hard thresholding method (AHTM). Numerical experiments of dynamic cardiac MRI illustrate the efficiency of our proposed algorithms.
The remainder of this paper is organised as follows. Section 2 reviews some existing related work and proposes our new models. Section 3 describes PADM and AHTM to solve our proposed models. Section 4 reports experimental results. Section 5 concludes this paper with some final remarks.
2. Model Analysis
We begin this section by introducing the low-rank and transformed sparse decomposition in [7]; that is,where is defined as the sum of all singular values and is defined as the sum of absolute values of all entries. is a real positive weighting parameter to balance the weights of rank and sparsity, and is the undersampled data. is a low-rank matrix, which can model the temporally correlated background. is a sparse matrix, which can model the dynamic information that lies on top of the background. Generally, the dynamic component has a sparse representation in some known sparsifying transformation (e.g., temporal frequency domain), hence the idea of minimizing and not itself. For acquisition with multiple receiver coils, is given by the frame-by-frame multicoil encoding operator, which performs a multiplication by coil sensitivities followed by a Fourier transform according to the sampling pattern.
Compared with a low-rank or sparse constraint, (1) can significantly increase compressibility of dynamic MRI data and thus enable high acceleration factors. However, (1) is the convex relaxation problem. The original optimization problem associated with (1) is formulated as follows:and its low-rank and sparse constrained optimization problem is as follows:where denotes the rank of , denotes the sparsity of , and is -norm, which counts the number of nonzero entries.
Note that if and reduce to identity matrix, (3) becomes a low-rank and sparse constrained optimization problem:where is given data. Zhou and Tao [15] developed GoDec method to estimate the low-rank part and the sparse part in (4). GoDec alternatively updated by SVD and by hard thresholding. To overcome the difficulty of SVD, they proposed bilateral random projections. Even so, it still costs much time.
We find sometimes that the background is exactly a rank-one matrix, instead of a general low-rank matrix, so Li et al. [16] proposed to replace low-rank matrix with rank-one matrix which avoids any SVD completely. Xiu and Kong [17] extended it to the case of tensor decomposition and showed that it performs well in surveillance video. Specifically speaking, we set , where and denotes the vector in whose entries are all , which leads to the following rank-one and transformed sparse matrix decomposition problem:and its unconstrained versionwhere is a regularized parameter.
3. Algorithm and Convergence
In this section, we present the projected alternative direction method (PADM) for (5) and alternative hard thresholding method (AHTM) for (6). Then we discuss their convergent properties.
3.1. Projected Alternative Direction Method for (5)
Let us look at the optimization problem (5). Inspired by ADM-based methods, it is easy to divide (5) into three subproblems:where is obtained by enforcing data consistency, where the aliasing artifacts corresponding to the residual in space data are subtracted from . Here, we denote these subproblems in (7) as -subproblem, -subproblem, and -subproblem, respectively.
The first -subproblem has the solutionwhere, for any , mean is defined as Before establishing the solution of the second -subproblem, we first give a proposition about sparse projection.
Proposition 1. Let . Then , where denotes all but largest absolute value components of .
Proof. Letting , we havewhere and the third equation is satisfied by introducing . Then, it holds that The proof is completed.
From Proposition 1, we derive the solution of -subproblem:Hence the scheme is summarized as the following projected alternative direction method (PADM).
In Algorithm 1, -subproblem is obtained by enforcing data consistency, so the convergence properties can be analyzed by considering the iterations of -subproblem and -subproblem. Now we will present the convergence result of PADM as follows.
| ||||||||||||||||||||||
Theorem 2. Let be the sequence generated by the above PADM, and is an accumulation point of . Then is a local minimizer of (5).
Proof. Let , for all . Indeed, one can observe thatIt follows that Hence, the sequence is nonincreasing. Since is an accumulation point of , there exists a subsequence such that . We then observe that is bounded, which together with the monotonicity of implies that is bounded below and hence exists. This observation and the continuity of yieldUsing these relations and taking limits on both sides of (13) as , we haveIn addition, from the definition of , we know that , which immediately implies . Thus, is a saddle point of (5). Let be a variant such that and let be a vector in . It then follows from (16) and the first-order optimality condition that and . Henceforth, we obtain thatwhich implies that is a local minimizer of (5).
3.2. Alternative Hard Thresholding Method for (6)
We can also divide the optimization problem (6) into the following -subproblem, -subproblem, and -subproblem:Compared with the above subsection, we find that the -subproblem and -subproblem are the same; hereafter we mainly explain how to solve the -subproblem. By applying the separability of the objective function and the technique of operator splitting, the -subproblem can be converted into corresponding single variable minimization problems. Therefore, the following proposition will give the solution of the corresponding single variable minimization problem.
Proposition 3. Suppose that is a solution of problem where is variable and is a parameter. Then,
Proof. Note that and when , .(i)If , , then .(ii)If , , then .(iii)If , , then or .Thus, we immediately obtain the proposition.
It is straightforward to establish the following result.
Proposition 4. For the general -regularized problem, where is some transform matrix and is a given matrix in . Then, the solution is described as follows:where is the hard thresholding operator defined as
Proof. The conclusion holds directly from Proposition 3 by separating the objective function.
From Proposition 4, the solution of -subproblem is characterized asBased on the above argument, we have derived the following alternative hard thresholding method (AHTM) for (6).
In the end of this subsection, we will establish the convergence result of Algorithm 2.
| ||||||||||||||||||||||
Theorem 5. Let be the sequence generated by the above AHTM, and is an accumulation point of . Then is a local minimizer of (6).
Proof. Let ; thenLike the proof of Theorem 2, we derive that is a saddle point of (6). We choose to be a vector in , and satisfies . Based on the convexity of , we have and . Using this relation, we obtain thatHence is a local minimizer of (6).
4. Numerical Experiments
In this section, we conduct numerical experiments to compare the performance of IST [7], PADM, and AHTM for dynamic cardiac MRI data. We apply all these methods on MR images: cardiac perfusion and cardiac cine, which are available at the website http://cai2r.net/resources/software/ls-reconstruction-matlab-code. All experiments are implemented in MATLAB (MathWorks, Natick, MA) on a desktop computer with Intel Core I5 2.60 GHz CPU and 8 GB of RAM.
We now address the initialization and the termination criteria for these methods. In particular, we set the stopping criterion as , whereWe choose in AHTM as suggested in [7]. However, we do not know the sparse parameter beforehand. Hence, in this paper, we initialize as one percent of pixels in each frame of .
4.1. Cardiac Perfusion
We next conduct numerical experiments to test the performance of IST [7], PADM, and AHTM for dynamic cardiac perfusion MRI data. The computational results of the three methods are presented in Figures 1–3. In detail, we use (the first row) to denote the original dynamic cardiac perfusion MRI frames, (the second row) to denote the reconstructed background parts, and (the third row) to denote the reconstructed moving parts. Comparing Figures 1–3, we can observe that PADM (see Figure 2) and AHTM (see Figure 3) are easier to recognize the pathological parts than IST (see Figure 1), which is due to the use of rather than its convex relaxation . This capacity may be useful to identify lesions that are difficult to visualize in the original images.
In addition, Table 1 reports the computational results for dynamic cardiac perfusion MRI data. Here, the number of iterations is given in the second column (i.e., Iteration). The third column (i.e., Time) lists the computing time in seconds. The relative error is reported in the last column (i.e., RelErr). It is evident to see that the “Time” and “Iteration” of PADM and AHTM are lower than these of IST. This is because we employ rank-one to avoid SVD.
4.2. Cardiac Cine
Similarly, Figures 4–6 and Table 2 demonstrate the performance of IST [7], PADM, and AHTM for dynamic cardiac cine MRI data. From the third row, we find that Figures 5 and 6 can get clearer moving parts, because we can tune the sparsity of and .
5. Conclusions
In this paper, we focused on the application of rank-one and transformed sparse decomposition for dynamic cardiac MRI. We established the projected alternative direction method and alternative hard thresholding method and showed their convergent results. Finally, numerical experiments demonstrated the efficiency and effectiveness of our proposed methods.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgment
The work was supported in part by the National Natural Science Foundation of China (11431002, 11171018).