# Compressed Sensing MRI Reconstruction with Multiple Sparsity Constraints on Radial Sampling

Authors: Jianping Huang; Lihui Wang; Yuemin Zhu
Journal: Mathematical Problems in Engineering
 (2019)

---

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


---

## Body

## 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 [1–3].One of the most promising approaches to deal with this issue is compressed sensing (CS) [1, 4–6] 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 l1 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 [9–12]. 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 [13–15]. 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 [16–21], 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 [27–29] 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. Section2 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 thatx is the desired MR image and y 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 x with a constraint of y=Φx by solving the following optimization problem:(1)minx∈RnJxsubject  toy=Φxwhere Jx is a regularizing function that intends to find a perfect sparse representation of the image x. 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:(2)x^=arg⁡minx⁡12Φx-y22+αxNLTV+βΨx1+∑g∈GΨxg2where the sensing matrix Φ is a partial Fourier transform, which can be expressed by Fu=P·F with F being the Fourier transform and P being the common undersampling pattern (mask), Ψ is a wavelet transform, G represents all the parent-child groups that encourage the tree sparsity and g is one of such groups, and xNLTV is the NLTV-based constraint defined as follows [13, 32]:(3)xNLTV=∑u∑vxu-xv2wu,vwith x(u) and x(v) indicating the image values at u and v and with w(u,v) representing the weight function formulated as(4)wu,v=1Zxexp-qxu-qxv222h2ifu-v2≤δ20otherwiseIn Eq. (4), qx(u) and qx(v) represent a small patch centering at the coordinates u and v, respectively, Zx signifies a normalization factor,h denotes a scale parameter that is related to the patch size m 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 variablez is introduced to solve Eq. (2):(5)z=BΨxwhere B is a binary matrix for wavelet coefficients group configuration. Accordingly, the optimization problem in Eq. (2) can be rewritten as follows:(6)x^=arg⁡minx⁡12Φx-y22+αxNLTV+βΨx1+∑i=1szgi2+λ2z-BΨx22where gi is the ith wavelet tree group and s is the total number of groups. The solution of Eq. (6) can be split into two subproblems:z subproblem and x subproblem, in which,z subproblem can be expressed as follows [29]:(7)zgi=arg⁡mingi⁡βzgi2+λ2zgi-BΦxgi22,i=1,2,…,sIt has a closed form solution by soft-thresholding:(8)zgi=maxri2-βλ,0riri2,i=1,2,…,swhere ri=BΦxgi. For convenience,z is represented as(9)z=shrinkgroupBΨx,βλThex subproblem can be expressed as(10)x^=arg⁡minx⁡12Φx-y22+αxNLTV+βΨx1+λ2z-BΨx22Letfx=(1/2)Φx-y22+(λ/2)z-BΨx22, which is a convex and smooth function with Lipschitz constant Lf, and g1x=αxNLTV,g2x=βΨx1, 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 l1-norm. Each subproblem is actually a convex function which can be solved by a proximal mapping operation proxgx [33]:(11)proxρgx=arg⁡minu⁡gu+12ρu-x22where ρ is a positive scalar.The subsequent algorithm is described as shown in Algorithm1.Algorithm 1
Input: ρ = 1 / L f , t 1 = 1 , x 0 = r 1 , α , β , λ For k=1 to K do z = s h r i n k g r o u p G Ψ x k - 1 , β λ x g = r k - ρ ∇ f r k x 1 = p r o x ρ 2 α x N L T V x g x 2 = p r o x ρ 2 β Ψ x 1 x g x k = x 1 + x 2 2 t k + 1 = 1 + 1 + 4 t k 2 2 r k + 1 = x k + t k - 1 t k + 1 x k - x k - 1 End for.In Algorithm1,∇frk=ΦTy-Φx+λψTBT(Bψrk-z), with ΦTand ψT representing, respectively, the inverse partial Fourier transform and inverse wavelet transform, and K 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 l2 norm error (RLNE) [36] were used. Denote that x and x^ 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 M×N. Accordingly, PSNR is defined as(12)PSNR=20log⁡10MAX1/MN∑i=1M∑j=1Nxi,j-x^i,j2Let p and q represent one block of the images x and x^, respectively, the local means of which are expressed by μp and μq with σp, σq, and σpq indicating the standard deviations and cross-covariance of images blocks p and q. Based on such notations, the structure similarity (SSIM) is calculated as follows:(13)SSIMp,q=2μpμq+C12σpq+C2μp2+μq2+C1σp2+σq2+C2where C1 and C2 are the tunable constants. If the total number of image blocks is Nb, the MSSIM can be formulated as(14)MSSIMx,x^=1Nb∑i=1NbSSIMpi,qiBoth PSNR and MSSIM are proportional to image quality. The better the image quality is, the higher the PSNR and MSSIM are.


## 2.1. CS-MRI Reconstruction Framework

Assume thatx is the desired MR image and y 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 x with a constraint of y=Φx by solving the following optimization problem:(1)minx∈RnJxsubject  toy=Φxwhere Jx is a regularizing function that intends to find a perfect sparse representation of the image x. 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:(2)x^=arg⁡minx⁡12Φx-y22+αxNLTV+βΨx1+∑g∈GΨxg2where the sensing matrix Φ is a partial Fourier transform, which can be expressed by Fu=P·F with F being the Fourier transform and P being the common undersampling pattern (mask), Ψ is a wavelet transform, G represents all the parent-child groups that encourage the tree sparsity and g is one of such groups, and xNLTV is the NLTV-based constraint defined as follows [13, 32]:(3)xNLTV=∑u∑vxu-xv2wu,vwith x(u) and x(v) indicating the image values at u and v and with w(u,v) representing the weight function formulated as(4)wu,v=1Zxexp-qxu-qxv222h2ifu-v2≤δ20otherwiseIn Eq. (4), qx(u) and qx(v) represent a small patch centering at the coordinates u and v, respectively, Zx signifies a normalization factor,h denotes a scale parameter that is related to the patch size m 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 variablez is introduced to solve Eq. (2):(5)z=BΨxwhere B is a binary matrix for wavelet coefficients group configuration. Accordingly, the optimization problem in Eq. (2) can be rewritten as follows:(6)x^=arg⁡minx⁡12Φx-y22+αxNLTV+βΨx1+∑i=1szgi2+λ2z-BΨx22where gi is the ith wavelet tree group and s is the total number of groups. The solution of Eq. (6) can be split into two subproblems:z subproblem and x subproblem, in which,z subproblem can be expressed as follows [29]:(7)zgi=arg⁡mingi⁡βzgi2+λ2zgi-BΦxgi22,i=1,2,…,sIt has a closed form solution by soft-thresholding:(8)zgi=maxri2-βλ,0riri2,i=1,2,…,swhere ri=BΦxgi. For convenience,z is represented as(9)z=shrinkgroupBΨx,βλThex subproblem can be expressed as(10)x^=arg⁡minx⁡12Φx-y22+αxNLTV+βΨx1+λ2z-BΨx22Letfx=(1/2)Φx-y22+(λ/2)z-BΨx22, which is a convex and smooth function with Lipschitz constant Lf, and g1x=αxNLTV,g2x=βΨx1, 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 l1-norm. Each subproblem is actually a convex function which can be solved by a proximal mapping operation proxgx [33]:(11)proxρgx=arg⁡minu⁡gu+12ρu-x22where ρ is a positive scalar.The subsequent algorithm is described as shown in Algorithm1.Algorithm 1
Input: ρ = 1 / L f , t 1 = 1 , x 0 = r 1 , α , β , λ For k=1 to K do z = s h r i n k g r o u p G Ψ x k - 1 , β λ x g = r k - ρ ∇ f r k x 1 = p r o x ρ 2 α x N L T V x g x 2 = p r o x ρ 2 β Ψ x 1 x g x k = x 1 + x 2 2 t k + 1 = 1 + 1 + 4 t k 2 2 r k + 1 = x k + t k - 1 t k + 1 x k - x k - 1 End for.In Algorithm1,∇frk=ΦTy-Φx+λψTBT(Bψrk-z), with ΦTand ψT representing, respectively, the inverse partial Fourier transform and inverse wavelet transform, and K 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 l2 norm error (RLNE) [36] were used. Denote that x and x^ 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 M×N. Accordingly, PSNR is defined as(12)PSNR=20log⁡10MAX1/MN∑i=1M∑j=1Nxi,j-x^i,j2Let p and q represent one block of the images x and x^, respectively, the local means of which are expressed by μp and μq with σp, σq, and σpq indicating the standard deviations and cross-covariance of images blocks p and q. Based on such notations, the structure similarity (SSIM) is calculated as follows:(13)SSIMp,q=2μpμq+C12σpq+C2μp2+μq2+C1σp2+σq2+C2where C1 and C2 are the tunable constants. If the total number of image blocks is Nb, the MSSIM can be formulated as(14)MSSIMx,x^=1Nb∑i=1NbSSIMpi,qiBoth 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 Figures1(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.Figure 1
The test images and radial undersampling masks.
(a)
MRI-Chest (b)
MRI-Renal_Arteries (c)
Radial uniform (d)
Radial golden (e)
Radial randomizedThe observation measurementy was modeled as y=Φx+ε, where ε represents complex Gaussian white noise with standard deviation σn. The associated input SNR (ISNR) [40] is defined as ISNR=20log10σx/σn, with σx 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.Table 1
Parameter settings for different reconstruction methods.
Parameters SparseMRI RecPF FCSA FCSA-NLTV WaTMRI Proposed α   0.001 0.001 0.001 0.001 β   0.035 0.035 0.035 0.035 λ   - -  0.007 0.007 h   - - 0.0625 - - 0.0625 δ   - - 5 × 5 - - 5 × 5 m   - - 11 × 11 - - 11 × 11


### 3.2. Visual Comparisons

Figures2 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.Figure 2
Results of reconstruction on chest MR image using different methods and different sampling masks with 20% sampling ratio. (a), (c), and (e) represent the reconstructed chest MR images with different methods, and (b), (d), and (f) represent the corresponding error images.
(a) (b) (c) (d) (e) (f)Figure 3
Results of reconstruction on renal arteries MR image using different methods and different sampling masks with 20% sampling ratio. (a), (c), and (e) represent the reconstructed renal arteries MR images with different methods, and (b), (d) and (f) represent the corresponding error images.
(a) (b) (c) (d) (e) (f)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 Tables2–5, 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.Table 2
PSNR of images reconstructed with 20% sampling (ISNR=20).
Dataset Sampling mask SparseMRI RecPF FCSA FCSA-NLTV WaTMRI Proposed Chest image Radial uniform 24.54 24.39 25.12 25.05 25.16 27.57 Radial golden 24.31 24.28 24.77 24.72 24.90 26.78 Radial randomized 22.66 22.65 23.44 23.41 23.48 25.77 Renal_Arteries Radial uniform 26.96 26.82 28.31 28.69 28.95 30.88 Radial golden 26.77 26.58 28.09 28.41 28.64 30.35 Radial randomized 25.03 24.86 26.26 26.66 26.84 28.27Table 3
MSSIM of images reconstructed with 20% sampling (ISNR=20).
Dataset Sampling mask SparseMRI RecPF FCSA FCSA-NLTV WaTMRI Proposed Chest image Radial uniform 0.64 0.62 0.65 0.64 0.65 0.76 Radial golden 0.62 0.60 0.62 0.62 0.63 0.75 Radial randomized 0.58 0.57 0.59 0.58 0.58 0.70 Renal_Arteries Radial uniform 0.60 0.54 0.78 0.79 0.80 0.84 Radial golden 0.58 0.52 0.77 0.77 0.79 0.82 Radial randomized 0.56 0.51 0.73 0.74 0.75 0.79Table 4
FSIM of images reconstructed with 20% sampling (ISNR=20).
Dataset Sampling mask SparseMRI RecPF FCSA FCSA-NLTV WaTMRI Proposed Chest image Radial uniform 0.78 0.79 0.81 0.81 0.81 0.85 Radial golden 0.77 0.79 0.80 0.80 0.80 0.84 Radial randomized 0.76 0.77 0.78 0.78 0.78 0.82 Renal_Arteries Radial uniform 0.83 0.82 0.85 0.86 0.87 0.89 Radial golden 0.81 0.81 0.84 0.85 0.85 0.88 Radial randomized 0.79 0.79 0.82 0.83 0.83 0.85Table 5
RLNE of images reconstructed with 20% sampling (ISNR=20).
Dataset Sampling mask SparseMRI RecPF FCSA FCSA-NLTV WaTMRI Proposed Chest image Radial uniform 0.18 0.19 0.17 0.17 0.17 0.13 Radial golden 0.19 0.19 0.18 0.18 0.18 0.14 Radial randomized 0.23 0.23 0.21 0.21 0.21 0.16 Renal_Arteries Radial uniform 0.16 0.17 0.14 0.13 0.13 0.10 Radial golden 0.17 0.17 0.14 0.14 0.13 0.11 Radial randomized 0.20 0.21 0.18 0.17 0.16 0.14


### 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 (10%~50% or 0.1~0.5) for the MR images reconstructed with different methods and undersampling schemes are shown in Figures 4–7. 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.Figure 4
Curve of PSNRs vs sampling ratios for both chest and renal arteries MR images reconstructed with different methods and sampling schemes.Figure 5
Curve of MSSIMs vs sampling ratios for both chest and renal arteries MR images reconstructed with different methods and sampling schemes.Figure 6
Curve of FSIMs vs sampling ratios for both chest and renal arteries MR images reconstructed with different methods and sampling schemes.Figure 7
Curve of RLNEs vs sampling ratios for both chest and renal arteries MR images reconstructed with different methods and sampling schemes.


### 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 Figure8. 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.Figure 8
Reconstruction performance of the proposed method with various ISNRs and sampling ratios (Rs) on different undersampling schemes for both chest and renal arteries MR images. (a) PSNR vs SNR. (b) MSSIM vs SNR. (c) FSIM vs SNR. (d) RLNE vs SNR.
(a) (b)  (c) (d)


### 3.6. Computation Time

Concerning computation time, the comparison between the proposed method and the other methods is given in Tables6 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.Table 6
Comparison of computation time (in seconds) of different methods on chest MR image (size 256 × 256) with different radial sampling schemes.
sampling scheme method SparseMRI RecPF FCSA FCSA-NLTV WaTMRI Proposed radial uniform 2.86 4.17 4.72 4.71 4.56 113.74 radial golden 2.92 4.39 4.85 4.73 4.81 115.77 radial randomized 2.83 4.14 4.30 4.25 4.34 104.95Table 7
Comparison of computation time (in seconds) of different methods on renal arteries MR image (size 256×256) with different radial sampling schemes.
sampling scheme method SparseMRI RecPF FCSA FCSA-NLTV WaTMRI Proposed radial uniform 3.02 4.33 2.74 5.21 5.41 116.74 radial golden 2.74 3.97 2.82 4.63 4.90 103.34 radial randomized 2.82 4.19 2.65 5.15 5.11 112.29To 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 Tables8 and 9. From Tables 6–9, 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.Table 8
Comparison of computation time (in seconds) of different methods on chest MR image (size 512×512) with different radial sampling schemes.
sampling scheme method SparseMRI RecPF FCSA FCSA-NLTV WaTMRI Proposed radial uniform 44.34 18.00 62.91 62.45 70.41 210.96 radial golden 44.57 18.06 65.82 64.27 69.82 210.49 radial randomized 45.38 17.94 74.83 63.22 73.45 212.51Table 9
Comparison of computation time (in seconds) of different methods on renal arteries MR image (size 512×512) with different radial sampling schemes.
sampling scheme method SparseMRI RecPF FCSA FCSA-NLTV WaTMRI Proposed radial uniform 45.36 17.49 34.36 94.62 104.90 212.48 radial golden 46.23 18.26 36.56 96.11 104.29 213.36 radial randomized 44.92 17.97 35.36 93.10 97.21 215.43


## 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 Figures1(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.Figure 1
The test images and radial undersampling masks.
(a)
MRI-Chest (b)
MRI-Renal_Arteries (c)
Radial uniform (d)
Radial golden (e)
Radial randomizedThe observation measurementy was modeled as y=Φx+ε, where ε represents complex Gaussian white noise with standard deviation σn. The associated input SNR (ISNR) [40] is defined as ISNR=20log10σx/σn, with σx 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.Table 1
Parameter settings for different reconstruction methods.
Parameters SparseMRI RecPF FCSA FCSA-NLTV WaTMRI Proposed α   0.001 0.001 0.001 0.001 β   0.035 0.035 0.035 0.035 λ   - -  0.007 0.007 h   - - 0.0625 - - 0.0625 δ   - - 5 × 5 - - 5 × 5 m   - - 11 × 11 - - 11 × 11


## 3.2. Visual Comparisons

Figures2 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.Figure 2
Results of reconstruction on chest MR image using different methods and different sampling masks with 20% sampling ratio. (a), (c), and (e) represent the reconstructed chest MR images with different methods, and (b), (d), and (f) represent the corresponding error images.
(a) (b) (c) (d) (e) (f)Figure 3
Results of reconstruction on renal arteries MR image using different methods and different sampling masks with 20% sampling ratio. (a), (c), and (e) represent the reconstructed renal arteries MR images with different methods, and (b), (d) and (f) represent the corresponding error images.
(a) (b) (c) (d) (e) (f)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 Tables2–5, 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.Table 2
PSNR of images reconstructed with 20% sampling (ISNR=20).
Dataset Sampling mask SparseMRI RecPF FCSA FCSA-NLTV WaTMRI Proposed Chest image Radial uniform 24.54 24.39 25.12 25.05 25.16 27.57 Radial golden 24.31 24.28 24.77 24.72 24.90 26.78 Radial randomized 22.66 22.65 23.44 23.41 23.48 25.77 Renal_Arteries Radial uniform 26.96 26.82 28.31 28.69 28.95 30.88 Radial golden 26.77 26.58 28.09 28.41 28.64 30.35 Radial randomized 25.03 24.86 26.26 26.66 26.84 28.27Table 3
MSSIM of images reconstructed with 20% sampling (ISNR=20).
Dataset Sampling mask SparseMRI RecPF FCSA FCSA-NLTV WaTMRI Proposed Chest image Radial uniform 0.64 0.62 0.65 0.64 0.65 0.76 Radial golden 0.62 0.60 0.62 0.62 0.63 0.75 Radial randomized 0.58 0.57 0.59 0.58 0.58 0.70 Renal_Arteries Radial uniform 0.60 0.54 0.78 0.79 0.80 0.84 Radial golden 0.58 0.52 0.77 0.77 0.79 0.82 Radial randomized 0.56 0.51 0.73 0.74 0.75 0.79Table 4
FSIM of images reconstructed with 20% sampling (ISNR=20).
Dataset Sampling mask SparseMRI RecPF FCSA FCSA-NLTV WaTMRI Proposed Chest image Radial uniform 0.78 0.79 0.81 0.81 0.81 0.85 Radial golden 0.77 0.79 0.80 0.80 0.80 0.84 Radial randomized 0.76 0.77 0.78 0.78 0.78 0.82 Renal_Arteries Radial uniform 0.83 0.82 0.85 0.86 0.87 0.89 Radial golden 0.81 0.81 0.84 0.85 0.85 0.88 Radial randomized 0.79 0.79 0.82 0.83 0.83 0.85Table 5
RLNE of images reconstructed with 20% sampling (ISNR=20).
Dataset Sampling mask SparseMRI RecPF FCSA FCSA-NLTV WaTMRI Proposed Chest image Radial uniform 0.18 0.19 0.17 0.17 0.17 0.13 Radial golden 0.19 0.19 0.18 0.18 0.18 0.14 Radial randomized 0.23 0.23 0.21 0.21 0.21 0.16 Renal_Arteries Radial uniform 0.16 0.17 0.14 0.13 0.13 0.10 Radial golden 0.17 0.17 0.14 0.14 0.13 0.11 Radial randomized 0.20 0.21 0.18 0.17 0.16 0.14


## 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 (10%~50% or 0.1~0.5) for the MR images reconstructed with different methods and undersampling schemes are shown in Figures 4–7. 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.Figure 4
Curve of PSNRs vs sampling ratios for both chest and renal arteries MR images reconstructed with different methods and sampling schemes.Figure 5
Curve of MSSIMs vs sampling ratios for both chest and renal arteries MR images reconstructed with different methods and sampling schemes.Figure 6
Curve of FSIMs vs sampling ratios for both chest and renal arteries MR images reconstructed with different methods and sampling schemes.Figure 7
Curve of RLNEs vs sampling ratios for both chest and renal arteries MR images reconstructed with different methods and sampling schemes.


## 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 Figure8. 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.Figure 8
Reconstruction performance of the proposed method with various ISNRs and sampling ratios (Rs) on different undersampling schemes for both chest and renal arteries MR images. (a) PSNR vs SNR. (b) MSSIM vs SNR. (c) FSIM vs SNR. (d) RLNE vs SNR.
(a) (b)  (c) (d)


## 3.6. Computation Time

Concerning computation time, the comparison between the proposed method and the other methods is given in Tables6 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.Table 6
Comparison of computation time (in seconds) of different methods on chest MR image (size 256 × 256) with different radial sampling schemes.
sampling scheme method SparseMRI RecPF FCSA FCSA-NLTV WaTMRI Proposed radial uniform 2.86 4.17 4.72 4.71 4.56 113.74 radial golden 2.92 4.39 4.85 4.73 4.81 115.77 radial randomized 2.83 4.14 4.30 4.25 4.34 104.95Table 7
Comparison of computation time (in seconds) of different methods on renal arteries MR image (size 256×256) with different radial sampling schemes.
sampling scheme method SparseMRI RecPF FCSA FCSA-NLTV WaTMRI Proposed radial uniform 3.02 4.33 2.74 5.21 5.41 116.74 radial golden 2.74 3.97 2.82 4.63 4.90 103.34 radial randomized 2.82 4.19 2.65 5.15 5.11 112.29To 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 Tables8 and 9. From Tables 6–9, 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.Table 8
Comparison of computation time (in seconds) of different methods on chest MR image (size 512×512) with different radial sampling schemes.
sampling scheme method SparseMRI RecPF FCSA FCSA-NLTV WaTMRI Proposed radial uniform 44.34 18.00 62.91 62.45 70.41 210.96 radial golden 44.57 18.06 65.82 64.27 69.82 210.49 radial randomized 45.38 17.94 74.83 63.22 73.45 212.51Table 9
Comparison of computation time (in seconds) of different methods on renal arteries MR image (size 512×512) with different radial sampling schemes.
sampling scheme method SparseMRI RecPF FCSA FCSA-NLTV WaTMRI Proposed radial uniform 45.36 17.49 34.36 94.62 104.90 212.48 radial golden 46.23 18.26 36.56 96.11 104.29 213.36 radial randomized 44.92 17.97 35.36 93.10 97.21 215.43


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



---

*Source: 3694604-2019-02-10.xml*