Combined Similarity to Reference Image with Joint Sparsifying Transform for Longitudinal Compressive Sensing MRI

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-basedmonotone 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.


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 [1][2][3] 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 [4][5][6][7].
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 [8][9][10][11], 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 twodimensional 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.

Proposed Model and Algorithm
2.1.Proposed Model.In the conventional CS-MRI, given the measurements  ∈ C  and the sampling matrix  ∈ C × , -pixel complex valued image  ∈ C  (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  0 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 becomes where Φ 1 and Φ 2 are wavelet and Contourlet transform, respectively.Combining the wavelet transform Φ 1 and the Contourlet transform Φ 2 can preserve the abundant geometric information and enforce the sparsity of MR images. 1,1 and  1,2 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 Φ 1  and Φ 2 . 2 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 x, by solving (5) with  1, = ,  = 1, 2, and  2 = 0. Next, we sample additional   samples and solve (5) where the elements of  1, , for  = 1, 2, and  2 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,  0 , and the variable density probability density function.

Optimization.
To solve the weighted ℓ 1 -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 ℓ 1 -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.

Iterations
Step : Output.Estimated image: x =   , where the notation ‖ ⋅ ‖ 2 for matrices denotes the largest singular value.The elementwise soft shrinkage operator Γ  () is defined as To simplify notations, we use the following expressions: Algorithm 1 aims to minimize (5), where the overall convergence is controlled by .We used  1 = 0.31,  2 = 0.69, and  = 0.0364 in our experiments.

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 followup 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 2 5 , 2 4 , 2 4 , and 2 3 directional subbands from coarse to fine scales are employed in all experiments.The parameter  is set to 0.1.All the compared algorithms are run on Matlab R2010a with the Intel(R) Core(TM) i3 CPU M 370 @2.40 GHz processor, the RAM is 2 GB, and the system is Windows 7 Ultimate.

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 1 of [18].Different values of  1 ,  2 , and  in the range of [0.001, 10] are taken for single transform and  1 and  2 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 6% 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.

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 T2weighted 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). 1 = 0.025,  = 0.03 and  2 = 0.025,  = 0.03 are for wavelet and Contourlet transform alone, respectively. 1 = 0.0077,  2 = 0.0172, and  = 0.03 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 as where  denotes the maximum possible pixel value in the image and   is the Mean Square Error (MSE) between the original image  and the reconstructed image x.
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.

Parameters Sensitivity Analysis for
where x is an estimation of original -sparse vector .(f) WaTMRI [22], (g) PANO [16], and (h) the proposed algorithm.Figures 3 and 4 show PSNR results.From two experiment scenarios, it can roughly be conclude that lower values of  1 ,  2 , and  provide better PSNR than higher ones whether using a single transform or using joint transforms.When  1 ,  2 , and  are in [1,10], 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.

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.

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.

Figure 1 :
Figure 1: Sampling pattern and reconstructed images based on similarity within repeated scans scenario.(a) Sampling pattern with 2% of K-space data sampled.(b) Sampling pattern with 6% of K-space data sampled.(c) The ground truth.(d) Reconstructed image of wavelet transform.(e) Reconstructed image of Contourlet transform.(f) Reconstructed image of joint sparsifying transform.

Figure 2 :
Figure 2: Sampling pattern and reconstructed images based on the similarity within same scans scenario.(a) Sampling pattern with 15% of Kspace data sampled.(b) The original T2 image.(c) The original FLAIR image.(d) Reconstructed image of wavelet transform.(e) Reconstructed image of Contourlet transform.(f) Reconstructed image of joint sparsifying transform.

FollowFigure 3 :Figure 4 :
Figure 3: Sensitivity analysis of follow-up experiment with 6% sampling: the PSNR results for various values of regular parameters  and  with single wavelet and  1 ,  2 , and  with joint transform.(a) PSNR for various values of  and  with wavelet transform.(b) PSNR for various values of  and  with Contourlet transform.(c) PSNR for various values of  1 and  with joint transform.(d) PSNR for various values of  2 and  with joint transform.(e) PSNR for various values of  1 ⋅  1 and  with joint transform.(f) PSNR for various values of  2 ⋅  2 and  with joint transform.

Figure 5 :
Figure 5: The PSNR and RE results with wavelet, Contourlet, and joint transform for various values of approximate parameter .(a) PSNR of follow-up experiment for various values of .(b) RE of follow-up experiment for various values of .(c) PSNR of T2-FLAIR experiment for various values of .(d) RE of T2-FLAIR experiment for various values of .It can be seen that  ∈ [1 − 7, 1] exhibits reliable results for all cases.

Figure 6 :
Figure 6: Reconstructed images with 20% sampling for follow-up similarity experiment.(a) PSNR.(b) Relative error.(c) CPU time.(d) CPU time for fourth algorithms with less time.(e) CPU time without PANO algorithm.(f) CPU time of PANO and the proposed algorithm.
Two Experiment Scenarios.There are two parameters involved in the proposed algorithm, that is, regularization parameters  with single transform or  1 and  2 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  1 ,  2 , , and .The relative error (RE) is defined as