# Compressed Sensing MRI Reconstruction from Highly Undersampledk-Space Data Using Nonsubsampled Shearlet Transform Sparsity Prior

Authors: Min Yuan; Bingxin Yang; Yide Ma; Jiuwen Zhang; Runpu Zhang; Caiyuan Zhang
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

Compressed sensing has shown great potential in speeding up MR imaging by undersamplingk-space data. Generally sparsity is used as a priori knowledge to improve the quality of reconstructed image. Compressed sensing MR image (CS-MRI) reconstruction methods have employed widely used sparsifying transforms such as wavelet or total variation, which are not preeminent in dealing with MR images containing distributed discontinuities and cannot provide a sufficient sparse representation and the decomposition at any direction. In this paper, we propose a novel CS-MRI reconstruction method from highly undersampled k-space data using nonsubsampled shearlet transform (NSST) sparsity prior. In particular, we have implemented a flexible decomposition with an arbitrary even number of directional subbands at each level using NSST for MR images. The highly directional sensitivity of NSST and its optimal approximation properties lead to improvement in CS-MRI reconstruction applications. The experimental results demonstrate that the proposed method results in the high quality reconstruction, which is highly effective at preserving the intrinsic anisotropic features of MRI meanwhile suppressing the artifacts and added noise. The objective evaluation indices outperform all compared CS-MRI methods. In summary, NSST with even number directional decomposition is very competitive in CS-MRI applications as sparsity prior in terms of performance and computational efficiency.


---

## Body

## 1. Introduction

Magnetic resonance imaging (MRI) is a widely used noninvasive imaging modality for clinical diagnosis. However, relatively slow imaging speed of MRI remains a great challenge for clinical application. An effective way to speed up MR imaging isk-space undersampling; in the course, only a small subset of k-space is measured. However, undersampling often violates the Nyquist sampling criterion, resulting in aliasing artifacts in linear reconstruction, inadequate image resolution, and excessive Gibbs ringing artifacts.Compressed sensing (CS) [1–3] proposed by Candès et al. and Donoho is a novel signal acquisition and compression theory. As a promising method, CS has been introduced to MR image reconstruction, so-called compressed sensing MRI (CS-MRI) [4–6]. CS-MRI allows high quality reconstruction from undersampled k-space data by solving a constrained minimization problem using nonlinear optimization algorithm by enforcing the sparsity of images in a certain predefined sparsifying transform, such as the traditional 2D separable wavelet transform [4], total variation (TV) [4, 7, 8], contourlet [9, 10], sharp frequency localization contourlet (SFLCT) [11, 12], dual-tree complex wavelet transform (DT-CWT) [13, 14], and complex double-density dual-tree DWT [15]. The structured sparsity such as Gaussian scale mixture (GSM) model [16, 17] and wavelet tree sparsity [18, 19] for exploiting the dependencies between wavelet coefficients has been introduced to CS-MRI reconstruction. In terms of the restricted isometry property (RIP) condition in CS, incorporating a prior knowledge to enhance the sparsity into the reconstruction can reduce the reconstruction error and improve the reconstruction of details effectively. In the past several years these analytical sparsifying transforms have been successfully applied to CS-MRI and demonstrate high quality reconstructions from undersampled data. Some of these transforms have been combined to further improve the reconstruction [20–23]. Dictionary learning enables adaptive sparser representation of MR images than the general sparsifying transforms, which has been applied in CS-MRI [24].The quality of reconstructed images largely depends on the performance of exploiting the sparsity prior in CS-MRI, which is key to accurate CS reconstruction. Therefore an outstanding sparsifying transform suitable for representing MR images should be adopted as sparsity prior in CS-MRI reconstruction. The motivation of our study comes from the morphology that MR images consist of different components (point-like and curve-like features) in various orientations, which cannot be sparsely represented by existing sparse representation sufficiently, and yet dictionary learning is at the expense of sacrificing time. In this paper we proposed a novel CS-MRI image reconstruction method from highly undersampledk-space data to improve the quality of reconstructed MR images by enhancing the sparsity in nonsubsampled shearlet transform (NSST) domain. With the ability to capture intrinsic geometrical features of multidimensional data efficiently and sparsely represent images containing edges optimally, the prominent sparse representation-NSST is adopted as prior knowledge for the regularization in CS-MRI. The numerical computation employs a corresponding iterative NSST thresholding algorithm to solve this inverse optimization problem.


## 2. Theory and Methods

### 2.1. Problem Formulation

The common model of data acquisition with incomplete measurements for CS-MRI is given by the following formulation:(1)y=Fux+n,where x∈Cn is the reconstructed image, y∈Cm is the acquired k-space measurement data corrupted with the noise n, and Fu is the undersampled Fourier transform operator which directly relies on the k-space undersampling scheme. Suppose that x is represented as x=Ψα, where Ψ represents the sparsity prior associated with the transform under which MR image x has a sparse representation or approximation. Then the measured data is given by y=FuΨα+n, where FuΨ represents the sensing matrix.According to CS theory, CS-MRI claims to reconstruct MR image from undersampledk-space data by enforcing the image sparsity [4]. That is, x can be accurately reconstructed from a small subset of k-space data by solving the following l0 norm minimization problem:(2)min⁡αα0s.t.y=FuΨα.However, thel0 norm is not convex and the computational complexity of the optimization is NP-hard [7]. It has been proven that, under certain condition, l0 norm problem is equivalent to l1 norm. Thus the reconstruction can be obtained by solving the following constrained convex optimization:(3)min⁡αα1s.t.y-FuΨα2<ε,where ε is a statistic describing the magnitude of the error, defined as the noise variance or the maximum allowable error in the approximation. Minimizing the objective function α1 promotes the sparsity of the images. The constraint y-FuΨα2<ε enforces the fidelity of the reconstruction to the measured k-space data. Besides, constrained l1 norm convex optimization problem in (3) can be written in unconstrained Lagrangian form:(4)min⁡α12y-FuΨα22+λα1.The second term of (4) is a regularization term that represents prior sparse information of original images. λ is a regularization parameter governing the tradeoff between the data fidelity and its sparsity.


### 2.2. Nonsubsampled Shearlet Transform

The shearlet transform, introduced by the authors and their collaborators in [25, 26], is a very recent sibling in the family of geometric image representations. The shearlet representation originally derived from the framework of affine systems with composite dilations [25–28]. The shearlet frame elements associated with shearlet transform are defined at various scales. The shearlet transform is a multiscale directional transform which is especially adapted to localize distributed discontinuities such as edges. Unlike conventional multiresolution analysis tool, the shearlet representation is theoretically optimal in representing directional and anisotropic features in images and has the ability to accurately and efficiently capture the geometric information of multidimensional data at various scales. As a result, this approach provides optimal approximation properties for 2D images.The shearlet representation forms an affine-like system and has a simpler mathematical construction. In fact, the elements of this system form a Parseval frame and are generated by applying dilations, shear transformations, and translations to a single well-localized window function. Shearlet transform extends naturally to higher dimensions and can be associated with a multiresolution analysis [28]. In addition, this approach has a fast algorithmic implementation and is very competitive for CS-MRI reconstruction.The discrete shearlet transform, obtained by discretizing the corresponding continuous transform, has different form in the numerical implementation. Both a frequency and time-domain based implementation of the discrete shearlet transform have been developed [25]. The features of each particular representation will have various advantages for specific applications. Refer to [25–28] for more details about the mathematical framework and the implement of discrete shearlet transform. The reason for using nonsubsampled shearlet transform (NSST) as sparsity prior for CS-MRI reconstruction will be discussed in Section 2.3.1 in detail.


### 2.3. Proposed Method

The aforementioned sparsifying transforms integrated with CS-MRI are acknowledged to have a limited capability in representing these directional information and important details. The properties of each particular representation will have different advantages for specific application. They may not provide sufficient sparse representation for sharp spatial gradients and intrinsic geometrical features contained in MR images. The contourlet transform is an efficient directional multiresolution image representation. However, nonideal filters used in the original contourlet result in substantial amount of aliasing components and blurring artifacts in representing smooth boundaries. To solve this problem, Lu and Do [11] propose a new construction of the contourlet, called sharp frequency localization contourlet (SFLCT). SFLCT has been integrated with CS-MRI reconstruction [12]. Since the combination of LP and 2D directional filter bank (DFB) makes the aliasing problem serious, new multiscale pyramid with different set of low pass and high pass filters for the first level and all other levels is employed. SFLCT alleviates the nonlocalization problem with the same redundancy of the original contourlet. Though SFLCT is sharply localized in the frequency domain, the downsampling of LP and DFBs stages makes it lack shift invariance, which could easily produce Gibbs ringing artifacts around the singularities, for example, edges. Although the contourlet basis is anisotropic, the directional subbands can be decomposed only at 2M directions in terms of directional selectivity, in which M denotes the cascade layers of DFBs.To overcome these limitations, on account of the MR images consisting of curve singularities and anisotropic directional features, a more appropriate sparsity prior with highly directional sensitivity and anisotropy should be applied to reconstruction. So we employ a special form of shearlet transform—NSST—as sparsity prior in CS-MRI reconstruction.


#### 2.3.1. Major Advantages of NSST Sparsity Prior

Taking account of measurement noise ink-space and the problems of aliasing and shift-variance caused by decimation, in our implementation for CS-MRI, we adopt the particular form of the time-domain based shearlet transform. This will be simply referred to as the nonsubsampled shearlet transform (NSST), which is to use the nonsubsampled Laplacian pyramid (LP) transform with several different combinations of the shearing filters. The idea of using NSST will be represented one by one below.(i) Essentially Optimal Sparsity and Approximation Property. The family of shearlet function forms a tight frame of L2(R2), and, thus, an image f can be represented by f=∑j,l,k〈f,ψj,l,k〉ψj,l,k. The coefficients αj,l,k(f)=f,ψj,l,k are called shearlet coefficients of the image f. The shearlet elements form a tight frame of well-localized waveforms, at various scales and directions, and are optimally sparse in representing images with edges. They are compactly supported in the frequency domain and have fast decay in the spatial domain.To make the statement of optimal sparsity and approximation property more rigorous, it is useful to quantify the approximation performance from the point of view of approximation theory. The asymptotic convergence rate is actually the correct optimal behavior for approximating general smooth objects having discontinuities along piecewiseC2 curves. Denoting as f^M the approximation of an image f by the M largest transform coefficients in the corresponding representation, the resulting approximation error (in l2-norm square) is f-f^M22. It is very helpful to achieve the best asymptotic decay rate for this error in application. Let C2 be the space of functions that are twice continuously differentiable. If the image f is C2 everywhere away from edge curves that are piecewise C2, the best M-term asymptotic approximation error using shearlets has a decay rate of O((log⁡M)3M-2) [28], which is essentially optimal in representing 2D images f which are piecewise C2 except for discontinuities along C2 curves and greatly outperforms that of wavelet approximations only with the decay rate of O(M-1) [29]. The shearlet transform has very similar asymptotic approximation properties with the curvelets [30] and the contourlets. The corresponding argument about optimal sparsity is proved in [9, 28, 30]. The error decay rate using shearlets is close to the theoretical optimal approximation, where the error decays as M-2 [31, 32]. In particular, the shearlet transform has some advantages over the contourlet transform [25]. In this sense, the shearlet representation provides optimally sparse representation of objects with singularities along piecewise C2 edges. Thus through shearlet transform, MR image can be decomposed into various frequency regions whose supports are contained in pair of trapezoidal regions symmetric with respect to the origin and directionally oriented.For CS-MRI, two most useful features of shearlets are their ability to efficiently approximate signals containing piecewise singularities and allowing for a much less redundant sparse tight frame representation. Consequently shearlet transform can represent images sparsely better than other representations. Moreover NSST can offer shift invariance. These properties are of vital importance for CS-MRI reconstruction because these properties results in MR images that are more compressible and hence more effective to be reconstructed. We can achieve better approximation performance and better sparse representation by shearlet coefficients compared to others. The reconstruction error in CS is proportional to approximation error. By using NSST as the sparse transform for some singularities in MR images, the reconstruction error decreases faster than that of other representations from undersampledk-space data. This is the reason why we use NSST as sparsity prior to better reconstruct some crucial features than using other representations in CS-MRI.(ii) Highly Directional Sensitivity and Anisotropy. The collection of shearlets φk,ψj,l,kd:j≥0,-2j≤l≤2j-1,k∈Z2,d=0,1 is a Parseval frame (tight frame) for L2(R2) [25]. This indicates that the decomposition is invertible and the transformation is numerically well-conditioned. Details about this construction can be found in [25, 28]. The tiling of the frequency plane R^2 induced by the shearlets is illustrated in Figure 1(a).Figure 1
Frequency tiling of the shearlets and an example of basis function. (a) The tiling of the frequency plane induced by the shearlets, (b) basis function, and (c) the frequency support of a shearlet on a pair of trapezoids.
(a)
(b)
(c)Shearlets are compactly supported in the frequency domain and have fast decay in the spatial domain with time-frequency localization. The frequency support of a shearletψj,l,k satisfies parabolic scaling. Each element ψ^j,l,k is supported on a pair of trapezoids, of approximate size 22j×2j (see Figure 1(c)), oriented along lines of slope l2-j. As a consequence, shearlets exhibit highly directional sensitivity. Unlike the isotropic elements of wavelet bases, a pair of trapezoids elements of shearlet transform possesses very high directional sensitivity and anisotropy. Display of a basis function for the shearlet transform is shown in Figure 1(b). As a result, the shearlet transform requires fewer coefficients to represent curve-like feature of MR images than the other transforms do.(iii) No Restrictions on the Number of Directions for the Shearing-Highly Directional Sensitivity and Flexibility. Such elements are very efficient in representing curve-like edges and texture features. Contourlet and SFLCT can capture the geometry of images only in 2M directions which is not optimal for MR images, because the feature information of organ, tissue, or blood vessel could be in all directions. It is worth being noted that an outstanding advantage of the shearlet transform is that there are no restrictions on the number of directions for the shearing compared with other “directional wavelets” including contourlet and SFLCT. It is impossible to obtain the flexibility of directional decomposition using shearlet transform for a fan filter implementation of DFBs. In addition, there are no constraints on the size of the supports for the shearing in shearlet transform, in which the shearing filter can have smaller support size than the directional filters used in contourlet transform. Note that the only restriction on the shearing filters is that the size of the filter should be more than or equal to the maximum number of directional subbands.(iv) The Small Support Sizes of the Shearing Filters. Another benefit of this implementation is that the shearing filtering can be performed “directly” in the time-domain using a convolution. In addition, the small support sizes of the filters reduce the Gibbs-type ringing phenomenon and improve the computational efficiency of the algorithm. The shearing filters constructed by using a Meyer window are shown in Figure 2.Figure 2
Examples of the shearing filters constructed using a Meyer wavelet as the window function.(v) Shift Invariance. In the time-domain implementation, data is undecimated at a specific scale using nonsubsampled LP. Consequently NSST exhibits excellent shift invariance at the same time, which can be helpful for CS-MRI reconstruction. Although this is a highly redundant decomposition, this version will be shown to be highly effective for the purpose of eliminating the effect of acquisition noise in k-space.NSST combines the power of multiscale methods with a unique ability to capture the geometry of multidimensional data and is optimally efficient in representing images containing edges. The highly directional sensitivity of the shearlet transform and its optimal approximation properties will lead to improvements in CS-MRI applications. So NSST is a good candidate for CS-MRI as sparsity prior. Incorporating NSST sparsity prior knowledge into the process of image reconstruction, we can get the sparse NSST coefficients and then get the reconstructed MR images by solving aforementioned constrained optimization problem.


#### 2.3.2. The Improvement on the Number of Directions and Parameters in Directional Decomposition

An important advantage of the shearlet transform over the contourlet transform is that there are no restrictions on the number of directions for the shearing. In the experiments, considering anisotropic characteristics existing in MR images and the accuracy of reconstruction we have implemented the decomposition with an arbitrary even number of directional subbands (2N, N=1,2,…) at each scale using NSST, which is more flexible than 2M directions of contourlet and SFLCT. In this implementation, we have large flexibility in the choice of the number of separated directional subbands for each decomposition level, not just restricted to 23, 24, and so forth directional subbands. We could achieve the decomposition at 6, 10, and 12 directions, which cannot be achieved by SFLCT. More flexible directional representation can provide the superior decomposing directions for anisotropic features in MR images. Thus it benefits us to capture the intrinsic characteristics information of MR images at different directions using NSST. Therefore, it only needs less transform domain coefficients than the other sparsifying transforms when the curve characteristics of MRI are expressed using NSST.After considerable quantity of experiments, it demonstrates that only three-level decomposition using NSST and 12 directions in each level for sparse representation of MR images has already been enough, which can achieve better reconstruction performance outperforming four-level and2M directions decomposition using SFLCT. Finer decomposition scales and much more directions cannot result in the improvement of reconstruction performance; instead they will increase the burden of computation. In practice, we use a Meyer wavelet window to construct the directional shearing filters, as shown in Figure 2. The shearing filters of sizes 12×12, 12×12, and 12×12 from finer to coarser were used with the number of shearing directions chosen to be 12, 12, and 12. NSST (12, 12, and 12) indicates that a Meyer-based shearing filter of size 12 with 12 directions was applied to the first, second, and third decomposition level.Figure3 illustrates the three-level shearlet decomposition of a T2 brain MR image with Ewing’s sarcoma from the Harvard University site (http://www.med.harvard.edu/AANLIB/home.html). The first level decomposition generates 6 directional subbands, the second level decomposition generates 8 directional subbands, and the third level decomposition generates 12 directional subbands.Figure 3
Decomposition coefficient of a T2 brain MR image with sarcoma based on NSST, of 3 resolution scales, 6, 8, and 12 directions for every scale, respectively. (a) T2 brain MR image with sarcoma. (b) Lowpass subband, the approximate shearlet coefficients. (c)~(e) The detail shearlet coefficients of decomposition level 1, level 2, and level 3 in every direction, respectively, with the full colormap scaling for better presentation.


#### 2.3.3. Compressed Sensing MRI Reconstruction Using NSST Sparsity Prior

The reconstruction problem from undersamplingk-space data actually belongs to an underdetermined system of linear equations. It is very helpful to incorporate a priori sparse information into the nonlinear reconstruction for further improvement of reconstruction performance. In this paper, we specially employ NSST with even number directional decomposition as sparsity prior in CS-MRI applications. For the reconstruction problem, we adopt the convex l1-norm regularization to promote sparsity in NSST domain. Different numerical reconstruction algorithms to solve (4) optimization problem [33–38] have been proposed in the past few decades. Among the existing algorithms, the iterative shrinkage thresholding algorithm (ISTA) is straightforward, often-used, and effective for solving linear inverse problems. The reason is that ISTA with good robustness is easy to implement and can be integrated easily with existing sparsifying transforms in practical application. These properties are significant for CS-MRI. Considering the effectiveness and the convenience of ISTA, we adopt improved iterative shearlet thresholding to solve CS-MRI.The interpretation and convergence of iterative thresholding algorithm for solving constraintl1 norm optimization were discussed in the literature [33, 35]. The core idea of iterative thresholding is that the objective function is guaranteed to decrease on each iteration of the algorithm. This algorithm iteratively performs soft thresholding to decrease the l1 norm of the coefficients α and a gradient descent to decrease the value of y-FuΨα2. Equation (4) can be simply solved by iterative shearlet thresholding:(5)αk+1=Sθk(1cFuΨH(y-FuΨαk)+αk)=Sθk1cFuΨHr+αk,where (FuΨ)H denotes the conjugate transpose operator of FuΨ and Ψ and ΨH denote the shearlet inverse and forward transform, respectively. The choice of c should satisfy c>(FuΨ)H(FuΨ)2=λmax⁡((FuΨ)H(FuΨ)) (the maximal eigenvalue of the matrix (FuΨ)H(FuΨ)). r=y-FuΨαk is the residual in k-space. The thresholding function can be used in the coefficient iteration. Here Sθk(α) is a soft thresholding operator to shrink each entry αi of vector α according to the following:(6)Sθk(αi)=0ifαi≤θkαi-αiαiθkotherwise.Since α is a complex vector for MR images, a complex thresholding operator should be defined. Daubechies et al. [33] extended the definition of the operator and derived a complex thresholding operator, which is defined as Sθk(αiejω)=Sθk(αi)ejω, ω∈[0,2π).In order to eliminate the effect of the magnitude of MR images on stop criteria, a relative error tolerance [12] Rk=y-FuΨα^2/y2 is adopted to replace absolute stop criteria related to ε of (3). The decreasing factor ρ (0<ρ<1) is adopted to decrease the threshold in each iteration. The computational parameters η and ρ in the algorithm are constants, and we set them to be the same in all the experiments. The selecting problem of the value of η and ρ was investigated in the paper [12].From empirical analysis,η=10-6 and ρ=0.8 are an appropriate choice to acquire good performance for CS-MRI. The improved CS-MRI reconstruction algorithm from undersampled k-space measurements using iterative NSST soft thresholding is summarized in Algorithm 1.Algorithm 1:NSST_IST CS-MRI.
(
1) Require:
y∈Cm, Ψ∈Cn×d, η, ρ,(0<ρ<1).
(2) Set θ0=max⁡((FuΨ)Hr0), α0=[(0,…,0)1×d]T, x0=Ψα0, r0=y
(3) Repeat
(4) Update coefficient vector according to αk+1=αk+Sθk((FuΨ)Hrk)
(5) Computer residual rk+1=y-FuΨαk+1
(6) Update the threshold θk+1=ρθk
(7) k←k+1
(8) until  Rk=y-FuΨαk2/y2≤η
(9) return  α^,x^=Ψα^.


## 2.1. Problem Formulation

The common model of data acquisition with incomplete measurements for CS-MRI is given by the following formulation:(1)y=Fux+n,where x∈Cn is the reconstructed image, y∈Cm is the acquired k-space measurement data corrupted with the noise n, and Fu is the undersampled Fourier transform operator which directly relies on the k-space undersampling scheme. Suppose that x is represented as x=Ψα, where Ψ represents the sparsity prior associated with the transform under which MR image x has a sparse representation or approximation. Then the measured data is given by y=FuΨα+n, where FuΨ represents the sensing matrix.According to CS theory, CS-MRI claims to reconstruct MR image from undersampledk-space data by enforcing the image sparsity [4]. That is, x can be accurately reconstructed from a small subset of k-space data by solving the following l0 norm minimization problem:(2)min⁡αα0s.t.y=FuΨα.However, thel0 norm is not convex and the computational complexity of the optimization is NP-hard [7]. It has been proven that, under certain condition, l0 norm problem is equivalent to l1 norm. Thus the reconstruction can be obtained by solving the following constrained convex optimization:(3)min⁡αα1s.t.y-FuΨα2<ε,where ε is a statistic describing the magnitude of the error, defined as the noise variance or the maximum allowable error in the approximation. Minimizing the objective function α1 promotes the sparsity of the images. The constraint y-FuΨα2<ε enforces the fidelity of the reconstruction to the measured k-space data. Besides, constrained l1 norm convex optimization problem in (3) can be written in unconstrained Lagrangian form:(4)min⁡α12y-FuΨα22+λα1.The second term of (4) is a regularization term that represents prior sparse information of original images. λ is a regularization parameter governing the tradeoff between the data fidelity and its sparsity.


## 2.2. Nonsubsampled Shearlet Transform

The shearlet transform, introduced by the authors and their collaborators in [25, 26], is a very recent sibling in the family of geometric image representations. The shearlet representation originally derived from the framework of affine systems with composite dilations [25–28]. The shearlet frame elements associated with shearlet transform are defined at various scales. The shearlet transform is a multiscale directional transform which is especially adapted to localize distributed discontinuities such as edges. Unlike conventional multiresolution analysis tool, the shearlet representation is theoretically optimal in representing directional and anisotropic features in images and has the ability to accurately and efficiently capture the geometric information of multidimensional data at various scales. As a result, this approach provides optimal approximation properties for 2D images.The shearlet representation forms an affine-like system and has a simpler mathematical construction. In fact, the elements of this system form a Parseval frame and are generated by applying dilations, shear transformations, and translations to a single well-localized window function. Shearlet transform extends naturally to higher dimensions and can be associated with a multiresolution analysis [28]. In addition, this approach has a fast algorithmic implementation and is very competitive for CS-MRI reconstruction.The discrete shearlet transform, obtained by discretizing the corresponding continuous transform, has different form in the numerical implementation. Both a frequency and time-domain based implementation of the discrete shearlet transform have been developed [25]. The features of each particular representation will have various advantages for specific applications. Refer to [25–28] for more details about the mathematical framework and the implement of discrete shearlet transform. The reason for using nonsubsampled shearlet transform (NSST) as sparsity prior for CS-MRI reconstruction will be discussed in Section 2.3.1 in detail.


## 2.3. Proposed Method

The aforementioned sparsifying transforms integrated with CS-MRI are acknowledged to have a limited capability in representing these directional information and important details. The properties of each particular representation will have different advantages for specific application. They may not provide sufficient sparse representation for sharp spatial gradients and intrinsic geometrical features contained in MR images. The contourlet transform is an efficient directional multiresolution image representation. However, nonideal filters used in the original contourlet result in substantial amount of aliasing components and blurring artifacts in representing smooth boundaries. To solve this problem, Lu and Do [11] propose a new construction of the contourlet, called sharp frequency localization contourlet (SFLCT). SFLCT has been integrated with CS-MRI reconstruction [12]. Since the combination of LP and 2D directional filter bank (DFB) makes the aliasing problem serious, new multiscale pyramid with different set of low pass and high pass filters for the first level and all other levels is employed. SFLCT alleviates the nonlocalization problem with the same redundancy of the original contourlet. Though SFLCT is sharply localized in the frequency domain, the downsampling of LP and DFBs stages makes it lack shift invariance, which could easily produce Gibbs ringing artifacts around the singularities, for example, edges. Although the contourlet basis is anisotropic, the directional subbands can be decomposed only at 2M directions in terms of directional selectivity, in which M denotes the cascade layers of DFBs.To overcome these limitations, on account of the MR images consisting of curve singularities and anisotropic directional features, a more appropriate sparsity prior with highly directional sensitivity and anisotropy should be applied to reconstruction. So we employ a special form of shearlet transform—NSST—as sparsity prior in CS-MRI reconstruction.


### 2.3.1. Major Advantages of NSST Sparsity Prior

Taking account of measurement noise ink-space and the problems of aliasing and shift-variance caused by decimation, in our implementation for CS-MRI, we adopt the particular form of the time-domain based shearlet transform. This will be simply referred to as the nonsubsampled shearlet transform (NSST), which is to use the nonsubsampled Laplacian pyramid (LP) transform with several different combinations of the shearing filters. The idea of using NSST will be represented one by one below.(i) Essentially Optimal Sparsity and Approximation Property. The family of shearlet function forms a tight frame of L2(R2), and, thus, an image f can be represented by f=∑j,l,k〈f,ψj,l,k〉ψj,l,k. The coefficients αj,l,k(f)=f,ψj,l,k are called shearlet coefficients of the image f. The shearlet elements form a tight frame of well-localized waveforms, at various scales and directions, and are optimally sparse in representing images with edges. They are compactly supported in the frequency domain and have fast decay in the spatial domain.To make the statement of optimal sparsity and approximation property more rigorous, it is useful to quantify the approximation performance from the point of view of approximation theory. The asymptotic convergence rate is actually the correct optimal behavior for approximating general smooth objects having discontinuities along piecewiseC2 curves. Denoting as f^M the approximation of an image f by the M largest transform coefficients in the corresponding representation, the resulting approximation error (in l2-norm square) is f-f^M22. It is very helpful to achieve the best asymptotic decay rate for this error in application. Let C2 be the space of functions that are twice continuously differentiable. If the image f is C2 everywhere away from edge curves that are piecewise C2, the best M-term asymptotic approximation error using shearlets has a decay rate of O((log⁡M)3M-2) [28], which is essentially optimal in representing 2D images f which are piecewise C2 except for discontinuities along C2 curves and greatly outperforms that of wavelet approximations only with the decay rate of O(M-1) [29]. The shearlet transform has very similar asymptotic approximation properties with the curvelets [30] and the contourlets. The corresponding argument about optimal sparsity is proved in [9, 28, 30]. The error decay rate using shearlets is close to the theoretical optimal approximation, where the error decays as M-2 [31, 32]. In particular, the shearlet transform has some advantages over the contourlet transform [25]. In this sense, the shearlet representation provides optimally sparse representation of objects with singularities along piecewise C2 edges. Thus through shearlet transform, MR image can be decomposed into various frequency regions whose supports are contained in pair of trapezoidal regions symmetric with respect to the origin and directionally oriented.For CS-MRI, two most useful features of shearlets are their ability to efficiently approximate signals containing piecewise singularities and allowing for a much less redundant sparse tight frame representation. Consequently shearlet transform can represent images sparsely better than other representations. Moreover NSST can offer shift invariance. These properties are of vital importance for CS-MRI reconstruction because these properties results in MR images that are more compressible and hence more effective to be reconstructed. We can achieve better approximation performance and better sparse representation by shearlet coefficients compared to others. The reconstruction error in CS is proportional to approximation error. By using NSST as the sparse transform for some singularities in MR images, the reconstruction error decreases faster than that of other representations from undersampledk-space data. This is the reason why we use NSST as sparsity prior to better reconstruct some crucial features than using other representations in CS-MRI.(ii) Highly Directional Sensitivity and Anisotropy. The collection of shearlets φk,ψj,l,kd:j≥0,-2j≤l≤2j-1,k∈Z2,d=0,1 is a Parseval frame (tight frame) for L2(R2) [25]. This indicates that the decomposition is invertible and the transformation is numerically well-conditioned. Details about this construction can be found in [25, 28]. The tiling of the frequency plane R^2 induced by the shearlets is illustrated in Figure 1(a).Figure 1
Frequency tiling of the shearlets and an example of basis function. (a) The tiling of the frequency plane induced by the shearlets, (b) basis function, and (c) the frequency support of a shearlet on a pair of trapezoids.
(a)
(b)
(c)Shearlets are compactly supported in the frequency domain and have fast decay in the spatial domain with time-frequency localization. The frequency support of a shearletψj,l,k satisfies parabolic scaling. Each element ψ^j,l,k is supported on a pair of trapezoids, of approximate size 22j×2j (see Figure 1(c)), oriented along lines of slope l2-j. As a consequence, shearlets exhibit highly directional sensitivity. Unlike the isotropic elements of wavelet bases, a pair of trapezoids elements of shearlet transform possesses very high directional sensitivity and anisotropy. Display of a basis function for the shearlet transform is shown in Figure 1(b). As a result, the shearlet transform requires fewer coefficients to represent curve-like feature of MR images than the other transforms do.(iii) No Restrictions on the Number of Directions for the Shearing-Highly Directional Sensitivity and Flexibility. Such elements are very efficient in representing curve-like edges and texture features. Contourlet and SFLCT can capture the geometry of images only in 2M directions which is not optimal for MR images, because the feature information of organ, tissue, or blood vessel could be in all directions. It is worth being noted that an outstanding advantage of the shearlet transform is that there are no restrictions on the number of directions for the shearing compared with other “directional wavelets” including contourlet and SFLCT. It is impossible to obtain the flexibility of directional decomposition using shearlet transform for a fan filter implementation of DFBs. In addition, there are no constraints on the size of the supports for the shearing in shearlet transform, in which the shearing filter can have smaller support size than the directional filters used in contourlet transform. Note that the only restriction on the shearing filters is that the size of the filter should be more than or equal to the maximum number of directional subbands.(iv) The Small Support Sizes of the Shearing Filters. Another benefit of this implementation is that the shearing filtering can be performed “directly” in the time-domain using a convolution. In addition, the small support sizes of the filters reduce the Gibbs-type ringing phenomenon and improve the computational efficiency of the algorithm. The shearing filters constructed by using a Meyer window are shown in Figure 2.Figure 2
Examples of the shearing filters constructed using a Meyer wavelet as the window function.(v) Shift Invariance. In the time-domain implementation, data is undecimated at a specific scale using nonsubsampled LP. Consequently NSST exhibits excellent shift invariance at the same time, which can be helpful for CS-MRI reconstruction. Although this is a highly redundant decomposition, this version will be shown to be highly effective for the purpose of eliminating the effect of acquisition noise in k-space.NSST combines the power of multiscale methods with a unique ability to capture the geometry of multidimensional data and is optimally efficient in representing images containing edges. The highly directional sensitivity of the shearlet transform and its optimal approximation properties will lead to improvements in CS-MRI applications. So NSST is a good candidate for CS-MRI as sparsity prior. Incorporating NSST sparsity prior knowledge into the process of image reconstruction, we can get the sparse NSST coefficients and then get the reconstructed MR images by solving aforementioned constrained optimization problem.


### 2.3.2. The Improvement on the Number of Directions and Parameters in Directional Decomposition

An important advantage of the shearlet transform over the contourlet transform is that there are no restrictions on the number of directions for the shearing. In the experiments, considering anisotropic characteristics existing in MR images and the accuracy of reconstruction we have implemented the decomposition with an arbitrary even number of directional subbands (2N, N=1,2,…) at each scale using NSST, which is more flexible than 2M directions of contourlet and SFLCT. In this implementation, we have large flexibility in the choice of the number of separated directional subbands for each decomposition level, not just restricted to 23, 24, and so forth directional subbands. We could achieve the decomposition at 6, 10, and 12 directions, which cannot be achieved by SFLCT. More flexible directional representation can provide the superior decomposing directions for anisotropic features in MR images. Thus it benefits us to capture the intrinsic characteristics information of MR images at different directions using NSST. Therefore, it only needs less transform domain coefficients than the other sparsifying transforms when the curve characteristics of MRI are expressed using NSST.After considerable quantity of experiments, it demonstrates that only three-level decomposition using NSST and 12 directions in each level for sparse representation of MR images has already been enough, which can achieve better reconstruction performance outperforming four-level and2M directions decomposition using SFLCT. Finer decomposition scales and much more directions cannot result in the improvement of reconstruction performance; instead they will increase the burden of computation. In practice, we use a Meyer wavelet window to construct the directional shearing filters, as shown in Figure 2. The shearing filters of sizes 12×12, 12×12, and 12×12 from finer to coarser were used with the number of shearing directions chosen to be 12, 12, and 12. NSST (12, 12, and 12) indicates that a Meyer-based shearing filter of size 12 with 12 directions was applied to the first, second, and third decomposition level.Figure3 illustrates the three-level shearlet decomposition of a T2 brain MR image with Ewing’s sarcoma from the Harvard University site (http://www.med.harvard.edu/AANLIB/home.html). The first level decomposition generates 6 directional subbands, the second level decomposition generates 8 directional subbands, and the third level decomposition generates 12 directional subbands.Figure 3
Decomposition coefficient of a T2 brain MR image with sarcoma based on NSST, of 3 resolution scales, 6, 8, and 12 directions for every scale, respectively. (a) T2 brain MR image with sarcoma. (b) Lowpass subband, the approximate shearlet coefficients. (c)~(e) The detail shearlet coefficients of decomposition level 1, level 2, and level 3 in every direction, respectively, with the full colormap scaling for better presentation.


### 2.3.3. Compressed Sensing MRI Reconstruction Using NSST Sparsity Prior

The reconstruction problem from undersamplingk-space data actually belongs to an underdetermined system of linear equations. It is very helpful to incorporate a priori sparse information into the nonlinear reconstruction for further improvement of reconstruction performance. In this paper, we specially employ NSST with even number directional decomposition as sparsity prior in CS-MRI applications. For the reconstruction problem, we adopt the convex l1-norm regularization to promote sparsity in NSST domain. Different numerical reconstruction algorithms to solve (4) optimization problem [33–38] have been proposed in the past few decades. Among the existing algorithms, the iterative shrinkage thresholding algorithm (ISTA) is straightforward, often-used, and effective for solving linear inverse problems. The reason is that ISTA with good robustness is easy to implement and can be integrated easily with existing sparsifying transforms in practical application. These properties are significant for CS-MRI. Considering the effectiveness and the convenience of ISTA, we adopt improved iterative shearlet thresholding to solve CS-MRI.The interpretation and convergence of iterative thresholding algorithm for solving constraintl1 norm optimization were discussed in the literature [33, 35]. The core idea of iterative thresholding is that the objective function is guaranteed to decrease on each iteration of the algorithm. This algorithm iteratively performs soft thresholding to decrease the l1 norm of the coefficients α and a gradient descent to decrease the value of y-FuΨα2. Equation (4) can be simply solved by iterative shearlet thresholding:(5)αk+1=Sθk(1cFuΨH(y-FuΨαk)+αk)=Sθk1cFuΨHr+αk,where (FuΨ)H denotes the conjugate transpose operator of FuΨ and Ψ and ΨH denote the shearlet inverse and forward transform, respectively. The choice of c should satisfy c>(FuΨ)H(FuΨ)2=λmax⁡((FuΨ)H(FuΨ)) (the maximal eigenvalue of the matrix (FuΨ)H(FuΨ)). r=y-FuΨαk is the residual in k-space. The thresholding function can be used in the coefficient iteration. Here Sθk(α) is a soft thresholding operator to shrink each entry αi of vector α according to the following:(6)Sθk(αi)=0ifαi≤θkαi-αiαiθkotherwise.Since α is a complex vector for MR images, a complex thresholding operator should be defined. Daubechies et al. [33] extended the definition of the operator and derived a complex thresholding operator, which is defined as Sθk(αiejω)=Sθk(αi)ejω, ω∈[0,2π).In order to eliminate the effect of the magnitude of MR images on stop criteria, a relative error tolerance [12] Rk=y-FuΨα^2/y2 is adopted to replace absolute stop criteria related to ε of (3). The decreasing factor ρ (0<ρ<1) is adopted to decrease the threshold in each iteration. The computational parameters η and ρ in the algorithm are constants, and we set them to be the same in all the experiments. The selecting problem of the value of η and ρ was investigated in the paper [12].From empirical analysis,η=10-6 and ρ=0.8 are an appropriate choice to acquire good performance for CS-MRI. The improved CS-MRI reconstruction algorithm from undersampled k-space measurements using iterative NSST soft thresholding is summarized in Algorithm 1.Algorithm 1:NSST_IST CS-MRI.
(
1) Require:
y∈Cm, Ψ∈Cn×d, η, ρ,(0<ρ<1).
(2) Set θ0=max⁡((FuΨ)Hr0), α0=[(0,…,0)1×d]T, x0=Ψα0, r0=y
(3) Repeat
(4) Update coefficient vector according to αk+1=αk+Sθk((FuΨ)Hrk)
(5) Computer residual rk+1=y-FuΨαk+1
(6) Update the threshold θk+1=ρθk
(7) k←k+1
(8) until  Rk=y-FuΨαk2/y2≤η
(9) return  α^,x^=Ψα^.


## 2.3.1. Major Advantages of NSST Sparsity Prior

Taking account of measurement noise ink-space and the problems of aliasing and shift-variance caused by decimation, in our implementation for CS-MRI, we adopt the particular form of the time-domain based shearlet transform. This will be simply referred to as the nonsubsampled shearlet transform (NSST), which is to use the nonsubsampled Laplacian pyramid (LP) transform with several different combinations of the shearing filters. The idea of using NSST will be represented one by one below.(i) Essentially Optimal Sparsity and Approximation Property. The family of shearlet function forms a tight frame of L2(R2), and, thus, an image f can be represented by f=∑j,l,k〈f,ψj,l,k〉ψj,l,k. The coefficients αj,l,k(f)=f,ψj,l,k are called shearlet coefficients of the image f. The shearlet elements form a tight frame of well-localized waveforms, at various scales and directions, and are optimally sparse in representing images with edges. They are compactly supported in the frequency domain and have fast decay in the spatial domain.To make the statement of optimal sparsity and approximation property more rigorous, it is useful to quantify the approximation performance from the point of view of approximation theory. The asymptotic convergence rate is actually the correct optimal behavior for approximating general smooth objects having discontinuities along piecewiseC2 curves. Denoting as f^M the approximation of an image f by the M largest transform coefficients in the corresponding representation, the resulting approximation error (in l2-norm square) is f-f^M22. It is very helpful to achieve the best asymptotic decay rate for this error in application. Let C2 be the space of functions that are twice continuously differentiable. If the image f is C2 everywhere away from edge curves that are piecewise C2, the best M-term asymptotic approximation error using shearlets has a decay rate of O((log⁡M)3M-2) [28], which is essentially optimal in representing 2D images f which are piecewise C2 except for discontinuities along C2 curves and greatly outperforms that of wavelet approximations only with the decay rate of O(M-1) [29]. The shearlet transform has very similar asymptotic approximation properties with the curvelets [30] and the contourlets. The corresponding argument about optimal sparsity is proved in [9, 28, 30]. The error decay rate using shearlets is close to the theoretical optimal approximation, where the error decays as M-2 [31, 32]. In particular, the shearlet transform has some advantages over the contourlet transform [25]. In this sense, the shearlet representation provides optimally sparse representation of objects with singularities along piecewise C2 edges. Thus through shearlet transform, MR image can be decomposed into various frequency regions whose supports are contained in pair of trapezoidal regions symmetric with respect to the origin and directionally oriented.For CS-MRI, two most useful features of shearlets are their ability to efficiently approximate signals containing piecewise singularities and allowing for a much less redundant sparse tight frame representation. Consequently shearlet transform can represent images sparsely better than other representations. Moreover NSST can offer shift invariance. These properties are of vital importance for CS-MRI reconstruction because these properties results in MR images that are more compressible and hence more effective to be reconstructed. We can achieve better approximation performance and better sparse representation by shearlet coefficients compared to others. The reconstruction error in CS is proportional to approximation error. By using NSST as the sparse transform for some singularities in MR images, the reconstruction error decreases faster than that of other representations from undersampledk-space data. This is the reason why we use NSST as sparsity prior to better reconstruct some crucial features than using other representations in CS-MRI.(ii) Highly Directional Sensitivity and Anisotropy. The collection of shearlets φk,ψj,l,kd:j≥0,-2j≤l≤2j-1,k∈Z2,d=0,1 is a Parseval frame (tight frame) for L2(R2) [25]. This indicates that the decomposition is invertible and the transformation is numerically well-conditioned. Details about this construction can be found in [25, 28]. The tiling of the frequency plane R^2 induced by the shearlets is illustrated in Figure 1(a).Figure 1
Frequency tiling of the shearlets and an example of basis function. (a) The tiling of the frequency plane induced by the shearlets, (b) basis function, and (c) the frequency support of a shearlet on a pair of trapezoids.
(a)
(b)
(c)Shearlets are compactly supported in the frequency domain and have fast decay in the spatial domain with time-frequency localization. The frequency support of a shearletψj,l,k satisfies parabolic scaling. Each element ψ^j,l,k is supported on a pair of trapezoids, of approximate size 22j×2j (see Figure 1(c)), oriented along lines of slope l2-j. As a consequence, shearlets exhibit highly directional sensitivity. Unlike the isotropic elements of wavelet bases, a pair of trapezoids elements of shearlet transform possesses very high directional sensitivity and anisotropy. Display of a basis function for the shearlet transform is shown in Figure 1(b). As a result, the shearlet transform requires fewer coefficients to represent curve-like feature of MR images than the other transforms do.(iii) No Restrictions on the Number of Directions for the Shearing-Highly Directional Sensitivity and Flexibility. Such elements are very efficient in representing curve-like edges and texture features. Contourlet and SFLCT can capture the geometry of images only in 2M directions which is not optimal for MR images, because the feature information of organ, tissue, or blood vessel could be in all directions. It is worth being noted that an outstanding advantage of the shearlet transform is that there are no restrictions on the number of directions for the shearing compared with other “directional wavelets” including contourlet and SFLCT. It is impossible to obtain the flexibility of directional decomposition using shearlet transform for a fan filter implementation of DFBs. In addition, there are no constraints on the size of the supports for the shearing in shearlet transform, in which the shearing filter can have smaller support size than the directional filters used in contourlet transform. Note that the only restriction on the shearing filters is that the size of the filter should be more than or equal to the maximum number of directional subbands.(iv) The Small Support Sizes of the Shearing Filters. Another benefit of this implementation is that the shearing filtering can be performed “directly” in the time-domain using a convolution. In addition, the small support sizes of the filters reduce the Gibbs-type ringing phenomenon and improve the computational efficiency of the algorithm. The shearing filters constructed by using a Meyer window are shown in Figure 2.Figure 2
Examples of the shearing filters constructed using a Meyer wavelet as the window function.(v) Shift Invariance. In the time-domain implementation, data is undecimated at a specific scale using nonsubsampled LP. Consequently NSST exhibits excellent shift invariance at the same time, which can be helpful for CS-MRI reconstruction. Although this is a highly redundant decomposition, this version will be shown to be highly effective for the purpose of eliminating the effect of acquisition noise in k-space.NSST combines the power of multiscale methods with a unique ability to capture the geometry of multidimensional data and is optimally efficient in representing images containing edges. The highly directional sensitivity of the shearlet transform and its optimal approximation properties will lead to improvements in CS-MRI applications. So NSST is a good candidate for CS-MRI as sparsity prior. Incorporating NSST sparsity prior knowledge into the process of image reconstruction, we can get the sparse NSST coefficients and then get the reconstructed MR images by solving aforementioned constrained optimization problem.


## 2.3.2. The Improvement on the Number of Directions and Parameters in Directional Decomposition

An important advantage of the shearlet transform over the contourlet transform is that there are no restrictions on the number of directions for the shearing. In the experiments, considering anisotropic characteristics existing in MR images and the accuracy of reconstruction we have implemented the decomposition with an arbitrary even number of directional subbands (2N, N=1,2,…) at each scale using NSST, which is more flexible than 2M directions of contourlet and SFLCT. In this implementation, we have large flexibility in the choice of the number of separated directional subbands for each decomposition level, not just restricted to 23, 24, and so forth directional subbands. We could achieve the decomposition at 6, 10, and 12 directions, which cannot be achieved by SFLCT. More flexible directional representation can provide the superior decomposing directions for anisotropic features in MR images. Thus it benefits us to capture the intrinsic characteristics information of MR images at different directions using NSST. Therefore, it only needs less transform domain coefficients than the other sparsifying transforms when the curve characteristics of MRI are expressed using NSST.After considerable quantity of experiments, it demonstrates that only three-level decomposition using NSST and 12 directions in each level for sparse representation of MR images has already been enough, which can achieve better reconstruction performance outperforming four-level and2M directions decomposition using SFLCT. Finer decomposition scales and much more directions cannot result in the improvement of reconstruction performance; instead they will increase the burden of computation. In practice, we use a Meyer wavelet window to construct the directional shearing filters, as shown in Figure 2. The shearing filters of sizes 12×12, 12×12, and 12×12 from finer to coarser were used with the number of shearing directions chosen to be 12, 12, and 12. NSST (12, 12, and 12) indicates that a Meyer-based shearing filter of size 12 with 12 directions was applied to the first, second, and third decomposition level.Figure3 illustrates the three-level shearlet decomposition of a T2 brain MR image with Ewing’s sarcoma from the Harvard University site (http://www.med.harvard.edu/AANLIB/home.html). The first level decomposition generates 6 directional subbands, the second level decomposition generates 8 directional subbands, and the third level decomposition generates 12 directional subbands.Figure 3
Decomposition coefficient of a T2 brain MR image with sarcoma based on NSST, of 3 resolution scales, 6, 8, and 12 directions for every scale, respectively. (a) T2 brain MR image with sarcoma. (b) Lowpass subband, the approximate shearlet coefficients. (c)~(e) The detail shearlet coefficients of decomposition level 1, level 2, and level 3 in every direction, respectively, with the full colormap scaling for better presentation.


## 2.3.3. Compressed Sensing MRI Reconstruction Using NSST Sparsity Prior

The reconstruction problem from undersamplingk-space data actually belongs to an underdetermined system of linear equations. It is very helpful to incorporate a priori sparse information into the nonlinear reconstruction for further improvement of reconstruction performance. In this paper, we specially employ NSST with even number directional decomposition as sparsity prior in CS-MRI applications. For the reconstruction problem, we adopt the convex l1-norm regularization to promote sparsity in NSST domain. Different numerical reconstruction algorithms to solve (4) optimization problem [33–38] have been proposed in the past few decades. Among the existing algorithms, the iterative shrinkage thresholding algorithm (ISTA) is straightforward, often-used, and effective for solving linear inverse problems. The reason is that ISTA with good robustness is easy to implement and can be integrated easily with existing sparsifying transforms in practical application. These properties are significant for CS-MRI. Considering the effectiveness and the convenience of ISTA, we adopt improved iterative shearlet thresholding to solve CS-MRI.The interpretation and convergence of iterative thresholding algorithm for solving constraintl1 norm optimization were discussed in the literature [33, 35]. The core idea of iterative thresholding is that the objective function is guaranteed to decrease on each iteration of the algorithm. This algorithm iteratively performs soft thresholding to decrease the l1 norm of the coefficients α and a gradient descent to decrease the value of y-FuΨα2. Equation (4) can be simply solved by iterative shearlet thresholding:(5)αk+1=Sθk(1cFuΨH(y-FuΨαk)+αk)=Sθk1cFuΨHr+αk,where (FuΨ)H denotes the conjugate transpose operator of FuΨ and Ψ and ΨH denote the shearlet inverse and forward transform, respectively. The choice of c should satisfy c>(FuΨ)H(FuΨ)2=λmax⁡((FuΨ)H(FuΨ)) (the maximal eigenvalue of the matrix (FuΨ)H(FuΨ)). r=y-FuΨαk is the residual in k-space. The thresholding function can be used in the coefficient iteration. Here Sθk(α) is a soft thresholding operator to shrink each entry αi of vector α according to the following:(6)Sθk(αi)=0ifαi≤θkαi-αiαiθkotherwise.Since α is a complex vector for MR images, a complex thresholding operator should be defined. Daubechies et al. [33] extended the definition of the operator and derived a complex thresholding operator, which is defined as Sθk(αiejω)=Sθk(αi)ejω, ω∈[0,2π).In order to eliminate the effect of the magnitude of MR images on stop criteria, a relative error tolerance [12] Rk=y-FuΨα^2/y2 is adopted to replace absolute stop criteria related to ε of (3). The decreasing factor ρ (0<ρ<1) is adopted to decrease the threshold in each iteration. The computational parameters η and ρ in the algorithm are constants, and we set them to be the same in all the experiments. The selecting problem of the value of η and ρ was investigated in the paper [12].From empirical analysis,η=10-6 and ρ=0.8 are an appropriate choice to acquire good performance for CS-MRI. The improved CS-MRI reconstruction algorithm from undersampled k-space measurements using iterative NSST soft thresholding is summarized in Algorithm 1.Algorithm 1:NSST_IST CS-MRI.
(
1) Require:
y∈Cm, Ψ∈Cn×d, η, ρ,(0<ρ<1).
(2) Set θ0=max⁡((FuΨ)Hr0), α0=[(0,…,0)1×d]T, x0=Ψα0, r0=y
(3) Repeat
(4) Update coefficient vector according to αk+1=αk+Sθk((FuΨ)Hrk)
(5) Computer residual rk+1=y-FuΨαk+1
(6) Update the threshold θk+1=ρθk
(7) k←k+1
(8) until  Rk=y-FuΨαk2/y2≤η
(9) return  α^,x^=Ψα^.


## 3. Experimental Results and Discussion

To evaluate the performance of the proposed approach, we performed a large number of experiments onin vivo MR scan and a standard phantom, which contain a representative T2-weighted axial MR complex image of human brain and water phantom complex data (from Computational Imaging Group at Xiamen University) [39–41]. Sampling schemes in the experiments include 2D variable density random sampling pattern and pseudoradial line sampling. All the experiments were implemented using a laptop equipped with an Intel Core i5-2450M CPU at 2.50 GHz and 6 GB memory, employing a 32-bit Windows 7 operating system. The routines were tested in MATLAB R2011b.The reconstruction results and performance of the proposed method are compared with four existing schemes: zero-filling, TV-CG [4], orthogonal discrete wavelet transform- (ODWT-) based [4], and SFLCT-based [13] CS-MRI reconstruction. For ODWT-based method, images are decomposed using the “db4” wavelet with 4 vanishing moments and four decomposition levels, denoted as ODWT. For implementing SFLCT method [13], the level of decomposition is set to 4. SFLCT is with decomposition level [5,4,4,3], which means four decomposition levels and 25, 24, 24, and 23 directional subbands from coarse to fine scales. The “pkva” quincunx/fan filter is employed as a decomposition filter in SFLCT. For NSST-based CS-MRI reconstruction, the level of decomposition, the size of shearing filters, and the maximum number of directional subbands have been mentioned in forenamed Section 2.3.2.


### 3.1. Performance without Noise

We first study the noiseless scenario to see the performance of NSST-based sparsity representation and its application in CS-MRI reconstruction.(i) T2-Weighted Image of the Brain. In Figure 4, comparison of NSST-based CS-MRI reconstruction versus existing four approaches from undersampled k-space data for T2-weighted MR image of the brain (256 × 256) of slice 10 [39–41] is illustrated. The complex raw data is acquired from a healthy volunteer at a 3T Siemens Trio Tim MRI scanner using the T2-weighted turbo spin echo sequence with sequence parameters (TR/TE = 6100/99 ms, 220 × 220 mm field of view, and 3 mm slice thickness). Variable density random sampling pattern with 24.96% sampling rate is employed.Figure 4
Comparison of reconstruction results for T2wBrain_slice_10. (a) Reconstructed image from fully sampled data (256 × 256), (b)k-space sampling mask of variable density random sampling with 24.96% sampling rate, (c)~(g) reconstruction using TV-CG, zero-filling, ODWT-based, SFLCT-based, and NSST-based CS-MRI, respectively, and (a1), (c1)~(g1) and (a2), (c2)~(g2) zoom-in of (a), (c)~(g) for two different regions, respectively. The PSNRs of them are 34.08, 28.65, 33.99, 34.95, and 37.45, respectively.Obviously, the reconstructed MRI images using other four methods suffer considerably from less contrast and less visibility in some tissue structures, which is clearly seen to have many undesirable artifacts and loss of features. TV-CG is unable to remove aliasing artifacts seen in the zero-filling result. On the other hand, as shown in Figure4(g), the reconstruction using NSST-based approach is relatively devoid of aliasing artifacts. Gibbs ringing artifacts have been drastically mitigated employing the proposed approach compared with the other four reconstructions due to NSST sparsity constraint in our method. The quality of reconstructed results has a large degree of improvement utilizing the proposed approach. Specifically, in comparison to other methods, the reconstructed images (g) using the proposed approach can provide good contrast between gray and white matter. The anatomic structure of bilateral basal ganglia was depicted well. Cerebral cisterns and sulci were present with clear border (see the zoom-in in Figure 4(g1)). (c1)~(f1) in Figure 4 show line-like and texture-like artifacts in the smooth region and the edges of tissue are blurred. The reconstructed MR images using NSST-based approach with higher contrast and spatial resolution can preserve integrity of boundary and texture of tissue perfectly. To better illustrate the reconstruction performance based on NSST for the edge details of MR image, the comparison visual appearance is given in (c1)~(g1) and (c2)~(g2) of Figure 4. The zoom-in (c2)~(f2) exists Gibbs ringing artifacts around the edges at different levels, while (g2) is observed to mitigate the artifacts considerably.The magnitude of reconstruction error for five methods is shown on the same scale[0,0.25] in Figure 5. The magnitudes of the reconstruction error of comparing four methods show much more regions of high error indicating loss of structured features, indicating the reconstruction error for the proposed approach is minimal among these methods. The error image of zero-filling and TV-CG is observed to have significantly more structured errors indicating loss of feature information in corresponding reconstructed images. The error of SFLCT-based is smaller than the previous three but still not as good as the result of the proposed approach. In other words, the reconstruction with NSST-based could preserve more details and produce an artifact-free reconstruction that looks close to the fully sampled reconstruction images.Figure 5
Magnitude of reconstruction error for T2wBrain_slice_10 with the gray scale of[0,0.25]. (a)~(e) The difference images between fully sampled MR image (a) and reconstruction images (c)~(g) in Figure 4, respectively. The RLNEs of them are 0.1316, 0.2461, 0.1331, 0.1192, and 0.0894, respectively.
(a)
(b)
(c)
(d)
(e)(ii) Water Phantom. Figure 6 shows the case of standard water phantom complex data [39–41]. The water phantom data is acquired at 7T Varian MRI system (Varian, Palo Alto, CA, USA) with the spin echo sequence (TR/TE = 2000/100 ms, 80 × 80 mm field of view, and 2 mm slice thickness). The phantom for MRI system is used for evaluation of image quality and testing a MRI machine performance. The phantom is well suited to assess the spatial resolution qualitatively and quantitatively over a wide range. It contains simple geometric objects which are the representation of anatomical structures. They are important in evaluating MR image quality. These objects possess different gray scale, contrast, and spatial resolution. Figure 6(a) shows fully sampled reconstruction phantom image (256 × 256). Pseudoradial sampling pattern (44 lines) with 16.35% undersampling rate (Figure 6(b)) is employed on a water phantom as shown in Figure 6.Figure 6
Comparison of NSST-based CS-MRI reconstruction versus other four schemes for a water phantom. (a) Fully sampled water phantom (256 × 256), (b) simulatedk-space trajectory (pseudoradial 44 lines mask with 16.35% sampling rate), (c)~(g) reconstruction using TV-CG, zero-filling, ODWT-based, SFLCT-based, and the proposed NSST-based CS-MRI, respectively, and (a1), (c1)~(g1) and (a2), (c2)~(g2) zoom-in of (a), (c)~(g) for two different regions, respectively. The PSNRs of them are 31.02, 23.61, 28.64, 30.34, and 32.52, respectively.From the experimental results, the proposed method can achieve better visual results than others. Considering the visual perception, the reconstructed result of NSST-based CS-MRI which looks much clearer on the whole is obviously superior to other four schemes. Reconstruction using zero-filling in Figure6(d) gives poor result which exhibits curve-like artifacts and appears obscure for the objects. The reconstructed results of TV-CG (c) and ODWT-based approach (e) with less contrast and loss of considerable details still show the presence of some Gibbs ringing artifacts and blurring around edges. The reconstruction with SFLCT-based CS approach (f) can suppress these artifacts around the objects but contains a number of random textures in the background region which do not exist in the fully sampled image. These textures could be mistaken as a certain feature to distort vision and affect the judgment of the doctor. By contrast, the proposed method performs better. Reconstructed phantom (g) of Figure 6 via the proposed method with improved visual quality looks much better compared to the other four images. It is clear enough to permit the visualization of small details with less graininess. More details are visible because this image has higher spatial resolution. Such MR images tend to be more diagnostic, more accurate, and generally desirable.The reconstruction results (a1), (c1)~(g1) and (a2), (c2)~(g2) of Figure6 are zoom-in of (a), (c)~(g) for two different regions, respectively. It is worth noting that the reconstructed image by the proposed algorithm with the uniform background just like the fully sampled image eliminates artifacts to some extent. The results demonstrate that NSST outperforms other compared sparsity priors in reconstructing the curve-like characteristics of image. Meanwhile the edge of objects is much clearer and sharper than four compared approaches. The homogeneity and image intensities of the NSST-based reconstruction appear well.The magnitude of reconstruction error for five methods is shown on the same scale[0, 0.25] in Figure 7. From Figure 7, it is found that the magnitude of reconstruction error for the proposed approach is the lowest among these five methods. The error image of zero-filling and TV-CG is observed to have significantly more structured errors indicating loss of feature information in corresponding reconstructed images. The error of SFLCT-based is smaller than the previous three, but still not as good as the result of the proposed approach. The magnitudes of the reconstruction error of comparing four methods show much more regions of high error indicating loss of structured features. In other words, the reconstruction with NSST-based could preserve more details and produce an artifact-free reconstruction that looks close to the fully sampled image.Figure 7
Magnitude of reconstruction error for water phantom with the scale of[0,0.25]. (a)~(e) The difference images between fully sampled MR image (a) and reconstruction images (c)~(g) in Figure 4, respectively. The RLNEs of them are 0.0591, 0.1388, 0.0778, 0.0639, and 0.0498, respectively.
(a)
(b)
(c)
(d)
(e)The above comparison demonstrates that NSST-based CS-MRI reconstruction approach is highly effective at preserving more intrinsic characteristics information and mitigating Gibbs ringing artifact. It can significantly improve the quality of reconstructed MR image which meets the requirements of clinical diagnostic applications, which outperforms compared methods.


### 3.2. Performance with Noise

To demonstrate the performance of the proposed method in the noisy case, zero-mean complex Gaussian white noise of the standard deviationσ=10.2 is added to the k-space data for the following examples.(i) T2wBrain_Slice_27. Figure 8 exhibits reconstructed results of the proposed method and other compared methods from k-space noisy data for T2-weighted brain image employing 2D variable density random sampling with 24.96% sampling rate. The reconstruction with ODWT, SFLCT is unable to sufficiently remove obvious artifacts and noise seen in the zero-filling result. The proposed method achieves relative high quality reconstruction with superior edge preserving and less noise. The noise is greatly reduced in the reconstruction of Figure 8(h), which shows the superior ability of noise suppression.Figure 8
Comparison of reconstruction results for T2wBrain_slice_27 from noise-added undersampledk-space data. (a) The noisy fully sampled reconstruction image, (b) variable density random sampling pattern at sampling rate 24.96%, (c) added noise magnitude, and (d)~(h) reconstructed images using zero-filling, TV-CG, ODWT-based, SFLCT-based, and the proposed NSST-based CS-MRI, respectively. The PSNRs of (d)~(h) are 27.12, 31.41, 32.85, 34.23, and 36.52, respectively.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)Figure9 demonstrates the magnitude of reconstruction error for five methods on the same scale [0,0.25] under noisy circumstances. The magnitude image of the reconstruction error for NSST shows much smaller error magnitude and loss of features than that of other methods. The PSNR of noisy fully sampled reconstruction (Figure 8(a)) with respect to the noise-free reference image is approximately 29.82 dB. The PSNR of NSST-based result is about 5.11 dB, 3.67 dB, and 2.29 dB higher than that of TV-CG, ODWT, and SFLCT reconstruction, respectively. It is also 6.7 dB higher than the PSNR of the fully sampled noisy image, in spite of only using 1/4 of the full data to reconstruct. The rest of the objective evaluation indices of the proposed method also outperform other methods.Figure 9
Reconstruction error magnitudes for T2wBrain_slice_27 with noise standard deviation 10.2 with the gray scale of[0,0.25]. (a)~(e) The difference images between fully sampled MR image (Figure 8(a)) and reconstruction images (d)~(h), respectively. The RLNEs of them are 0.27, 0.17, 0.14, 0.12, and 0.09, respectively.
(a)
(b)
(c)
(d)
(e)(ii) Water Phantom. Figure 10 demonstrates the performance of our method on water phantom of Figure 6 using pseudo radial line sampling with 16.35% sampling rate (as shown in Figure 6(b)) under the noisy circumstance.Figure 10
Performance with noise standard deviation 10.2 for a water phantom. (a) The noisy fully sampled reconstruction, (b) added noise magnitude, and (c)~(g) reconstructed images using zero-filling, TV-CG, ODWT-based, SFLCT-based, and the proposed NSST-based CS-MRI method, respectively. The PSNRs of (c)~(g) are 23.76, 31.31, 28.64, 30.34, and 32.52, respectively. (h)~(l) The difference images between fully sampled MR image and reconstruction images (c)~(g), respectively. The RLNEs of them are 0.1364, 0.0589, 0.0778, 0.0639, and 0.0498, respectively.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)
(k)
(l)The reconstruction with ODWT and SFLCT is unable to sufficiently remove the aliasing and noise seen in the zero-filled result. On the other hand, our method can suppress the aliasing and noise, meanwhile preserving edges and details. TV-CG provides a comparable reconstruction to NSST contributing to the superior denoising and edge-preserving property of total variation. The magnitude image of the reconstruction error for NSST-based method demonstrates errors of much smaller magnitude and less structure than that of zero-filling, ODWT-based, SFLCT-based method.The PSNR of the fully sampled noisy image with respect to the noise-free reference is about 30.76 dB. The PSNR of NSST-based reconstruction is about 2.2 dB higher than that of SFLCT and 1.2 dB higher than that of TV-CG. It is also higher than the PSNR of both the fully sampled noisy image and other method results indicating good denoising and aliasing removal. The RLNE of NSST-based reconstruction is also the lowest among compared methods.The results are observed to be highly effective in denoising the reconstructed image, which indicates the promising performance of our method in the presence of reasonable amount of noise. The reason is that the implementation of NSST applies a nonsubsampled Laplacian pyramid decomposition, which results in less effect on subband coefficients for noise compared with critical downsampling wavelet.


### 3.3. Contrast of Reconstruction Quality at Different Sampling Rates by Objective Assessment Indices

Besides the visual appearance, the objective criteria are also essential for the reconstructed Image Quality Assessment (IQA). The reconstructed MR images with excellent perceptual visual quality and high objective evaluation indices have a high diagnostic value clinically. We desire to obtain diagnostic-quality reconstruction from highly undersampledk-space data. In order to measure the performance of the proposed method fully, we need to evaluate reconstructions quality quantitatively and qualitatively. The quality of the reconstruction is quantified using four objective evaluation indices—peak signal-to-noise ratio (PSNR), structural similarity (SSIM) index [42], transferred edge information (TEI) [43], and relative l2 norm error (RLNE) [12]. To evaluate the reconstruction error, we use the relative l2 norm error (RLNE) defined as e(x^)=x^-x~2/x~2 to measure the difference between the reconstructed image x^ and the fully sampled image x~. A lower error implies that the reconstructed image is more consistent to the fully sampled image.The curves (a)~(d) of Figure11 show comparison of objective assessment indices PSNR, SSIM, RLNE, and TEI versus sampling rate, respectively, using zero-filling, TV-CG, ODWT-based, and SFLCT-based CS-MRI reconstruction with comparison to the proposed NSST-based algorithm at different undersampling rates for T2-weighted brain image (256 × 256) of slice 10 with variable density random sampling mask. These assessment parameters demonstrate NSST-based CS-MRI reconstruction algorithm obtains the highest objective criteria except SSIM. This also demonstrates that the proposed algorithm can give the best reconstructed image at considerable low sampling rate in k-space measurements.Figure 11
Objective assessment indices of five methods for T2-weighted brain image (256 × 256) of slice 10 with variable density random sampling mask with different sampling rates. (a)~(d) PSNR, SSIM, RLNE, and TEI versus sampling rate, respectively.
(a)
(b)
(c)
(d)Table1 depicts the comparison of objective evaluation indices for aforementioned two reconstructed images according to the proposed NSST-based CS-MRI reconstruction with existing methods using three different sampling patterns with different sampling rates vividly. The values of four typical evaluation indices indicate the proposed method can obtain a better reconstruction performance and significant improvements.Table 1
Comparison of the objective assessment indices of CS-MRI reconstruction via NSST-based to compared methods using three different sampling patterns with different sampling rates.
MR images and mask
Indices
Reconstruction methods
Zero-filling
TV-CG
ODWT
SFLCT
Proposed
T2wBrain_slice_10using mask_Carteian_040
PSNR
28.9486
33.1252
32.7337
34.1371
36.0823
TEI
0.4346
0.6417
0.6319
0.6889
0.7486
SSIM
0.7514
0.8815
0.8610
0.8673
0.9317
RLNE
0.2378
0.1470
0.1538
0.1308
0.1046
T2wBrain_slice_10using mask_VDS_2496
PSNR
28.6510
34.0847
33.9880
34.9453
37.4500
TEI
0.4385
0.6935
0.6722
0.7155
0.7830
SSIM
0.6166
0.8933
0.8714
0.8686
0.9485
RLNE
0.2461
0.1316
0.1331
0.1192
0.0894
Water phantomusing mask_VDS_2496
PSNR
26.5477
33.1360
32.5484
33.9574
35.0718
TEI
0.5916
0.6819
0.6664
0.6913
0.6984
SSIM
0.6150
0.6165
0.8005
0.8235
0.8296
RLNE
0.0990
0.0463
0.0496
0.0422
0.0371
Water phantomusing mask_PRS_44_1635
PSNR
23.6072
31.0209
28.6408
30.3416
32.5182
TEI
0.3537
0.6403
0.6023
0.6616
0.6799
SSIM
0.5233
0.6321
0.7086
0.7517
0.8090
RLNE
0.1388
0.0591
0.0778
0.0639
0.0498
Note. mask_Carteian_040 denotes with variable density Cartesian sampling mask in k-space with 40% sampling rate, mask_VDS_2496 denotes variable density random sampling mask with 24.96% sampling rate, and mask_PRS_44_1635 denotes pseudoradial 44 lines sampling mask with 16.35% sampling rate.From Figure11 and Table 1, it is obvious that the proposed method not only can obtain minor reconstruction error, good edge-preserving characteristics but also can improve the spatial detail information and preserve the structural similarity of image compared with the existing algorithms (especially in the zoom-in), which can also be justified by the obtained superior values of evaluation indices (see Table 1). These objective assessment findings agree with the visual assessment.Considering IQA from various comprehensive assessment standards the best overall reconstruction performance is achieved by the proposed NSST-based CS-MRI reconstruction method. Its performance is superior to other comparing methods. Furthermore, it has good stability of MR images reconstruction.Reasons of the superior reconstruction performance for the proposed method are obvious. This is due to the fact that NSST is an excellent multiresolution geometric analysis tool which possesses optimal approximation property and provides an optimal sparse representation in terms of better spatial localization, highly directional sensitivity, and anisotropy and shift invariance. This is helpful to extract all prominent information from MR images and provide more useful reconstruction clinically with improved visual quality. It is also verified by the fact that the superiority of the proposed approach over comparing approaches and shift-invariant decomposition based on NSST overcomes pseudo-Gibbs phenomena successfully and improves the quality of the reconstructed image around edges.


### 3.4. Computational Time

The computation time of these reconstruction methods is summarized in Table2. In the implementation of TV-CG by 8 iterations, the weight for TV penalty and for transform l1 penalty is set as 0.01 and 0.005, respectively. The iteration stopping criteria of ODWT-based, SFLCT-based, and proposed NSST method are all that relative error tolerance is less than a certain threshold.Table 2
Runtime of different reconstruction methods.
Methods
Zero-filling
TV-CG
ODWT
SFLCT
Proposed
Runtime (S)
0.003
13.339
5.650
36.004
163.026
Note. Computations were performed on a laptop equipped with 2.50 GHz Intel Core i5-2450 M CPU and 6 GB memory for Figure 6(b) with 16.35% k-space data.Table2 shows clearly that the runtime of proposed NSST-based reconstruction is about 163 seconds. It seems like the most time-consuming one. Indeed the computational time of TV-CG, SFLCT, and NSST belongs to the same order. This order of difference in time is not a vital problem in the current computational environment. So the computation time of proposed method is acceptable within a rational range.The same iterative reconstruction based on subsampled shearlet transform takes 81 seconds or so, which is much less than NSST method. But its reconstruction performance is similar with that of SFLCT-based method. The difference for PSNRs of the two methods is less than 1 dB, and other objective indices are close. Thus the advantage of nonsubsampled shearlet transform is obvious for CS-MRI reconstruction as sparsity prior. The relative time-consuming reason of proposed NSST method is undecimated at a specific scale using nonsubsampled LP; in contrast, the DFB used in the contourlets decimates 2D data along vertical and horizontal direction. NSST is still computational efficiency although it is highly redundant. In addition, the computational cost may be decreased substantially with the use of code optimization.


## 3.1. Performance without Noise

We first study the noiseless scenario to see the performance of NSST-based sparsity representation and its application in CS-MRI reconstruction.(i) T2-Weighted Image of the Brain. In Figure 4, comparison of NSST-based CS-MRI reconstruction versus existing four approaches from undersampled k-space data for T2-weighted MR image of the brain (256 × 256) of slice 10 [39–41] is illustrated. The complex raw data is acquired from a healthy volunteer at a 3T Siemens Trio Tim MRI scanner using the T2-weighted turbo spin echo sequence with sequence parameters (TR/TE = 6100/99 ms, 220 × 220 mm field of view, and 3 mm slice thickness). Variable density random sampling pattern with 24.96% sampling rate is employed.Figure 4
Comparison of reconstruction results for T2wBrain_slice_10. (a) Reconstructed image from fully sampled data (256 × 256), (b)k-space sampling mask of variable density random sampling with 24.96% sampling rate, (c)~(g) reconstruction using TV-CG, zero-filling, ODWT-based, SFLCT-based, and NSST-based CS-MRI, respectively, and (a1), (c1)~(g1) and (a2), (c2)~(g2) zoom-in of (a), (c)~(g) for two different regions, respectively. The PSNRs of them are 34.08, 28.65, 33.99, 34.95, and 37.45, respectively.Obviously, the reconstructed MRI images using other four methods suffer considerably from less contrast and less visibility in some tissue structures, which is clearly seen to have many undesirable artifacts and loss of features. TV-CG is unable to remove aliasing artifacts seen in the zero-filling result. On the other hand, as shown in Figure4(g), the reconstruction using NSST-based approach is relatively devoid of aliasing artifacts. Gibbs ringing artifacts have been drastically mitigated employing the proposed approach compared with the other four reconstructions due to NSST sparsity constraint in our method. The quality of reconstructed results has a large degree of improvement utilizing the proposed approach. Specifically, in comparison to other methods, the reconstructed images (g) using the proposed approach can provide good contrast between gray and white matter. The anatomic structure of bilateral basal ganglia was depicted well. Cerebral cisterns and sulci were present with clear border (see the zoom-in in Figure 4(g1)). (c1)~(f1) in Figure 4 show line-like and texture-like artifacts in the smooth region and the edges of tissue are blurred. The reconstructed MR images using NSST-based approach with higher contrast and spatial resolution can preserve integrity of boundary and texture of tissue perfectly. To better illustrate the reconstruction performance based on NSST for the edge details of MR image, the comparison visual appearance is given in (c1)~(g1) and (c2)~(g2) of Figure 4. The zoom-in (c2)~(f2) exists Gibbs ringing artifacts around the edges at different levels, while (g2) is observed to mitigate the artifacts considerably.The magnitude of reconstruction error for five methods is shown on the same scale[0,0.25] in Figure 5. The magnitudes of the reconstruction error of comparing four methods show much more regions of high error indicating loss of structured features, indicating the reconstruction error for the proposed approach is minimal among these methods. The error image of zero-filling and TV-CG is observed to have significantly more structured errors indicating loss of feature information in corresponding reconstructed images. The error of SFLCT-based is smaller than the previous three but still not as good as the result of the proposed approach. In other words, the reconstruction with NSST-based could preserve more details and produce an artifact-free reconstruction that looks close to the fully sampled reconstruction images.Figure 5
Magnitude of reconstruction error for T2wBrain_slice_10 with the gray scale of[0,0.25]. (a)~(e) The difference images between fully sampled MR image (a) and reconstruction images (c)~(g) in Figure 4, respectively. The RLNEs of them are 0.1316, 0.2461, 0.1331, 0.1192, and 0.0894, respectively.
(a)
(b)
(c)
(d)
(e)(ii) Water Phantom. Figure 6 shows the case of standard water phantom complex data [39–41]. The water phantom data is acquired at 7T Varian MRI system (Varian, Palo Alto, CA, USA) with the spin echo sequence (TR/TE = 2000/100 ms, 80 × 80 mm field of view, and 2 mm slice thickness). The phantom for MRI system is used for evaluation of image quality and testing a MRI machine performance. The phantom is well suited to assess the spatial resolution qualitatively and quantitatively over a wide range. It contains simple geometric objects which are the representation of anatomical structures. They are important in evaluating MR image quality. These objects possess different gray scale, contrast, and spatial resolution. Figure 6(a) shows fully sampled reconstruction phantom image (256 × 256). Pseudoradial sampling pattern (44 lines) with 16.35% undersampling rate (Figure 6(b)) is employed on a water phantom as shown in Figure 6.Figure 6
Comparison of NSST-based CS-MRI reconstruction versus other four schemes for a water phantom. (a) Fully sampled water phantom (256 × 256), (b) simulatedk-space trajectory (pseudoradial 44 lines mask with 16.35% sampling rate), (c)~(g) reconstruction using TV-CG, zero-filling, ODWT-based, SFLCT-based, and the proposed NSST-based CS-MRI, respectively, and (a1), (c1)~(g1) and (a2), (c2)~(g2) zoom-in of (a), (c)~(g) for two different regions, respectively. The PSNRs of them are 31.02, 23.61, 28.64, 30.34, and 32.52, respectively.From the experimental results, the proposed method can achieve better visual results than others. Considering the visual perception, the reconstructed result of NSST-based CS-MRI which looks much clearer on the whole is obviously superior to other four schemes. Reconstruction using zero-filling in Figure6(d) gives poor result which exhibits curve-like artifacts and appears obscure for the objects. The reconstructed results of TV-CG (c) and ODWT-based approach (e) with less contrast and loss of considerable details still show the presence of some Gibbs ringing artifacts and blurring around edges. The reconstruction with SFLCT-based CS approach (f) can suppress these artifacts around the objects but contains a number of random textures in the background region which do not exist in the fully sampled image. These textures could be mistaken as a certain feature to distort vision and affect the judgment of the doctor. By contrast, the proposed method performs better. Reconstructed phantom (g) of Figure 6 via the proposed method with improved visual quality looks much better compared to the other four images. It is clear enough to permit the visualization of small details with less graininess. More details are visible because this image has higher spatial resolution. Such MR images tend to be more diagnostic, more accurate, and generally desirable.The reconstruction results (a1), (c1)~(g1) and (a2), (c2)~(g2) of Figure6 are zoom-in of (a), (c)~(g) for two different regions, respectively. It is worth noting that the reconstructed image by the proposed algorithm with the uniform background just like the fully sampled image eliminates artifacts to some extent. The results demonstrate that NSST outperforms other compared sparsity priors in reconstructing the curve-like characteristics of image. Meanwhile the edge of objects is much clearer and sharper than four compared approaches. The homogeneity and image intensities of the NSST-based reconstruction appear well.The magnitude of reconstruction error for five methods is shown on the same scale[0, 0.25] in Figure 7. From Figure 7, it is found that the magnitude of reconstruction error for the proposed approach is the lowest among these five methods. The error image of zero-filling and TV-CG is observed to have significantly more structured errors indicating loss of feature information in corresponding reconstructed images. The error of SFLCT-based is smaller than the previous three, but still not as good as the result of the proposed approach. The magnitudes of the reconstruction error of comparing four methods show much more regions of high error indicating loss of structured features. In other words, the reconstruction with NSST-based could preserve more details and produce an artifact-free reconstruction that looks close to the fully sampled image.Figure 7
Magnitude of reconstruction error for water phantom with the scale of[0,0.25]. (a)~(e) The difference images between fully sampled MR image (a) and reconstruction images (c)~(g) in Figure 4, respectively. The RLNEs of them are 0.0591, 0.1388, 0.0778, 0.0639, and 0.0498, respectively.
(a)
(b)
(c)
(d)
(e)The above comparison demonstrates that NSST-based CS-MRI reconstruction approach is highly effective at preserving more intrinsic characteristics information and mitigating Gibbs ringing artifact. It can significantly improve the quality of reconstructed MR image which meets the requirements of clinical diagnostic applications, which outperforms compared methods.


## 3.2. Performance with Noise

To demonstrate the performance of the proposed method in the noisy case, zero-mean complex Gaussian white noise of the standard deviationσ=10.2 is added to the k-space data for the following examples.(i) T2wBrain_Slice_27. Figure 8 exhibits reconstructed results of the proposed method and other compared methods from k-space noisy data for T2-weighted brain image employing 2D variable density random sampling with 24.96% sampling rate. The reconstruction with ODWT, SFLCT is unable to sufficiently remove obvious artifacts and noise seen in the zero-filling result. The proposed method achieves relative high quality reconstruction with superior edge preserving and less noise. The noise is greatly reduced in the reconstruction of Figure 8(h), which shows the superior ability of noise suppression.Figure 8
Comparison of reconstruction results for T2wBrain_slice_27 from noise-added undersampledk-space data. (a) The noisy fully sampled reconstruction image, (b) variable density random sampling pattern at sampling rate 24.96%, (c) added noise magnitude, and (d)~(h) reconstructed images using zero-filling, TV-CG, ODWT-based, SFLCT-based, and the proposed NSST-based CS-MRI, respectively. The PSNRs of (d)~(h) are 27.12, 31.41, 32.85, 34.23, and 36.52, respectively.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)Figure9 demonstrates the magnitude of reconstruction error for five methods on the same scale [0,0.25] under noisy circumstances. The magnitude image of the reconstruction error for NSST shows much smaller error magnitude and loss of features than that of other methods. The PSNR of noisy fully sampled reconstruction (Figure 8(a)) with respect to the noise-free reference image is approximately 29.82 dB. The PSNR of NSST-based result is about 5.11 dB, 3.67 dB, and 2.29 dB higher than that of TV-CG, ODWT, and SFLCT reconstruction, respectively. It is also 6.7 dB higher than the PSNR of the fully sampled noisy image, in spite of only using 1/4 of the full data to reconstruct. The rest of the objective evaluation indices of the proposed method also outperform other methods.Figure 9
Reconstruction error magnitudes for T2wBrain_slice_27 with noise standard deviation 10.2 with the gray scale of[0,0.25]. (a)~(e) The difference images between fully sampled MR image (Figure 8(a)) and reconstruction images (d)~(h), respectively. The RLNEs of them are 0.27, 0.17, 0.14, 0.12, and 0.09, respectively.
(a)
(b)
(c)
(d)
(e)(ii) Water Phantom. Figure 10 demonstrates the performance of our method on water phantom of Figure 6 using pseudo radial line sampling with 16.35% sampling rate (as shown in Figure 6(b)) under the noisy circumstance.Figure 10
Performance with noise standard deviation 10.2 for a water phantom. (a) The noisy fully sampled reconstruction, (b) added noise magnitude, and (c)~(g) reconstructed images using zero-filling, TV-CG, ODWT-based, SFLCT-based, and the proposed NSST-based CS-MRI method, respectively. The PSNRs of (c)~(g) are 23.76, 31.31, 28.64, 30.34, and 32.52, respectively. (h)~(l) The difference images between fully sampled MR image and reconstruction images (c)~(g), respectively. The RLNEs of them are 0.1364, 0.0589, 0.0778, 0.0639, and 0.0498, respectively.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)
(k)
(l)The reconstruction with ODWT and SFLCT is unable to sufficiently remove the aliasing and noise seen in the zero-filled result. On the other hand, our method can suppress the aliasing and noise, meanwhile preserving edges and details. TV-CG provides a comparable reconstruction to NSST contributing to the superior denoising and edge-preserving property of total variation. The magnitude image of the reconstruction error for NSST-based method demonstrates errors of much smaller magnitude and less structure than that of zero-filling, ODWT-based, SFLCT-based method.The PSNR of the fully sampled noisy image with respect to the noise-free reference is about 30.76 dB. The PSNR of NSST-based reconstruction is about 2.2 dB higher than that of SFLCT and 1.2 dB higher than that of TV-CG. It is also higher than the PSNR of both the fully sampled noisy image and other method results indicating good denoising and aliasing removal. The RLNE of NSST-based reconstruction is also the lowest among compared methods.The results are observed to be highly effective in denoising the reconstructed image, which indicates the promising performance of our method in the presence of reasonable amount of noise. The reason is that the implementation of NSST applies a nonsubsampled Laplacian pyramid decomposition, which results in less effect on subband coefficients for noise compared with critical downsampling wavelet.


## 3.3. Contrast of Reconstruction Quality at Different Sampling Rates by Objective Assessment Indices

Besides the visual appearance, the objective criteria are also essential for the reconstructed Image Quality Assessment (IQA). The reconstructed MR images with excellent perceptual visual quality and high objective evaluation indices have a high diagnostic value clinically. We desire to obtain diagnostic-quality reconstruction from highly undersampledk-space data. In order to measure the performance of the proposed method fully, we need to evaluate reconstructions quality quantitatively and qualitatively. The quality of the reconstruction is quantified using four objective evaluation indices—peak signal-to-noise ratio (PSNR), structural similarity (SSIM) index [42], transferred edge information (TEI) [43], and relative l2 norm error (RLNE) [12]. To evaluate the reconstruction error, we use the relative l2 norm error (RLNE) defined as e(x^)=x^-x~2/x~2 to measure the difference between the reconstructed image x^ and the fully sampled image x~. A lower error implies that the reconstructed image is more consistent to the fully sampled image.The curves (a)~(d) of Figure11 show comparison of objective assessment indices PSNR, SSIM, RLNE, and TEI versus sampling rate, respectively, using zero-filling, TV-CG, ODWT-based, and SFLCT-based CS-MRI reconstruction with comparison to the proposed NSST-based algorithm at different undersampling rates for T2-weighted brain image (256 × 256) of slice 10 with variable density random sampling mask. These assessment parameters demonstrate NSST-based CS-MRI reconstruction algorithm obtains the highest objective criteria except SSIM. This also demonstrates that the proposed algorithm can give the best reconstructed image at considerable low sampling rate in k-space measurements.Figure 11
Objective assessment indices of five methods for T2-weighted brain image (256 × 256) of slice 10 with variable density random sampling mask with different sampling rates. (a)~(d) PSNR, SSIM, RLNE, and TEI versus sampling rate, respectively.
(a)
(b)
(c)
(d)Table1 depicts the comparison of objective evaluation indices for aforementioned two reconstructed images according to the proposed NSST-based CS-MRI reconstruction with existing methods using three different sampling patterns with different sampling rates vividly. The values of four typical evaluation indices indicate the proposed method can obtain a better reconstruction performance and significant improvements.Table 1
Comparison of the objective assessment indices of CS-MRI reconstruction via NSST-based to compared methods using three different sampling patterns with different sampling rates.
MR images and mask
Indices
Reconstruction methods
Zero-filling
TV-CG
ODWT
SFLCT
Proposed
T2wBrain_slice_10using mask_Carteian_040
PSNR
28.9486
33.1252
32.7337
34.1371
36.0823
TEI
0.4346
0.6417
0.6319
0.6889
0.7486
SSIM
0.7514
0.8815
0.8610
0.8673
0.9317
RLNE
0.2378
0.1470
0.1538
0.1308
0.1046
T2wBrain_slice_10using mask_VDS_2496
PSNR
28.6510
34.0847
33.9880
34.9453
37.4500
TEI
0.4385
0.6935
0.6722
0.7155
0.7830
SSIM
0.6166
0.8933
0.8714
0.8686
0.9485
RLNE
0.2461
0.1316
0.1331
0.1192
0.0894
Water phantomusing mask_VDS_2496
PSNR
26.5477
33.1360
32.5484
33.9574
35.0718
TEI
0.5916
0.6819
0.6664
0.6913
0.6984
SSIM
0.6150
0.6165
0.8005
0.8235
0.8296
RLNE
0.0990
0.0463
0.0496
0.0422
0.0371
Water phantomusing mask_PRS_44_1635
PSNR
23.6072
31.0209
28.6408
30.3416
32.5182
TEI
0.3537
0.6403
0.6023
0.6616
0.6799
SSIM
0.5233
0.6321
0.7086
0.7517
0.8090
RLNE
0.1388
0.0591
0.0778
0.0639
0.0498
Note. mask_Carteian_040 denotes with variable density Cartesian sampling mask in k-space with 40% sampling rate, mask_VDS_2496 denotes variable density random sampling mask with 24.96% sampling rate, and mask_PRS_44_1635 denotes pseudoradial 44 lines sampling mask with 16.35% sampling rate.From Figure11 and Table 1, it is obvious that the proposed method not only can obtain minor reconstruction error, good edge-preserving characteristics but also can improve the spatial detail information and preserve the structural similarity of image compared with the existing algorithms (especially in the zoom-in), which can also be justified by the obtained superior values of evaluation indices (see Table 1). These objective assessment findings agree with the visual assessment.Considering IQA from various comprehensive assessment standards the best overall reconstruction performance is achieved by the proposed NSST-based CS-MRI reconstruction method. Its performance is superior to other comparing methods. Furthermore, it has good stability of MR images reconstruction.Reasons of the superior reconstruction performance for the proposed method are obvious. This is due to the fact that NSST is an excellent multiresolution geometric analysis tool which possesses optimal approximation property and provides an optimal sparse representation in terms of better spatial localization, highly directional sensitivity, and anisotropy and shift invariance. This is helpful to extract all prominent information from MR images and provide more useful reconstruction clinically with improved visual quality. It is also verified by the fact that the superiority of the proposed approach over comparing approaches and shift-invariant decomposition based on NSST overcomes pseudo-Gibbs phenomena successfully and improves the quality of the reconstructed image around edges.


## 3.4. Computational Time

The computation time of these reconstruction methods is summarized in Table2. In the implementation of TV-CG by 8 iterations, the weight for TV penalty and for transform l1 penalty is set as 0.01 and 0.005, respectively. The iteration stopping criteria of ODWT-based, SFLCT-based, and proposed NSST method are all that relative error tolerance is less than a certain threshold.Table 2
Runtime of different reconstruction methods.
Methods
Zero-filling
TV-CG
ODWT
SFLCT
Proposed
Runtime (S)
0.003
13.339
5.650
36.004
163.026
Note. Computations were performed on a laptop equipped with 2.50 GHz Intel Core i5-2450 M CPU and 6 GB memory for Figure 6(b) with 16.35% k-space data.Table2 shows clearly that the runtime of proposed NSST-based reconstruction is about 163 seconds. It seems like the most time-consuming one. Indeed the computational time of TV-CG, SFLCT, and NSST belongs to the same order. This order of difference in time is not a vital problem in the current computational environment. So the computation time of proposed method is acceptable within a rational range.The same iterative reconstruction based on subsampled shearlet transform takes 81 seconds or so, which is much less than NSST method. But its reconstruction performance is similar with that of SFLCT-based method. The difference for PSNRs of the two methods is less than 1 dB, and other objective indices are close. Thus the advantage of nonsubsampled shearlet transform is obvious for CS-MRI reconstruction as sparsity prior. The relative time-consuming reason of proposed NSST method is undecimated at a specific scale using nonsubsampled LP; in contrast, the DFB used in the contourlets decimates 2D data along vertical and horizontal direction. NSST is still computational efficiency although it is highly redundant. In addition, the computational cost may be decreased substantially with the use of code optimization.


## 4. Conclusions

In this paper, a novel CS-MRI reconstruction based on NSST is proposed. NSST provides good localization and improved directional selectivity. More intrinsic information can be preserved in reconstructed MR image with improved visual quality. In order to show the practical applicability of the proposed method, two representative MR images and a standard phantom are considered. The experiments demonstrate that the performance of CS-MRI reconstruction using NSST sparsity prior outperforms that of existing zero-filling, TV-CG, ODWT-based, and SFLCT-based methods. From the subjective visual and the comparisons of objective evaluation indices, it is not hard to see the proposed method’s superiority, which can preserve more prominent details of MR image and can greatly improve visual quality of reconstruction image with much less graininess and less information distortion than the others. The usage of NSST sparsity prior with optimal approximation and highly directional sensitivity can reduce aliasing and noise, by means of which reducing unexpected artifacts like in other reconstructed images while also capturing local intrinsic features effectively. Efficient implementation of the proposed method, for example, by using the computation power of Graphics Processing Unit is expected. Moreover, the application of some powerful adaptive sparse representations would be explored, such as multiscale dictionary learning in CS-MRI in a future work.



---

*Source: 615439-2015-03-25.xml*