Dual-domain denoising in three dimensional magnetic resonance imaging.

Denoising is a crucial preprocessing procedure for three dimensional magnetic resonance imaging (3D MRI). Existing denoising methods are predominantly implemented in a single domain, ignoring information in other domains. However, denoising methods are becoming increasingly complex, making analysis and implementation challenging. The present study aimed to develop a dual-domain image denoising (DDID) algorithm for 3D MRI that encapsulates information from the spatial and transform domains. In the present study, the DDID method was used to distinguish signal from noise in the spatial and frequency domains, after which robust accurate noise estimation was introduced for iterative filtering, which is simple and beneficial for computation. In addition, the proposed method was compared quantitatively and qualitatively with existing methods for synthetic and in vivo MRI datasets. The results of the present study suggested that the novel DDID algorithm performed well and provided competitive results, as compared with existing MRI denoising filters.

Introduction

Magnetic resonance imaging (MRI) is a medical imaging tool that can provide highly detailed images of tissues and organs in the human body (1,2). It is primarily used to demonstrate pathological and physiological alterations in living tissues (1,2), since it is non-invasive and can produce accurate three dimensional (3D) representations of internal structures of the human body (3). However, MRI images are easily contaminated by random noise, which limits the accuracy of quantitative measurements and the validity of diagnoses (4). Commonly used image denoising methods deal with Gaussian noise; however, denoising methods for MRI must be adapted for Rician noise MRI images provides a more accurate representation of information, which may benefit the denoising procedure (5).

Numerous MRI denoising approaches have been proposed by researchers. Sijbers et al (6,7) applied the maximum likelihood (ML) approach to reduce bias associated with the conventional denoising method, estimate the Rician noise and perform signal reconstruction. Conversely, in 2008, Aja-Fernandez et al (8,9) used the linear minimum mean square error estimator, in which the true value for each noisy pixel was estimated by a set of pixels selected from a local neighborhood. In addition, Wiener filtering is an adaptive optimal estimation method that uses the minimum mean square error as a criterion for denoising, and has been used to remove Rician noise in MRI (10).

Denoising in the spatial domain is the most direct approach. In 1985, Mcveigh et al (11) developed the spatial and temporal filters for reducing Gaussian noise in MRI images; however, both the noise and signal were reduced by the same factor and the signal-to-noise ratio (SNR) was unaffected. In 1990, Perona and Malik (12) established a multiscale smoothing scheme termed the anisotropic diffusion filter to overcome this defect; they smoothed markedly in the slow changing areas and suppressed at the boundaries by detecting the local gradient strengths in different directions in order to preserve the edges. In 2005, Buades et al (13) proposed the nonlocal means (NLM) filter, which uses the redundancy of patterns within the images to remove noise, and then adapted it for Rician noise correction in MRI (14,15). These local methods preserve the main structure and eliminate image blurring caused by subtle structures (15). However, as a result of the complexity of 3D MRI data, computational burden has emerged as a major challenge.

As well as denoising methods in the spatial domain, a number of methods, such as wavelet-based filers (16) and discrete cosine transform (17), are based on the processing of images in the transform domain. These approaches attempt to separate edge information from noise in the transform domain, and then remove noise while simultaneously preserving the edges and fine details (4). Denoising methods that are based on denoising patches, such as block matching with 3D filtering (BM3D) (18), which stacks the collected patches on top of one another and denoises them by 3D-wavelet shrinkage, have emerged as the most state-of-the-art image denoising methods in recent years (17). However, such methods are complex and thus prohibit a thorough analysis.

Knaus and Zwicker (19) developed a dual-domain image denoising (DDID) method that can produce high-quality results. In the present study, the DDID algorithm was adapted for a 3D MRI dataset with robust Rician noise estimation and 3D implementation in both the spatial and transform domains.

Data and methods

Algorithm

The present method for 3D MRI denoising was adapted from a recently proposed DDID algorithm for natural and synthetic images (20). Initially, a gradient descent was performed by progressively estimating Rician noise differentials and subtracting them iteratively from a noisy dataset. Subsequently, the noise differentials were estimated using robust kernels in both the spatial and frequency domains.

Noise removal by gradient descent

According to previous studies, a noise-free MRI image (x) will be corrupted by Rician noise (n) with zero mean and σ standard deviation (6,7). Hence, the noisy MRI dataset (y) may be represented as follows:

The denoising procedure may be regarded as an energy minimization problem for the optimal estimation of n̂. Thus, the denoising procedure was formulated as a gradient descent, as follows:

Where i is the number of iterations and λ is the scale factor that controls the step size in the direction of the gradient descent. The algorithm was initialed at x0=y. The energy term E(xi) was then estimated as a noise estimate (λni) for iteration i, as follows:

This permitted the reinterpretation of the gradient descent as a progressive removal of noise differentials (λni), which was integrated over time (i) to estimate the total noise, as follows:

Robust noise estimation in dual-domains

To accomplish progressive noise removal, the noise for each iteration should be robustly estimated. In the DDID method, the noise was estimated via wavelet shrinkage, and denoising was performed by shrinkage of the wavelet coefficients in the wavelet domain to preserve the signal and to discard the noise (20). In the present study, the signal was considered the outlier and the noise as the inlier being estimated; the denoised signal was not estimated directly, but determined by subtracting the estimated noise. Noisy MRI images may be divided into three categories: Large amplitude signals, medium amplitude signals and small amplitude noise (19,20). Large amplitude signals, like step signals, have large gradients and are easily recognized and detected in the spatial domain. Signal and noise are uncorrelated when the signal and noise amplitudes are similar (i.e. medium amplitude signal). In the present study, medium amplitude signals were detected as large amplitudes in the frequency domain, thereby allowing the small amplitude noise to be estimated by rejecting the large amplitudes.

Signals should be distinguished from noise so as to compute the noise estimation (ni) for iteration (i). Initially, the robust noise estimator for i was implemented in the spatial domain, where large amplitude signals were removed. Considering pixel p using pixels q in a neighborhood window Np with window radius r, the gradient (gi,p,q) was computed by subtracting the center pixel value (xi,p) from all the neighboring pixels (xi,p), as follows:

In this framework, the gradient was used to define a smoothly decaying range kernel (kr), which masked large signals from being denoised. If those signals could not be rejected, the noise estimation would have been biased in the frequency domain. In addition, a smoothly decaying spatial kernel (ks) was implemented to limit bias from spatially distant pixels. Subsequently, the two kernels (kr and ks) were combined to an unnormalized bilateral kernel. Finally, a discrete Fourier transform was performed to obtain the masked signal in the frequency domain (Fp), yielding the Fourier coefficients (Gi,p,f) for frequency f as follows:

The imaginary number j=√-1 was used to avoid confusion with the iteration time i. Scale parameters Ti and Si were obtained by deterministic annealing (DA) (21,22), as described previously:

Where a−1 is the rate of an exponential decay, γ is a large initial scale factor, σ2 is the noise variance of the noisy input y and σ2 defines a reference standard deviation for the spatial kernel. Similar to the parameter γr of the range kernel, γs is a small initial scale factor for the scale σs of the spatial kernel.

After removing large amplitude signals, the medium amplitude signals were similarly masked in the frequency domain in order to obtain the remaining small amplitude noise. This procedure involved the use of another range kernel K to remove large Fourier coefficients. Finally, the noise was estimated by taking the center pixel following inverse Fourier transformation of the signal. To obtain this value, the Fourier slicing theorem was applied and averaged for all Fourier coefficients, as follows:

Where Vi is another scale parameter similar to the spatial domain. It is in this case the variance of the Fourier coefficients Gi,p,f, defined using the noise variance σ2, as follows:

Comparison of DDID method with NLM and BM3D. The DDID method proposed in the present study was compared with unbiased NLM (13) and BM3D (17). The parameters in the proposed method were empirically identified as mentioned above and were the same for all noise levels. A total of 30 iterations (N=30) and a window radius of 15 (r=15) was selected. The parameters for the NLM and BM3D methods were the optimal ones, as suggested in previous studies (13,17). All algorithms were processed on a desktop computer with an Intel Core™ i3 Processor and 4 GB random-access memory, and scripted in Matlab 2013a (http://uk.mathworks.com/products/new_products/release2013a.html).

Simulated and real datasets were used to evaluate the proposed method. The 3D simulated brain database BrainWeb (http://brainweb.bic.mni.mcgill.ca/brainweb/) was used for quantitative and qualitative comparisons in different modalities, including T1-weighted (T1w), proton density (PD)-weighted and T2w, according to previous studies (23,24). All of the datasets were isotropic voxel with a resolution of 1 mm3 and a spatial resolution of 181×217×181. With this phantom, the denoising results in the various iterations are presented in Fig. 1, and are qualitatively compared with NLM and BM3D in Figs. 2–4. The quantitative evaluation of these methods in the three modalities with varying levels of intensity (1–15% of maximum intensity) is presented in Table I. Three quality measures were used to evaluate the results: i) The root mean squared error (RMSE) metric, which is commonly used in image processing to analyze the differences between original and restored data; ii) the peak-to-noise ratio (PSNR), which measures the extent to which the noise has been suppressed; and iii) the structural similarity index (SSIM), which is a measure more consistent with the human visual system (25), as follows:

Figure 1.

Denoising results following various iterations (0, 10, 20, 30 and 30). The Rician noise level was 15%. The second number denotes the increasing peak-signal-to-noise ratio.

Figure 2.

Example filtering results for an axial slice of the T1-weighted BrainWeb phantom (Rician noise level of 15%). The third row shows the absolute value of the image residuals for the various methods. NLM, nonlocal means filter; BM3D, block matching with 3D filtering.

Figure 4.

Example filtering results for an axial slice of the T2-weighted BrainWeb phantom (Rician noise level of 15%). The third row shows the absolute value of the image residuals for the various methods. NLM, nonlocal means filter; BM3D, block matching with 3D filtering.

Table I.

RMSE, PSNR and SSIM of the compared methods for different MRI modalities and Rician noise levels.

	Noise level (%)
Modality	1	3	5	7	9	11	13	15
T1w
RMSE
NLM	3.54	6.96	8.56	9.88	11.30	12.89	14.61	16.72
BM3D	1.68	3.66	5.27	6.81	8.58	10.54	12.53	15.00
Proposed	1.75	3.63	5.25	6.82	8.67	10.68	12.76	15.34
PSNR
NLM	36.74	31.07	29.27	28.03	26.86	25.72	24.63	23.46
BM3D	43.40	36.65	33.48	31.26	29.26	27.47	25.97	24.40
Proposed	43.06	36.72	33.52	31.25	29.16	27.35	25.81	24.21
SSIM
NLM	0.9700	0.9042	0.8666	0.8319	0.7968	0.7617	0.7275	0.6907
BM3D	0.9904	0.9632	0.9320	0.8993	0.8661	0.8326	0.8029	0.7697
Proposed	0.9908	0.9669	0.9356	0.9008	0.8638	0.8278	0.7967	0.7606
PDw
RMSE
NLM	4.63	7.96	10.70	13.00	14.77	16.17	17.34	18.58
BM3D	1.72	3.86	5.53	6.97	8.42	9.87	11.26	13.00
Proposed	1.78	3.94	5.61	7.01	8.41	9.76	11.15	12.80
PSNR
NLM	34.81	30.11	27.55	25.85	24.74	23.95	23.35	22.75
BM3D	43.43	36.40	33.28	31.27	29.62	28.25	27.10	25.86
Proposed	43.14	36.22	33.16	31.22	29.64	28.35	27.19	25.98
SSIM
NLM	0.9677	0.9051	0.8637	0.8259	0.7850	0.7467	0.7078	0.6731
BM3D	0.9898	0.9627	0.9334	0.9084	0.8768	0.8497	0.8201	0.7918
Proposed	0.9899	0.9633	0.9323	0.9038	0.8736	0.8472	0.8189	0.7949
T2w
RMSE
NLM	5.47	1073	15.31	18.07	19.89	21.32	22.74	24.31
BM3D	1.90	4.55	6.64	8.54	10.42	12.28	14.55	17.04
Proposed	1.94	4.62	6.62	8.40	10.29	12.20	14.45	16.90
PSNR
NLM	33.37	27.52	24.43	22.99	22.16	21.55	20.99	20.42
BM3D	42.54	34.97	31.68	29.50	27.77	26.34	24.87	23.50
Proposed	42.36	34.84	31.71	29.64	27.88	26.40	24.94	23.57
SSIM
NLM	0.9728	0.9075	0.8401	0.7856	0.7435	0.7088	0.6742	0.6463
BM3D	0.9915	0.9665	0.9403	0.9144	0.8885	0.8678	0.8406	0.8163
Proposed	0.9922	0.9694	0.9443	0.9187	0.8925	0.8707	0.8441	0.8201

RMSE, root mean squared error; PSNR, peak-to-noise ratio; SSIM, structural similarity index; MRI, magnetic resonance imaging; T1w, T1-weighted; PDw, proton density-weighted; T2-w, T2-weighted; NLM, nonlocal means filter; BM3D, block matching with 3D filtering.

Where µx and µy are the mean values of images x and y, σx and σy are the standard deviation of images x and y, σxy is the covariance of x and y, x and C1=(k1L)2 and C2=(k2L)2 (where L is the dynamic range, k1=0.01 and k2=0.03). As suggested by Wang et al (25) in 2004, the SSIM was locally estimated using the Gaussian kernel of 3×3×3 voxels. Finally, the mean value of all estimations was used as a quality metric. For the sake of clarity, both measures were estimated only in the region of interest (ROI; head tissues) obtained by removing the background, as previously described (26).

The real 3D MRI dataset (spinal cord; 5 male and 5 female; age range, 20–26 years) was acquired using a Siemens 3T scanner (echo time=4.7 msec; repetition time=2,040 msec; inversion time=900 msec; voxel size=1.0×1.0×1.0 nm; image size=180×192×192 voxels; flip angle=8°; Siemens AG, Munich, Germany). The Rician noise level was estimated to be ~3% of the maximum gray level intensity using the method proposed by Aja-Fernandez et al (8) in 2008. Since the original dataset was noisy already, only the restored images for the various methods are presented in Fig. 5. The second row is the enlarged images of the corresponding block in the first row. In addition to visual comparisons, the results of segmentation with a region growing algorithm (27) are also provided in the third row. The ROI was defined by dilating the seed for segmentation five times, and the resulting image was used as an ROI for quantitative analysis.

Figure 5.

Example of in vivo MRI data (spinal cord). The first row is the original image, the second row is the enlarged part of block and the third row is the segmentation result (demonstrated in white). The black asterisk is the seed point. NLM, nonlocal means filter; BM3D, block matching with 3D filtering.

Results and Discussion

Evaluation of simulated datasets

Fig. 1 demonstrates that there were improvements in the visual quality and PSNR over time. During the initial 20 iterations, the PSNR markedly increased; however, it subsequently slowed down until convergence. This is not surprising, since the proposed method was able to effectively distinguish noise and large amplitude signals in the spatial domain, and then moderate the medium amplitude signal and noise in the frequency domain until convergence.

Figs. 2–4 present the denoised MRI images using the NLM, BM3D and proposed DDID methods in the T1w, PDw and T2w modalities, respectively. Only images with the highest noise intensities (15%) are presented in order to emphasize the differences among the algorithms. As can be observed, the proposed DDID method produced a visually more pleasant image, as compared with the NLM method, showing less grainy and more definite restored images. In addition, the denoised images generated by the proposed algorithm appeared clearer and less blurred in homogeneous areas (for example, T1w images in Fig. 2), and the boundaries were more clearly defined (for example, PD and T2w images in Figs. 3 and 4). These results may be due to the fact that the proposed method utilized information in both the spatial and transform domains to denoise MRI datasets, whereas the other methods were predominantly implemented in a single domain. Furthermore, the residual images in the third row demonstrated that the proposed and BM3D methods were superior, as compared with the NLM method, since there only minor anatomical information could be observed in their residual images, whereas the opposite occurred for NLM.

Figure 3.

Example filtering results for an axial slice of the proton density-weighted BrainWeb phantom (Rician noise level of 15%). The third row shows the absolute value of the image residuals for the various methods. NLM, nonlocal means filter; BM3D, block matching with 3D filtering.

Table I presents a quantitative comparison of the three methods. The proposed method out-performed the NLM and BM3D methods and exhibited improvements over the NLM method. The NLM filter suffered from its lack of adaptation to particular properties of Rician noise (28). Conversely, BM3D is a state-of-the-art denoising method that has recently been developed (29); however, the proposed method was superior to the BM3D method in numerous cases. For example, the proposed method achieved better results at low noise levels for T1w data, and at high noise levels for PDw data. In addition, the best performance across the majority of noise levels was observed for the T2w data, which may be due to the high contrast on the T2w data.

Evaluation of clinical data

The consistency of the proposed method on spinal cord data was evaluated. As the original datasets already had noise, neither ideal residual images nor quantitative results could be obtained. As shown in Fig. 5, the proposed and BM3D methods exceeded NLM in terms of noise removal and sharpness. As compared with the other methods, the proposed method achieved a superior smoothing result, which was likely due to the robust estimation of noise and effective differentiation of signal. Furthermore, upon close inspection of the central area of the ROI, the boundary of the spinal cord appeared to be clearer, which may have been due to the edge protection ability afforded by the DDID method. These results support the effects observed in the simulated images. The segmentation results in the proposed denoised MRI image exhibited smoother spinal cord contour, which appeared more physiologically realistic and may be beneficial to further applications.

In conclusion, the present study developed a novel DDID algorithm for denoising 3D MRI data that is assumed to be corrupted by Rician noise. This DDID method was implemented using a simple iterative filtering scheme that is beneficial for computation. In addition, Rician noise was robustly estimated using the algorithm, which unified the spatial and frequency domains. This accurate noise estimation was introduced by filtering iteratively, which ensured the plausible denoising result. Furthermore, the proposed method was implemented for both synthetic and in vivo MRI datasets, and was shown to achieve competitive results quantitatively and qualitatively.