Magnetic Resonance Imaging 32 (2014) 372–378
Contents lists available at ScienceDirect
Magnetic Resonance Imaging journal homepage: www.mrijournal.com
Adaptive fixed-point iterative shrinkage/thresholding algorithm for MR imaging reconstruction using compressed sensing Geming Wu, Shuqian Luo ⁎ School of Biomedical Engineering, Capital Medical University, Beijing, China
a r t i c l e
i n f o
Article history: Received 15 June 2013 Revised 29 August 2013 Accepted 1 December 2013 Keywords: Magnetic resonance imaging Compressed sensing Non-linear reconstruction Regularization parameter
a b s t r a c t Recently compressed sensing (CS) has been applied to under-sampling MR image reconstruction for significantly reducing signal acquisition time. To guarantee the accuracy and efficiency of the CS-based MR image reconstruction, it necessitates determining several regularization and algorithm-introduced parameters properly in practical implementations. The regularization parameter is used to control the trade-off between the sparsity of MR image and the fidelity measures of k-space data, and thus has an important effect on the reconstructed image quality. The algorithm-introduced parameters determine the global convergence rate of the algorithm itself. These parameters make CS-based MR image reconstruction a more difficult scheme than traditional Fourier-based method while implemented on a clinical MR scanner. In this paper, we propose a new approach that reveals that the regularization parameter can be taken as a threshold in a fixed-point iterative shrinkage/thresholding algorithm (FPIST) and chosen by employing minimax threshold selection method. No extra parameter is introduced by FPIST. The simulation results on synthetic and real complex-valued MRI data show that the proposed method can adaptively choose the regularization parameter and effectively achieve high reconstruction quality. The proposed method should prove very useful for practical CS-based MRI applications. © 2014 Elsevier Inc. All rights reserved.
1. Introduction In recent years, compressed sensing (CS) has attracted much attention in many areas, especially in signal processing, by offering the possibility of accurately recovering a sparse or compressible signal from fewer measurements than that suggested by the conventional Nyquist sampling theory [1–3]. Due to the implicit sparsity of MR images and the data acquisition mode of MRI meet the requirements of CS theory [4–7], CS-based MR image reconstruction has the potential to reduce scan time considerably while keeping the images of high quality. Applying CS to MR image reconstruction from under-sampled kspace data, one seeks to solve the following ℓ1-norm minimization problem [4]: 2
min k F u x−y k2 þ ΛkΨxk1 ;
x∈ℂn
ð1Þ
where x is the reconstructed image, y is the measured k-space data acquired by an MR scanner, n is the number of pixels of x and λ is the
⁎ Corresponding author at: School of Biomedical Engineering, Capital Medical University, No. 10 Xitoutiao, Youanmen, Beijing 100069, P.R. China. Tel.: +86 108 3911 566. E-mail addresses: [email protected] (G. Wu), [email protected] (S. Luo). 0730-725X/$ – see front matter © 2014 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.mri.2013.12.009
regularization parameter. Fu and Ψ denote the under-sampling Fourier transform and the sparsifying transform respectively. There are several algorithms that have been proposed in the literature to find the optimal solution for the minimization problem (1), such as non-linear conjugate gradient algorithm (NLCG) [4], two-step iterative shrinkage/thresholding algorithm (TwIST) [8], and fast iterative shrinkage/thresholding algorithm (FISTA) [9], etc. In order to implement these algorithms for CS-based MR image reconstruction, several parameters should be determined properly in advance. One of them is the regularization parameter λ, which controls the trade-off between the data fidelity of the reconstruction to the measurements and the sparsity of MR images [4]. Proper selection of λ is important in guaranteeing the reconstructed image quality. Others are imposed by the algorithms themselves, such as the parameter L in FISTA which plays the role of a step size related to the convergence rate and should be larger than an unknown Lipschitz constant. These parameters make it much difficult to implement reliable and fast CS-based MR image reconstruction on clinical MR scanners. There are two ways to determine the value of λ. One is to solve (1) with different values of λ and select the one which satisfies ‖Fux − y‖2 ≈ e, where e controls the upper bound of the recovery error [4]. Unfortunately, e is usually unknown previously. The other is to extend existing parameter choice methods [10], proposed for Tikhonov regularization problems, to (1) for finding out the optimal
G. Wu, S. Luo / Magnetic Resonance Imaging 32 (2014) 372–378
373
Fig. 1. Data distribution in the highest sub-band HH of wavelet transform of 80% under-sampled MR image: (a) MR image reconstructed from fully sampled data by inverse FFT; (b) 2D radial under-sampling scheme; (c) MR image reconstructed from 80% under-sampled data by zero-padded inverse FFT; (d) wavelet transform of under-sampled MR image, the ROI region contains the pixels of highest-frequency sub-band; (e) probability density distribution of real parts of pixels in highest-frequency sub-band; (f) probability density distribution of imaginary parts of pixels in highest-frequency sub-band; (g) probability density distribution of absolute values of pixels in highestfrequency sub-band.
value of the regularization parameter as suggested in Ref. [3]. However, most popular parameter choice approaches used in practical problems, such as Generalized Cross-Validation (GCV) [11] and L-curve (LC) [12], need to solve the target regularization problem many times to find the best regularization parameter that meets their criteria. In CS-based MR image reconstruction, almost all e x i s t i n g a l go r i t h m s f o r s o l vi n g ( 1 ) e m p l o y i t e r a t i v e schemes [4,5,8,9,13,14], which make slow processes to find the final solution. Thus it will be very numerically costly while employing these parameter choice methods to determine the regularization parameter in (1). Algorithm-introduced parameters are usually estimated empirically in advance. So the decision of the proper selection of these parameters is very important for practical CSbased MRI applications. In this paper, we present a fixed-point iterative shrinkage/ thresholding algorithm (FPIST) to solve the optimization problem (1) for under-sampling MR image reconstruction. No extra parameters are introduced in FPIST. Based on this algorithm, the regularization parameter λ is regarded as a threshold and chosen by estimating background noise level in wavelet domain. We evaluate the performance of the proposed method by simulations based on synthetic complex-valued data. The remaining parts of this article are organized as follows: Section 2 describes the proposed algorithm with adaptive regularization parameter selection method in details, and gives the framework of the proposed approach. Section 3 shows the simulation results from a numerical phantom and in vivo MRI data, and followed by the study conclusion in Section 4.
2. Methods 2.1. Fixed-point iterative shrinkage/thresholding algorithm Take the orthogonal wavelets as the sparsifying transform Ψ, and define z ≡ Ψx;
ð2Þ
F ≡ Fu þ Fu;
ð3Þ −1
Rðu;uÞ ≡ F ðu;uÞ Ψ
;
ð4Þ
where the subscript u indicates that the under-sampling operation is applied, while the subscript u denotes its complementary operation. F is the normalized full-sampling Fourier transform. Rewrite (1) as follows: 2
min kRu z−y k2 þ Λkz k1 :
z∈ℂn
ð5Þ
Since that Ru ¼ R−Ru and y is the under-sampled k-space data, we obtain 2
2
min kRz−y k2 −kRu z k2 þ Λkz k1 :
z∈ℂn
ð6Þ
The unique minimizer of (6) satisfies H
H
H
R Rz−R y−Ru Ru z þ
Λ γðz Þ ¼ 0; 2
ð7Þ
374
G. Wu, S. Luo / Magnetic Resonance Imaging 32 (2014) 372–378
where γðz Þ ¼ ∂kz k1 ¼
Table 1 MRI scan setup for in vivo experiments.
e
i argðz Þ
0
z≠0 : z¼0
ð8Þ
Note the fact that both Fourier transform and wavelet transform H are unitary, and F H u F u ¼ F F u . Eq. (7) reduces to H
z−R ðy þ Ru z Þ þ
Λ γðz Þ ¼ 0: 2
ð9Þ
The above equation can be rewritten into the following fixedpoint form H z ¼ Γτ R ðy þ Ru z Þ ;
ð10Þ
where τ ¼ Λ=2 and Γα : ℂ n → ℂ n is a complex shrinkage operator defined by Γα ðz Þ ¼ e
i argðz Þ
maxðjz j−α; 0Þ:
ð11Þ
We solve Eq. (10) by using fixed-point iterative scheme, and the basic step is given by the following iteration H z k ¼ Γτ R ðy þ Ru z k−1 Þ :
ð12Þ
FPIST is similar to iterative shrinkage/thresholding algorithm (ISTA) [14] but without extra parameters introduced. So the convergence rate of FPIST is comparable to that of ISTA, and the speed-up schemes that are used in fast IST algorithms, such as TwIST [8] and FISTA [9,15], can also be employed by FPIST. 2.2. Adaptive regularization parameter selection Eq. (12) shows that the regularization parameter λ appears only as a part of the threshold τ. Thus we can determine λ by choosing τ properly. The iterative scheme given by Eq. (12) can be rewritten as follows 2
z k ¼ arg min kz−bk−1 k2 þ 2τkz k1 ; z∈ℂn
ð13Þ
Image
Pulse sequence
Effective echo Repetition time Field of view time (ms) (ms) (mm)
Brain Fast spin echo (FSE) 130 Abdomen Fast spin echo (FSE) 110 Knee Spin echo (SE) 12
4000 4000 450
240 340 240
where bk−1 ¼ RH ðy þ Ru z k−1 Þ . The above equation indicates that each iteration of FPIST is a wavelet shrinkage denoising (WSD) process. So the optimal threshold selection schemes of WSDs can be employed to choose τ. Suppose we are given the observations bk − 1 according to the following model bk−1 ¼ w k þ k−1 ;
ð14Þ
where wk is the ideal result of Eq. (13) and ξk − 1 is the noise to be removed. It is well known that the real and imaginary parts of raw MR data are corrupted by additive Gaussian white noise ε. By taking the magnitude of complex-valued MR image, the noise can be described by a Rician distribution [16] which leads to a Rayleigh distribution in background [17]. Recall that the under-sampling Fourier transform Fu is incoherent as required in CS [3], and Fourier transform and wavelet transform will not change the characteristics of noise due to the linearity and orthogonality of these transforms. Thus ξk − 1 has the same distribution as ε. Fig. 1 illustrates the behaviors of coefficients in highest-frequency sub-band of wavelet transform of an in vivo head image, which is reconstructed from radial under-sampled MR data. It is shown that the real and imaginary parts of the high-level wavelet coefficients are Gaussian distributed and the coefficients in magnitude can be well described by Rayleigh distribution. Note that the border pixels are excluded to avoid the filtering effect of AD conversion in practical. The conditions of CS and orthogonal transforms guarantee the distribution characteristics of noise unchanged. Denote σk − 1 the standard deviation of the real or imaginary part of ξk − 1. Considering MR signal is complex-valued and defining the function ρτk ðw k Þ ¼ E kΓτk ðbk−1 Þ‐w k k22 , we employ the minimax threshold selection method proposed in Ref. [18] to derive τ at the kth iteration by defining the following minimax quantities
Λk;n
ρτk w ik ; ¼ inf sup τk ‐1 2 i 2 w ik 2σ2 k−1 n þ min w k ; 2σk−1
ð15Þ
∗ above. Suppose zk⁎ is the optimizer of and choosing τk∗ to attain Λk,n Eq. (13) with threshold τk∗. Then the overall risk satisfies
( ) n n X 2 X i 2 i 2 2 E z k −w k 2 ¼ ρτk w k ≤k;n 2σ k−1 þ min w k ; 2σ k−1 ; i¼1
i¼1
ð16Þ for all wk ∈ ℂ n. Thus the optimal selection of λ at the kth iteration can be calculated with the following equation
Λk ¼ 2τ k :
Fig. 2. 2D radial sampling trajectory with 64 projections (white lines indicate sample locations).
ð17Þ
To compute τk∗, we should determine σk − 1. In practice, σk − 1 can be estimated by taking the median absolute deviation (MAD) of the real or imaginary part of high-level coefficients in wavelet domain as
G. Wu, S. Luo / Magnetic Resonance Imaging 32 (2014) 372–378
375
Fig. 3. Complex-valued Shepp Logan phantoms corrupted by noise with different levels. From left to right, the noise levels range from 3% to 9% with the step of 2%.
HH ^ : σ k−1 ¼ 1:4826 MAD real bk−1
ð18Þ
2.3. Numerical algorithm framework Considering adaptive regularization parameter selection as discussed above and the acceleration strategy presented in FISTA [9], we propose an adaptive approach for CS-based MR image reconstruction. The framework of the proposed algorithm is summarized as follows: Algorithm. aFPIST for CS-MRI reconstruction Initialization Input under-sampled k-space data y, orthogonal wavelet transform Ψ, under-sampling Fourier transform Fu, maximum iteration number kmax, and tolerance of main loop tol. Initialize k = 0, t0 = 1, obj0 = ‖y‖22, v0 = θ, z0 = θ, where θ is null matrix. Main Repeat 1. Increase the iteration counter k by 1 2. Compute bk−1 ¼ RH ðy þRu v k−1 Þ 3. Estimate σk − 1 with Eq. (18). Let τk⁎ be the largest τk ∗ attaining Λ k,n in Eq. (15) and compute λk⁎ with Eq. (17). 4. Compute 4.1 z k ¼ Γτk ðbk −1Þ pffiffiffiffiffiffiffiffiffiffiffiffiffi 1þ 1þ4t 2k−1 4.2 t k ¼ 2 4.3 v k ¼ z k þ
t k−1 −1 tk
ðz k −z k−1 Þ
5. Calculate the objective-based consecutive error 5.1 objk = ‖Ruzk − y‖22 + λk⁎‖zk‖1 k−1 j 5.2 cerr ¼ jobjk −obj obj
3. Results To assess the efficiency of the proposed method, we carried out simulations on both synthetic and real MR data. For the synthetic case, we used the popular Shepp Logan (S-L) phantom [21,22] of size 256 × 256 as test image. Complex Gaussian white noise of zero mean and appropriate variance was added to the Fourier transform of S-L phantom to simulate noisy data. For experiments with real MRI data, three volunteers participated in this study. Fully sampled data sets of in vivo brain, abdomen and knee were acquired using a 0.35-T permanent MRI system, which was manufactured by JHMed of China, with scanning parameter setups listed in Table 1. In all experiments, radial sampling trajectory as illustrated in Fig. 2 was employed and measured data were simulated by applying a mask, which was confined to the Cartesian grids, on the Fourier transform of test images. NLCG, FISTA and aFPIST had been implemented using MATLAB (The MathWorks, Inc., Natick, MA, USA), and were employed to solve (5). The initial value of parameter L was set to 1 for FISTA. The backtracking scheme [9] was used to guarantee the convergence of FISTA, and its increasing ratio of parameter L was set to 2. Parameters introduced by NLCG were same as those used in Ref. [4]. Daubechies wavelet with four scales of orthogonal decomposition was taken as the sparsifying transform, and maximum number of iterations was limited to 100 in all algorithms for fair comparisons. For quantitative comparison, the
34 32
NL=3% NL=5% NL=7% NL=9%
30
PSNR(dB)
proposed in Refs. [19,20]. Denote bkHH− 1 as the coefficients of highestfrequency sub-band of bk − 1. MAD estimator to σk − 1 calculated with real part of wavelet coefficients is as follows
28 26
k
Until k ≤ kmax or cerr ≤ tol −1
^ x¼Ψ
zk
24 22 20
Output: reconstructed MR image ^ x In the iterative loop, τk⁎ relies on the estimation of σk − 1. After a few iterations, σk − 1 will be stable. Then the third step of iterative loop in the above algorithm can be skipped to reduce the computational cost for selecting the regularization parameter.
−6
−5
−4
−3
−2
log10λ Fig. 4. Plots of PSNR as a function of regularization parameter λ for complex-valued Shepp Logan phantom reconstruction with fixed under-sampling ratio of 74% and different noise level (NL) which ranges from 3% to 9% with the step of 2%. Upwardpointing triangles makers indicate the selected λ-values achieved by aFPIST.
376
G. Wu, S. Luo / Magnetic Resonance Imaging 32 (2014) 372–378
Fig. 5. Reconstruction results of complex-valued Shepp Logan phantom with fixed under-sampling ratio of 74% and different noise level by aFPIST. Upper row are inverse FFT reconstructions after zero padding from under-sampled data with noise level of 3%, 5%, 7%, and 9% (from left to right). Lower row are the reference and reconstructions from under-sampled data with noise level of 3%, 5%, 7%, and 9% (from left to right).
quality of a reconstructed image was measured by peak signal-tonoise ratio (PSNR), which is computed as the ratio of the peak intensity value of the reference image to the root mean square reconstruction error relative to the reference. All experiments were carried out on a laptop with a 2.4-GHz Intel Core i5 processor and 2GB memory. 3.1. Reconstruction from noisy data with different noise levels In order to study the stability of the proposed method to noise intensity, we used S-L phantom data set with complex Gaussian noise of different levels (3% to 9% with the step of 2%, as shown in Fig. 3.). The noisy data of level r were built by adding Gaussian white noise with zero mean and standard variance of r ⋅ mmax to real and imaginary parts of the Fourier transform of noise-free S-L phantom image, where mmax is the amplitude value of the brightest pixel in the simulated phantom image. In this simulation, the undersampling ratio was fixed to 74%.
34
USR=64% USR=74% USR=84%
32
The regularization parameter λ varied over a large range for FISTA, and the noise-free S-L phantom image was taken as the reference for computing PSNR. Fig. 4 gives the plots of PSNR as a function of the regularization parameter λ while employing FISTA to reconstruct S-L phantom image with noise of different levels, and upward-pointing triangles mark out the chosen λ-values by the proposed method. Reconstruction results with the proposed method are illustrated in Fig. 5. Fig. 4 illustrates that the profile of PSNR can be divided into three stages. PSNR rises slowly at the beginning, then undulates in the middle, and then drops rapidly. It presents the relationship of noise level and the threshold. Small threshold related to the noise level has no effect on improving image quality and large threshold will destroy the image structure by shrinking small signal into zero. Only when the threshold matches the noise level, the reconstructed image quality will be improved and the best image quality is achieved at the end of the second stage. Fig. 4 also shows that the regularization parameter achieving the best image quality increases with noise level. It is shown that reconstructed images with optimal regularization parameters chosen by the proposed method achieve high image quality in all four noise level cases. 3.2. Reconstruction from noisy data with different under-sampling ratios
PSNR(dB)
30 28 26 24 22
−6
−5
−4
log10 λ
−3
−2
Fig. 6. Plots of PSNR as a function of regularization parameter λ for complex-valued Shepp Logan phantom reconstruction with fixed noise level of 5% and different undersampling ratio (USR) of 64%, 74%, and 84%. Upward-pointing triangles makers indicate the selected λ-values achieved by aFPIST.
In order to study the stability of the proposed method to data with different under-sampling ratios, we repeat the previous experiment using different under-sampling ratios (64% to 84% with the step of 10%) with fixed noise level of 5%. Fig. 6 gives the plots of PSNRs of reconstructed images from different under-sampled data versus the regularization parameter λ while employing FISTA. The selected λ-values by the proposed method are also marked out with upward-pointing triangles on the plots. Reconstruction results with the proposed method are illustrated in Fig. 7. Fig. 6 shows that the regularization parameter achieving the best image quality decreases with increased under-sampling ratio. This is due to the artifact caused by using radial under-sampling trajectory that changes the distribution of noise. However, it is observed that reconstructed images with optimal regularization parameters chosen by the proposed method achieve high image quality and are robust to under-sampling ratio.
G. Wu, S. Luo / Magnetic Resonance Imaging 32 (2014) 372–378
377
Fig. 7. Reconstruction results of complex-valued Shepp Logan phantom with fixed noise level of 5% and different under-sampling ratio by aFPIST. Upper row are inverse FFT reconstructions after zero padding from 64%, 74%, and 84% under-sampled data (from left to right). Lower row are the reference and reconstructions from 64%, 74%, and 84% under-sampled data (from left to right).
3.3. Real MR image reconstruction To evaluate the practical effectiveness of aFPIST, the fully sampled MR data sets (MR scan setup listed in Table 1) were used to obtain reference images for quantitative assessment by using FFT reconstruction and modified multi-scale wavelet thresholding [23]. We used 74% down-sampled k-space data as the measurement.
Fig. 8 shows the reconstruction results by employing NLCG, FISTA and aFPIST. NLCG and FISTA used the optimal regularization parameter λ which was selected by the proposed method. The reconstruction accuracies and elapsed times are summarized in Table 2. The quantitative results show that aFPIST is faster than FISTA and NLCG even if it employs an adaptive approach to properly choose the regularization parameter and achieves the almost same
Fig. 8. Experiment with in vivo MR data. From left to right are the references, inverse FFT reconstructions from under-sampled k-space data after zero padding, NLCG reconstructions, FISTA reconstructions with backtracking scheme, and aFPIST reconstructions.
378
G. Wu, S. Luo / Magnetic Resonance Imaging 32 (2014) 372–378
References
Table 2 Comparison of reconstruction accuracies and elapsed times. Image
Method
PSNR (dB)
Elapsed time (sec)
Brain
NLCG FISTA aFPIST NLCG FISTA aFPIST NLCG FISTA aFPIST
29.98 29.90 29.91 30.94 31.04 31.04 34.61 34.84 34.84
10.27 2.30 1.50 13.22 2.16 1.66 7.23 2.59 2.03
Abdomen
Knee
reconstruction accuracy as other algorithms. Choosing the regularization parameter by using noise properties in wavelet domain contributes little computational cost to image reconstruction approach. Also, no external parameter introduced makes the proposed method simpler than existing algorithms.
4. Conclusion Compressed sensing-based MR image reconstruction offers a new potential method to significantly reduce the acquisition time for MR scanner. In this article, we present an adaptive fixed-point iterative shrinkage/thresholding algorithm (aFPIST) to make CSbased method practical in clinical MRI applications. No extra parameter is introduced by the algorithm itself and the regularization parameter is properly selected by a data-driven approach. Simulation results reveal that the aFPIST algorithm is able to obtain favorable image quality without choosing any parameter empirically and is easy to implement. In this article, the proposed method to estimate the regularization parameter is based on the assumption that Rician distribution can well describe the noise in MR images. In some cases, such as multichannel signal acquisition with sum-of-squares reconstruction for phasedarray MRI, the noise is modeled by non-central χ distribution [24,25]. For these applications, the method to estimating the regularization parameter should be changed according to the noise characteristics. Also, this data-driven approach can also be employed by other ℓ1 norm-based MR image reconstruction algorithms.
Acknowledgment This study was supported by the National Natural Science Foundation of China, Grant Nos. 61227802, 60532090 and 30770593, and by the 7th Framework Programme of the European Community, Grant Agreement Number PIRSES-GA-2009-269124.
[1] Donoho DL. Compressed sensing. IEEE Trans Inf Theory 2006;52(4):1289–306. [2] Candès EJ, Romberg J, Tao T. Robust uncertainty principles: exact signal reconstruction from highly incomplete frequency information. IEEE Trans Inf Theory 2006;52(2):489–509. [3] Eldar Yonina C, Kutyniok Gitta. Compressed sensing: theory and applications. Cambridge University Press; 2012. [4] Lustig M, Donoho DL, Pauly JM. Sparse MRI: the application of compressed sensing for rapid MR imaging. Magn Reson Med 2007;58:1182–95. [5] Guerquin-Kern M, Haberlin M, Pruessmann KP, Unser M. A fast wavelet-based reconstruction method for magnetic resonance imaging. IEEE Trans Med Imaging 2011;30(9):1649–60. [6] Chartrand R. Fast algorithms for nonconvex compressive sensing: MRI reconstruction from very few data. IEEE International Symposium on Biomedical Imaging (ISBI); 2009. p. 262–5. [7] Lustig M, Donoho DL, Santos JM, Pauly JM. Compressed sensing MRI. IEEE Trans Signal Process Mag 2008;25(2):72–82. [8] Bioucas-Dias JM, Figueiredo MAT. A new TwIST: two-step iterative shrinkage/ thresholding algorithms for image restoration. IEEE Trans Image Process 2007;16(12):2992–3004. [9] Beck A, Teboulle M. A fast iterative shrinkage-thresholding algorithm for linear inverse problems. SIAM J Imaging Sci 2009;2(1):183–202. [10] Bauer F, Lukas MA. Comparing parameter choice methods for regularization of ill-posed problems. Math Comput Simul 2011;81(9):1795–841. [11] Golub GH, Heath M, Wahba G. Generalized cross-validation as a method for choosing a good ridge parameter. Technometrics 1979;21(2):215–23. [12] Hansen PC. The L-curve and its use in the numerical treatment of inverse problems. Computational Inverse Problems in Electrocardiology. WIT Press; 2000. [13] Majumdar A, Ward RK. An algorithm for sparse MRI reconstruction by Schattern p-norm minimization. Magn Reson Imaging 2011;29(3):408–17. [14] Ye JC, Tak S, Han Y, Park HW. Projection reconstruction MR imaging using FOCUSS. Magn Reson Med 2007;57:764–75. [15] Nesterov YE. Gradient methods for minimizing composite objective function. CORE Discussion Papers 2007076. Université catholique de Louvain, Center for Operations Research and Econometrics (CORE); 2007 [Available at http://ideas. repec.org/p/cor/louvco/2007076.html]. [16] Aja-Fernandez S, Alberola-Lopez C, Westin CF. Noise and signal estimation in magnitude MRI and Rician distributed images: a LMMSE approach. IEEE Trans Image Process 2008;17(8):1383–98. [17] Koay CG, Basser PJ. Analytically exact correction scheme for signal extraction from noisy magnitude MR signals. J Magn Reson 2006;179(2):317–22. [18] Sardy S. Minimax threshold for denoising complex signals with Waveshrink. IEEE Trans Signal Process 2000;48(4):1023–8. [19] Donoho DL. De-noising by soft-thresholding. IEEE Trans Inf Theory 1995;41 (3):613–27. [20] Donoho DL, Johnstone IM. Ideal spatial adaptation by wavelet shrinkage. Biometrika 1994;81(3):425–55. [21] Van de Walle R, Barnett HH, Myers KJ, Aitbach MI, Desplanques B, Gmitro AF, et al. Reconstruction of MR images from data acquired on a general non-regular grid by pseudoinverse calculation. IEEE Trans Image Process 2000;19(12):1160–7. [22] Ouwekerk Ronald, Hopkins Johns. http://www.mathworks.com/matlabcentral/ fileexchange/1759-mriphantom. [23] Jansen M, Bultheel A. Multiple wavelet threshold estimation by generalized cross validation for images with correlated noise. IEEE Trans Image Process July 1999;8(7):947–53. [24] Larsson Erik G, Erdogmus Deniz, Yan Rui, Principe Jose C, Fitzsimmons Jeffrey R. SNR-optimality of sum-of-squares reconstruction for phased-array magnetic resonance imaging. J Magn Reson 2003;163(1):121–3. [25] Dietrich O, Raya JG, Reeder SB, Ingrisch M, Reiser MF, Schoenberg SO. Influence of multichannel combination, parallel imaging and other reconstruction techniques on MRI noise characteristics. Magn Reson Imaging 2008;26(6): 754–62.
Contents lists available at ScienceDirect
Magnetic Resonance Imaging journal homepage: www.mrijournal.com
Adaptive fixed-point iterative shrinkage/thresholding algorithm for MR imaging reconstruction using compressed sensing Geming Wu, Shuqian Luo ⁎ School of Biomedical Engineering, Capital Medical University, Beijing, China
a r t i c l e
i n f o
Article history: Received 15 June 2013 Revised 29 August 2013 Accepted 1 December 2013 Keywords: Magnetic resonance imaging Compressed sensing Non-linear reconstruction Regularization parameter
a b s t r a c t Recently compressed sensing (CS) has been applied to under-sampling MR image reconstruction for significantly reducing signal acquisition time. To guarantee the accuracy and efficiency of the CS-based MR image reconstruction, it necessitates determining several regularization and algorithm-introduced parameters properly in practical implementations. The regularization parameter is used to control the trade-off between the sparsity of MR image and the fidelity measures of k-space data, and thus has an important effect on the reconstructed image quality. The algorithm-introduced parameters determine the global convergence rate of the algorithm itself. These parameters make CS-based MR image reconstruction a more difficult scheme than traditional Fourier-based method while implemented on a clinical MR scanner. In this paper, we propose a new approach that reveals that the regularization parameter can be taken as a threshold in a fixed-point iterative shrinkage/thresholding algorithm (FPIST) and chosen by employing minimax threshold selection method. No extra parameter is introduced by FPIST. The simulation results on synthetic and real complex-valued MRI data show that the proposed method can adaptively choose the regularization parameter and effectively achieve high reconstruction quality. The proposed method should prove very useful for practical CS-based MRI applications. © 2014 Elsevier Inc. All rights reserved.
1. Introduction In recent years, compressed sensing (CS) has attracted much attention in many areas, especially in signal processing, by offering the possibility of accurately recovering a sparse or compressible signal from fewer measurements than that suggested by the conventional Nyquist sampling theory [1–3]. Due to the implicit sparsity of MR images and the data acquisition mode of MRI meet the requirements of CS theory [4–7], CS-based MR image reconstruction has the potential to reduce scan time considerably while keeping the images of high quality. Applying CS to MR image reconstruction from under-sampled kspace data, one seeks to solve the following ℓ1-norm minimization problem [4]: 2
min k F u x−y k2 þ ΛkΨxk1 ;
x∈ℂn
ð1Þ
where x is the reconstructed image, y is the measured k-space data acquired by an MR scanner, n is the number of pixels of x and λ is the
⁎ Corresponding author at: School of Biomedical Engineering, Capital Medical University, No. 10 Xitoutiao, Youanmen, Beijing 100069, P.R. China. Tel.: +86 108 3911 566. E-mail addresses: [email protected] (G. Wu), [email protected] (S. Luo). 0730-725X/$ – see front matter © 2014 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.mri.2013.12.009
regularization parameter. Fu and Ψ denote the under-sampling Fourier transform and the sparsifying transform respectively. There are several algorithms that have been proposed in the literature to find the optimal solution for the minimization problem (1), such as non-linear conjugate gradient algorithm (NLCG) [4], two-step iterative shrinkage/thresholding algorithm (TwIST) [8], and fast iterative shrinkage/thresholding algorithm (FISTA) [9], etc. In order to implement these algorithms for CS-based MR image reconstruction, several parameters should be determined properly in advance. One of them is the regularization parameter λ, which controls the trade-off between the data fidelity of the reconstruction to the measurements and the sparsity of MR images [4]. Proper selection of λ is important in guaranteeing the reconstructed image quality. Others are imposed by the algorithms themselves, such as the parameter L in FISTA which plays the role of a step size related to the convergence rate and should be larger than an unknown Lipschitz constant. These parameters make it much difficult to implement reliable and fast CS-based MR image reconstruction on clinical MR scanners. There are two ways to determine the value of λ. One is to solve (1) with different values of λ and select the one which satisfies ‖Fux − y‖2 ≈ e, where e controls the upper bound of the recovery error [4]. Unfortunately, e is usually unknown previously. The other is to extend existing parameter choice methods [10], proposed for Tikhonov regularization problems, to (1) for finding out the optimal
G. Wu, S. Luo / Magnetic Resonance Imaging 32 (2014) 372–378
373
Fig. 1. Data distribution in the highest sub-band HH of wavelet transform of 80% under-sampled MR image: (a) MR image reconstructed from fully sampled data by inverse FFT; (b) 2D radial under-sampling scheme; (c) MR image reconstructed from 80% under-sampled data by zero-padded inverse FFT; (d) wavelet transform of under-sampled MR image, the ROI region contains the pixels of highest-frequency sub-band; (e) probability density distribution of real parts of pixels in highest-frequency sub-band; (f) probability density distribution of imaginary parts of pixels in highest-frequency sub-band; (g) probability density distribution of absolute values of pixels in highestfrequency sub-band.
value of the regularization parameter as suggested in Ref. [3]. However, most popular parameter choice approaches used in practical problems, such as Generalized Cross-Validation (GCV) [11] and L-curve (LC) [12], need to solve the target regularization problem many times to find the best regularization parameter that meets their criteria. In CS-based MR image reconstruction, almost all e x i s t i n g a l go r i t h m s f o r s o l vi n g ( 1 ) e m p l o y i t e r a t i v e schemes [4,5,8,9,13,14], which make slow processes to find the final solution. Thus it will be very numerically costly while employing these parameter choice methods to determine the regularization parameter in (1). Algorithm-introduced parameters are usually estimated empirically in advance. So the decision of the proper selection of these parameters is very important for practical CSbased MRI applications. In this paper, we present a fixed-point iterative shrinkage/ thresholding algorithm (FPIST) to solve the optimization problem (1) for under-sampling MR image reconstruction. No extra parameters are introduced in FPIST. Based on this algorithm, the regularization parameter λ is regarded as a threshold and chosen by estimating background noise level in wavelet domain. We evaluate the performance of the proposed method by simulations based on synthetic complex-valued data. The remaining parts of this article are organized as follows: Section 2 describes the proposed algorithm with adaptive regularization parameter selection method in details, and gives the framework of the proposed approach. Section 3 shows the simulation results from a numerical phantom and in vivo MRI data, and followed by the study conclusion in Section 4.
2. Methods 2.1. Fixed-point iterative shrinkage/thresholding algorithm Take the orthogonal wavelets as the sparsifying transform Ψ, and define z ≡ Ψx;
ð2Þ
F ≡ Fu þ Fu;
ð3Þ −1
Rðu;uÞ ≡ F ðu;uÞ Ψ
;
ð4Þ
where the subscript u indicates that the under-sampling operation is applied, while the subscript u denotes its complementary operation. F is the normalized full-sampling Fourier transform. Rewrite (1) as follows: 2
min kRu z−y k2 þ Λkz k1 :
z∈ℂn
ð5Þ
Since that Ru ¼ R−Ru and y is the under-sampled k-space data, we obtain 2
2
min kRz−y k2 −kRu z k2 þ Λkz k1 :
z∈ℂn
ð6Þ
The unique minimizer of (6) satisfies H
H
H
R Rz−R y−Ru Ru z þ
Λ γðz Þ ¼ 0; 2
ð7Þ
374
G. Wu, S. Luo / Magnetic Resonance Imaging 32 (2014) 372–378
where γðz Þ ¼ ∂kz k1 ¼
Table 1 MRI scan setup for in vivo experiments.
e
i argðz Þ
0
z≠0 : z¼0
ð8Þ
Note the fact that both Fourier transform and wavelet transform H are unitary, and F H u F u ¼ F F u . Eq. (7) reduces to H
z−R ðy þ Ru z Þ þ
Λ γðz Þ ¼ 0: 2
ð9Þ
The above equation can be rewritten into the following fixedpoint form H z ¼ Γτ R ðy þ Ru z Þ ;
ð10Þ
where τ ¼ Λ=2 and Γα : ℂ n → ℂ n is a complex shrinkage operator defined by Γα ðz Þ ¼ e
i argðz Þ
maxðjz j−α; 0Þ:
ð11Þ
We solve Eq. (10) by using fixed-point iterative scheme, and the basic step is given by the following iteration H z k ¼ Γτ R ðy þ Ru z k−1 Þ :
ð12Þ
FPIST is similar to iterative shrinkage/thresholding algorithm (ISTA) [14] but without extra parameters introduced. So the convergence rate of FPIST is comparable to that of ISTA, and the speed-up schemes that are used in fast IST algorithms, such as TwIST [8] and FISTA [9,15], can also be employed by FPIST. 2.2. Adaptive regularization parameter selection Eq. (12) shows that the regularization parameter λ appears only as a part of the threshold τ. Thus we can determine λ by choosing τ properly. The iterative scheme given by Eq. (12) can be rewritten as follows 2
z k ¼ arg min kz−bk−1 k2 þ 2τkz k1 ; z∈ℂn
ð13Þ
Image
Pulse sequence
Effective echo Repetition time Field of view time (ms) (ms) (mm)
Brain Fast spin echo (FSE) 130 Abdomen Fast spin echo (FSE) 110 Knee Spin echo (SE) 12
4000 4000 450
240 340 240
where bk−1 ¼ RH ðy þ Ru z k−1 Þ . The above equation indicates that each iteration of FPIST is a wavelet shrinkage denoising (WSD) process. So the optimal threshold selection schemes of WSDs can be employed to choose τ. Suppose we are given the observations bk − 1 according to the following model bk−1 ¼ w k þ k−1 ;
ð14Þ
where wk is the ideal result of Eq. (13) and ξk − 1 is the noise to be removed. It is well known that the real and imaginary parts of raw MR data are corrupted by additive Gaussian white noise ε. By taking the magnitude of complex-valued MR image, the noise can be described by a Rician distribution [16] which leads to a Rayleigh distribution in background [17]. Recall that the under-sampling Fourier transform Fu is incoherent as required in CS [3], and Fourier transform and wavelet transform will not change the characteristics of noise due to the linearity and orthogonality of these transforms. Thus ξk − 1 has the same distribution as ε. Fig. 1 illustrates the behaviors of coefficients in highest-frequency sub-band of wavelet transform of an in vivo head image, which is reconstructed from radial under-sampled MR data. It is shown that the real and imaginary parts of the high-level wavelet coefficients are Gaussian distributed and the coefficients in magnitude can be well described by Rayleigh distribution. Note that the border pixels are excluded to avoid the filtering effect of AD conversion in practical. The conditions of CS and orthogonal transforms guarantee the distribution characteristics of noise unchanged. Denote σk − 1 the standard deviation of the real or imaginary part of ξk − 1. Considering MR signal is complex-valued and defining the function ρτk ðw k Þ ¼ E kΓτk ðbk−1 Þ‐w k k22 , we employ the minimax threshold selection method proposed in Ref. [18] to derive τ at the kth iteration by defining the following minimax quantities
Λk;n
ρτk w ik ; ¼ inf sup τk ‐1 2 i 2 w ik 2σ2 k−1 n þ min w k ; 2σk−1
ð15Þ
∗ above. Suppose zk⁎ is the optimizer of and choosing τk∗ to attain Λk,n Eq. (13) with threshold τk∗. Then the overall risk satisfies
( ) n n X 2 X i 2 i 2 2 E z k −w k 2 ¼ ρτk w k ≤k;n 2σ k−1 þ min w k ; 2σ k−1 ; i¼1
i¼1
ð16Þ for all wk ∈ ℂ n. Thus the optimal selection of λ at the kth iteration can be calculated with the following equation
Λk ¼ 2τ k :
Fig. 2. 2D radial sampling trajectory with 64 projections (white lines indicate sample locations).
ð17Þ
To compute τk∗, we should determine σk − 1. In practice, σk − 1 can be estimated by taking the median absolute deviation (MAD) of the real or imaginary part of high-level coefficients in wavelet domain as
G. Wu, S. Luo / Magnetic Resonance Imaging 32 (2014) 372–378
375
Fig. 3. Complex-valued Shepp Logan phantoms corrupted by noise with different levels. From left to right, the noise levels range from 3% to 9% with the step of 2%.
HH ^ : σ k−1 ¼ 1:4826 MAD real bk−1
ð18Þ
2.3. Numerical algorithm framework Considering adaptive regularization parameter selection as discussed above and the acceleration strategy presented in FISTA [9], we propose an adaptive approach for CS-based MR image reconstruction. The framework of the proposed algorithm is summarized as follows: Algorithm. aFPIST for CS-MRI reconstruction Initialization Input under-sampled k-space data y, orthogonal wavelet transform Ψ, under-sampling Fourier transform Fu, maximum iteration number kmax, and tolerance of main loop tol. Initialize k = 0, t0 = 1, obj0 = ‖y‖22, v0 = θ, z0 = θ, where θ is null matrix. Main Repeat 1. Increase the iteration counter k by 1 2. Compute bk−1 ¼ RH ðy þRu v k−1 Þ 3. Estimate σk − 1 with Eq. (18). Let τk⁎ be the largest τk ∗ attaining Λ k,n in Eq. (15) and compute λk⁎ with Eq. (17). 4. Compute 4.1 z k ¼ Γτk ðbk −1Þ pffiffiffiffiffiffiffiffiffiffiffiffiffi 1þ 1þ4t 2k−1 4.2 t k ¼ 2 4.3 v k ¼ z k þ
t k−1 −1 tk
ðz k −z k−1 Þ
5. Calculate the objective-based consecutive error 5.1 objk = ‖Ruzk − y‖22 + λk⁎‖zk‖1 k−1 j 5.2 cerr ¼ jobjk −obj obj
3. Results To assess the efficiency of the proposed method, we carried out simulations on both synthetic and real MR data. For the synthetic case, we used the popular Shepp Logan (S-L) phantom [21,22] of size 256 × 256 as test image. Complex Gaussian white noise of zero mean and appropriate variance was added to the Fourier transform of S-L phantom to simulate noisy data. For experiments with real MRI data, three volunteers participated in this study. Fully sampled data sets of in vivo brain, abdomen and knee were acquired using a 0.35-T permanent MRI system, which was manufactured by JHMed of China, with scanning parameter setups listed in Table 1. In all experiments, radial sampling trajectory as illustrated in Fig. 2 was employed and measured data were simulated by applying a mask, which was confined to the Cartesian grids, on the Fourier transform of test images. NLCG, FISTA and aFPIST had been implemented using MATLAB (The MathWorks, Inc., Natick, MA, USA), and were employed to solve (5). The initial value of parameter L was set to 1 for FISTA. The backtracking scheme [9] was used to guarantee the convergence of FISTA, and its increasing ratio of parameter L was set to 2. Parameters introduced by NLCG were same as those used in Ref. [4]. Daubechies wavelet with four scales of orthogonal decomposition was taken as the sparsifying transform, and maximum number of iterations was limited to 100 in all algorithms for fair comparisons. For quantitative comparison, the
34 32
NL=3% NL=5% NL=7% NL=9%
30
PSNR(dB)
proposed in Refs. [19,20]. Denote bkHH− 1 as the coefficients of highestfrequency sub-band of bk − 1. MAD estimator to σk − 1 calculated with real part of wavelet coefficients is as follows
28 26
k
Until k ≤ kmax or cerr ≤ tol −1
^ x¼Ψ
zk
24 22 20
Output: reconstructed MR image ^ x In the iterative loop, τk⁎ relies on the estimation of σk − 1. After a few iterations, σk − 1 will be stable. Then the third step of iterative loop in the above algorithm can be skipped to reduce the computational cost for selecting the regularization parameter.
−6
−5
−4
−3
−2
log10λ Fig. 4. Plots of PSNR as a function of regularization parameter λ for complex-valued Shepp Logan phantom reconstruction with fixed under-sampling ratio of 74% and different noise level (NL) which ranges from 3% to 9% with the step of 2%. Upwardpointing triangles makers indicate the selected λ-values achieved by aFPIST.
376
G. Wu, S. Luo / Magnetic Resonance Imaging 32 (2014) 372–378
Fig. 5. Reconstruction results of complex-valued Shepp Logan phantom with fixed under-sampling ratio of 74% and different noise level by aFPIST. Upper row are inverse FFT reconstructions after zero padding from under-sampled data with noise level of 3%, 5%, 7%, and 9% (from left to right). Lower row are the reference and reconstructions from under-sampled data with noise level of 3%, 5%, 7%, and 9% (from left to right).
quality of a reconstructed image was measured by peak signal-tonoise ratio (PSNR), which is computed as the ratio of the peak intensity value of the reference image to the root mean square reconstruction error relative to the reference. All experiments were carried out on a laptop with a 2.4-GHz Intel Core i5 processor and 2GB memory. 3.1. Reconstruction from noisy data with different noise levels In order to study the stability of the proposed method to noise intensity, we used S-L phantom data set with complex Gaussian noise of different levels (3% to 9% with the step of 2%, as shown in Fig. 3.). The noisy data of level r were built by adding Gaussian white noise with zero mean and standard variance of r ⋅ mmax to real and imaginary parts of the Fourier transform of noise-free S-L phantom image, where mmax is the amplitude value of the brightest pixel in the simulated phantom image. In this simulation, the undersampling ratio was fixed to 74%.
34
USR=64% USR=74% USR=84%
32
The regularization parameter λ varied over a large range for FISTA, and the noise-free S-L phantom image was taken as the reference for computing PSNR. Fig. 4 gives the plots of PSNR as a function of the regularization parameter λ while employing FISTA to reconstruct S-L phantom image with noise of different levels, and upward-pointing triangles mark out the chosen λ-values by the proposed method. Reconstruction results with the proposed method are illustrated in Fig. 5. Fig. 4 illustrates that the profile of PSNR can be divided into three stages. PSNR rises slowly at the beginning, then undulates in the middle, and then drops rapidly. It presents the relationship of noise level and the threshold. Small threshold related to the noise level has no effect on improving image quality and large threshold will destroy the image structure by shrinking small signal into zero. Only when the threshold matches the noise level, the reconstructed image quality will be improved and the best image quality is achieved at the end of the second stage. Fig. 4 also shows that the regularization parameter achieving the best image quality increases with noise level. It is shown that reconstructed images with optimal regularization parameters chosen by the proposed method achieve high image quality in all four noise level cases. 3.2. Reconstruction from noisy data with different under-sampling ratios
PSNR(dB)
30 28 26 24 22
−6
−5
−4
log10 λ
−3
−2
Fig. 6. Plots of PSNR as a function of regularization parameter λ for complex-valued Shepp Logan phantom reconstruction with fixed noise level of 5% and different undersampling ratio (USR) of 64%, 74%, and 84%. Upward-pointing triangles makers indicate the selected λ-values achieved by aFPIST.
In order to study the stability of the proposed method to data with different under-sampling ratios, we repeat the previous experiment using different under-sampling ratios (64% to 84% with the step of 10%) with fixed noise level of 5%. Fig. 6 gives the plots of PSNRs of reconstructed images from different under-sampled data versus the regularization parameter λ while employing FISTA. The selected λ-values by the proposed method are also marked out with upward-pointing triangles on the plots. Reconstruction results with the proposed method are illustrated in Fig. 7. Fig. 6 shows that the regularization parameter achieving the best image quality decreases with increased under-sampling ratio. This is due to the artifact caused by using radial under-sampling trajectory that changes the distribution of noise. However, it is observed that reconstructed images with optimal regularization parameters chosen by the proposed method achieve high image quality and are robust to under-sampling ratio.
G. Wu, S. Luo / Magnetic Resonance Imaging 32 (2014) 372–378
377
Fig. 7. Reconstruction results of complex-valued Shepp Logan phantom with fixed noise level of 5% and different under-sampling ratio by aFPIST. Upper row are inverse FFT reconstructions after zero padding from 64%, 74%, and 84% under-sampled data (from left to right). Lower row are the reference and reconstructions from 64%, 74%, and 84% under-sampled data (from left to right).
3.3. Real MR image reconstruction To evaluate the practical effectiveness of aFPIST, the fully sampled MR data sets (MR scan setup listed in Table 1) were used to obtain reference images for quantitative assessment by using FFT reconstruction and modified multi-scale wavelet thresholding [23]. We used 74% down-sampled k-space data as the measurement.
Fig. 8 shows the reconstruction results by employing NLCG, FISTA and aFPIST. NLCG and FISTA used the optimal regularization parameter λ which was selected by the proposed method. The reconstruction accuracies and elapsed times are summarized in Table 2. The quantitative results show that aFPIST is faster than FISTA and NLCG even if it employs an adaptive approach to properly choose the regularization parameter and achieves the almost same
Fig. 8. Experiment with in vivo MR data. From left to right are the references, inverse FFT reconstructions from under-sampled k-space data after zero padding, NLCG reconstructions, FISTA reconstructions with backtracking scheme, and aFPIST reconstructions.
378
G. Wu, S. Luo / Magnetic Resonance Imaging 32 (2014) 372–378
References
Table 2 Comparison of reconstruction accuracies and elapsed times. Image
Method
PSNR (dB)
Elapsed time (sec)
Brain
NLCG FISTA aFPIST NLCG FISTA aFPIST NLCG FISTA aFPIST
29.98 29.90 29.91 30.94 31.04 31.04 34.61 34.84 34.84
10.27 2.30 1.50 13.22 2.16 1.66 7.23 2.59 2.03
Abdomen
Knee
reconstruction accuracy as other algorithms. Choosing the regularization parameter by using noise properties in wavelet domain contributes little computational cost to image reconstruction approach. Also, no external parameter introduced makes the proposed method simpler than existing algorithms.
4. Conclusion Compressed sensing-based MR image reconstruction offers a new potential method to significantly reduce the acquisition time for MR scanner. In this article, we present an adaptive fixed-point iterative shrinkage/thresholding algorithm (aFPIST) to make CSbased method practical in clinical MRI applications. No extra parameter is introduced by the algorithm itself and the regularization parameter is properly selected by a data-driven approach. Simulation results reveal that the aFPIST algorithm is able to obtain favorable image quality without choosing any parameter empirically and is easy to implement. In this article, the proposed method to estimate the regularization parameter is based on the assumption that Rician distribution can well describe the noise in MR images. In some cases, such as multichannel signal acquisition with sum-of-squares reconstruction for phasedarray MRI, the noise is modeled by non-central χ distribution [24,25]. For these applications, the method to estimating the regularization parameter should be changed according to the noise characteristics. Also, this data-driven approach can also be employed by other ℓ1 norm-based MR image reconstruction algorithms.
Acknowledgment This study was supported by the National Natural Science Foundation of China, Grant Nos. 61227802, 60532090 and 30770593, and by the 7th Framework Programme of the European Community, Grant Agreement Number PIRSES-GA-2009-269124.
[1] Donoho DL. Compressed sensing. IEEE Trans Inf Theory 2006;52(4):1289–306. [2] Candès EJ, Romberg J, Tao T. Robust uncertainty principles: exact signal reconstruction from highly incomplete frequency information. IEEE Trans Inf Theory 2006;52(2):489–509. [3] Eldar Yonina C, Kutyniok Gitta. Compressed sensing: theory and applications. Cambridge University Press; 2012. [4] Lustig M, Donoho DL, Pauly JM. Sparse MRI: the application of compressed sensing for rapid MR imaging. Magn Reson Med 2007;58:1182–95. [5] Guerquin-Kern M, Haberlin M, Pruessmann KP, Unser M. A fast wavelet-based reconstruction method for magnetic resonance imaging. IEEE Trans Med Imaging 2011;30(9):1649–60. [6] Chartrand R. Fast algorithms for nonconvex compressive sensing: MRI reconstruction from very few data. IEEE International Symposium on Biomedical Imaging (ISBI); 2009. p. 262–5. [7] Lustig M, Donoho DL, Santos JM, Pauly JM. Compressed sensing MRI. IEEE Trans Signal Process Mag 2008;25(2):72–82. [8] Bioucas-Dias JM, Figueiredo MAT. A new TwIST: two-step iterative shrinkage/ thresholding algorithms for image restoration. IEEE Trans Image Process 2007;16(12):2992–3004. [9] Beck A, Teboulle M. A fast iterative shrinkage-thresholding algorithm for linear inverse problems. SIAM J Imaging Sci 2009;2(1):183–202. [10] Bauer F, Lukas MA. Comparing parameter choice methods for regularization of ill-posed problems. Math Comput Simul 2011;81(9):1795–841. [11] Golub GH, Heath M, Wahba G. Generalized cross-validation as a method for choosing a good ridge parameter. Technometrics 1979;21(2):215–23. [12] Hansen PC. The L-curve and its use in the numerical treatment of inverse problems. Computational Inverse Problems in Electrocardiology. WIT Press; 2000. [13] Majumdar A, Ward RK. An algorithm for sparse MRI reconstruction by Schattern p-norm minimization. Magn Reson Imaging 2011;29(3):408–17. [14] Ye JC, Tak S, Han Y, Park HW. Projection reconstruction MR imaging using FOCUSS. Magn Reson Med 2007;57:764–75. [15] Nesterov YE. Gradient methods for minimizing composite objective function. CORE Discussion Papers 2007076. Université catholique de Louvain, Center for Operations Research and Econometrics (CORE); 2007 [Available at http://ideas. repec.org/p/cor/louvco/2007076.html]. [16] Aja-Fernandez S, Alberola-Lopez C, Westin CF. Noise and signal estimation in magnitude MRI and Rician distributed images: a LMMSE approach. IEEE Trans Image Process 2008;17(8):1383–98. [17] Koay CG, Basser PJ. Analytically exact correction scheme for signal extraction from noisy magnitude MR signals. J Magn Reson 2006;179(2):317–22. [18] Sardy S. Minimax threshold for denoising complex signals with Waveshrink. IEEE Trans Signal Process 2000;48(4):1023–8. [19] Donoho DL. De-noising by soft-thresholding. IEEE Trans Inf Theory 1995;41 (3):613–27. [20] Donoho DL, Johnstone IM. Ideal spatial adaptation by wavelet shrinkage. Biometrika 1994;81(3):425–55. [21] Van de Walle R, Barnett HH, Myers KJ, Aitbach MI, Desplanques B, Gmitro AF, et al. Reconstruction of MR images from data acquired on a general non-regular grid by pseudoinverse calculation. IEEE Trans Image Process 2000;19(12):1160–7. [22] Ouwekerk Ronald, Hopkins Johns. http://www.mathworks.com/matlabcentral/ fileexchange/1759-mriphantom. [23] Jansen M, Bultheel A. Multiple wavelet threshold estimation by generalized cross validation for images with correlated noise. IEEE Trans Image Process July 1999;8(7):947–53. [24] Larsson Erik G, Erdogmus Deniz, Yan Rui, Principe Jose C, Fitzsimmons Jeffrey R. SNR-optimality of sum-of-squares reconstruction for phased-array magnetic resonance imaging. J Magn Reson 2003;163(1):121–3. [25] Dietrich O, Raya JG, Reeder SB, Ingrisch M, Reiser MF, Schoenberg SO. Influence of multichannel combination, parallel imaging and other reconstruction techniques on MRI noise characteristics. Magn Reson Imaging 2008;26(6): 754–62.
