Undersampled MR Image Reconstruction with Data-Driven Tight Frame

Undersampled magnetic resonance image reconstruction employing sparsity regularization has fascinated many researchers in recent years under the support of compressed sensing theory. Nevertheless, most existing sparsity-regularized reconstruction methods either lack adaptability to capture the structure information or suffer from high computational load. With the aim of further improving image reconstruction accuracy without introducing too much computation, this paper proposes a data-driven tight frame magnetic image reconstruction (DDTF-MRI) method. By taking advantage of the efficiency and effectiveness of data-driven tight frame, DDTF-MRI trains an adaptive tight frame to sparsify the to-be-reconstructed MR image. Furthermore, a two-level Bregman iteration algorithm has been developed to solve the proposed model. The proposed method has been compared to two state-of-the-art methods on four datasets and encouraging performances have been achieved by DDTF-MRI.


Introduction
Magnetic Resonance Imaging (MRI) is one of the major diagnostic imaging modalities with noninvasive and nonionizing radiation nature. However, relatively low imaging speed limits its wide application in clinic [1]. To accelerate MRI, one popular way is to reduce the number of acquired data [2]. With ∈ C and ∈ C , respectively, denoting the -space measurement and the original image, the imaging model could be described as follows: where is the complex additive white Gaussian noise with standard deviation , is the undersampling operator, is the Fourier operator, and = .
Recovering from the undersampled measurement is an ill-posed inverse problem. To address this ill-posed nature, it is necessary to utilize some prior knowledge to regularize the MR image so as to make up the missing information. With the popularity of compressed sensing (CS) theory [3,4], the sparsity-promoting regularization for MR image reconstruction has attracted many researchers [2,5,6]. Specifically, the CS theory has shown that if an image has a sparse representation under certain transform, we can precisely restore the original image from its partial measurements under the RIP condition [3,4]. With such a transform, the MR image reconstruction from its undersampledspace data can be achieved by nonlinear algorithms, like ℓ 1 minimization or orthogonal match pursuit (OMP) algorithm [2,5,6].
The commonly used transforms in MR image reconstruction include total variation (TV) and wavelet transform, both of which can be regarded as tight frames [7][8][9][10][11][12][13]. As the tight frame satisfies the perfect reconstruction property, which ensures that the given signal can be perfectly represented by its canonical representation, it has been leveraged in diverse inverse problems [7,11,12]. However, TV prior assumes that the image consists of piecewise constant areas which may not be valid in many real MR images. When the measurements are highly undersampled, compressed sensing based MRI using TV regularization (CSMRI-TV) could lead to severe blocky artifacts [2,14]. To improve the image quality, 2 Computational and Mathematical Methods in Medicine Liang et al. [15] applied the nonlocal total variation (NLTV) regularization in parallel MR imaging by replacing the gradient functional in conventional TV with a weighted nonlocal gradient function. However, although NLTV reduces the blocky effects, it is still based on fixed transform, which does not adapt to the target image. Furthermore, other analytically designed X-lets [8][9][10]16], such as wavelets and shearlets, also suffer from their intrinsic deficiencies [17], such as the dependency between parent and child wavelet coefficients and the lack of adaptability. Qu et al. [18] introduced the patch-based directional wavelet (PBDW) in MR image reconstruction by exploiting the geometric direction of image patches, which has shown encouraging performances on edge preservation and noise removal. Nevertheless, PBDW is still a simplified form of bandlets and the adaptability can still be explored [19].
Since the fixed tight frame/transform might not be universally optimal for all images, the data-driven tight frame/transform has been developed [5,6,20,21] to adaptively capture the structure information. One most popular direction is the incorporation of the dictionary learning (DL) into MRI [6,21,22], among which Ravishankar and Bresler proposed a benchmark work named DLMRI [5] with outstanding reconstruction results achieved. However, dictionary learning based optimization is a large-scale and highly nonconvex problem, which requires high computational complexity [6,20]. More recently, Cai et al. [20] presented a scheme to learn a data-driven tight frame from the measurements and applied it to solve the image denoising problem. As the data-driven tight frame satisfies the perfect reconstruction property, it is very efficient to obtain the results with less artifacts. Therefore, Wang et al. [23] tried to incorporate data-driven tight frame into dynamic MRI. However, the proposed method relies on a perfect reference image, which is quite hard to be obtained from its undersampledspace data, and the data-driven tight frame is learnt from the reference image instead of the target image.
Based on these observations and motivated by the efficiency and effectiveness of the data-driven tight frame (DDTF), in this work, we propose an undersampled MR imaging method by incorporating DDTF into the reconstruction model with the aim of further improving the image reconstruction accuracy without introducing too much computation. Specifically, a tight frame has been adaptively trained for each to-be-reconstructed MR image. To solve the proposed model, a two-level Bregman iteration algorithm has been developed. We name the proposed approach as DDTF-MRI and compare it to two state-of-the-art MR image reconstruction algorithms, including CSMRI-TV and DLMRI, on one simulated MR image and three in vivo datasets under different undersampling scenarios. The experiments have shown encouraging performances of the proposed method.
The rest of the paper is organized as follows. Section 2 presents the preliminary and previous work. Our proposed method DDTF-MRI is illustrated in Section 3. The experimental results are provided in Section 4 demonstrating that the proposed algorithm has potential to improve the MR image reconstruction from its undersampled -space data. Finally, conclusions are given in Section 5.

Preliminary and Previous Work
where ⟨⋅, ⋅⟩ is the inner product operator and the space 2 (C) represents the set of all the functions ( ) satisfying ‖ ‖ 2 (C) := (∫ C | ( )| ) 1/2 . For a given tight frame, we can define the corresponding analysis operator as and the synthesis operator as The system is a tight frame if and only if = , where is the identity operator. Moreover, a tight frame can be generated from a filter bank satisfying the Unitary Extension Principle (UEP) condition [7]. Let { } = 2 =1 denote the set of 2D filters and the size of is × and then the synthesis operator ∈ C × 2 and the analysis operator where ∈ C × denotes the convolution matrix associated with the filter , which is a block-wise Toeplitz matrix under Neumann boundary conditions [7,24].

Sparsity-Regularized MR Image
Reconstruction. This subsection briefly introduces two representative sparsityregularized MR image reconstruction methods, namely, CSMRI and DLMRI. CSMRI is a classical approach exploiting sparsity via fixed transform and its mathematical model can be defined as follows: where is the analysis operator of its corresponding frame system, that is, wavelet or total variation. DLMRI proposed by Ravishankar and Bresler [5], on the other hand, employs the -singular value decomposition ( -SVD) algorithm [21] to learn an adaptive dictionary to sparsely represent the patches extracted from the target image. The DLMRI model can be written as Computational and Mathematical Methods in Medicine   3 where represents the dictionary, denotes the operator that extracts the th patch from the image , and Γ denotes the set of all sparse coefficients { }.

Image Reconstruction Model.
To enhance the sparsity in the frame domain, the proposed DDTF-MRI method is devoted to learning an adaptive data-driven tight frame to effectively sparsify the target image. The proposed model can be described as where Λ = { | = } and the tight frame is constructed by the corresponding 2D filters { } = 2 =1 in the manner of (5).

The Proposed Algorithm.
A two-level Bregman iteration algorithm has been developed to attack the target model.

The Inner Bregman Iteration.
With introduction of an assistant variable, V = , we can obtain the constrained version of the first subproblem in (9): Employing the Bregman iteration technique again under the assistance of a dual variable, , in the inner-level Bregman iteration, we target solving where > 0 and 0 < ≤ 1 [26]. (11). For the -subproblem in (11), as = , its least squares solution can be obtained as

Subproblems of
Multiplying both sides of (12) by and letting 1 = − and 2 = [ (V − )], we obtain the interpolation formulation in the frequency domain: where Ω stands for the subset of -space that has been sampled. The image can then be updated via the inverse Fourier transform. Next, we apply the alternating-minimization strategy to solve the {V, }-subproblems: With temporarily fixed, the V-subproblem becomes Following the iterative shrinkage/thresholding algorithm (ISTA) [27], it yields where shrink( , ) = sign( ) max{0, | | − }. Now let V be kept fixed; the -subproblem turns to be Instead of directly optimizing , we apply the technique of [20] to solve this subproblem by using SVD to obtain the corresponding filters { }.
Furthermore, for a better restoration result, we also update the variable V in the updated tight frame domain and obtain To help the readers grasp the overall picture better, we summarize the entire process in Algorithm 1. (1) Initialization: upda te +1 according (13)  (5) upda teV +1/2 according (16)  (6) upda te +1 according (17)

Experimental Setup.
We evaluated the proposed DDTF-MRI on four datasets and compared it to two state-of-theart methods, including CSMRI-TV (the CSMRI-TV code is available in http://www.eecs.berkeley.edu/∼mlustig/CS.html) and DLMRI (the DLMRI code is available in http://www .ifp.illinois.edu/∼yoram/DLMRI-Lab/DLMRI.html) [5]. The four datasets consist of one simulated data and three in vivo datasets. The simulated data is from the paper in [5]. For the in vivo data, informed consent was obtained from the volunteer in accordance with the institutional review board policy. The first in vivo data is a sagittal brain collected on a GE 3T scanner (GE Healthcare, Waukesha, WI, USA) with 32-channel head coil and 3D T1-weighted spoiled gradient echo sequence: TE = minimum full, TR = 7.5 ms, FOV = 24 × 24 cm, matrix = 256 × 256, and slice thickness = 1.7 mm. The heart dataset was acquired on a 1.5 T Philips scanner using the steady-state free precession (SSFP) sequence with a flip angle of 50 degree and TR = 3.45 msec. The field of view (FOV) was 345 mm × 270 mm and the slice thickness was 10 mm. Retrospective cardiac gating was used with a heart rate of 66 bpm. The third in vivo data was acquired on a 3 T commercial scanner (GE Healthcare, Waukesha, WI, USA) and eight-channel head coil (Invivo, Gainesville, FL, USA) with a 2D T1-weighted spin-echo protocol (axial plane, TE/TR = 11/700 ms, 22 cm FOV, 10 slices, matrix size = 256 × 256). The adaptive combination method [28] is applied to integrate the multichannel data into a singlechannel complex-valued image before we apply the three methods. Then the full -space data corresponding to the single channel image is retrospectively undersampled under different sampling schemes.
As shown in Figure 1(a), we employed the 3-level shiftinvariant Haar wavelet filters (the size of each filter is 8 × 8) as the initialization of the tight frame in DDTF-MRI. As for the parameter settings, both CSMRI-TV and DLMRI were implemented with their default settings. For DDTF-MRI, we set = 3, = 1, = 1, = 10, and = 10. The outer-level iteration of DDTF-MRI continues until > 25. Furthermore, we used both the peak signal-to-noise ratio (PSNR) (the PSNR is defined as PSNR = 20log 10 255/RMSE, where the RMSE is the root mean error estimated between the ground truth and the reconstructed image), high-frequency error norm (HFEN) [5], and structural similarity (SSIM) index [29] for a quantitative comparison of recovered results.

Test on the Simulated Data.
We firstly applied each algorithm to reconstruct the simulated T2-weighted image under the 1D Cartesian sampling schemes with accelerating factor = 6.7 (sampling ratio 15%). The filters learned by DDTF-MRI are shown in Figure 1(b). Compared with the initialization filter presented in Figure 1(a), we can see that the learnt filters have captured more directional information.  Moreover, a region-of-interest (ROI) analysis was performed by calculating means and standard deviations in selected ROIs. Table 1 shows the numerical result calculated from three selected ROIs in this experiment. The numbers of pixels in ROI 1, ROI 2, and ROI 3 are 1944, 1248, and 1944, respectively. From Table 1, we can find the values of our reconstruction result are closest to the values of reference image, which illustrates that our method can give a more precise reconstruction result.

Test on In Vivo Datasets.
We then employed these three methods to reconstruct the sagittal brain and test their sensitivity to the acceleration factors. Figure 3(c) plots the PSNR values versus different accelerating factors under the random sampling trajectory. It shows that DDTF-MRI performs better than the other two methods at all the acceleration factors. For a visual comparison, Figure 4 provides the reconstructed results at accelerating factor 3. For a close-up comparison, we have enlarged two edge parts, from which we can see DDTF-MRI provides clearer details.
To investigate the sensitivity of various methods to noise, CSMRI-TV, DLMRI, and DDTF-MRI were applied to reconstruct the heart image under pseudoradial sampling scheme with accelerating factor = 3. Figure 5(c) provides the PSNR values of the MR images reconstructed by CSMRI-TV (red curves), DLMRI (green curves), and DDTF-MRI (blue curves) versus diverse noise levels ( = 2, 4, 6, 8).  We can see that our proposed DDTF-MRI exhibits the best performance for these noise levels. Figures 5(d), 5(e), and 5(f) show the reconstructed error magnitudes from the 3-fold radial undersampled data with noise standard deviation as 4. We can find that our method can preserve more detail and recover the structures more precisely. In Figure 6, we evaluated the proposed method with another in vivo brain data, which contains more fine-detailed structures, using the radial sampling trajectory. The reconstructed results using CSMRI-TV, DLMRI, and DDTF-MRI with a higher acceleration factor = 4 are displayed in Figures 6(b), 6(c), and 6(d), respectively. The zoom-in results are also provided in Figure 6. Compared to the reference image shown in Figure 6(a), it can be observed that the result obtained by CSMRI-TV suffers from blocky artifacts. Meanwhile, the reconstructed image obtained by DDTF-MRI is clearer and sharper than those reconstructed by CSMRI-TV and DLMRI. This reveals that our proposed method can provide a more accurate recovered image.

Computation Time.
In our experiments, all algorithms were implemented in MATLAB and performed on a laptop equipped with Intel 2.60 GHz CPU and 4 GB RAM. Table 2 illustrates the cost time for these experiments. These computational times were obtained by averaging 10 times for each experiment. The size of the simulated image (Figure 2 on fixed transform, it is the fastest. However, this efficiency comes at the cost of accuracy. DLMRI takes the longest time in all cases. On the other hand, DDTF-MRI has achieved the best reconstruction result without taking the longest reconstruction time, which is due to three aspects: (1) ℓ 1 norm is used in DDTF-MRI instead of ℓ 0 norm. While ℓ 0 norm minimization is a NP-hard problem, ℓ 1 norm minimization is convex and can be solved with much less time while promoting sparsity; (2) the tight frame used in this paper ensures that the given signal can be perfectly represented by its canonical expansion as = ; therefore it can be efficiently implemented; (3) two-level Bregman iteration technique is used in the proposed method, which can further promote fast and accurate MR image reconstruction. So the proposed method has improved the accuracy without sacrificing too much efficiency.

Conclusions
A DDTF-MRI method has been proposed in this paper to effectively reconstruct MR image from undersampled -space data. DDTF-MRI trains an adaptive tight frame for each to-be-reconstructed image. Furthermore, a twolevel Bregman iteration algorithm has been developed to solve the proposed model. Results on both simulated and in vivo datasets demonstrate the superior performance of DDTF-MRI over the other two state-of-the-art MR image reconstruction methods, including CSMRI-TV and DLMRI, in artifacts suppression and edge preservation. Although DDTF-MRI is slower than the fixed transform based method CSMRI-TV, its accuracy is much higher than CSMRI-TV. Furthermore, as an adaptive training method, DDTF-MRI is much faster than the typical training method, DLMRI. This indicates that our method can improve the reconstruction accuracy without introducing too much computation load. In the future, we may further optimize the implementation and consider sparser representations. also like to thank Lustig and Ravishankar for sharing their experiment materials and source codes. This work was partly supported by China NSFC 61102043, 11301508, 81120108012, 61471350, and KQCX20120816155710259 and Shenzhen Basical Research Project (JCYJ20140610152828678).