110
IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 62, NO. 1, JANUARY 2015
Compressed Sensing MRI via Two-stage Reconstruction Yang Yang, Feng Liu∗ , Wenlong Xu, and Stuart Crozier, Member, IEEE
Abstract—Compressed sensing (CS) has been applied to magnetic resonance imaging for the acceleration of data collection. However, existing CS techniques usually produce images with residual artifacts, particularly at high reduction factors. In this paper, we propose a novel, two-stage reconstruction scheme, which takes advantage of the properties of k-space data and undersampling patterns that are useful in CS. In this algorithm, the under-sampled k-space data is segmented into low-frequency and high-frequency domains. Then, in stage one, using dense measurements, the low-frequency region of k-space data is faithfully reconstructed. The fully reconstituted low-frequency k-space data from the first stage is then combined with the high-frequency k-space data to complete the second stage reconstruction of the whole of k-space. With this two-stage approach, each reconstruction inherently incorporates a lower data under-sampling rate and more homogeneous signal magnitudes than conventional approaches. Because the restricted isometric property is easier to satisfy, the reconstruction consequently produces lower residual errors at each step. Compared with a conventional CS reconstruction, for the cases of cardiac cine, brain and angiogram imaging, the proposed method achieves a more accurate reconstruction with an improvement of 2–4 dB in peak signal-to-noise ratio respectively, using reduction factors of up to 6. Index Terms—Compressed sensing (CS), K-space segmentation, stationary magnetic resonance imaging (MRI), two-stage reconstruction.
I. INTRODUCTION OMPRESSED sensing (CS), a relatively new signal processing theory, has been successfully applied to accelerate the scan speed of magnetic resonance imaging (MRI) [1]. MR images can be represented by a relatively small number of significant coefficients with nonzero values under an appropriate transform domain. CS takes advantage of the MR image sparsity and reconstructs the image with an under-sampling rate far lower than that required by the Nyquist sampling theorem [2]–[4]. However, in stationary two-dimensional (2D) imaging,
C
Manuscript received October 20, 2013; accepted July 10, 2014. Date of publication July 23, 2014; date of current version December 18, 2014. This work was supported in part by the Australian Research Council, Quarantine of China under Grant 201210079 and Science Technology Department of Zhejiang Province of China under Grant 2013C24019. Asterisk indicates corresponding author. Y. Yang and S. Crozier are with the School of Information Technology and Electrical Engineering, The University of Queensland, Brisbane, QLD 4072, Australia (e-mail: [email protected]; [email protected]). ∗ F. Liu is with the School of Information Technology and Electrical Engineering, The University of Queensland, Brisbane, Qld. 4072, Australia (e-mail: [email protected]). W. Xu is with the School of Information Engineering, China Jiliang University, Hangzhou 310018, China (e-mail: [email protected]). Color versions of one or more of the figures in this paper are available online at http://ieeexplore.ieee.org. Digital Object Identifier 10.1109/TBME.2014.2341621
CS is still limited in providing diagnostically accurate images because it usually results in MR images with large residual artifacts [5]. Over the past few years, various methods have been proposed to improve the reconstruction quality of CS MRI. These methods can be generally categorized in three ways. First, designing irregular k-space sampling, such as variable density k-space sampling and random sampling, is an effective method. It makes the aliasing interference less visually prominent and incoherent in the sparse transform domain, which is an essential aspect for CS [5]–[7]. However, artifacts are difficult to distinguish from real signals at larger reduction factors (the ratio between the amount of subsampled data (D) and fully sampled data (N), D ) [2]. This is because the final reconstructed image i.e., R = N quality not only relies on the limited randomness of the k-space sampling, but it is also largely determined by another major factor, that is the total amount and property of the collected signal [5], [8]–[10]. A second approach attempts to use a range of sparsity bases in spatial and temporal dimensions to provide sufficient sparsity for faithful reconstruction using a subset of the largest transform coefficients. For example, the discrete wavelet transform [11], curvelet transform [12], Walsh transform [13], and singular value decomposition-based transforms [14], [15] are commonly used as sparsity transforms. However, without sufficient k-space data, sparse coefficients transformed from the sampling signals cannot faithfully represent the original object, and thus, the quality of the reconstructed images is degraded. The third category comprises optimization methods for nonlinear signal recovery [11], [16], [17]. These have been a focus of research interest since CS was introduced. These algorithms, such as the focal underdetermined system solution, the iterative reweighted least squares and the conjugate gradients, enforce both the image sparsity and data consistency between the reconstruction and the acquisition. In these methods, a simple regularization scheme is typically used to balance the sparsity of the whole image, relative to the required data consistency [11], [16]. However, one threshold may not satisfy both conditions simultaneously, thus affecting the final reconstruction accuracy, especially with high reduction factors [11], [18], [19]. Typically, the sparsity properties of the low-frequency and highfrequency k-space data can be significantly different [5], [19]. Frequency-specific regularization schemes could be used to resolve this issue. In this paper, we introduce a novel reconstruction scheme to relieve some of the shortcomings of previous methods. The two-part CS framework proposed in [20] treats the whole reconstruction with two algorithms, while our method focuses on
0018-9294 © 2014 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications standards/publications/rights/index.html for more information.
YANG et al.: COMPRESSED SENSING MRI VIA TWO-STAGE RECONSTRUCTION
111
the study of the underlying signal properties of under-sampled k-space data. It consists of a two-stage reconstruction procedure to treat the low-frequency and high-frequency parts progressively. In this algorithm, an image outline will first be reconstructed based on the densely sampled k-space center, and then a full image will be recaptured using the recovered low-frequency k-space data and the sparsely sampled high-frequency data. The feasibility and accuracy of the proposed method will first be theoretically analyzed and then validated through typical cardiac cine, brain, and angiography MRI imaging examples.
part of image] is nonzero (except in very few situations) [5]. The leaves of the trees (the high frequency subbands) are typically the most sparse. During the reconstruction process, the sparsity is constrained by the l1 norm, which conventionally includes both the coarser scale and high-frequency subbands. It is difficult to balance the different sparsities of these parts during optimization [11], [19]. Therefore, it is not an optimal method for sparsifying the images by global sparisty constraint. In addition, the global sparsity constraint in conventional CS methods tends to select the signals with large magnitude in the k-space [19]. The k-space signal magnitudes corresponding to these two parts of the wavelet coefficients also show different features. As shown in Fig. 1(c), the signal magnitudes in the central part of the k-space are much higher than the peripheral part. Thus, the global sparsity constraint leads the CS reconstruction to focus on the central region of k-space data. However, the l2 norm operation in (2) is restricted by this in a conventional CS algorithm. The l2 norm of the entire k-space is minimized to ensure the data fidelity. Although the signal magnitude of the k-space is typically smaller in the peripheral data, if it is not reconstructed faithfully, the error can accumulate over a population much larger than in the central region, and can strongly affect the effectiveness of the l2 norm. To solve this problem, different weights are applied to balance l1 and l2 norms. However, this is a somewhat difficult implementation and the accuracy of reconstruction cannot be guaranteed [23]. Focusing on nonuniform sampling, [19] introduced a method with local sparsity constraints instead of a global one. Here, we will propose a new method based on k-space segmentation.
II. THEORY A. CS MRI CS MRI can be briefly described as follows. Suppose an N × N object MR image, which is known to have a sparse representation in some basisΨ. Its k-space data are collected in a random manner. The CS MRI reconstruction model can be described as the following convex optimization problem [11] Minimize : ||Ψ(m)||0 , subject to ||P ϕ (m) − y||2 < ε
(1)
where || · ||0 counts the number of nonzero components, Ψ denotes the sparse transform, m is the reconstructed image, ϕ the Fourier transform, P the random under-sampling pattern, y consists of the measurements, ε is the tolerance that controls the reconstruction data fidelity and relates to the expected noise level. In this optimization problem, the first l0 norm term promotes the sparsity of the reconstructed image m and the constraint ||P ϕ (m) − y||2 < ε enforces the data fidelity [11]. In CS MRI, the restricted isometry property (RIP) is usually used to analyze the performance and robustness of the CS algorithms. If the RIP condition holds [17], the reconstructed image m will be close to the original image. Equation (1) is a well-known nondeterministic polynomialtime (NP) hard problem [21]. By replacing the l0 norm with the l1 norm to enforce the sparsity, it can be solved efficiently as a convex optimization, as follows: Minimize : ||Ψ(m)||1 , subject to || P ϕ (m) − y||2 < ε |xi |. Here, the l1 norm is defined as || x ||1 =
(2)
i
B. Nonuniform Signal Sparsity It is known that the k-space distribution concentrates in the region of origin while it has low levels at the borders [see Fig. 1(c), which shows the k-space data of a phantom]. This reveals that the signal sparsity varies in different sections of the k-space. The nonuniform sparsity distribution can also be clearly seen in 2-D Fourier transform [see Fig. 1(c)]. Naturally indexed on dyadic cubes, the wavelet has a tree structure under inclusion [22]. The wavelet coefficients are monotonically nonincreasing along the tree branches, outwards from the root. In the example shown in Fig. 1, the tree structure [see Fig. 1(a)] illustrates that a wavelet coefficient is nonzero only if the wavelet coefficient in the coarser scale [the area in the red box in Fig. 1(b), corresponding to the lower frequency
C. k-Space Segmentation As shown in Section II-B, the energy of the k-space is concentrated close to the central region for the majority of MR images. The variable density-sampling pattern generally collects more data near the k-space center (low-frequency parts) and less in the periphery of the k-space (high-frequency parts). One example of a random sampling pattern (256 × 256 pixels) with the ∗ 256 reduction factor R = 256 16395 = 3.99 is illustrated in Fig. 2(a). Its central region (the red square box with 64 × 64 pixels) has a ∗ 64 reduction factor Rcenter = 64 2296 = 1.36, which is much lower than the reduction factor of the whole k-space. Hence, the central region of the k-space can be reconstructed more faithfully with a low reduction factor [24]. In addition, the aforementioned wavelet coefficients in the coarser-scale part can be obtained from the faithfully recovered k-space data in the central region. Using the tree structure of the wavelet coefficients, the nonzero coefficients in the high-frequency subband can be constrained based on the reconstructed coarser-scale data [5]. Because the sparsities of both the coarser-scale and high-frequency subbands in the wavelet domain have been improved, the reconstruction accuracy of the whole image is meliorated. In this section, we will describe how to segment the k-space in accordance with different sampling patterns. Suppose the size of the k-space is N × N, a square box with size n × n is chosen in
112
IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 62, NO. 1, JANUARY 2015
Fig. 1. Relationship between the wavelet domain and k-space domain. (a) The tree structure of the Wavelet. Wavelet coefficients flow from the coarsest scale to the finest as the arrow shows. (b) The Wavelet coefficients of a phantom image. The nonzero coefficients focus on the coarser scale as in the red box. (c) The k-space data are transformed from the Wavelet coefficients by k = ϕΨ −1 w. The signal magnitudes of the central part in the red box [corresponding to the wavelet coefficients in the red box of (b)] are much higher than the peripheral part.
Fig. 2. k-space segmentation. (a) Selected central part (the area in the red square box with n × n pixels) in a random sampling pattern (256 × 256 pixels) when the reduction factor R = 4. (b) Sampling density of the selected central part (n × n pixels) in the sampling pattern as in (a) with n = 2 m, m ε [1, 128]. (c) Under-sampled k-space of a cardiac cine MR image with the sampling pattern as in (a). (d) Power ratio of the selected central part (n × n pixels) in (c) when n = 2m, m ε [1, 128]. (e) Average power ratio in the selected central part of the under-sampled k-space data as in (c).
the k-space center, as shown in the previous example in Fig. 2(a). The sampling density of the selected part can be calculated by D = k/n2 , where k is the number of the sampled points in the square box. With the increase of n as n = 2 m, m ε [1, N/2], the curve of D can be obtained. Fig. 2(b) illustrates the sampling density curve of the random sampling pattern in Fig. 2(a). It is obvious that the sampling density of the selected central part is much higher than that in the whole of the k-space. Not only the sampling density, but also the k-space power of the selection is essential for the faithful reconstruction of the
representative outline. If the central part is too small to obtain a high sampling density, it will contain too little k-space information. Thus, the image outline corresponding to the selected part cannot be representative of the whole image. The determination of the size of the central part is important to the reconstruction of the image outline and also the final CS MRI solution. To find the proper size of a central k-space region, we define the k-space power ratio of the selected part p(n) , n = 2 m, m ε [1, N/2], where the function p(n) by P = p(N) means the sum of the power in the selected square area with size
YANG et al.: COMPRESSED SENSING MRI VIA TWO-STAGE RECONSTRUCTION
113
Fig. 4. Catesian k-space under-sampling patterns used for reconstructions with R = 6. (a) 2-D random sampling pattern. (b) 1-D random sampling pattern.
higher sampling density ensures more accurate information in the central k-space region. A larger power ratio indicates that the k-space data in the selected region contains more information for the whole image. However, with the increase of n (size of the selected area), the sampling density and the power ratio show the opposite trends [see Fig. 2(b) and 2(d)]. Therefore, a proper size for segmenting the k-space should be considered. To balance the two contrary elements of the segmentation and to designate the segmentation, we define the average power ratio (AP) of the selected part as AP = D × P.
(3)
The average power ratio can be thought of as a metric of the k-space power/information density of the selected region. One example is shown in Fig. 2(e) by combining the sampling density (D) data from Fig. 2(b) and the power ratio (P) data from Fig. 2(d). The peak point of the AP curve in Fig. 2(e) at n = 42 is obviously the optimal choice of selected region size. For the consideration of 2-D wavelets, in the sparse base used here, the n value will be rounded to be the power of 2. However, the peak point of AP n = 42 does not satisfy the condition n = 2j , j ∈ N, which is a fundamental property of the 2-D wavelet transformation. Therefore, the proper choice n of the selected central part should satisfy n = 2j , j ∈ N and is closest to the peak point of the average power ratio curve. In the example shown in Fig. 2(e), n1 = 32 and n2 = 64 are the two satisfied choices. Because AP (n1 ) < AP (n2 ), n2 = 64 is the more suitable choice. D. Proposed Two-Stage Reconstruction
Fig. 3.
Flowchart of the proposed two-stage reconstruction algorithm.
n × n. Fig. 2(d) shows the power ratio curve of under-sampled cardiac cine MRI data with the sampling pattern from Fig. 2(a). As Fig. 2(d) shows, it is clear that the power ratio P increases with the extension of the selected area. A proper balance of the sampling density and power ratio is considered for the selection of the central region size. A
We now illustrate the flowchart of two-stage process in Fig. 3. The image reconstruction procedure of the proposed method consists of the following four steps: Step 1: Initialization. Suppose m is a MR image with N × N pixels, p is the under-sampling pattern, y consists of the undersampled measurements with pattern p, ϕ is the Fourier transform, Ψ is the 2-D wavelet transform in advance, and ε is the tolerance in (1). Step 2: k-space segmentation. Using the strategy explained in Section II-C, we segment the k-space into two parts. Step 3: The reconstruction of the image outline (Level 1 reconstruction). According to the segmentation, the central part of
114
IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 62, NO. 1, JANUARY 2015
Fig. 5. Performance comparison between conventional CS and two-stage CS for cardiac cine MR data under different reduction factors: (a) with 2-D random sampling pattern and (b) with 1-D random sampling pattern.
the k-space measurements and sampling pattern can be obtained as yC and pC . According to the truncated-Fourier transform, the central part of the k-space data corresponds to the image outline [1]. The outline of the image mC can be reconstructed by the CS algorithm for solving the following problem: minx ||Ψ(mC )||1 , subject to ||pC ϕ (mC ) − yC ||2 < ε.
(4)
Step 4: k-space combination. Transformed from the reconstructed image outline mc , the recovered k-space central part kC = ϕ(mc ) will replace the responding part yC in the whole k-space measurement y to obtain the new k-space data kR . Step 5: The reconstruction of the whole image (Level 2 reconstruction). Because the central region of the k-space has been re-
Fig. 6. Performance comparison between conventional CS and two-stage CS for brain MR data under different reduction factors: (a) with 2-D random sampling pattern and (b) with 1-D random sampling pattern.
constructed, for the Level 2 reconstruction, we use pR = p + pac (pac represents the central region to be all 1s) as the new sampling pattern. Using mR = ϕ−1 (kR ) as the initial solution, the whole image mR can be reconstructed by the CS as follows: minx ||Ψ(mR )||1 , subject to ||pR ϕ (mR ) − y||2 < ε.
(5)
III. MATERIALS AND METHODS A. Data Acquisition The performance of the proposed method was tested with three typical MR datasets: cardiac cine MR data, brain MR data, and angiography MR data [25]–[27]. The cardiac cine MR data was a frame of full k-space data acquired using a 1.5 T Philips system, as a typical test case. The brain MR image, as a complex form with large contrast in
YANG et al.: COMPRESSED SENSING MRI VIA TWO-STAGE RECONSTRUCTION
115
transform domain [3]. Whilst 2-D random sampling patterns cannot be physically realized in an MRI scanner, it is mainly considered as an ideal case for algorithm evaluation. The 1-D random sampling pattern [see Fig. 4(b)] was implemented only in the phase-encoded direction and by using a variable-density sampling scheme with a denser sampling at a low-frequency. It was easily implemented by randomly changing the amplitudes of the phase-encoding [7]. However, the incoherence of a 1-D random sampling pattern is weaker than a 2-D random sampling pattern, which can cause regular aliasing interference along the phase-encoded direction (data sampled randomly in this direction) in the reconstructed image. We will mainly evaluate the robustness and sensitivity of the proposed algorithm using 1-D random sampling scheme. B. Experimental Methods
Fig. 7. Performance comparison between conventional CS and two-stage CS for angiography MR data under different reduction factors: (a) with 2-D random sampling pattern and (b) with 1-D random sampling pattern.
the pixel domain, was employed to test the algorithm performance for the cases with general MR images. The angiography MR images were chosen to represent MR images already being sparse in the pixel domain. The k-space was fully acquired by a Siemens (Siemens Magnetom TIM Trio, Erlangen, Germany) 1.5 T MRI system (CCAI, 2012). For all these cases, the image size was 256 × 256 pixels. Two different Cartesian under-sampling schemes: 2-D random sampling pattern and 1-D (only in phase direction) random sampling pattern were simulated based on the fully sampled k-space data. For all the cases, the reduction factor R = 6, which means only 16.67% of the k-space data were sampled (see Fig. 4). The 2-D random sampling pattern [see Fig. 4(a)] is guaranteed to have incoherent aliasing interference in the sparse
We utilized nonlinear conjugate gradients and a backtracking line-search [15] for solving the optimization in (4) and (5). For all the reconstruction methods, the initial estimations are directly set as the under-sampled k-space data. With an adjustment to the dataset, the parameter λ, as the weighting of l1 and l2 norms, are experimentally set to be 0.005, 0.0044, and 0.006 for cardiac cine, brain, and angiography. All reconstructions were implemented on a desktop computer with Intel Core i7–3770 eight-core Processor and 12GB of RAM utilizing MATLAB (R2013b: The Mathworks, Natick, MA, USA). The peak signal-to-noise ratio (PSNR) is used as the metric for the quality assessment. PSNR is commonly used to measure the quality of the processed image, which includes the mean square error of the reconstructed image in comparison with the reference image [28]. The PSNR was calculated as ⎛ ⎞ MAX I ⎠ PSNR = 20log10 ⎝ m n 2 1 [I (i, j) − I (i, j)] rec i=1 j =1 mn (6) where I is the benchmark (fully sampled image), MAXI is the maximum possible pixel value of the image, and Irec is the reconstructed image. In addition, the intensity difference between the reconstructed images and the fully sampled images are also used to visually assess the reconstruction quality. IV. RESULTS A. Reconstruction Performance With Different Reduction Factors The PSNR of the reconstructed images by a conventional CS method and the proposed method was recorded with respect to the different reduction factors. The results of cardiac cine, brain, and angiography MR datasets are shown as Figs. 5–7, respectively. The circles in the figure represent the conventional method, while the squares denote the two-stage method. As these figures illustrate, the improvements of the proposed method were very obvious in various cases and sampling patterns. They yielded an enhancement in quality of about 1.5 to 4 dB.
116
IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 62, NO. 1, JANUARY 2015
Fig. 8. Reconstructed images of cardiac cine MR data (top row images) and difference images (bottom row images) under 2-D random sampling pattern at R = 6: (left) with conventional CS (right) with the proposed two-stage CS.
Fig. 9. Reconstructed images of cardiac cine MR data (top row images) and difference images (bottom row images) under 1-D random sampling pattern at R = 6: (left) with conventional CS, (right) with the proposed two-stage CS.
With the comparison of three different cases, the proposed method shows its capability in handling the cases with larger reduction factors. Using the 2-D random sampling pattern, less k-space data in the peripheral part were obtained with a higher reduction factor. The resulting images were blurred both at the edges and in the details. The proposed method can more
Fig. 10. Reconstructed images of brain MR data (top row images) and difference images (bottom row images) under 2-D random sampling pattern at R = 6: (left) with conventional CS, (right) with the proposed two-stage CS.
faithfully reconstruct an outline of the image with a higher sampling ratio and use the tree structure of the 2-D wavelet to control the sparsity of the whole image. The sparsity control provided better reconstruction quality compared with the conventional method. For detailed comparison, see Figs. 5(a), 6(a) and 7(a). As Fig. 4 shows, the 1-D random sampling pattern with Cartesian trajectories is more coherent compared with the 2D random sampling pattern. Comparing the PSNR curves of the conventional CS method under these two random sampling patterns (see Figs. 5–7), we can see the aliasing artifacts accompanied by 1-D random sampling pattern obviously degrade the quality of the reconstructed images. Under the 1-D random sampling pattern, the PSNR curve of a conventional CS method decreases steeply with the increase of reduction factor R. The proposed method uses the higher sampling ratio in the selected region to suppress the aliasing artifacts in the reconstructed image outlines. For the tree structure of a 2-D wavelet, the more faithfully recovered image outline leads to an accurate sparsity control of the whole image. As shown in Figs. 5(b), 6(b) and 7(b), the proposed method showed substantial improvement compared with the conventional ones.
B. Comparison of Image Reconstruction The reconstruction results of the conventional method and the proposed method are presented in Figs. 8–13. Three imaging cases (cardiac cine, brain, and angiography) are compared and R = 6 in each case. The differences between the reconstructed images and the originals are also shown to facilitate comparison.
YANG et al.: COMPRESSED SENSING MRI VIA TWO-STAGE RECONSTRUCTION
117
Fig. 11. Reconstructed images of brain MR data (top row images) and difference images (bottom row images) under 2-D random sampling pattern at R = 6: (left) with conventional CS, (right) with the proposed two-stage CS.
Fig. 13. Reconstructed images of angiography MR data (top row images) and difference images (bottom row images) under 2-D random sampling pattern at R = 6: (left) with conventional CS, (right) with the proposed two-stage CS.
reconstruction of image details than the conventional method. For the examples, note the improvements marked by red arrows in Figs. 8–13. V. DISCUSSION
Fig. 12. Reconstructed images of angiography MR data (top row images) and difference images (bottom row images) under 2-D random sampling pattern at R = 6: (left) with conventional CS, (right) with the proposed two-stage CS.
For better visualization, the difference image in each case is normalized by its individual maximum. The imaging studies indicate that, the proposed method effectively reduces the residual artifacts in 2-D random sampling patterns and aliasing artifacts in 1-D random sampling patterns. In most cases, the proposed method offers higher fidelity in the
Compared with the conventional method, the reconstruction of a two-stage method localizes different parts of the k-space to exploit the redundancy of the sampled data. To balance the reduction factor and reconstruction quality, the two-stage CS method takes advantage of the inhomogeneous distribution of the k-space power and the variable-density of the sampling pattern. The higher sampling ratio of the selected central k-space data leads to a more faithful reconstruction. The reconstructed image outline in stage one can conduce a better initial solution for the whole k-space reconstruction of stage two. The corresponding wavelet coefficients in the coarser scale can approach the sparsity control of the whole image. The solution also avoids most incoherent artifacts caused by the estimation errors at the low-frequencies. For the 2-D random sampling pattern with a high level of incoherence, the faithful reconstruction of the image outline will lead to less residual errors in the whole image reconstruction, because it provides more general features. For the 1-D random sampling pattern, the reconstruction of the conventional CS schemes usually has more aliasing artifacts, because the sampling pattern is only random in one direction. Using the proposed method, the low-frequency k-space information was recovered with a lower reduction factor first that could reduce the aliasing effects. With a more faithful image outline (coarser scale of wavelet coefficient), the whole image
118
IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 62, NO. 1, JANUARY 2015
reconstruction was more accurate and with less aliasing artifacts by the sparsity control of wavelet tree structure. As mentioned, the sparsity of the wavelet coefficients in the coarser scale and the high-frequency subband are different. It is useful to adjust the weights of l1 and l2 in the two-stage optimization process of the CS algorithm. From previous case studies, the results manifested the effects of weighting adjustment in different k-space parts. We also note that, the 2-D wavelet sparsity in the two-stage reconstruction can be replaced by other proper sparsifying transformations [13], [15]. In the implementation of the proposed two-stage reconstruction, the algorithm convergence of each iteration is better owing to smaller problem size (in stage one) and better initial condition (in stage two). Therefore, although the total iteration times (including two stages) of the proposed method are similar to those of the conventional method, the two-stage CS has less total computational time. VI. CONCLUSION To reduce the large reconstruction artifacts that occur when using existing CS schemes with high reduction factors, we have proposed a two-stage CS method that recovers the central and peripheral k-space data in a sequential manner. From the relatively dense sampled inner k-space part, an image outline is first reconstructed. Using image outline from stage one and the sampled peripheral k-space data, a full image recovery is then implemented. This stepped reconstruction naturally solves a number of technical difficulties in conventional CS methods, including the processing of a large magnitude contrast between the low-frequency and high-frequency k-space data. The new method has been applied to cardiac cine, brain, and angiography MRI datasets and compared with conventional CS schemes. We have demonstrated that the new method improves the reconstruction accuracy and reduces the reconstruction time. Because the presence of aliasing artifacts also limits the application of the CS in dynamic MRI, in our future work, we plan to extend the developed CS algorithm to dynamic MRI studies. REFERENCES [1] A. Haase, “Snapshot flash mri. applications to t1, t2, and chemical-shift imaging,” Magn. Reson. Med., vol. 13, pp. 77–89, 1990. [2] M. Stehling, R. Turner, and P. Mansfield, “Echo-planar imaging: Magnetic resonance imaging in a fraction of a second,” Science, vol. 254, pp. 43–50, 1991. [3] M. Lustig, D. Donoho, J. M. Santos, and J. M. Pauly, “A look at how CS can improve on current imaging techniques,” IEEE Signal Process. Mag., vol. 25, no. 2, pp. 72–82, Mar. 2008. [4] E. Cand`es, J. Romberg, and T. Tao, “Robust uncertainty principles: Exact signal reconstruction from highly incomplete frequency information,” IEEE Trans. Inf. Theory, vol. 52, no. 2, pp. 489–509, Feb. 2006. [5] K. Sung and B. A. Hargreaves. (2012). High-frequency subband compressed sensing MRI using quadruplet sampling. Magn. Reson. Med., vol. 70, no. 5, pp. 1306–1318. [Online]. Available: http://onlinelibrary. wiley.com/doi/10.1002/mrm.24592/abstract;jsessionid = 403CC7AB962031BE9C94F0D97A960083.f02t01 [6] E. Cand`es and T. Tao, “Near optimal signal recovery from random projections: Universal encoding strategies,” IEEE Trans. Inf. Theory, vol. 52, no. 12, pp. 5406–5425, Dec. 2006.
[7] H. Wang, D. Liang, and L. Ying, “Pseudo 2D random sampling for compressed sensing MRI,” in Proc. IEEE Eng. Med. Biol. Soc. Annu. Int. Conf., 2009, pp. 2672–2675. [8] P. M. Margosian, G. DeMeester, and H. Liu, “Partial fourier acquisition in MRI,” in Encyclopaedia of Magnetic Resonance. New York, NY, USA: Wiley, 2007. [9] N. Chauffert, P. Ciuciu, J. Kahn, and P. Weiss, “Variable density sampling with continuous sampling trajectories,” arXiv:1311.6039, 2013. [10] R. W. Chan, E. A. Ramsay, E. Y. Cheung, and D. B. Plewes, “The influence of radial undersampling schemes on compressed sensing reconstruction in breast MRI,” Magn. Reson. Med., vol. 67, pp. 363–377, 2012. [11] M. Lustig, D. Donoho, and J. M. Pauly, “Sparse MRI: The application of compressed sensing for rapid MR imaging,” Magn. Reson. Med., vol. 58, pp. 1182-1195, 2007. [12] E. Cande‘s and D. Donoho, Curvelets—A Surprisingly Effectivenonadaptive Representation for Objects With Edges, in A. Cohen, C. Rabut, and L. Schumaker, Eds. Nashville, TN, USA: Vanderbilt Univ. Press, 2000, pp. 105–120. [13] Z. Feng, F. Liu, M. Jiang, S. Crozier, H. Guo. and Y. Wang, “Improved l1-SPIRiT using 3D Walsh transform-based sparsity basis,” Magn. Reson. Imag., vol. 32, no. 7, pp. 924–933, 2014. [14] M. Hong, Y. Yu, H. Wang, F. Liu, and S. Crozier, “Compressed sensing MRI with singular value decomposition based sparsity basis,” Phys. Med. Biol., vol. 56, no. 19, pp. 6311–6325, 2011. [15] Y. Yu, J. Jin, F. Liu, and S. Crozier, “Multidimensional compressed sensing MRI using tensor decomposition-based sparsifying transform,” PLoS One, vol. 9, no. 6, p. e98441, Jun. 5, 2014. [16] D. S. Taubman and M. W. Marcellin, JPEG 2000: Image Compression Fundamentals, Standards and Practice (International Series in Engineering and Computer Science). Norwell, MA, USA: Kluwer, 2002. [17] D. Donoho, “Compressed sensing,” IEEE Trans. Inf. Theory, vol. 52, no. 4, pp. 1289–1306, Apr. 2006. [18] R. Chartrand, “Fast algorithms for nonconvex compressive sensing: MRI reconstruction from very few data,” in Proc. IEEE Int. Symp. Biomed. Imag.: From Nano to Macro, Jun. 28/Jul. 1, 2009, pp. 262–265. [19] F. A. Razzaq, S. Mohamed, A. Bhatti, and S. Nahavandi, “Non-uniform sparsity in rapid compressive sensing MRI,” in Proc. IEEE Int. Conf. Syst. Man, Cybern., 2012, pp. 2253–2258. [20] Y. Ma, D. Baron, and D. Needell, “Two-part reconstruction in compressed sensing,” in Proc. IEEE Global Conf. Signal Inf. Process., 2013, pp. 1041–1044. [21] E. J. Cand`es, J. K. Romberg, and T. Tao, “Stable signal recovery from incomplete and inaccurate measurements,” Commun. Pure Appl. Math., vol. 59, pp. 1207–1223, 2006. [22] R. G. Baraniuk, R. A. DeVore, G. Kyriazis, and X. M. Yu, “Near best tree approximation,” Adv. Comput. Math., vol. 16, pp. 357–373, 2012. [23] H. Oh and S. Lee, “Visually weighted reconstruction of compressive sensing MRI,” Magn. Reson. Imag., vol. 32, pp. 270–280, 2014. [24] E. J. Cand`es, “Compressive sampling,” presented at the Int. Congr. Mathematicians, Madrid, Spain, 2006. [25] K. T. Block, M. Uecker, and J. Frahm, “Undersampled radial MRI with multiple coils. Iterative image reconstruction using a total variation constraint,” Magn. Reson. Med., vol. 57, no. 6, pp. 1086–1098, 2007. [26] 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,” Magn. Reson. Med., vol. 61, no. 1, pp. 103–116, 2009. [27] D. Liang, E. V. DiBella, R. R. Chen, and L. Ying, “k-t ISD: dynamic cardiac MR imaging using compressed sensing with iterative support detection,” Magn. Reson. Med., vol. 68, no. 1, pp. 41–53, 2012. [28] S. G. Lingala, H. Yue, E. DiBella, and M. Jacob, “Accelerated dynamic MRI exploiting sparsity and low-rank structure: k-t SLR,” IEEE Trans. Med. Imag., vol. 30, no. 5, pp. 1042–1054, May 2011.
Authors’ photographs and biographies not available at the time of publication.
IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 62, NO. 1, JANUARY 2015
Compressed Sensing MRI via Two-stage Reconstruction Yang Yang, Feng Liu∗ , Wenlong Xu, and Stuart Crozier, Member, IEEE
Abstract—Compressed sensing (CS) has been applied to magnetic resonance imaging for the acceleration of data collection. However, existing CS techniques usually produce images with residual artifacts, particularly at high reduction factors. In this paper, we propose a novel, two-stage reconstruction scheme, which takes advantage of the properties of k-space data and undersampling patterns that are useful in CS. In this algorithm, the under-sampled k-space data is segmented into low-frequency and high-frequency domains. Then, in stage one, using dense measurements, the low-frequency region of k-space data is faithfully reconstructed. The fully reconstituted low-frequency k-space data from the first stage is then combined with the high-frequency k-space data to complete the second stage reconstruction of the whole of k-space. With this two-stage approach, each reconstruction inherently incorporates a lower data under-sampling rate and more homogeneous signal magnitudes than conventional approaches. Because the restricted isometric property is easier to satisfy, the reconstruction consequently produces lower residual errors at each step. Compared with a conventional CS reconstruction, for the cases of cardiac cine, brain and angiogram imaging, the proposed method achieves a more accurate reconstruction with an improvement of 2–4 dB in peak signal-to-noise ratio respectively, using reduction factors of up to 6. Index Terms—Compressed sensing (CS), K-space segmentation, stationary magnetic resonance imaging (MRI), two-stage reconstruction.
I. INTRODUCTION OMPRESSED sensing (CS), a relatively new signal processing theory, has been successfully applied to accelerate the scan speed of magnetic resonance imaging (MRI) [1]. MR images can be represented by a relatively small number of significant coefficients with nonzero values under an appropriate transform domain. CS takes advantage of the MR image sparsity and reconstructs the image with an under-sampling rate far lower than that required by the Nyquist sampling theorem [2]–[4]. However, in stationary two-dimensional (2D) imaging,
C
Manuscript received October 20, 2013; accepted July 10, 2014. Date of publication July 23, 2014; date of current version December 18, 2014. This work was supported in part by the Australian Research Council, Quarantine of China under Grant 201210079 and Science Technology Department of Zhejiang Province of China under Grant 2013C24019. Asterisk indicates corresponding author. Y. Yang and S. Crozier are with the School of Information Technology and Electrical Engineering, The University of Queensland, Brisbane, QLD 4072, Australia (e-mail: [email protected]; [email protected]). ∗ F. Liu is with the School of Information Technology and Electrical Engineering, The University of Queensland, Brisbane, Qld. 4072, Australia (e-mail: [email protected]). W. Xu is with the School of Information Engineering, China Jiliang University, Hangzhou 310018, China (e-mail: [email protected]). Color versions of one or more of the figures in this paper are available online at http://ieeexplore.ieee.org. Digital Object Identifier 10.1109/TBME.2014.2341621
CS is still limited in providing diagnostically accurate images because it usually results in MR images with large residual artifacts [5]. Over the past few years, various methods have been proposed to improve the reconstruction quality of CS MRI. These methods can be generally categorized in three ways. First, designing irregular k-space sampling, such as variable density k-space sampling and random sampling, is an effective method. It makes the aliasing interference less visually prominent and incoherent in the sparse transform domain, which is an essential aspect for CS [5]–[7]. However, artifacts are difficult to distinguish from real signals at larger reduction factors (the ratio between the amount of subsampled data (D) and fully sampled data (N), D ) [2]. This is because the final reconstructed image i.e., R = N quality not only relies on the limited randomness of the k-space sampling, but it is also largely determined by another major factor, that is the total amount and property of the collected signal [5], [8]–[10]. A second approach attempts to use a range of sparsity bases in spatial and temporal dimensions to provide sufficient sparsity for faithful reconstruction using a subset of the largest transform coefficients. For example, the discrete wavelet transform [11], curvelet transform [12], Walsh transform [13], and singular value decomposition-based transforms [14], [15] are commonly used as sparsity transforms. However, without sufficient k-space data, sparse coefficients transformed from the sampling signals cannot faithfully represent the original object, and thus, the quality of the reconstructed images is degraded. The third category comprises optimization methods for nonlinear signal recovery [11], [16], [17]. These have been a focus of research interest since CS was introduced. These algorithms, such as the focal underdetermined system solution, the iterative reweighted least squares and the conjugate gradients, enforce both the image sparsity and data consistency between the reconstruction and the acquisition. In these methods, a simple regularization scheme is typically used to balance the sparsity of the whole image, relative to the required data consistency [11], [16]. However, one threshold may not satisfy both conditions simultaneously, thus affecting the final reconstruction accuracy, especially with high reduction factors [11], [18], [19]. Typically, the sparsity properties of the low-frequency and highfrequency k-space data can be significantly different [5], [19]. Frequency-specific regularization schemes could be used to resolve this issue. In this paper, we introduce a novel reconstruction scheme to relieve some of the shortcomings of previous methods. The two-part CS framework proposed in [20] treats the whole reconstruction with two algorithms, while our method focuses on
0018-9294 © 2014 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications standards/publications/rights/index.html for more information.
YANG et al.: COMPRESSED SENSING MRI VIA TWO-STAGE RECONSTRUCTION
111
the study of the underlying signal properties of under-sampled k-space data. It consists of a two-stage reconstruction procedure to treat the low-frequency and high-frequency parts progressively. In this algorithm, an image outline will first be reconstructed based on the densely sampled k-space center, and then a full image will be recaptured using the recovered low-frequency k-space data and the sparsely sampled high-frequency data. The feasibility and accuracy of the proposed method will first be theoretically analyzed and then validated through typical cardiac cine, brain, and angiography MRI imaging examples.
part of image] is nonzero (except in very few situations) [5]. The leaves of the trees (the high frequency subbands) are typically the most sparse. During the reconstruction process, the sparsity is constrained by the l1 norm, which conventionally includes both the coarser scale and high-frequency subbands. It is difficult to balance the different sparsities of these parts during optimization [11], [19]. Therefore, it is not an optimal method for sparsifying the images by global sparisty constraint. In addition, the global sparsity constraint in conventional CS methods tends to select the signals with large magnitude in the k-space [19]. The k-space signal magnitudes corresponding to these two parts of the wavelet coefficients also show different features. As shown in Fig. 1(c), the signal magnitudes in the central part of the k-space are much higher than the peripheral part. Thus, the global sparsity constraint leads the CS reconstruction to focus on the central region of k-space data. However, the l2 norm operation in (2) is restricted by this in a conventional CS algorithm. The l2 norm of the entire k-space is minimized to ensure the data fidelity. Although the signal magnitude of the k-space is typically smaller in the peripheral data, if it is not reconstructed faithfully, the error can accumulate over a population much larger than in the central region, and can strongly affect the effectiveness of the l2 norm. To solve this problem, different weights are applied to balance l1 and l2 norms. However, this is a somewhat difficult implementation and the accuracy of reconstruction cannot be guaranteed [23]. Focusing on nonuniform sampling, [19] introduced a method with local sparsity constraints instead of a global one. Here, we will propose a new method based on k-space segmentation.
II. THEORY A. CS MRI CS MRI can be briefly described as follows. Suppose an N × N object MR image, which is known to have a sparse representation in some basisΨ. Its k-space data are collected in a random manner. The CS MRI reconstruction model can be described as the following convex optimization problem [11] Minimize : ||Ψ(m)||0 , subject to ||P ϕ (m) − y||2 < ε
(1)
where || · ||0 counts the number of nonzero components, Ψ denotes the sparse transform, m is the reconstructed image, ϕ the Fourier transform, P the random under-sampling pattern, y consists of the measurements, ε is the tolerance that controls the reconstruction data fidelity and relates to the expected noise level. In this optimization problem, the first l0 norm term promotes the sparsity of the reconstructed image m and the constraint ||P ϕ (m) − y||2 < ε enforces the data fidelity [11]. In CS MRI, the restricted isometry property (RIP) is usually used to analyze the performance and robustness of the CS algorithms. If the RIP condition holds [17], the reconstructed image m will be close to the original image. Equation (1) is a well-known nondeterministic polynomialtime (NP) hard problem [21]. By replacing the l0 norm with the l1 norm to enforce the sparsity, it can be solved efficiently as a convex optimization, as follows: Minimize : ||Ψ(m)||1 , subject to || P ϕ (m) − y||2 < ε |xi |. Here, the l1 norm is defined as || x ||1 =
(2)
i
B. Nonuniform Signal Sparsity It is known that the k-space distribution concentrates in the region of origin while it has low levels at the borders [see Fig. 1(c), which shows the k-space data of a phantom]. This reveals that the signal sparsity varies in different sections of the k-space. The nonuniform sparsity distribution can also be clearly seen in 2-D Fourier transform [see Fig. 1(c)]. Naturally indexed on dyadic cubes, the wavelet has a tree structure under inclusion [22]. The wavelet coefficients are monotonically nonincreasing along the tree branches, outwards from the root. In the example shown in Fig. 1, the tree structure [see Fig. 1(a)] illustrates that a wavelet coefficient is nonzero only if the wavelet coefficient in the coarser scale [the area in the red box in Fig. 1(b), corresponding to the lower frequency
C. k-Space Segmentation As shown in Section II-B, the energy of the k-space is concentrated close to the central region for the majority of MR images. The variable density-sampling pattern generally collects more data near the k-space center (low-frequency parts) and less in the periphery of the k-space (high-frequency parts). One example of a random sampling pattern (256 × 256 pixels) with the ∗ 256 reduction factor R = 256 16395 = 3.99 is illustrated in Fig. 2(a). Its central region (the red square box with 64 × 64 pixels) has a ∗ 64 reduction factor Rcenter = 64 2296 = 1.36, which is much lower than the reduction factor of the whole k-space. Hence, the central region of the k-space can be reconstructed more faithfully with a low reduction factor [24]. In addition, the aforementioned wavelet coefficients in the coarser-scale part can be obtained from the faithfully recovered k-space data in the central region. Using the tree structure of the wavelet coefficients, the nonzero coefficients in the high-frequency subband can be constrained based on the reconstructed coarser-scale data [5]. Because the sparsities of both the coarser-scale and high-frequency subbands in the wavelet domain have been improved, the reconstruction accuracy of the whole image is meliorated. In this section, we will describe how to segment the k-space in accordance with different sampling patterns. Suppose the size of the k-space is N × N, a square box with size n × n is chosen in
112
IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 62, NO. 1, JANUARY 2015
Fig. 1. Relationship between the wavelet domain and k-space domain. (a) The tree structure of the Wavelet. Wavelet coefficients flow from the coarsest scale to the finest as the arrow shows. (b) The Wavelet coefficients of a phantom image. The nonzero coefficients focus on the coarser scale as in the red box. (c) The k-space data are transformed from the Wavelet coefficients by k = ϕΨ −1 w. The signal magnitudes of the central part in the red box [corresponding to the wavelet coefficients in the red box of (b)] are much higher than the peripheral part.
Fig. 2. k-space segmentation. (a) Selected central part (the area in the red square box with n × n pixels) in a random sampling pattern (256 × 256 pixels) when the reduction factor R = 4. (b) Sampling density of the selected central part (n × n pixels) in the sampling pattern as in (a) with n = 2 m, m ε [1, 128]. (c) Under-sampled k-space of a cardiac cine MR image with the sampling pattern as in (a). (d) Power ratio of the selected central part (n × n pixels) in (c) when n = 2m, m ε [1, 128]. (e) Average power ratio in the selected central part of the under-sampled k-space data as in (c).
the k-space center, as shown in the previous example in Fig. 2(a). The sampling density of the selected part can be calculated by D = k/n2 , where k is the number of the sampled points in the square box. With the increase of n as n = 2 m, m ε [1, N/2], the curve of D can be obtained. Fig. 2(b) illustrates the sampling density curve of the random sampling pattern in Fig. 2(a). It is obvious that the sampling density of the selected central part is much higher than that in the whole of the k-space. Not only the sampling density, but also the k-space power of the selection is essential for the faithful reconstruction of the
representative outline. If the central part is too small to obtain a high sampling density, it will contain too little k-space information. Thus, the image outline corresponding to the selected part cannot be representative of the whole image. The determination of the size of the central part is important to the reconstruction of the image outline and also the final CS MRI solution. To find the proper size of a central k-space region, we define the k-space power ratio of the selected part p(n) , n = 2 m, m ε [1, N/2], where the function p(n) by P = p(N) means the sum of the power in the selected square area with size
YANG et al.: COMPRESSED SENSING MRI VIA TWO-STAGE RECONSTRUCTION
113
Fig. 4. Catesian k-space under-sampling patterns used for reconstructions with R = 6. (a) 2-D random sampling pattern. (b) 1-D random sampling pattern.
higher sampling density ensures more accurate information in the central k-space region. A larger power ratio indicates that the k-space data in the selected region contains more information for the whole image. However, with the increase of n (size of the selected area), the sampling density and the power ratio show the opposite trends [see Fig. 2(b) and 2(d)]. Therefore, a proper size for segmenting the k-space should be considered. To balance the two contrary elements of the segmentation and to designate the segmentation, we define the average power ratio (AP) of the selected part as AP = D × P.
(3)
The average power ratio can be thought of as a metric of the k-space power/information density of the selected region. One example is shown in Fig. 2(e) by combining the sampling density (D) data from Fig. 2(b) and the power ratio (P) data from Fig. 2(d). The peak point of the AP curve in Fig. 2(e) at n = 42 is obviously the optimal choice of selected region size. For the consideration of 2-D wavelets, in the sparse base used here, the n value will be rounded to be the power of 2. However, the peak point of AP n = 42 does not satisfy the condition n = 2j , j ∈ N, which is a fundamental property of the 2-D wavelet transformation. Therefore, the proper choice n of the selected central part should satisfy n = 2j , j ∈ N and is closest to the peak point of the average power ratio curve. In the example shown in Fig. 2(e), n1 = 32 and n2 = 64 are the two satisfied choices. Because AP (n1 ) < AP (n2 ), n2 = 64 is the more suitable choice. D. Proposed Two-Stage Reconstruction
Fig. 3.
Flowchart of the proposed two-stage reconstruction algorithm.
n × n. Fig. 2(d) shows the power ratio curve of under-sampled cardiac cine MRI data with the sampling pattern from Fig. 2(a). As Fig. 2(d) shows, it is clear that the power ratio P increases with the extension of the selected area. A proper balance of the sampling density and power ratio is considered for the selection of the central region size. A
We now illustrate the flowchart of two-stage process in Fig. 3. The image reconstruction procedure of the proposed method consists of the following four steps: Step 1: Initialization. Suppose m is a MR image with N × N pixels, p is the under-sampling pattern, y consists of the undersampled measurements with pattern p, ϕ is the Fourier transform, Ψ is the 2-D wavelet transform in advance, and ε is the tolerance in (1). Step 2: k-space segmentation. Using the strategy explained in Section II-C, we segment the k-space into two parts. Step 3: The reconstruction of the image outline (Level 1 reconstruction). According to the segmentation, the central part of
114
IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 62, NO. 1, JANUARY 2015
Fig. 5. Performance comparison between conventional CS and two-stage CS for cardiac cine MR data under different reduction factors: (a) with 2-D random sampling pattern and (b) with 1-D random sampling pattern.
the k-space measurements and sampling pattern can be obtained as yC and pC . According to the truncated-Fourier transform, the central part of the k-space data corresponds to the image outline [1]. The outline of the image mC can be reconstructed by the CS algorithm for solving the following problem: minx ||Ψ(mC )||1 , subject to ||pC ϕ (mC ) − yC ||2 < ε.
(4)
Step 4: k-space combination. Transformed from the reconstructed image outline mc , the recovered k-space central part kC = ϕ(mc ) will replace the responding part yC in the whole k-space measurement y to obtain the new k-space data kR . Step 5: The reconstruction of the whole image (Level 2 reconstruction). Because the central region of the k-space has been re-
Fig. 6. Performance comparison between conventional CS and two-stage CS for brain MR data under different reduction factors: (a) with 2-D random sampling pattern and (b) with 1-D random sampling pattern.
constructed, for the Level 2 reconstruction, we use pR = p + pac (pac represents the central region to be all 1s) as the new sampling pattern. Using mR = ϕ−1 (kR ) as the initial solution, the whole image mR can be reconstructed by the CS as follows: minx ||Ψ(mR )||1 , subject to ||pR ϕ (mR ) − y||2 < ε.
(5)
III. MATERIALS AND METHODS A. Data Acquisition The performance of the proposed method was tested with three typical MR datasets: cardiac cine MR data, brain MR data, and angiography MR data [25]–[27]. The cardiac cine MR data was a frame of full k-space data acquired using a 1.5 T Philips system, as a typical test case. The brain MR image, as a complex form with large contrast in
YANG et al.: COMPRESSED SENSING MRI VIA TWO-STAGE RECONSTRUCTION
115
transform domain [3]. Whilst 2-D random sampling patterns cannot be physically realized in an MRI scanner, it is mainly considered as an ideal case for algorithm evaluation. The 1-D random sampling pattern [see Fig. 4(b)] was implemented only in the phase-encoded direction and by using a variable-density sampling scheme with a denser sampling at a low-frequency. It was easily implemented by randomly changing the amplitudes of the phase-encoding [7]. However, the incoherence of a 1-D random sampling pattern is weaker than a 2-D random sampling pattern, which can cause regular aliasing interference along the phase-encoded direction (data sampled randomly in this direction) in the reconstructed image. We will mainly evaluate the robustness and sensitivity of the proposed algorithm using 1-D random sampling scheme. B. Experimental Methods
Fig. 7. Performance comparison between conventional CS and two-stage CS for angiography MR data under different reduction factors: (a) with 2-D random sampling pattern and (b) with 1-D random sampling pattern.
the pixel domain, was employed to test the algorithm performance for the cases with general MR images. The angiography MR images were chosen to represent MR images already being sparse in the pixel domain. The k-space was fully acquired by a Siemens (Siemens Magnetom TIM Trio, Erlangen, Germany) 1.5 T MRI system (CCAI, 2012). For all these cases, the image size was 256 × 256 pixels. Two different Cartesian under-sampling schemes: 2-D random sampling pattern and 1-D (only in phase direction) random sampling pattern were simulated based on the fully sampled k-space data. For all the cases, the reduction factor R = 6, which means only 16.67% of the k-space data were sampled (see Fig. 4). The 2-D random sampling pattern [see Fig. 4(a)] is guaranteed to have incoherent aliasing interference in the sparse
We utilized nonlinear conjugate gradients and a backtracking line-search [15] for solving the optimization in (4) and (5). For all the reconstruction methods, the initial estimations are directly set as the under-sampled k-space data. With an adjustment to the dataset, the parameter λ, as the weighting of l1 and l2 norms, are experimentally set to be 0.005, 0.0044, and 0.006 for cardiac cine, brain, and angiography. All reconstructions were implemented on a desktop computer with Intel Core i7–3770 eight-core Processor and 12GB of RAM utilizing MATLAB (R2013b: The Mathworks, Natick, MA, USA). The peak signal-to-noise ratio (PSNR) is used as the metric for the quality assessment. PSNR is commonly used to measure the quality of the processed image, which includes the mean square error of the reconstructed image in comparison with the reference image [28]. The PSNR was calculated as ⎛ ⎞ MAX I ⎠ PSNR = 20log10 ⎝ m n 2 1 [I (i, j) − I (i, j)] rec i=1 j =1 mn (6) where I is the benchmark (fully sampled image), MAXI is the maximum possible pixel value of the image, and Irec is the reconstructed image. In addition, the intensity difference between the reconstructed images and the fully sampled images are also used to visually assess the reconstruction quality. IV. RESULTS A. Reconstruction Performance With Different Reduction Factors The PSNR of the reconstructed images by a conventional CS method and the proposed method was recorded with respect to the different reduction factors. The results of cardiac cine, brain, and angiography MR datasets are shown as Figs. 5–7, respectively. The circles in the figure represent the conventional method, while the squares denote the two-stage method. As these figures illustrate, the improvements of the proposed method were very obvious in various cases and sampling patterns. They yielded an enhancement in quality of about 1.5 to 4 dB.
116
IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 62, NO. 1, JANUARY 2015
Fig. 8. Reconstructed images of cardiac cine MR data (top row images) and difference images (bottom row images) under 2-D random sampling pattern at R = 6: (left) with conventional CS (right) with the proposed two-stage CS.
Fig. 9. Reconstructed images of cardiac cine MR data (top row images) and difference images (bottom row images) under 1-D random sampling pattern at R = 6: (left) with conventional CS, (right) with the proposed two-stage CS.
With the comparison of three different cases, the proposed method shows its capability in handling the cases with larger reduction factors. Using the 2-D random sampling pattern, less k-space data in the peripheral part were obtained with a higher reduction factor. The resulting images were blurred both at the edges and in the details. The proposed method can more
Fig. 10. Reconstructed images of brain MR data (top row images) and difference images (bottom row images) under 2-D random sampling pattern at R = 6: (left) with conventional CS, (right) with the proposed two-stage CS.
faithfully reconstruct an outline of the image with a higher sampling ratio and use the tree structure of the 2-D wavelet to control the sparsity of the whole image. The sparsity control provided better reconstruction quality compared with the conventional method. For detailed comparison, see Figs. 5(a), 6(a) and 7(a). As Fig. 4 shows, the 1-D random sampling pattern with Cartesian trajectories is more coherent compared with the 2D random sampling pattern. Comparing the PSNR curves of the conventional CS method under these two random sampling patterns (see Figs. 5–7), we can see the aliasing artifacts accompanied by 1-D random sampling pattern obviously degrade the quality of the reconstructed images. Under the 1-D random sampling pattern, the PSNR curve of a conventional CS method decreases steeply with the increase of reduction factor R. The proposed method uses the higher sampling ratio in the selected region to suppress the aliasing artifacts in the reconstructed image outlines. For the tree structure of a 2-D wavelet, the more faithfully recovered image outline leads to an accurate sparsity control of the whole image. As shown in Figs. 5(b), 6(b) and 7(b), the proposed method showed substantial improvement compared with the conventional ones.
B. Comparison of Image Reconstruction The reconstruction results of the conventional method and the proposed method are presented in Figs. 8–13. Three imaging cases (cardiac cine, brain, and angiography) are compared and R = 6 in each case. The differences between the reconstructed images and the originals are also shown to facilitate comparison.
YANG et al.: COMPRESSED SENSING MRI VIA TWO-STAGE RECONSTRUCTION
117
Fig. 11. Reconstructed images of brain MR data (top row images) and difference images (bottom row images) under 2-D random sampling pattern at R = 6: (left) with conventional CS, (right) with the proposed two-stage CS.
Fig. 13. Reconstructed images of angiography MR data (top row images) and difference images (bottom row images) under 2-D random sampling pattern at R = 6: (left) with conventional CS, (right) with the proposed two-stage CS.
reconstruction of image details than the conventional method. For the examples, note the improvements marked by red arrows in Figs. 8–13. V. DISCUSSION
Fig. 12. Reconstructed images of angiography MR data (top row images) and difference images (bottom row images) under 2-D random sampling pattern at R = 6: (left) with conventional CS, (right) with the proposed two-stage CS.
For better visualization, the difference image in each case is normalized by its individual maximum. The imaging studies indicate that, the proposed method effectively reduces the residual artifacts in 2-D random sampling patterns and aliasing artifacts in 1-D random sampling patterns. In most cases, the proposed method offers higher fidelity in the
Compared with the conventional method, the reconstruction of a two-stage method localizes different parts of the k-space to exploit the redundancy of the sampled data. To balance the reduction factor and reconstruction quality, the two-stage CS method takes advantage of the inhomogeneous distribution of the k-space power and the variable-density of the sampling pattern. The higher sampling ratio of the selected central k-space data leads to a more faithful reconstruction. The reconstructed image outline in stage one can conduce a better initial solution for the whole k-space reconstruction of stage two. The corresponding wavelet coefficients in the coarser scale can approach the sparsity control of the whole image. The solution also avoids most incoherent artifacts caused by the estimation errors at the low-frequencies. For the 2-D random sampling pattern with a high level of incoherence, the faithful reconstruction of the image outline will lead to less residual errors in the whole image reconstruction, because it provides more general features. For the 1-D random sampling pattern, the reconstruction of the conventional CS schemes usually has more aliasing artifacts, because the sampling pattern is only random in one direction. Using the proposed method, the low-frequency k-space information was recovered with a lower reduction factor first that could reduce the aliasing effects. With a more faithful image outline (coarser scale of wavelet coefficient), the whole image
118
IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 62, NO. 1, JANUARY 2015
reconstruction was more accurate and with less aliasing artifacts by the sparsity control of wavelet tree structure. As mentioned, the sparsity of the wavelet coefficients in the coarser scale and the high-frequency subband are different. It is useful to adjust the weights of l1 and l2 in the two-stage optimization process of the CS algorithm. From previous case studies, the results manifested the effects of weighting adjustment in different k-space parts. We also note that, the 2-D wavelet sparsity in the two-stage reconstruction can be replaced by other proper sparsifying transformations [13], [15]. In the implementation of the proposed two-stage reconstruction, the algorithm convergence of each iteration is better owing to smaller problem size (in stage one) and better initial condition (in stage two). Therefore, although the total iteration times (including two stages) of the proposed method are similar to those of the conventional method, the two-stage CS has less total computational time. VI. CONCLUSION To reduce the large reconstruction artifacts that occur when using existing CS schemes with high reduction factors, we have proposed a two-stage CS method that recovers the central and peripheral k-space data in a sequential manner. From the relatively dense sampled inner k-space part, an image outline is first reconstructed. Using image outline from stage one and the sampled peripheral k-space data, a full image recovery is then implemented. This stepped reconstruction naturally solves a number of technical difficulties in conventional CS methods, including the processing of a large magnitude contrast between the low-frequency and high-frequency k-space data. The new method has been applied to cardiac cine, brain, and angiography MRI datasets and compared with conventional CS schemes. We have demonstrated that the new method improves the reconstruction accuracy and reduces the reconstruction time. Because the presence of aliasing artifacts also limits the application of the CS in dynamic MRI, in our future work, we plan to extend the developed CS algorithm to dynamic MRI studies. REFERENCES [1] A. Haase, “Snapshot flash mri. applications to t1, t2, and chemical-shift imaging,” Magn. Reson. Med., vol. 13, pp. 77–89, 1990. [2] M. Stehling, R. Turner, and P. Mansfield, “Echo-planar imaging: Magnetic resonance imaging in a fraction of a second,” Science, vol. 254, pp. 43–50, 1991. [3] M. Lustig, D. Donoho, J. M. Santos, and J. M. Pauly, “A look at how CS can improve on current imaging techniques,” IEEE Signal Process. Mag., vol. 25, no. 2, pp. 72–82, Mar. 2008. [4] E. Cand`es, J. Romberg, and T. Tao, “Robust uncertainty principles: Exact signal reconstruction from highly incomplete frequency information,” IEEE Trans. Inf. Theory, vol. 52, no. 2, pp. 489–509, Feb. 2006. [5] K. Sung and B. A. Hargreaves. (2012). High-frequency subband compressed sensing MRI using quadruplet sampling. Magn. Reson. Med., vol. 70, no. 5, pp. 1306–1318. [Online]. Available: http://onlinelibrary. wiley.com/doi/10.1002/mrm.24592/abstract;jsessionid = 403CC7AB962031BE9C94F0D97A960083.f02t01 [6] E. Cand`es and T. Tao, “Near optimal signal recovery from random projections: Universal encoding strategies,” IEEE Trans. Inf. Theory, vol. 52, no. 12, pp. 5406–5425, Dec. 2006.
[7] H. Wang, D. Liang, and L. Ying, “Pseudo 2D random sampling for compressed sensing MRI,” in Proc. IEEE Eng. Med. Biol. Soc. Annu. Int. Conf., 2009, pp. 2672–2675. [8] P. M. Margosian, G. DeMeester, and H. Liu, “Partial fourier acquisition in MRI,” in Encyclopaedia of Magnetic Resonance. New York, NY, USA: Wiley, 2007. [9] N. Chauffert, P. Ciuciu, J. Kahn, and P. Weiss, “Variable density sampling with continuous sampling trajectories,” arXiv:1311.6039, 2013. [10] R. W. Chan, E. A. Ramsay, E. Y. Cheung, and D. B. Plewes, “The influence of radial undersampling schemes on compressed sensing reconstruction in breast MRI,” Magn. Reson. Med., vol. 67, pp. 363–377, 2012. [11] M. Lustig, D. Donoho, and J. M. Pauly, “Sparse MRI: The application of compressed sensing for rapid MR imaging,” Magn. Reson. Med., vol. 58, pp. 1182-1195, 2007. [12] E. Cande‘s and D. Donoho, Curvelets—A Surprisingly Effectivenonadaptive Representation for Objects With Edges, in A. Cohen, C. Rabut, and L. Schumaker, Eds. Nashville, TN, USA: Vanderbilt Univ. Press, 2000, pp. 105–120. [13] Z. Feng, F. Liu, M. Jiang, S. Crozier, H. Guo. and Y. Wang, “Improved l1-SPIRiT using 3D Walsh transform-based sparsity basis,” Magn. Reson. Imag., vol. 32, no. 7, pp. 924–933, 2014. [14] M. Hong, Y. Yu, H. Wang, F. Liu, and S. Crozier, “Compressed sensing MRI with singular value decomposition based sparsity basis,” Phys. Med. Biol., vol. 56, no. 19, pp. 6311–6325, 2011. [15] Y. Yu, J. Jin, F. Liu, and S. Crozier, “Multidimensional compressed sensing MRI using tensor decomposition-based sparsifying transform,” PLoS One, vol. 9, no. 6, p. e98441, Jun. 5, 2014. [16] D. S. Taubman and M. W. Marcellin, JPEG 2000: Image Compression Fundamentals, Standards and Practice (International Series in Engineering and Computer Science). Norwell, MA, USA: Kluwer, 2002. [17] D. Donoho, “Compressed sensing,” IEEE Trans. Inf. Theory, vol. 52, no. 4, pp. 1289–1306, Apr. 2006. [18] R. Chartrand, “Fast algorithms for nonconvex compressive sensing: MRI reconstruction from very few data,” in Proc. IEEE Int. Symp. Biomed. Imag.: From Nano to Macro, Jun. 28/Jul. 1, 2009, pp. 262–265. [19] F. A. Razzaq, S. Mohamed, A. Bhatti, and S. Nahavandi, “Non-uniform sparsity in rapid compressive sensing MRI,” in Proc. IEEE Int. Conf. Syst. Man, Cybern., 2012, pp. 2253–2258. [20] Y. Ma, D. Baron, and D. Needell, “Two-part reconstruction in compressed sensing,” in Proc. IEEE Global Conf. Signal Inf. Process., 2013, pp. 1041–1044. [21] E. J. Cand`es, J. K. Romberg, and T. Tao, “Stable signal recovery from incomplete and inaccurate measurements,” Commun. Pure Appl. Math., vol. 59, pp. 1207–1223, 2006. [22] R. G. Baraniuk, R. A. DeVore, G. Kyriazis, and X. M. Yu, “Near best tree approximation,” Adv. Comput. Math., vol. 16, pp. 357–373, 2012. [23] H. Oh and S. Lee, “Visually weighted reconstruction of compressive sensing MRI,” Magn. Reson. Imag., vol. 32, pp. 270–280, 2014. [24] E. J. Cand`es, “Compressive sampling,” presented at the Int. Congr. Mathematicians, Madrid, Spain, 2006. [25] K. T. Block, M. Uecker, and J. Frahm, “Undersampled radial MRI with multiple coils. Iterative image reconstruction using a total variation constraint,” Magn. Reson. Med., vol. 57, no. 6, pp. 1086–1098, 2007. [26] 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,” Magn. Reson. Med., vol. 61, no. 1, pp. 103–116, 2009. [27] D. Liang, E. V. DiBella, R. R. Chen, and L. Ying, “k-t ISD: dynamic cardiac MR imaging using compressed sensing with iterative support detection,” Magn. Reson. Med., vol. 68, no. 1, pp. 41–53, 2012. [28] S. G. Lingala, H. Yue, E. DiBella, and M. Jacob, “Accelerated dynamic MRI exploiting sparsity and low-rank structure: k-t SLR,” IEEE Trans. Med. Imag., vol. 30, no. 5, pp. 1042–1054, May 2011.
Authors’ photographs and biographies not available at the time of publication.
