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 CS-based reconstruction scheme, which combines complex double-density dual-tree discrete wavelet transform (CDDDT-DWT) with fast iterative shrinkage/soft thresholding algorithm (FISTA) to efficiently reduce such visual artifacts. The CDDDT-DWT 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 CS-based reconstruction methods, the experimental results demonstrate that this novel approach achieves higher peak signal-to-noise ratio (PSNR), larger signal-to-noise 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 soft-tissue 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, CS-MRI 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 high-frequency part contains less image information than the low-frequency 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 CS-MRI, 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 pseudo-Gibbs 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 double-density dual-tree discrete wavelet transform (CDDDT-DWT) 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], two-step 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 time-complexity. ISTA is quite sensitive to the step size and its convergence speed may be rather slow especially when the measurement matrix is seriously ill-conditioned. 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 large-scale 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 CDDDT-DWT with FISTA. Although dual-tree complex wavelet transform has also been exploited in the literature [25], its directional selectivity is inferior to CDDDT-DWT. It may suffer from artifacts as well, especially when the original image contains information in several directions. The CS-MRI combining with CDDDT-DWT 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 FISTA-CDDDT method. The experimental results of the proposed approach and its comparison with other state-of-the-art 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 Double-Density Dual-Tree DWT
The CDDDT-DWT is an overcompleted discrete wavelet transform that combines double-density DWT [19] with dual-tree CWT [26]. It consists of two scale functions and and four distinct wavelets and , where is an offset from by one-half and is an offset from by one-half. One pair of wavelets and form an approximation of the Hilbert transform pair [18]; namely, , .
Two-dimensional (2D) double-density dual-tree DWT includes 2D real double-density dual-tree DWT and 2D complex double-density dual-tree DWT. The former is constructed from two oversampled 2D double-density 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 low-frequency subbands (i.e., , ) and 16 high-frequency subbands (i.e., , , , , , , , , , , , , , , , and ) to describe the details of the recovered image.
The latter is formed by utilizing four oversampled 2D double-density 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 first-level decomposition, where represents a scale filter and depicts eight wavelet filters, while and denote the filter bank structures of the second-level decomposition. Each level generates four low-frequency subbands (, , and 32 high-frequency subbands (, , ) through 2D CDDDT-DWT transform. Similar to other wavelet transforms, the redundant transform is achieved by recursively applying low-frequency subbands to complete the decomposition of each level. For each pair of subbands, CDDDT-DWT 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 complex-valued 2D wavelet function and the imaginary part, respectively [27].
2.2. Proposed FISTA-CDDDT Algorithm
The CS-MRI 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 high-frequency subbands of CDDDT-DWT, 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 CDDDT-DWT to a local optimal image , we can getwhere is the new wavelet coefficients, which can be adjusted by the proximal forward-backward iterative scheme [28] to catch the accurate coefficients. Although -norm is nonsmooth, it is separable and CDDDT-DWT 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 CDDDT-DWT () 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 double-density dual-tree and fast iterative shrinkage thresholding algorithm (FISTA-CDDDT) 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 double-density dual-tree wavelet and conventional wavelet using the software in [27, 29]. The experiments are conducted on three typical MR datasets: Shepp-Logan phantom [17], axial brain MR data, and spine MR data, as shown in Figures 3(a)–3(c). The first Shepp-Logan 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 i5-4460 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 low-frequency signals in the central region of -space and less high-frequency 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 (FISTA-DWT), the second algorithm integrates the complex dual-tree wavelet transform with FISTA (FISTA-CDT), and the third approach incorporates the complex double-density dual-tree wavelet transform with FISTA (FISTA-CDDDT). For both the simulation and experiments on in vivo data, FISTA-DWT uses a Daubechies wavelet frame with four decomposition levels as a sparsity basis. FISTA-CDT utilizes a biorthogonal Daubechies wavelet with the 9/7 filters in the first stage and then exploits the Q-filter by Kingsbury in the second stage [25]. FISTA-CDDDT applies the finite impulse response (FIR) to perfectly reconstruct the filter banks.
We first conduct the experiment on the Shepp-Logan 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 signal-to-noise ratio (PSNR), signal-to-noise ratio (SNR), structural similarity (SSIM) index, and relative error (Rel.Err) are used to evaluate the FISTA-CDDDT 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 FISTA-DWT, the pink dotted lines denote FISTA-CDT, and the blue lines represent the proposed FISTA-CDDDT.
Figure 5 gives the reconstructions by FISTA-DWT, FISTA-CDT, and FISTA-CDDDT using different sampling schemes at the same sampling ratio. The Gaussian random sampling mask shown in Figure 5(a) is applied on Shepp-Logan phantom. According to the reconstruction results presented in Figure 5, it can be seen that the image produced by FISTA-DWT contains serious artifacts due to undersampling and the one recovered by FISTA-CDT shows visible artifacts which are visually better than FISTA-DWT, whereas the artifacts in FISTA-CDDDT recovery are much less noticeable than FISTA-DWT and FISTA-CDT under the same conditions. Using the radial sampling mask illustrated in Figure 5(e), both FISTA-DWT and FISTA-CDT 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 CDDDT-DWT adopts more wavelets and performs better in directional selectivity, the image reconstructed by this novel algorithm is much better than those by FISTA-CDT and FISTA-DWT, regardless of the adopted sampling patterns. Figures 6(b)–6(d) demonstrate that FISTA-CDDDT has less reconstruction errors than FISTA-CDT and FISTA-DWT.
(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 FISTA-DWT, FISTA-CDT, and the proposed FISTA-CDDDT. 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 sum-of-square 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 FISTA-DWT, FISTA-CDT, and FISTA-CDDDT 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 FISTA-CDDDT 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 FISTA-CDDDT 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 FISTA-DWT, FISTA-CDT, and FISTA-CDDDT. 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 FISTA-CDT are more effective than FISTA-DWT 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 FISTA-CDT.
(a)
(b)
(a)
(b)
4. Discussion
Considering the superiority of CDDDT-DWT in preserving edges and maintaining higher directional selectivity, the proposed reconstruction approach combines the CDDDT-DWT with FISTA to produce better recovery results with a faster convergence rate. Although the ISTA and TwIST can be integrated with CDDDT-DWT as well, both of them were designed for simple regularization problems. Besides, they have some drawbacks that cannot be ignored. ISTA based on the operator-splitting strategy is a promising method, which has been successfully used in signal recovery. However, it belongs to the first-order 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 CDDDT-DWT transform at the same time. The absolute value in iCS is approximated by a smooth function. In addition, the searching step size of backtracking line-search 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, FISTA-CDDDT first uses proximal forward-back 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 CS-MRI based on complex double-density dual-tree wavelet transform. The filter bank structure of the CDDDT-DWT is explored. This novel approach has been applied to Shepp-Logan phantom and axial brain and spine image reconstruction and compared with two popular methods, namely, FISTA-DWT and FISTA-CDT. 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.