Abstract

Compressed Sensing Magnetic Resonance Imaging (CS-MRI) is a promising technique for accelerating MRI acquisitions by using fewer k-space data. Exploiting more sparsity is an important approach to improving the CS-MRI reconstruction quality. We propose a novel CS-MRI framework based on multiple sparse priors to increase reconstruction accuracy. The wavelet sparsity, wavelet tree structured sparsity, and nonlocal total variation (NLTV) regularizations were integrated in the CS-MRI framework, and the optimization problem was solved using a fast composite splitting algorithm (FCSA). The proposed method was evaluated on different types of MR images with different radial sampling schemes and different sampling ratios and compared with the state-of-the-art CS-MRI reconstruction methods in terms of peak signal-to-noise ratio (PSNR), feature similarity (FSIM), relative l2 norm error (RLNE), and mean structural similarity (MSSIM). The results demonstrated that the proposed method outperforms the traditional CS-MRI algorithms in both visual and quantitative comparisons.

1. Introduction

Magnetic resonance imaging (MRI) has been widely used in the radiological diagnosis due to its high spatial resolution, noninvasive, and nonionizing radiation merits. It has become a powerful and promising technique to visualize and investigate the anatomical and physiological properties of many organs. Generally, to obtain a high spatial resolution image with a high signal-to-noise ratio (SNR), the image should be acquired with more k-space data and for several times to average the noise. This however will undoubtedly increase the acquisition time. Long acquisition time may introduce several problems, such as the increase of motion artifacts, especially for special patients with motion control problems and moving organs. Therefore, how to reduce the acquisition time of MRI while maintaining image quality is currently a great challenge [13].

One of the most promising approaches to deal with this issue is compressed sensing (CS) [1, 46] which reduces the acquisition time by acquiring only a small fraction of k-space data whereas enabling us to reconstruct the image from such highly undersampled measurements with high quality. It is well known that, to apply CS to MRI, finding proper k-space undersampling patterns, sparsifying transformations, and the corresponding reconstruction algorithms is fundamental for guaranteeing the reconstructed image quality [1]. In the pioneering work of CS-MRI, Lusting et al. [7] reconstructed magnetic resonance (MR) images from the Cartesian undersampled k-space data by solving a norm minimization equation with wavelet transform and total variation (TV) sparsity constraints. According to the theory of CS, the sparser the representation is, the better the reconstruction quality is. Therefore, numerous efforts have been devoted to investigating sparse representations and sparsity regularizations. For instance, sparse transformations based on some fixed basis were widely used, including discrete cosine transform (DCT), contourlet transform [8], and Shearlet transform [912]. To further explore the degree of sparsity, several sparsity regularizations and the corresponding combinations were developed, such as the nonlocal self-similarity constraint, wavelet transform regularization, TV regularization, and nonlocal TV (NLTV) regularization [1315]. The fixed sparsifying transformation is very fast for implementation, but it does not guarantee getting the best sparse representation for a specific image. Accordingly, data-driven learning was proposed which learns adaptively the sparse representations from the image itself [1621], including the dictionary learning such as k-singular value decomposition (K-SVD) [22] and the transform learning based on adaptive tight frame algorithms [23, 24]. With the success of deep learning in solving inverse problems such as image denoising and reconstruction, data-driven sparsity learning using autoencoder (AE) [25] and restricted Boltzmann machine (RBM) [26] have been proposed. Although data-driven based approaches indeed improved image reconstruction quality, compared to traditional CS-MRI, the learning process requires longer time since the sparse basis training is very computation expensive, which is not practical for clinical use. Considering the complexity of optimization problem and the speed of computation, combining several simpler regularization terms to further improve the reconstruction performance of traditional CS-MRI algorithms may be an alternative. Therefore, novel CS-MRI algorithms using more sparse priors to accelerate reconstruction process while maintaining reconstruction accuracy are still greatly desired.

In this paper, we propose to use multiple sparse constraints to improve MR image reconstruction precision. More precisely, wavelet sparsity, wavelet tree sparsity, and NLTV regularization are all integrated in the CS-MRI framework. The wavelet coefficients of MR images not only are approximately sparse, but also have a tree structure. The latter provides better image reconstruction quality than standard wavelet sparsity prior [2729] and enables us to exploit the correlation between the parent and child of wavelet coefficients and reduce the required number of k-space data for MR image reconstruction. Meanwhile, NLTV model extends the conventional TV term to a nonlocal version, which can effectively avoid staircase artifacts caused by TV regularization while better preserving image edges and fine details. It is therefore expected that the combination of these three sparse constraints can give more accurate reconstruction result. In such framework, the optimization problem is solved by means of a fast composite splitting algorithm (FCSA) [30], which is an iterative algorithm in which each iteration involves a series shrinkage step. To validate the performance of the method proposed, the experiments are conducted on both chest and renal arteries MR images with several radial sampling schemes, including radial uniform, radial golden, and radial randomized sampling. The results are compared to previous CS-MRI reconstruction methods, including SparseMRI [7], RecPF [31], FCSA [30], FCSA-NLTV [32], and WaTMRI [29].

The remainder of this manuscript is organized as follows. Section 2 describes the multiple sparsity constraints-based CS-MRI reconstruction framework as well as the performance evaluation indices. The experimental results and discussions are given in Section 3, including the experimental setup, visual and quantitative comparison, and the influence of sampling ratio on the performance of the proposed method. Section 4 concludes the work.

2. Methods

In this section, we introduce the CS-MRI framework based on multiple sparse constraints and the corresponding evaluation criteria.

2.1. CS-MRI Reconstruction Framework

Assume that is the desired MR image and is an undersampled measurement in k-space domain with the sensing matrix denoted as . The aim of CS-MRI is to reconstruct the MR image with a constraint of by solving the following optimization problem:where is a regularizing function that intends to find a perfect sparse representation of the image . To exploit more sparse priors, wavelet sparsity, wavelet tree sparsity, and nonlocal TV regularization items are combined in the CS-MRI framework proposed in this paper. This multiple regularizations-based MR image reconstruction problem can be formulated as follows:where the sensing matrix is a partial Fourier transform, which can be expressed by with being the Fourier transform and being the common undersampling pattern (mask), is a wavelet transform, represents all the parent-child groups that encourage the tree sparsity and is one of such groups, and is the NLTV-based constraint defined as follows [13, 32]:with and indicating the image values at and and with representing the weight function formulated as

In Eq. (4), and represent a small patch centering at the coordinates and , respectively, signifies a normalization factor, h denotes a scale parameter that is related to the patch size and the standard deviation of noise and controls to what extent similarity between patches is enforced, and is a positive value used to control the nonlocality of the method and to speed up computation, which means that only the neighbors in a window around the target pixel are considered when calculating the nonlocal image gradient [13].

An auxiliary variable is introduced to solve Eq. (2):where is a binary matrix for wavelet coefficients group configuration. Accordingly, the optimization problem in Eq. (2) can be rewritten as follows:where is the wavelet tree group and is the total number of groups. The solution of Eq. (6) can be split into two subproblems: z subproblem and subproblem, in which, z subproblem can be expressed as follows [29]:It has a closed form solution by soft-thresholding:where . For convenience, z is represented as

The subproblem can be expressed as

Let , which is a convex and smooth function with Lipschitz constant , and , which are both convex and nonsmooth functions. Eq. (10) can be solved using the fast composite splitting algorithm (FCSA) [30], which means that the problem can be divided into multiple subproblems: the regularization of NLTV-norm and wavelet -norm. Each subproblem is actually a convex function which can be solved by a proximal mapping operation [33]:where is a positive scalar.

The subsequent algorithm is described as shown in Algorithm 1.

Input:
For k=1 to K do
End for.

In Algorithm 1,, with and representing, respectively, the inverse partial Fourier transform and inverse wavelet transform, and is the maximum iteration times.

2.2. Quantitative Evaluation Criteria

To evaluate quantitatively the performance of the proposed CS-MRI reconstruction method, the indices including peak signal-to-noise ratio (PSNR), mean structural similarity (MSSIM) [34], feature similarity (FSIM) [35], and relative norm error (RLNE) [36] were used. Denote that and represent, respectively, the original image and that reconstructed from the undersampled k-space data with CS-MRI algorithms, with a maximum value noted as MAX and the dimension of . Accordingly, PSNR is defined asLet and represent one block of the images and , respectively, the local means of which are expressed by and with , , and indicating the standard deviations and cross-covariance of images blocks and . Based on such notations, the structure similarity (SSIM) is calculated as follows:where and are the tunable constants. If the total number of image blocks is , the MSSIM can be formulated asBoth PSNR and MSSIM are proportional to image quality. The better the image quality is, the higher the PSNR and MSSIM are.

3. Experimental Results and Discussion

3.1. Experimental Setup

The experiments were performed with two real-value images, chest and renal arteries MR images with size of 256 × 256, respectively, as illustrated in Figures 1(a) and 1(b). We observed that one has rich texture and the other is relative smooth (Note: data from Ref. [29]). Since the radial sampling schemes have been demonstrated to be more feasible in practice and better than Cartesian sampling [37], we used several radial sampling masks, including radial uniform, radial golden [38, 39], and radial randomized, as shown in Figures 1(c)1(e), where 10% k-space data is kept.

The observation measurement y was modeled as , where represents complex Gaussian white noise with standard deviation . The associated input SNR (ISNR) [40] is defined as , with denoting the standard deviation of the original image. The ISNR was set to 20 dB.

In order to evaluate the benefits of multiple sparse constraints proposed in this paper, we compared quantitatively the proposed method with several reconstruction approaches, including SparseMRI [7], RecPF [31], FCSA [30], FCSA-NLTV [32], and WaTMRI [29], in terms of PSNR, MSSIM, FSIM, and RLNE indices. SparseMRI is the first work of CS-MRI, which models MR image reconstruction as a linear combination of least square fitting, wavelet sparsity, and TV regularization. The optimization problem is solved by conjugate gradient (CG) method. RecPF uses a variable splitting method to solve the same model as that in SparseMRI. FCSA decomposes the above problem into two easier subproblems and separately solves each of them with FISTA. To overcome the intrinsic drawback of the TV model, NLTV reconstruction method was proposed in FCSA-NLTV that replaces the TV in classical CS-MRI model with NLTV. WaTMRI is based on using a multiple regularization method, which combines wavelet sparsity, gradient sparsity, and tree sparsity in one model, and is formulated as the group sparsity problem for validating the benefit of wavelet tree structure in CS-MRI. For fair comparisons, all the experiments were performed with MATLAB 2012a on a desktop equipped with 3.10 GHz Intel core i5 2400 CPU and 8 GB RAM. In the experiments, Daubechies wavelets with four decomposition levels were used to form a quadtree. Daubechies wavelets were chosen because they are the most commonly used ones. They are orthogonal and compactly supported wavelets, which have excellent edge detection characteristics and fast implementation. Daubechies wavelets were also used as sparse bases in previously reported CS-MRI works and were shown to give good image reconstruction performance. Nevertheless, the use of many kinds of wavelets (Haar, Daubechies, Symmlets, etc.) other than Daubechies wavelets is possible. Concerning the number of decomposition levels, since the image size is 256 × 256, we have chosen 4 levels of decomposition. Too many decomposition levels will increase too much computation cost; too few decomposition levels will weaken tree structure benefit. Four decomposition levels in the present study seemed the best compromise. The parameter settings for different reconstruction methods were given in Table 1.

3.2. Visual Comparisons

Figures 2 and 3 give the visual comparison of the chest and renal arteries MR images reconstructed using different methods and different radial undersampling masks. The first raw shows the reconstructed images, and the second row shows the error maps between the original image and the reconstructed ones. The first column titled by “Zero-filling” represents the figures of the undersampled measurements corrupted by noise. The undersampled ratio is set to be approximately 20%. We observe that, for the chest MR image with rich texture, the reconstructions with SparseMRI, RecPF, FCSA, FCSA-NLTV, and WaTMRI contain obvious steaking artifacts, especially for the radial randomized sampling patterns, whereas the images reconstructed by the proposed method present much less artifacts than the others for all the sampling masks. For the renal arteries MR image with less edges, the reconstruction quality of all the methods is better than that of chest image. The results verify that, using our method, we can achieve more satisfying results with clear contours, sharp edges, and fine image details for various radial sampling masks and for both types of MR images. From the error maps, we see that large errors occur mainly in the zones of edges.

In the proposed method, the wavelet and wavelet tree sparsity constraints mainly contribute to reconstructing the smooth regions, and the NLTV enables us to recover the texture and edges due to the use of nonlocal similarity. The combination of these constraints makes them complementary to each other and therefore results in relatively pleasant reconstruction results.

3.3. Quantitative Comparison

In order to provide more quantitative comparison for different methods, the PSNR, MSSIM, FSIM, and RLNE are calculated for the MR images reconstructed above with different methods and various sampling schemes (undersampled ratio is about 20%), as given in Tables 25, respectively. It is clearly seen that, compared with the WaTMRI, combining NLTV, wavelet, and wavelet tree sparsity constraints improved PSNR of about 2 dB for different sampling schemes. Moreover, the proposed method achieves the best performance in the radial uniform sampling scheme for both kinds of MR images, which has the highest PSNR value.

3.4. Influence of Sampling Ratio

In order to investigate the influence of different sampling ratios on the performance of the reconstruction methods, the curve of PSNR, MSSIM, FSIM, and RLNE versus different undersampled ratios ( or 0.1~0.5) for the MR images reconstructed with different methods and undersampling schemes are shown in Figures 47. It can be seen that, for both the chest and renal arteries MR images, PSNR, MSSIM, and FSIM of the images reconstructed by the proposed method (WaTMR_NLTV) are higher than the other methods for all the sampling schemes with all the sampling ratios. As to RLNE, for both chest and renal arteries images, the RLNE obtained by our method is lower than the others for all the sampling ratios. Especially, the proposed method is much better than the others with low sampling ratio.

3.5. Influence of Noise

The noise robustness of the proposed method on different undersampling pattern MR images for different values of ISNR and sampling ratio (R) is shown in Figure 8. It is seen that the reconstruction performance changes with the sampling ratio and ISNR. Increasing sampling ratio reduces the reconstruction error (except in the cases ISNR=5dB and ISNR=10dB) and reducing ISNR increases the reconstruction error. On the other hand, the reconstruction performance changes slightly when ISNR is larger than 10 dB when fixing sampling ratios; for example, R=0.3 in Figure 8.

3.6. Computation Time

Concerning computation time, the comparison between the proposed method and the other methods is given in Tables 6 and 7. The computation time changed in fact slightly for different sampling ratios. Therefore, we have chosen to show in Tables 6 and 7 the average time over different sampling ratios for assessing the computation time. It is seen that the SparseMRI method is the fastest among the six methods and that the proposed method takes the longest time.

To assess computational cost in case larger image sizes are involved, experiments on an image size of 512 × 512 were also performed and the computation time is compared in Tables 8 and 9. From Tables 69, it can be observed that when image size increases, the computational cost rises for all methods. However, compared with other methods, the computation time for the proposed method has a relatively small increase.

4. Conclusion

In this work, by combining the wavelet sparsity, wavelet tree structured sparsity, and NLTV, a novel method to accurately reconstruct MR images from highly undersampled k-space data has been developed. Several experiments were performed by changing the test images, sampling schemes, sampling ratios, and reconstruction methods. The results showed that the proposed WaTMRI-NLTV outperforms the conventional CS-MRI algorithms with less artifacts, richer texture, higher PSNR, MSSIM, and FSIM, and lower reconstruction RLNE, especially for the case of low sampling ratio. In the future, according to image properties, exploring more accurate sparse priors for further improving reconstruction quality would be desired.

Data Availability

The test images (a) and (b) in Figure 1 can be obtained by downloading the code attached to Ref. [29].

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

This work was supported in part by the National Natural Science Foundation of China (nos. 61661010, 61701105), the Funds for Talents of Guizhou University (No. 2013(33)), the Nature Science Foundation of Guizhou Province (Qiankehe J No.20152044), the Natural Science Foundation of Heilongjiang Province of China (no. QC2017066), French ANR (under MOSIFAH ANR-13-MONU-0009-01), and the Program PHC-Cai Yuanpei 2018 N° 41400TC.