Abstract

It is challenging to save acquisition time and reconstruct a medical magnetic resonance (MR) image with important details and features from its compressive measurements. In this paper, a novel method is proposed for longitudinal compressive sensing (LCS) MR imaging (MRI), where the similarity between reference and acquired image is combined with joint sparsifying transform. Furthermore, the joint sparsifying transform with the wavelet and the Contourlet can efficiently represent both isotropic and anisotropic features and the objective function is solved by extended smooth-based monotone version of the fast iterative shrinkage thresholding algorithm (SFISTA). The experiment results demonstrate that the existing regularization model obtains better performance with less acquisition time and recovers both edges and fine details of MR images, much better than the existing regularization model based on the similarity and the wavelet transform for LCS-MRI.

1. Introduction

Magnetic resonance imaging (MRI) is a noninvasive and nonionizing imaging modality that offers a variety of contrast mechanisms and enables excellent visualization of anatomical structures and physiological functions. However, the MRI data sampled in the spatial Fourier transform or K-space of an object are acquired sequentially in time. Long waiting time and slow scanning may result in patients’ annoyance and affect both clinical throughput and image quality due to local motion such as breathing and heart beating.

Recently the developments in compressive sensing (CS) theories [13] show that the number of measurements required by Nyquist theorem can be significantly reduced under the condition that the underlying signal is sparse or approximately sparse in some transform or dictionary and that the measurement acquisition procedure is incoherent, in an appropriate sense, with the transform. According to the CS theories, MRI images with a sparse representation can be recovered from randomly undersampled K-space data which provided an appropriate nonlinear recovery scheme. Compressive sensing magnetic resonance imaging (CS-MRI) becomes one of the most successful applications of compressive sensing, since MR scanning time can be significantly reduced [47].

An appropriate choice of a sparsifying transform is of vital importance to the quality of the solution. The existing sparsifying transforms, for example, wavelet, total variation (TV), and TV-like transforms [811], can exploit the image structure with piecewise constant intensities; however, they are not suitable for recovering images with details and fine structure. Although dictionaries trained with samples [12] perform much better, they still suffer from artifacts and expensive computation.

There is another method that exploits similarity within a series of MRI images. In many clinical MRI scenarios, the existing imaging information can be used to significantly reduce acquisition time and improve image quality. In some cases, a previously acquired image (denote by reference image) may exhibit similarity to the image being acquired. In fact, the similarity between different contrasts in the same scan exists in multiple-contrast MRI [13] and between different scans of the same patient in dynamic MRI [14]. These similarities can be exploited by enforcing sparsity along the temporal dimension, which are called prior information (reconstruct a target CS signal with the aid of a similar signal that is known beforehand) in [15]. It is shown in [15] that if priors were employed, the bound on the number of measurements that is required for perfect reconstruction can be much smaller than the bound for the classical CS. PANO [16] makes use of the similarity of image patches to provide sparse representation and achieves lower reconstruction error and higher visual quality than conventional CS-MRI methods.

Recently, for LCS-MRI, Weizman et al. [17, 18] proposed a general sampling and reconstruction scheme, where the degree of similarity to the reference image is considered and the sparsity is exploited in both image and transform domain via the similarity and the wavelet, respectively. But their method still suffers from point-like artifacts due to the inherent shortcoming of wavelet transform.

The wavelet transform has been successfully used in signal processing as it can approximate piecewise smooth one-dimension signals very efficiently. However, as the two-dimensional wavelet transform is constructed by the tensor product of one-dimension wavelet transform, it has only 2.5 dimensions, which makes it hard to capture complex geometric properties which widely existed in images. To cope with this problem, a series of multiscale geometric transforms [23] were proposed to provide better approximation of image with fewer atoms than the wavelet. Among various multiscale geometric transforms, it is regarded that the Contourlet has the best sparse representation of smooth contours, which common exist in MRI images [24, 25]. However, the main drawback of the Contourlet is that it is not suitable for representing point-like features, which can be sparsely represented by the wavelet. It is natural to combine wavelet and Contourlet as the joint sparsifying transform to produce more physically plausible solutions instead of using either of them individually [26].

Inspired by this idea, in this paper, a novel reconstruction method is proposed for the LCS-MRI, where the similarity between the reference image and the acquired image is employed to save acquisition time and improve peak signal to noise ratio (PSNR). Furthermore, the wavelet and the Contourlet are combined to sparsely represent the underlying image. In addition, to deal with the nondifferentiable terms in our model and to achieve fast convergence of the reconstruction, SFISTA [27] is extended to solve the proposed model. Experiment results indicate that the proposed algorithm improves the quality of reconstructed images in terms of both visible quality and the PSNR. Moreover, it outperforms state-of-the-art MR image reconstruction methods including SparseMRI [4], TVCMRI [19], RecPF [11], FCSA [20, 21], and WaTMRI [22] in reducing artifacts and improving both the PSNR and the visual quality. Although patch-based nonlocal operator (PANO) [16] is the best, it consumes very long time which should be tried to avoid for MRI. Therefore, our proposed method is valuable and competitive for LCS-MRI.

The rest of the paper is organized as follows. The proposed model and an extension of the reconstruction algorithm are presented in Section 2. In Section 3, numerical results are illustrated to show the consistent performance of the proposed method. Finally, a conclusion is given in Section 4.

2. Proposed Model and Algorithm

2.1. Proposed Model

In the conventional CS-MRI, given the measurements and the sampling matrix , -pixel complex valued image (the vectorization of a two-dimension image) can be reconstructed by solving the following problem:where is a sparsifying transform operator and controls the fidelity of the reconstruction to the measured data. Conventionally, problem (1) is solved via the Lagrangian method:where is a regularization parameter. In many clinical MRI scenarios, a previously acquired image can be served as a reference image, which may exhibit similarity to the image being acquired. Thus, the similarity can be exploited to enforce sparsity in image domain and reduce sampling rates.

It is assumed that is the previously acquired image of the same patient, which is similar to in most image regions. This temporal similarity is embedded into model (2) as a regularization term and leads to the following LCS-MRI reconstruction problem:where and are regularization parameters that trade the sparsity between wavelet and spatial similarity domain.

Note that the similarity between the reference image and the acquiring image may be low or even invalid; therefore, it is necessary to embed reweight scheme into (3) to improve the performance. The resulted adaptive-weighted LCS-MRI recovery problem is as follows:where is a diagonal matrix; that is, The wavelet transform represents isotropic image features much more sparsely than anisotropic ones. In comparison, the Contourlet transform can efficiently represent anisotropic features but not be suitable for representing point-like features. To sparsify more different types of image features, we take joint transform into account for problem (4). Thus the problem becomeswhere and are wavelet and Contourlet transform, respectively. Combining the wavelet transform and the Contourlet transform can preserve the abundant geometric information and enforce the sparsity of MR images. and of the second term in (5) are used to weight specific wavelet and Contourlet atoms in the reconstruction process, which would relax the demand for sparsity on elements in the support of and . of the third term in (5) is used to weight image regions according to their similarity with the reference scan, which allows exploitation of temporal similarity when it exists and prevents degradation of image quality if the consecutive scans are significantly different. To achieve these goals, similar to [17], we first sample K-space samples randomly and reconstruct , by solving (5) with , and . Next, we sample additional samples and solve (5) where the elements of , for , and are chosen as follows:where denotes the th element of the vector in brackets. The rest of the samples are taken iteratively, taking into account both the difference between the K-spaces of the reconstructed image, , and the previous image in the time series, , and the variable density probability density function.

2.2. Optimization

To solve the weighted -minimization problem (5) in the weighted reconstruction phase, an extension of the SFISTA [27] is employed. Recently, the SFISTA is a fast iterative thresholding algorithm- (FISTA-) based algorithm, which allows for solving the -norm minimization problem based on analysis sparse model:where is an analysis dictionary under which is sparse. The extended SFISTA algorithm is summarized as follows.

Algorithm 1 (extended SFISTA algorithm for solving (5)).
Input. K-space measurements: Sparsifying transform operator: An K-space undersampling operator: Reference image: Tuning constants: , , , , , ,An upper bound: InitializeIterationsStep : () Compute + Output. Estimated image: ,
where the notation for matrices denotes the largest singular value. The elementwise soft shrinkage operator is defined asTo simplify notations, we use the following expressions:

Algorithm 1 aims to minimize (5), where the overall convergence is controlled by . We used , , and in our experiments.

3. Experiment Results

To verify the performance of the proposed algorithm for LCS-MRI, the experiments are designed in two different clinical MRI scenarios. The first scenario utilizes similarity between different scans of the same patient for fast scanning of follow-up scans. The second scenario utilizes similarity between the T2-weighted and the fluid-attenuated inversion recovery (FLAIR) for the fast FLAIR scanning. In this scenario, the reference image is a different imaging contrast of the patient from the same scan. In these applications, the performance of joint sparsifying transform is compared to the performance of single transform for LCS-MRI. The source data and partial code are available for download at the link listed in http://webee.technion.ac.il/people/YoninaEldar/software_det13.php.

The Daubechies 4 wavelet transform and the Contourlet with , , , and directional subbands from coarse to fine scales are employed in all experiments. The parameter is set to . All the compared algorithms are run on Matlab R2010a with the Intel(R) Core(TM) i CPU M @GHz processor, the RAM is GB, and the system is Windows 7 Ultimate.

3.1. Combination with Similarity between Baseline of Follow-Up Scans

Repeated brain MRI scans of the same patient every few weeks or months are very common for follow-up of brain tumors. This subsection examines the performance of proposed algorithm in the case that uses a previous scan in the time series as a reference scan for reconstruction of a follow-up scan. According to a power of distance from the origin, a variable density random undersampling scheme [4] is adopted in this paper. We begin with variable density random undersampling of 2% of the K-space and acquire additional 2% of the K-space at each iteration as described in Algorithm of [18]. Different values of , , and in the range of [0.001, 10] are taken for single transform and and in the range of [3.1 × 10−4, 10] and in the range of [0.001, 10] are taken for joint transform.

The reconstruction results of a follow-up brain scan under different sparsifying transform are shown in Figure 1. Results are obtained using only of K-space data. It is shown in Figure 1(d) that the wavelet is efficient representations for disparity maps, while it can not capture smooth contours. As comparison, the Contourlet has better representation for lines and curves, while it is not suitable for point feature as shown in Figure 1(e). Figure 1(f) presents the result by jointly utilizing wavelet and Contourlet as sparsifying transform. It has better performance by combining wavelet and Contourlet than using either of them alone. It can be seen that the results obtained by using joint transform exhibit imaging features that are hardly visible when using the wavelet or Contourlet transform alone.

3.2. Combination with Similarity between T2-Weighted and FLAIR

In the scenario, the similarity between T2-weighted scan and FLAIR scan is employed to further validate the proposed method. MRI consists of multiple imaging contrasts among various soft tissues and organs. The FLAIR and the T2-weighted scans are known to exhibit high similarity in regions with low fluid concentration [17, 18]. Based on this similarity, the FLAIR image is reconstructed from 15% of K-space data and the T2 image is served as the reference. Sampling pattern is shown in Figure 2(a). , and , are for wavelet and Contourlet transform alone, respectively. , , and are for the joint transform.

Figure 2 shows the original and reconstructed images with different transforms. It can be observed that the joint transform based FLAIR reconstruction outperforms the traditional wavelet or the Contourlet based methods, as either wavelet or Contourlet is limited in which it can just suppress certain kinds of artifact, while producing inferior results for other kinds of artifact. The superiority of the proposed approach is achieved due to the joint transform that can effectively suppress square-like and curve-like artifact.

To evaluate the quality of reconstructed results, the PSNR is formulated aswhere denotes the maximum possible pixel value in the image and is the Mean Square Error (MSE) between the original image and the reconstructed image .

Table 1 shows the PSNR values of the reconstruction results based on different sparsifying transforms for two similarity scenarios. It indicated that joint sparsifying transform outperforms the single sparsifying transform which is consistent with the visual performance in Figures 1 and 2. In addition, it is shown that the T2-FLAIR experiment provides a lower range of PSNR values in comparison to the follow-up experiment. This can be explained by the fact that similarity is not enforced over the entire image in this case, due to the difference of several regions between FLAIR and T2 weighted scans.

3.3. Parameters Sensitivity Analysis for Two Experiment Scenarios

There are two parameters involved in the proposed algorithm, that is, regularization parameters with single transform or and with joint transform and and approximate parameter . The sensitivity of these parameters is examined, and PSNR and RE results of each experiment are also investigated for various values of , , , and . The relative error (RE) is defined aswhere is an estimation of original -sparse vector .

For two clinical scenarios, the single transform with and values of [0.001 0.01 0.05 0.1 1 10 20 30 50 80 100] and the joint transform with , , and values of [0.001 0.01 0.05 0.1 1 10 20 30 50 80 100] are tested, when and .

Figures 3 and 4 show PSNR results. From two experiment scenarios, it can roughly be conclude that lower values of , , and provide better PSNR than higher ones whether using a single transform or using joint transforms. When , , and are in [], high PSNR is obtained. In other words, it shows the fact that overpromoting sparsity versus consistency to measurements degrades the reconstruction quality.

As it was argued in [28], the SFISTA is sensitive to the approximate parameter . How the approximate parameter affects the empirical convergence of LCS-MRI is also investigated and two MRI scenarios with values of [0.0000001 0.000001 0.00001 0.0001 0.001 0.01 0.1 1 10 100] are tested for the extended SFISTA, which is shown in Figure 5.

It can be seen that smaller leads to a more accurate reconstruction, which is consistent with [27, 28]. In our method, we inherit the continuation technique of [27] to accelerate the convergence while obtaining a desired reconstruction accuracy.

3.4. Comparison with the State-of-the-Art Algorithms for Two Experiment Scenarios

In this subsection, the proposed method is compared with the state-of-the-art methods, that is SparseMRI [4], TVCMRI [19], RecPF [11], FCSA [20, 21], WaTMRI [22], and PANO [16]. For fair comparisons, all codes are downloaded from the authors’ websites and their default settings are used. It is worth emphasizing that, among the 7 comparison methods, only PANO and the proposed method use the similarity prior information for MR image reconstruction. The zero-filling image and the baseline scan are used as the reference image in PANO and the proposed method, respectively. All algorithms terminate after 30 iterations. To provide quantitative measures for our experiential results, PSNR, RE, and CPU time performance are examined.

The yellow line is the results of the proposed method in Figure 6. Both Figures 6 and 7 show that the proposed method outperforms the state-of-the-art methods in terms of PSNR, RE, and visual quality except the PANO. This result is reasonable because the proposed method and PANO exploit the prior information of similarity, which can reduce requirement for the number of measurements or increase the accuracy of the solution with the same measurements. It can be seen that the PANO method is quite well; however it needs too much time. As is well known, the core problem of MRI is to reduce the time required for sampling and reconstruction. Long waiting time may result in patients’ annoyance and blur on images due to local motion such as breathing and heart beating. Therefore, our proposed method is valuable and competitive for LCS-MRI. The result of the T2-FLAIR similarity experiment is consistent with Figures 6 and 7; here it is not presented.

4. Conclusion

In this paper, a new method has been proposed for LCS-MRI, where the similarity of longitudinal MR images can be exploited to significantly reduce scan time and improve the performance. The combination of the wavelet and the Contourlet transform is employed instead of the single transform to preserve the abundant geometric information. To deal with the nondifferentiable terms in the proposed model and achieve fast convergence of the reconstruction, the SFISTA is extended to solve the optimization problem. Numerical results show the improvement of the reconstruction quality over the compared models in terms of both visible quality and the PSNR. Future work will consider the more efficient method to enforce the sparsity of MRI image with plentiful features and to apply it into a wider range of medical imaging applications.

Competing Interests

The authors declare that there are no competing interests regarding the publication of this paper.

Acknowledgments

This work is partially supported by the National Natural Science Foundation of China (61271012; 61671004).