Abstract
Compressed sensing (CS) has been applied to accelerate magnetic resonance imaging (MRI) for many years. Due to the lack of translation invariance of the wavelet basis, undersampled MRI reconstruction based on discrete wavelet transform may result in serious artifacts. In this paper, we propose a CSbased reconstruction scheme, which combines complex doubledensity dualtree discrete wavelet transform (CDDDTDWT) with fast iterative shrinkage/soft thresholding algorithm (FISTA) to efficiently reduce such visual artifacts. The CDDDTDWT has the characteristics of shift invariance, high degree, and a good directional selectivity. In addition, FISTA has an excellent convergence rate, and the design of FISTA is simple. Compared with conventional CSbased reconstruction methods, the experimental results demonstrate that this novel approach achieves higher peak signaltonoise ratio (PSNR), larger signaltonoise ratio (SNR), better structural similarity index (SSIM), and lower relative error.
1. Introduction
Magnetic resonance imaging (MRI) is a powerful noninvasive imaging modality, which is ubiquitously used in modern medical diagnosis [1]. MRI provides comparable spatial resolution with ultrasound and yields superior performance than CT in softtissue imaging. Nevertheless, the long scanning time limits its applications. Compressed sensing (CS) can exploit the sparsity of MR images in the transform domain and perfectly recover images from fewer measurements than those suggested by the traditional Nyquist sampling theory [2–4]. Furthermore, CSMRI can reduce the number of samples, effectively shorten the scanning time, and then obtain successful recovery if two prerequisites are satisfied: (i) the raw imaging data must have a sparse representation in a known transform domain and (ii) the undersampling artifacts appear sufficiently incoherent in the sparsifying transform domain [3]. However, the quality of the reconstructed images is poor when the space data are highly undersampled and the representation is not sparse enough.
In recent years, a variety of techniques have been proposed to enhance the quality of MRI, which can be roughly classified into three categories [5]: incoherent undersampling pattern [6], sparse representation [7], and nonlinear reconstruction algorithms [8–11]. The first strategy (e.g., variable density random space sampling [6], spirals sampling [12], radial sampling [13], and Gaussian random sampling [14]) takes advantage of designing the space sampling pattern to shorten the sampling time, increase the imaging speed, and reduce the motion artifacts. However, aliasing artifacts may occur in these sampling methods. The highfrequency part contains less image information than the lowfrequency part. Hence, by using undersampling patterns, the information about details is lost in the reconstructed images. Besides, the substantial aliasing artifacts appear incoherent. In case that the sampling ratio is extremely low, it is almost impossible to remove the significant aliasing artifacts from real signals [6, 8, 9]. For the second approach, it is essential to find a suitable sparsifying transform to recover images from highly undersampling space data. The discrete wavelet transform (DWT) is widely applied in CSMRI, but it is sensitive to shift, lacks information about phase, and has poor directionality [15]. Wavelets cannot sparsely represent curves and may lead to visible artifacts. The contourlet sparse transform is another popular alternative that can efficiently capture the contour information. This transform exhibits superior performance in representing curves, but it may fail in representing singular points [16]. The stationary wavelet transform (SWT) can noticeably reduce pseudoGibbs artifacts [17]. Similar to DWT, SWT can only possess three spatial directions. Thus, when the original image involves rich directional information, the recovered images may become blurred. The complex doubledensity dualtree discrete wavelet transform (CDDDTDWT) has the characteristics of antialiasing properties and shift invariance and is approximate to continuous wavelet transform. Moreover, it has excellent directional selectivity that can better describe the direction of the original image [18, 19]. The third technique explores an effective nonlinear reconstruction algorithm to solve the optimization problem, which is usually a combination of least square fitting and norm regularization. These approaches, such as conjugate gradient [8], iterative shrinkage/soft thresholding algorithm (ISTA) [20], twostep ISTA (TwIST) [21, 22], and fast iterative shrinkage/soft thresholding algorithm (FISTA) [23], have been investigated intensively in the literature. However, each of them has limitations. For instance, the convergence speed of conjugate gradient is very slow due to the high timecomplexity. ISTA is quite sensitive to the step size and its convergence speed may be rather slow especially when the measurement matrix is seriously illconditioned. For TwIST and FISTA, their estimates are not only dependent on the previous one, but also related to two or more previous estimates. Moreover, the global convergence rate of TwIST has not been thoroughly studied, while FISTA inspired by Nesterov’s optimal algorithm [24] can be easily implemented and is sufficient to solve largescale convex problems. It has been proved that the convergence rate of FISTA is , where is the number of iterations.
To enhance the image reconstruction quality and reduce the reconstruction artifacts, in this paper, we propose a novel reconstruction scheme, which combines CDDDTDWT with FISTA. Although dualtree complex wavelet transform has also been exploited in the literature [25], its directional selectivity is inferior to CDDDTDWT. It may suffer from artifacts as well, especially when the original image contains information in several directions. The CSMRI combining with CDDDTDWT was first introduced in [15]. In comparison with [15], the FISTA algorithm [23] with faster convergence rate was utilized to replace conventional conjugate gradient algorithm that costs more computational time.
The remainder of this paper is organized as follows. Section 2 presents the new sparsity transform and briefly describes the basics of CS as well as the proposed FISTACDDDT method. The experimental results of the proposed approach and its comparison with other stateoftheart techniques are illustrated in Section 3. In Section 4, the discussion for our algorithm is presented. Section 5 concludes the paper.
2. Materials and Methods
2.1. Complex DoubleDensity DualTree DWT
The CDDDTDWT is an overcompleted discrete wavelet transform that combines doubledensity DWT [19] with dualtree CWT [26]. It consists of two scale functions and and four distinct wavelets and , where is an offset from by onehalf and is an offset from by onehalf. One pair of wavelets and form an approximation of the Hilbert transform pair [18]; namely, , .
Twodimensional (2D) doubledensity dualtree DWT includes 2D real doubledensity dualtree DWT and 2D complex doubledensity dualtree DWT. The former is constructed from two oversampled 2D doubledensity DWT in parallel, which is redundant by a factor of two. Figure 1 shows its filter bank structure, where the row and the column filters produce two lowfrequency subbands (i.e., , ) and 16 highfrequency subbands (i.e., , , , , , , , , , , , , , , , and ) to describe the details of the recovered image.
The latter is formed by utilizing four oversampled 2D doubledensity DWT in parallel to the input image. The filter bank structure of this transform can be obtained by extending the one illustrated in Figure 1. As shown in Figure 2, and make up the filter banks of the firstlevel decomposition, where represents a scale filter and depicts eight wavelet filters, while and denote the filter bank structures of the secondlevel decomposition. Each level generates four lowfrequency subbands (, , and 32 highfrequency subbands (, , ) through 2D CDDDTDWT transform. Similar to other wavelet transforms, the redundant transform is achieved by recursively applying lowfrequency subbands to complete the decomposition of each level. For each pair of subbands, CDDDTDWT takes their summation and difference to produce the 32 oriented wavelets, describing a total of 16 main directions. Besides, each main direction contains two distinct wavelet representations, which indicate the real part of a complexvalued 2D wavelet function and the imaginary part, respectively [27].
2.2. Proposed FISTACDDDT Algorithm
The CSMRI image reconstruction problem is defined as follows:where denotes the fully sampled image, is the space data acquired from a MR scanner, and () is a parameter appropriately chosen based on the noise level, which controls the difference between the object image and the reconstructed one. is the undersampled Fourier transform in MRI. is called the regularization function in the transform domain, which is generally nonsmooth. This optimization problem can be potentially solved by total variation (TV) based approaches [10], but we will not discuss them in this paper. Here, the constrained optimization problem in (1) can be transformed into the following unconstrained one by using Lagrangian function:where is a positive regularization parameter. To solve (2), “norm” is chosen as the regularization function. especially provides the simplest way to enforce the sparsity:where can be sparsely represented in this selected domain . Here, () is the 16 highfrequency subbands of CDDDTDWT, which serves as a new sparse basis. However, the solution of (3) is a NP hard problem, which means that a solution within polynomial time is not guaranteed [11].
As an alternative formulation, applying norm directly to the regularization function produces the result formally defined as
Since norm is nonsmooth and convex, (4) can be considered as the convex relaxation of (3) to effectively solve the quadratic convex problem. In the underdetermined problem (4), represents the quadratic term, which is a convex function with Lipschitz continuous gradient, and is a nonsmooth convex regularizer.
The FISTA algorithm is applied to solve the optimization problem of (4). For a given point , we can get the gradient of at bywhere denotes the gradient of at the given point , which is a specific combination of the previous estimate values , . The original FISTA based on wavelet transform has been well studied in the literature with a backtracking step size or a constant step size , both of which can provide an improved global convergence rate of [23]. For simplicity, most algorithms adopt a constant step size in the direction of the negative gradient of the convex function.
Applying the sparsity transform CDDDTDWT to a local optimal image , we can getwhere is the new wavelet coefficients, which can be adjusted by the proximal forwardbackward iterative scheme [28] to catch the accurate coefficients. Although norm is nonsmooth, it is separable and CDDDTDWT has the characteristics of tight frame. It is known that the shrinkage thresholding function with threshold is utilized to obtain the modified wavelet coefficients :
The recovered image is updated by
In (8), the inverse CDDDTDWT () is applied by the synthesis filter bank structure, constituted by inverse order of the analysis filter bank [27]. is a threshold relaxation factor to adjust , which optimizes and reduces the calculation time. When the stop condition is satisfied, we obtain the optimal solution of (4).
The proposed algorithm combining the complex doubledensity dualtree and fast iterative shrinkage thresholding algorithm (FISTACDDDT) for solving (4) is depicted as in Algorithm 1.

3. Experiments
3.1. Experimental Setup
To evaluate the performance of the proposed reconstruction algorithm, we implement the complex doubledensity dualtree wavelet and conventional wavelet using the software in [27, 29]. The experiments are conducted on three typical MR datasets: SheppLogan phantom [17], axial brain MR data, and spine MR data, as shown in Figures 3(a)–3(c). The first SheppLogan phantom is piecewise smooth and strictly sparse, which involves the directional curves and thus can be used for testing the proposed algorithm. The complex space data of the axial brain are acquired by a 3T GE MR750 scanner using fast spin echo sequence (TR/TE = 500/12.9 ms, field of view = 240 240 mm, and slice thickness = 5 mm). The spine MR data are a fully sampled space data obtained by a 3T GE MR750 system with FRFSE sequence (TR/TE = 2500/110 ms, field of view = 240 240 mm). For the sake of brevity, the size of all testing images is scaled to 256 256. All experiments are performed using MATLAB 2014b on a desktop computer with a 3.2 GHz Intel core i54460 CPU.
(a)
(b)
(c)
Gaussian random space pattern and radial undersampling pattern are used to undersample the fully sampled space raw data. For most of MR images, the space signal with a large magnitude is generally localized in the central part. Figure 4(a) shows a Gaussian random sampling pattern, which randomly collects more lowfrequency signals in the central region of space and less highfrequency signals in the peripheral region of space. The radial undersampling pattern displayed in Figure 4(b) contains 22 radial lines with a sampling ratio of 9% in the Fourier transform domain. The sampling ratio, defined as the number of sampled points divided by the total size of original image, depends on the number of radial lines. The more the radial lines are, the higher the sampling ratio will be. It is worth noting that all the experiments can use the spiral or Cartesian sampling pattern as well.
(a)
(b)
3.2. Experimental Methods
In this work, FISTA based on three different sparsity transforms is utilized to solve the optimization problem of (4). These three techniques are implemented under the same conditions. The first method combines the discrete wavelet transform with FISTA (FISTADWT), the second algorithm integrates the complex dualtree wavelet transform with FISTA (FISTACDT), and the third approach incorporates the complex doubledensity dualtree wavelet transform with FISTA (FISTACDDDT). For both the simulation and experiments on in vivo data, FISTADWT uses a Daubechies wavelet frame with four decomposition levels as a sparsity basis. FISTACDT utilizes a biorthogonal Daubechies wavelet with the 9/7 filters in the first stage and then exploits the Qfilter by Kingsbury in the second stage [25]. FISTACDDDT applies the finite impulse response (FIR) to perfectly reconstruct the filter banks.
We first conduct the experiment on the SheppLogan phantom image shown in Figure 3(a). The optimal values are experimentally set and all the three methods terminate after 100 iterations with 20% undersampling space data. In the axial brain experiment, we set the optimal parameters , , , and , while parameters in the spine experiment are set as , , , and . For different testing datasets, the tuning parameter is set to different values and the total number of iterations in all the cases is set to 120. Additionally, in both simulation and experiments on in vivo data, we add Gaussian white noise to simulate a realistic environment. For simplicity, the same standard derivation is used.
The peak signaltonoise ratio (PSNR), signaltonoise ratio (SNR), structural similarity (SSIM) index, and relative error (Rel.Err) are used to evaluate the FISTACDDDT recovery performance. The PSNR is calculated using the following equation:where denotes images reconstructed from fully sampled data, and are the number of rows and columns in the input image, respectively. MAX means the maximum possible pixel value in the input image data, and is the reconstructed image.
The SNR is defined as
The definition of the SSIM index is given bywhere and are the various local windows from the same local window in the two different images to be compared. and are the mean of and , respectively, while and represent their variance. is the covariance of and . , , and are three variables used to increase the stability of the results. Both SSIM index and SNR have the same criterion as PSNR; that is, the reconstruction quality is directly proportional to the value of the metrics.
The Rel.Err is defined asA smaller Rel.Err indicates a higher similarity between the original image and the reconstructed one.
3.3. Experimental Results
Note that, for all the figures in this part, various approaches are labeled by different colors below the images. The green dotted lines mean FISTADWT, the pink dotted lines denote FISTACDT, and the blue lines represent the proposed FISTACDDDT.
Figure 5 gives the reconstructions by FISTADWT, FISTACDT, and FISTACDDDT using different sampling schemes at the same sampling ratio. The Gaussian random sampling mask shown in Figure 5(a) is applied on SheppLogan phantom. According to the reconstruction results presented in Figure 5, it can be seen that the image produced by FISTADWT contains serious artifacts due to undersampling and the one recovered by FISTACDT shows visible artifacts which are visually better than FISTADWT, whereas the artifacts in FISTACDDDT recovery are much less noticeable than FISTADWT and FISTACDT under the same conditions. Using the radial sampling mask illustrated in Figure 5(e), both FISTADWT and FISTACDT have the streaking artifacts, while the image recovered by the proposed method is the most similar one to the original image with no streaking artifacts. These streaking artifacts may be caused by low sampling ratio. Under such a sampling ratio, classic norm based CS techniques may result in poor performance with substantial artifacts. Consequently, CS methods cannot guarantee the quality of the reconstructed image at low sampling ratios.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
Figure 6 illustrates that the proposed algorithm yields the best result with the highest PSNR and the lowest Rel.Err, where (a) and (c) describe the change of PSNR with the increasing number of iterations. All the three methods based on FISTA reconstruction have the same convergence rate. Since CDDDTDWT adopts more wavelets and performs better in directional selectivity, the image reconstructed by this novel algorithm is much better than those by FISTACDT and FISTADWT, regardless of the adopted sampling patterns. Figures 6(b)–6(d) demonstrate that FISTACDDDT has less reconstruction errors than FISTACDT and FISTADWT.
(a)
(b)
(c)
(d)
For further analysis, the PSNRs of the reconstructed images using different methods are plotted at various sampling ratios. It is clear from Figure 7 that, with the increase of the sampling ratio, the PSNR of all algorithms grows as well. Compared with traditional wavelet, the improvement of PNSR is approximately 6 dB applying radial sampling scheme. Therefore, the proposed scheme outperforms the other two in simulation.
(a)
(b)
Figures 8 and 9 present all the reconstruction results of FISTADWT, FISTACDT, and the proposed FISTACDDDT. Additionally, images in the first row are recovered by using the same sampling mask and images in the second row are magnified images of the marked regions in the first row. Two imaging cases (axial brain image and spine MR image) are compared with each other at a 20% sampling ratio. Figure 8(a) especially is reconstructed from full space data. The space data from each coil are reconstructed separately, and the final image is generated by the sumofsquare method [30].
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
Note that the PSNRs of reconstructed axial brain images using radial sampling mask by FISTADWT, FISTACDT, and FISTACDDDT are 33.99 dB, 37.58 dB, and 38.87 dB, respectively. The magnified images are shown in Figures 8(e)–8(h). From these figures, we can see that the brain structure in the local area becomes more and more distinct. The significant artifacts existing in Figure 8(b) may be caused by imperfect filter bank in the traditional wavelet. The synthesis and analysis filter banks adopted by our FISTACDDDT are more appropriate to obtain the curve details, especially in curve processing. The spine experiments (see Figure 9) demonstrate very similar results to the axial brain. Once again, this proves that the proposed method is more accurate and effective.
Table 1 gives the SNR, Rel.Err, and SSIM index for the reconstruction on an axial brain MR image at different radial sampling ratios with . These results further demonstrate that the proposed FISTACDDDT is superior to the other two methods, because it exhibits the highest SNR, best SSIM index, and lowest Rel.Err.
Figures 10 and 11 show the PSNR of the reconstructed images by FISTADWT, FISTACDT, and FISTACDDDT. In most cases, the proposed method has better performance than the other two approaches. However, when the Gaussian random sampling ratio is lower than 15%, the PSNR of our method is slightly higher than the other two. This is because when the sampling ratio is very low, the useful information about the main feature of the image is missing. The different evaluation criteria presented in Table 2 indicate that our method and FISTACDT are more effective than FISTADWT for performing reconstruction on a spine image at different sampling ratios. As the tissue structure of the cervical spine is too complicated, there is no obvious difference between the proposed algorithm and FISTACDT.
(a)
(b)
(a)
(b)
4. Discussion
Considering the superiority of CDDDTDWT in preserving edges and maintaining higher directional selectivity, the proposed reconstruction approach combines the CDDDTDWT with FISTA to produce better recovery results with a faster convergence rate. Although the ISTA and TwIST can be integrated with CDDDTDWT as well, both of them were designed for simple regularization problems. Besides, they have some drawbacks that cannot be ignored. ISTA based on the operatorsplitting strategy is a promising method, which has been successfully used in signal recovery. However, it belongs to the firstorder algorithm that converges quite slow. As a variant of ISTA, TwIST is also an iterative thresholding algorithm, which is not guaranteed to converge globally. In contrast, FISTA has a faster convergence rate and better reconstruction accuracy, as proved in [23].
It is worth noting that Zhu et al. [15] designed an improved compressed sensing MRI algorithm (iCS), a variant of nonlinear conjugate gradient descent approach, to minimize the traditional CS model in which the image should be sparse in both the total variation and the specific CDDDTDWT transform at the same time. The absolute value in iCS is approximated by a smooth function. In addition, the searching step size of backtracking linesearch in iCS was set as 5, which may be too large, leading to an inexact solution and more computational time. Unlike iCS associated with composite regularization, the simple optimization model with norm regularization is studied in this work. Furthermore, FISTACDDDT first uses proximal forwardback optimization to approximate the linearized function , then applies shrinkage thresholding function to solve the minimization problem due to the separable characteristics of norm, and finally adopts the specific linear combination of and to smartly select the search points. Besides, the threshold relaxation technique can further reduce the computational cost. Consequently, our method can gain a more accurate solution with a dramatically improved complexity of .
5. Conclusion
In this paper, we develop a new image reconstruction method for CSMRI based on complex doubledensity dualtree wavelet transform. The filter bank structure of the CDDDTDWT is explored. This novel approach has been applied to SheppLogan phantom and axial brain and spine image reconstruction and compared with two popular methods, namely, FISTADWT and FISTACDT. The reconstructed results demonstrate that our scheme improves the PSNR and SNR as well as SSIM index and reduces the reconstructed artifacts significantly. In both simulation and experiments on in vivo data, we use the FISTA as the reconstruction algorithm. However, it can only solve the unconstrained minimization problems. An algorithm that can solve both unconstrained and constrained convex optimization problems will be studied in the future work.
Competing Interests
The authors declare that they do not have any commercial or associative interest that represents a conflict of interests in connection with the work submitted.
Acknowledgments
This work was supported by National Science Foundation of China (Grants nos. 81371537, 91432301, and 81527802), Major State Basic Research Development Program of China (973 Program) (Grant no. 2013CB733803), and Fundamental Research Funds for the Central Universities of China (Grant no. WK2070000033). The authors also thank the two healthy volunteers from whom they obtained the original MR datasets.