Printer Friendly

Accelerating dynamic cardiac MR imaging using structured sparse representation.

1. Introduction

Dynamic cardiac cine MR imaging aims at simultaneously providing a series of dynamic magnetic resonance image in spatial and temporal domains (x-t space) at a high frame rate. It usually acquires the k-space at each time frame and collects the raw data in the spatial frequency and temporal domain, the so called k-t space. Therefore, it is necessary to reconstruct each time frame and get a series of dynamic images. However, the relatively low acquisition speed of the dynamic MRI is an important factor to limit its application in clinics. How to accelerate k-space sampling for each time frame and reconstruct them without sacrificing spatial resolution is a challenging problem.

In recent years, many advanced techniques [1-10] were proposed to effectively address this issue and can be divided into two categories. One is based on compressed sensing (CS) theory [11,12] utilizing the sparsity in dynamic images to be reconstructed, and the other is based on the partial separable theory [13] exploiting the low-rank property of images in x-t space. The application of CS in dynamic MRI has drawn a lot of attention, since this theory demonstrates that the signal can be accurately reconstructed from a small amount of linear undersampled measurements by exploiting the inherent sparsity in signal. For example, Jung et al. [7, 9] uncovered an intriguing link between the compressed sensing and k-t BLAST/SENSE and proposed the k-t FOCUSS algorithm to achieve high spatiotemporal resolution in cardiac cine imaging. Liang et al. [5] developed k-t iterative support detection (k-t ISD) method to further utilize the detected partial support information besides the sparsity in cardiac cine images.

Recently, image restoration with patch-based sparse representations has attracted a lot of attention. The similarity of works in this topic is seeking for a more appropriate way to sparsify the image patches than conventional fixed transform. One approach is to provide additional information when using fixed transform on patches. For example, Qu et al. [14] presented to provide the sparsest representation for each image patch by estimating geometric directions. In [15], the nonlocal patches with intensity similarity, instead of those from neighbors, are grouped and then transformed by using a 3D Haar wavelet to produce sparser representation. These two patch-based methods exhibited consistent improvements in reconstruction accuracy over conventional CS-MRI methods. Another approach is based on the dictionary learning technique which aims to learn an adaptive basis from image patches and has shown impressive image restoration results [16]. The essential difference between these two approaches is that the latter uses adaptive learned dictionary instead of fixed basis, such as temporal Fourier transform, as the sparsifying transform. Ravishankar and Bresler [17] have applied the dictionary learning technique to static MR image reconstruction and obtained better reconstruction results than the state-of-the-art methods using fixed sparsifying transform. Liu et al. [18] presented to train dictionaries from the patches of the horizontal and vertical gradient images instead of the pixel domain image. The sparser training samples from the gradient images that are already sparsified by gradient operators result in sparser representation. In [19], a two-level Bregman method with dictionary learning updating is developed by applying the outer-level and inner-level Bregman iterative procedures to update the whole image and image patches, respectively. Experimental results on static MR images demonstrate the superiority of the presented algorithm to the state-of-the-art reconstruction methods.

However, dictionary learning based optimization is a large-scale and highly nonconvex problem, which requires high computational complexity. The coherence of the dictionary and the large degree of freedom may become sources of instability and errors. Structured sparse representation model was proposed to reduce the degree of freedom in the estimations and was thus more stable than conventional sparse representation model. The structured learned overcomplete dictionary, composed of a union of bases of principal component analysis (PCA), was widely used in image restoration [20,21]. Recently, Dong et al. [22] proposed nonlocally centralized sparse representation (NCSR) model for single natural image restoration, specifically, clustering image patches by K-means algorithm at first and then learning PCA subdictionary of each cluster to sparsely represent image patches. Finally, the so-called sparse coding noise (SCN) was minimized to improve the performance of sparsity-constrained image restoration. This model has gotten the satisfactory results on image deblurring, image denoising, and image super resolution.

In this work, motivated by the effective representation ability of NCSR, a novel method based on the NCSR model is proposed to accelerate dynamic cardiac MRI applications. The method utilizes structured sparse dictionary learning to adaptively represent image sequence and reduces the error between the sparse coding coefficients learned by such dictionary and true sparse coding. Improvement of the proposed method over the basic CS approach is demonstrated using retrospectively undersampled in vivo cardiac cine MR datasets.

The rest of the paper is organized as follows. In Section 2, the NCSR model is briefly described and a detailed account of structured sparse representation-based dynamic cardiac MR imaging method is provided. We present experimental validation of our method and compare it to previous state-of-the-art method in Section 3. Conclusions are drawn in Section 4.

2. Materials and Methods

2.1. Nonlocally Centralized Sparse Representation (NCSR). Image restoration often requires solving an inverse problem. It amounts to estimate original image vector x from a vector of measurements y; that is, we have

y = Hx + v, (1)

which is obtained through the noninvertible linear degradation operator H and is contaminated by the additive noise v.

Mathematically, image vector x [member of] [C.sup.N] can be represented as x [approximately equal to] [PHI] [alpha] under the sparse representation framework, where [PHI] [member of] [c.sub.N x M], n < M is a dictionary, and [alpha] [member of] [C.sup.M] represents the sparse coefficients. The sparse decomposition of x can be obtained by solving a convex 11-minimization problem:


In the scenario of image restoration, to recover x from the degraded image, y is first sparsely coded with respect to O by solving the following optimization problem:


And then x is reconstructed by [??] = [PHI][[alpha].sub.y]. In order to achieve an effective image restoration, ay are expected to be as close as possible to approach the true sparse codes ax of the original image x. Dong et al. defined the sparse coding noise (SCN) as the difference between [a.sub.y] and [a.sub.x]:

[V.sub.[alpha]] = [[alpha].sub.y] - [[alpha].sub.x]. (4)

Thus, the quality of image restoration can be improved by suppressing SCN. However, [[alpha].sub.x] is unknown so SCN cannot be directly measured. To address this issue, a good estimation [beta] of [[alpha].sub.x] is necessary. There are various ways to obtain an accurate estimate of [[alpha].sub.x]. Dong et al. tried to learn the estimate [beta] by computing the weighted average of the sparse codes of nonlocal similar patches.

The NCSR model was proposed as follows:


where dictionary [PHI] can be designed as a PCA-based structured sparse dictionary, at denotes the sparse coding vector of ith image patch on a certain subdictionary, and [lambda] is the regularization parameter controlling the tradeoff between data consistency and sparse coding noise. To solve this problem, firstly, the training patches extracted from the given image are clustered into K clusters, and a PCA subdictionary is learned for each corresponding cluster. Then one PCA subdictionary is adaptively selected to code a given patch. Finally, an iterative shrinkage algorithm [23] can be used to solve the NCSR objective function in (5).

Algorithm 1: NCSR-based dynamic cardiac cine MR imaging.

(1) Initialization:
    Set the initial estimate [[??].sup.(1)] = Y, [[??].sup.(0)] = 0 and
      the regularization parameter [lambda];
(2) Outer loop (dictionary learning and clustering): iterate
        on l = 1,2, ..., L
    (a) Update the dictionaries {[[PHI].sub.k] via fc-means and PCA;
    (b) Inner loop (clustering): iterate on j = 1,2, ..., J
        (b.1) a (j) = 1 + (j - 1) / (j - 2) + b (j) = -(j - 1)/(j + 2)
        (b.2) Compute [[g.sup.(i).sub.s] = a (j) + b
                 [(x).sup.(j).sub.s] + a (j) + b
                 [(x).sup.(j - 1).sub.s] and
                 for s = 1,2, ..., S;
        (b.3) Compute v(i) = [[MATHEMATICAL EXPRESSION NOT
                 REPRODUCIBLE IN ASCII]]], where
                 ASCII]. is the dictionary assigned to patch;
        (b.4) Compute the [[alpha].sup.(j+1).sub.i] using the
                 shrinkage operator given in (8);
        (b.5) Update the estimate {[[beta].sub.i]};
        (b.6) Image matrix estimate update: [X.sup.(j+1)] =

2.2. NCSR-Based Dynamic Cardiac MR Imaging. When the degradation operator H is the under-sampled Fourier encoding operator and y is the acquired k-space data, we can modify the above model to MR image reconstruction. Based on the NCSR model, we propose a new method to reconstruct a time series of dynamic cardiac cine MR images which have high correlations in the spatial- and temporal-domain.

We define a matrix of image series X = [[x.sub.1], ... , [x.sub.s]] [epsilon] [C.sup.NxS], whose columns are the image vectors {[x.sub.s]}, s = 1, ..., S. S denotes the number of image frames. In order to reconstruct the image series, we propose the following cardiac cine MRI reconstruction model:


where X = [PHI][alpha] and Y [member of] [C.sup.NxS] is the acquired k-space data matrix, whose columns are the vector form of k-space data of images {xs}. Matrix operator [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] performs the under-sampled Fourier encoding, and set Q = |1,2,..., S} indicates that the undersampling masks are different for each frame to enforce incoherence. For the selection of the dictionary O, we adopt the PCA-based structured sparse dictionary like in NCSR. However, the way to get the patches from the image series {[x.sub.s]} is different from NCSR. We first transform the image series {[x.sub.s]} to the image matrix X defined above and then regard the transpose vector of each row [x.sub.n], n = 1,...,N, of X as a patch vector. In other words, we learn the PCA subdictionaries along the temporal dimension to exploit the inherent correlation in dynamic image series. After getting the patches, we use the K-means algorithm to partition the patch set into K clusters {[C.sub.1],[C.sub.2], .... [C.sub.K]} and then compute the covariance matrix [[SIGMA].sub.k] of each cluster [C.sub.k]. An orthogonal transformation matrix [P.sub.k] can be obtained by applying PCA to [[SIGMA].sub.k]. We set [P.sub.k] as the PCA subdictionary which can constitute the dictionary [PHI].

The iterative shrinkage algorithm is used to solve this problem and the final image series can be obtained from the solved sparse coding vector. Specifically, at each iteration, we use the same method as in NCSR to compute [[beta].sub.i]. For each local patch, the Euclidean distance was used to search for the first P (P = 13 in our experiments) closest patches. We applied the corresponding subdictionary to these nonlocal similar patches to obtain their sparse codes. [[beta].sub.i] was estimated by computing the weighted average of these sparse codes. This nonlocal method can produce accurate enough estimates of true sparse codes. Finally, the following minimization problem can be solved for a given [[beta].sub.i]:


where [[alpha].sub.i] (j) and [[beta].sub.i] (j) are the jth elements of [alpha].sub.i], [[beta].sub.i]. We adopted the surrogate algorithm in [23] to solve (7). In the (l + 1)th iteration, the proposed shrinkage operator for the jth elements of [alpha].sub.i] is

[[alpha].sup.(l+1).sub.i](j) = [S.sub.[tau]]([v.sup.(l).sub.i,j] - [[beta].sub.j](j)) + [[beta].sub.i](j)), (8)

where [S.sub.[tau]] (x) is the classic soft thresholding operator. [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] and c is a parameter guaranteeing the convexity of the surrogate function (c = 1 in our experiments).

Since one drawback of this iterative framework is the slow convergence rate of O(1/n), we introduce an accelerated method described in [24] to achieve a fast O(1/[n.sup.2]) convergence rate. In this method, two prior iterates are used to obtain the next updated solution in the soft thresholding framework. The detailed design is described in Algorithm 1 (b.1) and (b.2).

3. Experimental Results

The feasibility of the proposed method was validated on two sets of in vivo dynamic cardiac cine data. Informed consent was obtained from the volunteer in accordance with the institutional review board policy. The full k-space data with the size of 256 x 256 x 25 (number of frequency encoding x number of phase encoding x number of frames) of the first dataset was acquired using a steady-state free precession (SSFP) sequence on a 1.5 T Philips scanner. The flip angle was 50 degrees and TE/TR = 1.7/3.45 msec. The field of view (FOV) was 345 mm x 270 mm and the slice thickness was 10 mm. Retrospective cardiac gating was used with a heart rate of 66 bpm. The second dataset was acquired on a 3 T Siemens Trio scanner (Siemens Medical Solutions, Erlangen, Germany) with a flip angle of 44 degrees and TE/TR = 42.5/1.22 msec. The fully acquired k-t measurements were of size 304 x 165 x 26 x 12. The FOV was 340 mm x 276 mm and the slice thickness was 6 mm.

The image series reconstructed from the full k-t data was used as the reference for comparison, while, for the dataset acquired using multiple coils, the image from each channel was reconstructed from full samples and combined using square root of sum-of-squares (SOS) as the reference. To simulate the undersampled k-space data, the sampling masks corresponding to reduction factors, R = 3 and 4, were generated using the function provided in the k-t FOCUSS toolbox, where the central 8 phase encoding lines were fully sampled. The fully sampled data were then retrospectively undersampled using the designed sampling masks.

The proposed NCSR-based method and k-t FOCUSS were used to reconstruct the desired image series with the same sampling patterns for a given undersampled dataset. All methods were implemented in Matlab and the code for k-t FOCUSS was obtained from Simulations run on a dual core 2.6 GHz CPU laptop with 4 GB RAM. The running time of our program is about 30 minutes. This time is relatively long due to slow speed in dictionary learning and Thmeans clustering and can be reduced by optimizing the code and utilizing GPU for acceleration.

The reconstructions with reduction factors of R = 3 and 4 and the corresponding difference images at the fourth frame of the first dynamic cardiac dataset are shown in Figures 1(a) and 1(b), respectively. In Figure 1, the first row shows the results of k-t FOCUSS (the left two) and proposed method (the right two) in R = 3, and the second row is the results of two methods in R = 4. It can be seen from the reconstructed images that the k-t FOCUSS presents more undersampling artifacts along the phase encoding direction. Our method greatly suppresses the artifacts and obtains high quality reconstructions. The superiority of structured sparse representation based method is also clearly seen in the difference images as shown by the red arrows. Figures 2(b) and 2(c) show the region-of-interest (ROI) reconstructions using k-t FOCUSS and the proposed method (from left to right) with R = 3 and 4 (from top to left) of the second dynamic cardiac dataset. We can find some artifacts appearing in the k-t FOCUSS result especially with a reduction factor of 4.

To quantify the improvement of the proposed method over k-t FOCUSS, the normalized mean-squared error (NMSE) between the reconstruction and the reference at R = 3 and 4 was calculated and plotted as a function of time frame in Figures 3(a) and 3(b), respectively. The nRMSE was calculated using the following formula:


where [x.sub.rec] is the reconstructed images from the undersampled data, x is the reference, and N is the image size. The dotted lines are for our method and dashed lines for k-t FOCUSS, respectively. Our method is seen to have a lower MSE than k-t FOCUSS for all frames at specified reduction factors.

The ability to catch the dynamic motion along temporal direction is a key factor for comparing different dynamic reconstruction methods. To evaluate the temporal fidelity, we show in Figure 4 the reconstructions in x-t plane of the first dataset with R = 3 and 4 for a fixed position in the frequency-encoding direction. It is seen that k-t FOCUSS shows some loss of contrast. In comparison, the proposed method preserves more temporal variations especially in regions indicated by red arrows.

In our algorithm, a regularization parameter [lambda] was introduced. This parameter controls the tradeoff between the data fidelity and the accuracy of the sparse codes, and it also affects the thresholds [tau] in (8). In this work, [lambda] was elaborately tuned in a parameter range. To show the effects of this parameter on final reconstructions, the curves of NMSE with respect to parameter [lambda] for the 10th frame of the first dynamic cardiac dataset at R = 3 and 4 were plotted in Figure 5. We can find that the reconstructions are relatively robust to this parameter. Results with least NMSE are obtained when A = 0.0015 with R = 3 and [lambda] = 0.002 with R = 4. In our experiments, we empirically set [lambda] = 0.002.

The convergence behavior is an important factor in evaluating the performance of the proposed method. The corresponding NMSE-iteration plots are shown in Figure 6 when R = 3 and 4 for the first dataset. It can be seen that the NMSE decreases fast at the first few iterations and then becomes flatter and reaches the convergence zone after 6 outer iterations.

From the above experimental results, we can find that our method produces more accurate reconstruction on image sequence than k-t FOCUSS. It is because we force the sparse coefficients of dictionary learning to approach the true sparse coding, which is estimated through the nonlocal similarity technique. This technique was proved to be an effective method using image redundancy and therefore the accurate sparse representation promotes the quality of reconstruction.

4. Conclusions

In this work, we propose a novel dynamic cardiac MR imaging method based on the NCSR model. This method sparsely codes the image sequence by adaptively learning PCA-based structured sparse dictionary and recovers the true sparse coding coefficients with a centralized sparse constraint, which effectively exploits the image nonlocal redundancy. An accelerated iterative shrinkage method was presented for solving the proposed model. From the experimental results from in vivo dynamic cardiac cine MR imaging, it is proved that the proposed method could produce fewer artifacts and preserve contrast than the state-of-the-art method.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.


The authors would like to thank Ye et al. for sharing their experiment materials and source codes of k-t FOCUSS. This work was partly supported by the National Natural Science Foundation of China (Grant nos. 61001179, 61102043, 61372173, and 81120108012), the Shenzhen Peacock Plan (KQCX20120816155710259), and the Basic Research Program of Shenzhen (JC201104220219A).


[1] A. Majumdar, R. K. Ward, and T. Aboulnasr, "Compressed sensing based Real-Time dynamic MRI reconstruction," IEEE Transactions on Medical Imaging, vol. 31, no. 12, pp. 2253-2266, 2012.

[2] B. Zhao, J. P. Haldar, C. Brinegar, and Z.-P. Liang, "Low rank matrix recovery for real-time cardiac MRI," in Proceedings of the 7th IEEE International Symposium on Biomedical Imaging: From Nano to Macro (ISBI '10), pp. 996-999, Rotterdam, The Netherlands, April 2010.

[3] J. P Haldar and Z.-P. Liang, "Low-rank approximations for dynamic imaging," in Proceedings of the 8th IEEE International Symposium on Biomedical Imaging: From Nano to Macro (ISBI '11), pp. 1052-1055, Chicago, Ill, USA, April 2011.

[4] S. G. Lingala, Y. Hu, E. Dibella, and M. Jacob, "Accelerated dynamic MRI exploiting sparsity and low-rank structure: k-t SLR," IEEE Transactions on Medical Imaging, vol. 30, no. 5, pp. 1042-1054, 2011.

[5] D. Liang, E. V. R. Dibella, R.-R. Chen, and L. Ying, "k-t ISD: dynamic cardiac MR imaging using compressed sensing with iterative support detection," Magnetic Resonance in Medicine, vol. 68, no. 1, pp. 41-53, 2012.

[6] M. Lustig, J. M. Santos, D. L. Donoho, and J. M. Pauly, "k-t SPARSE: high frame rate dynamic MRI exploiting spatiotemporal sparsity," in Proceedings of the 14th Annual Meeting of ISMRM, p. 2420, Seattle, Wash, USA, 2006.

[7] H. Jung, J. C. Ye, and E. Y. Kim, "Improved k-t BLAST and k-t SENSE using FOCUSS," Physics in Medicine and Biology, vol. 52, no. 11, article 018, pp. 3201-3226, 2007.

[8] U. Gamper, P Boesiger, and S. Kozerke, "Compressed sensing in dynamic MRI," Magnetic Resonance in Medicine, vol. 59, no. 2, pp. 365-373, 2008.

[9] H. Jung, K. Sung, K. S. Nayak, E. Y. Kim, and J. C. Ye, "K-t FOCUSS: a general compressed sensing framework for high resolution dynamic MRI," Magnetic Resonance in Medicine, vol. 61, no. 1, pp. 103-116, 2009.

[10] L. Feng, M. B. Srichai, R. P Lim et al., "Highly accelerated realtime cardiac cine MRI using kt SPARSE-SENSE," Magnetic Resonance in Medicine, vol. 70, no. 1, pp. 64-74, 2013.

[11] D. L. Donoho, "Compressed sensing," IEEE Transactions on Information Theory, vol. 52, no. 4, pp. 1289-1306, 2006.

[12] E. J. Candes, J. Romberg, and T. Tao, "Robust uncertainty principles: exact signal reconstruction from highly incomplete frequency information," IEEE Transactions on Information Theory, vol. 52, no. 2, pp. 489-509, 2006.

[13] Z.-P. Liang, "Spatiotemporal imaging with partially separable functions," in Proceedings of the 4th IEEE International Symposium on Biomedical Imaging: From Nano to Macro (ISBI '07), pp. 988-991, Arlington, Va, USA, April 2007.

[14] X. Qu, D. Guo, B. Ning et al., "Undersampled MRI reconstruction with patch-based directional wavelets," Magnetic Resonance Imaging, vol. 30, no. 7, pp. 964-977, 2012.

[15] X. B. Qu, Y. K. Hou, F. Lam, D. Guo, J. H. Zhong, and Z. Chen, "Magnetic resonance image reconstruction from undersampled measurements using a patch-based nonlocal operator," Medical Image Analysis, 2013.

[16] M. Aharon, M. Elad, and A. Bruckstein, "K-SVD: an algorithm for designing overcomplete dictionaries for sparse representation," IEEE Transactions on Signal Processing, vol. 54, no. 11, pp. 4311-4322, 2006.

[17] S. Ravishankar and Y. Bresler, "MR image reconstruction from highly undersampled k-space data by dictionary learning," IEEE Transactions on Medical Imaging, vol. 30, no. 5, pp. 1028-1041, 2011.

[18] Q. G. Liu, S. S. Wang, L. Ying, X. Peng, Y. J. Zhu, and D. Liang, "Adaptive dictionary learning in sparse gradient domain for image recovery," IEEE Transactions on Image Processing, vol. 22, no. 12, pp. 4652-4663, 2013.

[19] Q. G. Liu, S. S. Wang, K. Yang, J. H. Luo, Y. M. Zhu, and D. Liang, "Highly undersampled magnetic resonance imaging reconstruction using two-level Bregman method with dictionary updating," IEEE Transactions on Medical Imaging, vol. 32, no. 7, pp. 1290-1301, 2013.

[20] G. Yu, G. Sapiro, and S. Mallat, "Solving inverse problems with piecewise linear estimators: from gaussian mixture models to structured sparsity," IEEE Transactions on Image Processing, vol. 21, no. 5, pp. 2481-2499, 2012.

[21] D. Zoran and Y. Weiss, "From learning models of natural image patches to whole image restoration," in Proceedings of the IEEE International Conference on Computer Vision (ICCV '11), pp. 479-486, Barcelona, Spain, November 2011.

[22] W. Dong, L. Zhang, G. Shi, and X. Li, "Nonlocally centralized sparse representation for image restoration," IEEE Transactions on Image Processing, vol. 22, no. 4, pp. 1620-1630, 2013.

[23] I. Daubechies, M. Defrise, and C. de Mol, "An iterative thresholding algorithm for linear inverse problems with a sparsity constraint," Communications on Pure and Applied Mathematics, vol. 57, no. 11, pp. 1413-1457, 2004.

[24] K. Khare, C. J. Hardy, K. F. King, P. A. Turski, and L. Marinelli, "Accelerated MR imaging using compressive sensing with no free parameters," Magnetic Resonance in Medicine, vol. 68, no. 5, pp. 1450-1457, 2012.

Nian Cai, (1) Shengru Wang, (1,2,3) Shasha Zhu, (1) and Dong Liang (2,3)

(1) School of Information Engineering, Guangdong University of Technology, Guangzhou 510006, China

(2) Paul C. Lauterbur Research Centre for Biomedical Imaging, Shenzhen Institutes of Advanced Technology, Shenzhen 518055, China

(3) Shenzhen Key Laboratory for MRI, Guangdong, Shenzhen 518055, China

Correspondence should be addressed to Dong Liang;

Received 20 October 2013; Accepted 21 November 2013

Academic Editor: Linwei Wang
COPYRIGHT 2014 Hindawi Limited
No portion of this article can be reproduced without the express written permission from the copyright holder.
Copyright 2014 Gale, Cengage Learning. All rights reserved.

Article Details
Printer friendly Cite/link Email Feedback
Title Annotation:Research Article
Author:Cai, Nian; Wang, Shengru; Zhu, Shasha; Liang, Dong
Publication:Computational and Mathematical Methods in Medicine
Article Type:Report
Geographic Code:9CHIN
Date:Jan 1, 2014
Previous Article:Coarse-grained multifractality analysis based on structure function measurements to discriminate healthy from distressed foetuses.
Next Article:Large-scale modeling of epileptic seizures: scaling properties of two parallel neuronal network simulation algorithms.

Terms of use | Privacy policy | Copyright © 2021 Farlex, Inc. | Feedback | For webmasters |