A Bayesian Multiresolution Approach for Noise Removal in Medical Magnetic Resonance Images

1 Introduction

During the acquisition process, noise affects the magnetic resonance (MR) image and de-noising is required because clinical diagnosis accuracy depends on the visual quality of the MR image. MR imaging (MRI) provides detailed images of organs and tissues in the human body. MRI is generally affected by random noise during the reconstruction process. Artifacts due to scanned object, magnetic susceptibility, radiofrequency coil, eddy current, pulse sequence design, and rigid and non-rigid motions are the sources of noise in MRI [26, 33, 34, 41]. The scanned object in the acquisition process causes thermal noise in MRI [31]. The noise in MRI can be Gaussian distributed or Rician distributed, depending on the reconstruction process of MRI [37]. The acquisition system of MR images can be a series or parallel acquisition system. A series acquisition system works in single-coil technology, while a parallel acquisition system works in multiple-coil technology. Single-coil technology introduces Rician distribution of noise, and multiple-coil technology introduces zero mean complex additive Gaussian noise with equal variance (ση2) in each coil [19].

Noise in MRI may be reduced by acquiring the data several times and averaging them. This method takes more acquisition time, and for that reason de-noising methods are preferred. A variety of studies on de-noising MR images have been published in the literature. The de-noising methods can be classified as a filtering approach, transform domain approach, and statistical approach [31]. Filters like spatial filters [30], non-local mean filters [28], bilateral filters [20], and anisotropic diffusion filters [24] are filtering approach filters. Filtering using wavelet transform [22], curvelet transform [7], and contourlet transform [3] are transform domain approach filtering methods. Filters based on statistical modeling and estimation methods like the maximum a posteriori (MAP) approach, maximum likelihood estimator (MLE) approach, and minimum mean square error (MMSE) approach come under statistical approaches.

De-noising using wavelet filters based on the multiresolution approach is more popular than the spatial domain and frequency domain filters. Wavelet filters decompose the noisy image into sub-bands at various levels. The signal components are estimated using thresholding (hard and soft) techniques. Neigh Shrink, Visu Shrink, Bayes Shrink, and Sure Shrink are popular wavelet de-noising methods. Wu et al. [39] proposed a wavelet-based de-noising filter for MR images. Anand and Sahambi [2] presented a de-noising method based on wavelet-based bilateral filtering for MR images. De-noising with complex wavelet transform was proposed by Zaroubi and Goelman [40]. An adaptive multiscale thresholding was proposed by Bao and Zhang [5] to de-noise MR images. Gai et al. [16] proposed a new multiscale-based de-noising algorithm. The de-noising algorithm is based on hidden Markov tree model utilizing the quaternion wavelet transform.

De-noising based on the statistical approach has been promising, as the de-noising efficiency depends on the correct estimation of noise and signal variances by use of estimators. The MLE-based de-noising approach was proposed by Sijbers and Den Dekker [36] for estimation of Rician noise level and de-noising MR images. A linear MMSE-based filter was proposed by Krissian and Aja-Fernandez [24] for filtering MRI for Rician noise. A wavelet-based median absolute deviation (MAD) estimator was developed by Coupe et al. [13] for estimating and de-noising MR images. He and Green Shields [21] proposed a post-acquisition de-noising method for MRI, using non-local maximum likelihood estimation (NLME). The authors proved that NLME performs better than MLE. Manjon et al. [29] proposed a de-noising method for removing both Rician and Gaussian distributed noise from MR images. They implemented a non-local noise estimation method, which adjusts the de-noising strength of the filter to remove the spatially varying noise present in the MR image. A homomorphic approach was proposed by Aja-Fernandez et al. [1] to reduce spatially variant noise in MR images. Rajan et al. [33] designed an NLME based on the Kolmogorov-Smirnov test to remove Rician noise from MR images. The method is an improvement over the NLME estimated method, based on the Euclidian distance. Awate and Whitaker [4] proposed a Bayes estimator for de-noising MR images. Gai et al. [17] developed a color image de-noising method based on color monogenic wavelet transform. The authors applied trivariate Gaussian distribution to capture the statistical dependencies between the wavelet coefficients and further used the MAP estimator to derive the shrinkage filter. In Ref. [14], the authors proposed an image de-noising algorithm using an improved sparse representation in three-dimensional (3D) transform domain. The performance of the algorithm highly depends on the effective patch size/block. The algorithm estimates the codes of the overlapping patches and averages the estimates. Reconstruction of the patches by finding similar ones in the image is known as block matching. Further, the patches are stacked into 3D blocks and de-noised using hard thresholding and a Wiener filter.

In this study, the proposed method combined the features of wavelet transform and Bayesian estimator to remove additive and signal-independent noise. This method utilizes the properties of wavelet coefficient: that they are independent and identically distributed [12]. It is assumed that the wavelet coefficients are random variables and distributed by a probability density function (PDF). An appropriate PDF is required for effective modeling of the wavelet coefficients. The Bayesian approach estimator utilizes the PDF to obtain noise-free wavelet coefficients. Further, the proposed method utilizes the normal inverse Gaussian PDF for modeling wavelet coefficients. Our method is based on statistical modeling of wavelet coefficients, and it is applied for medical images. The effectiveness of the proposed method highly depends on the correct choice of PDFs.

Therefore, the proposed de-noising method is totally different from the method proposed in Ref. [14].

2 Background

2.1 Wavelet Transform and Statistical Modeling of Coefficients

Wavelet transform effectively performs image de-noising by converting image information into transform coefficients. Wavelet transform based on thresholding performs de-noising by retaining the large coefficients and setting others to zero [23]. The de-noising accuracy of the technique depends on the correct choice of the threshold value. In this paper, the authors attempt to find the threshold value by modeling the wavelet coefficients.

The application of wavelet transform converts the image into four sub-bands, namely approximation, horizontal, vertical, and diagonal sub-bands.

Let x(i, j) be the (i, j)th pixel in an MR image, corrupted by additive Gaussian noise n(i, j), resulting in noisy image y(i, j).

The noisy image can be expressed as

(1) y(i,j)=x(i,j)+n(i,j).

The PDF of additive noise is given by

(2) pn(n)=(12πσn2)e(−n22σn2),

where, σn2 is the variance of noise.

Wavelet transform up to m level converts the image into the LLm, LHk, HLk, and HHk sub-bands [27]. The LLm sub-band contains the low-frequency information of the image, which possesses most of the information. Discrete wavelet transform (DWT) has the capability to describe the local features either spectrally or spatially. This feature makes DWT perform de-noising while preserving corners and edges. The wavelet-transformed image can be written as

(3) ykl(i,j)=xkl(i,j)+nkl(i,j),

where l=1, 2, 3; k=1, 2, … m; l=1; corresponds to horizontal orientation; l=2 corresponds to vertical orientation; l=3 corresponds to diagonal orientation; and y_k(i, j), x_k(i, j) and n_k(i, j) are the (i, j)th wavelet coefficients of y(i, j), x(i, j), and n(i, j), respectively, at level k. Figure 1 shows the wavelet decomposition of the MRI into level 3.

Figure 1:

Three-Level Wavelet Decomposition of MR Image.

The distribution of wavelet coefficients is peaked at zero and has a heavy tailed structure [27], as shown in Figure 2. NIG PDF fits heavy-tailed data and is analytically tractable [11]. NIG PDF has been implemented in modeling of wavelet coefficients and results in effective noise reduction [8, 9]. The NIG distribution is the mixed distribution of normal distribution and inverse Gaussian distribution. NIG distribution defines a homogeneous Levy-type process [6]. The proposed NIG prior is defined as

Figure 2:

Distribution of Wavelet Coefficients.

(4) px(x)=αδeαδkj(αδ2+x2)/πδ2+x2,

where k is the modified Bessel function with j (index)=1 and δ, α are the parameters. α defines tail heaviness and the steepness of the distribution can be controlled by it. δ defines the scale parameter [9].

The steepness of the NIG distribution curve increases with α. For heavier tails, a small value of α is preferred. When α→∞ and δ→∞, the NIG distribution results in Gaussian distribution. When α→0, the NIG PDF converts to a Cauchy distribution [9].

Figure 3 shows that the NIG distribution fits better than the Gaussian (normal) distribution. The empirical distribution of wavelet coefficients in sub-band LL at level 1 and the corresponding fitting of NIG distribution are shown in Figure 4. The HH sub-band contains the most noise information. Thus, a Gaussian distribution can perfectly model the wavelet coefficients in the HH sub-band. The empirical cumulative distribution of wavelet coefficients in the HH sub-band in level 1 and the corresponding fitting of Gaussian distribution are shown in Figure 5.

Figure 3:

Probability Plot for MRI Data.

Figure 4:

Probability Plot for LL Sub-band.

Figure 5:

Cumulative Distribution Plot for HH Sub-band.

3 Proposed De-noising Algorithm

This section discusses the de-noising algorithm for removing additive noise in MR images. In the proposed method, the input noisy image is first converted to wavelet sub-bands by applying wavelet transform. Wavelet sub-bands consist of wavelet coefficients. Further wavelet coefficients are statistically modeled by NIG and Gaussian PDFs, and processed by MLE and MAD estimators. By utilizing the dependency between the wavelet coefficients, the signal and noise information are collected. The MMSE estimator utilizes the variance information of signal and noise, and calculates the shrinkage factor. The wavelet coefficients are shrunk by the shrinkage factor. The resulting wavelet coefficients are then inverse wavelet transformed to obtain the desired de-noised image. Figure 6 shows the block diagram of the proposed Bayesian multiresolution approach for de-noising MR images. A detailed explanation of block diagram and summary of the proposed de-noising algorithm are discussed through the following steps.

Figure 6:

Block Diagram of Bayesian Multiresolution Approach for De-noising MR Image.

1. Perform wavelet transform of the input noisy image, degraded by additive noise.

2. Model the wavelet coefficients using NIG PDF, generated in step 1.

3. By applying MLE, find the signal parameters α^ and δ^ using Eqs. (7) and (8), respectively.

4. Find the signal and noise variances using Eqs. (9) and (10), respectively.

5. Obtain the modified wavelet coefficients applying the MMSE estimator given in Eq. (6).

6. Carry out inverse wavelet transform to reconstruct the noise-free image.

4 Simulation Results

Simulation is carried out using an MR image of size 400×400 [32]. A three-level DWT using Daubechies 8 (db8) wavelet is carried out to decompose the MRI. Qualitative and quantitative comparisons are performed, to evaluate the effectiveness of the proposed method. The performance and quality parameters of the proposed method are compared with Donoho’s soft thresholding [15], Bayes Shrink [10], and method in Ref. [8]. Figure 7 shows a qualitative comparison of the proposed method and state-of-the-art methods. It can be noticed that the proposed algorithm has the best visual quality compared with the other compared methods. The comparison of the proposed method and state-of-the-art methods in terms of peak signal-to-noise ratio (PSNR), structure similarity index (SSIM), Pratt’s figure of merit (pratt-FOM) [18], Bhattacharyya coefficient (BC) [38], mean squared difference (MSD) [35], and signal-to-noise ratio (SNR) is given in Tables 1–6, respectively. PSNR is the performance evaluation parameter. SSIM and BC are the similarity measures. The maximum value in SSIM and BC shows high similarity between the measured quantities. pratt-FOM is the edge preservation parameter. It shows the edge preservation accuracy by obtaining the displacement between the detected edge location and ideal edge location. A high value of pratt-FOM indicates less difference between the detected and actual edge points. The value varies between 0 and 1, for SSIM, pratt-FOM, and BC parameters. MSD shows the degree of dissimilarity. A high value of MSD shows less similarity, and less value shows high similarity between the measured quantities and pixel intensities. SNR is a performance parameter and is defined as the ratio of signal power to noise power. A high value of SNR is required for better results. The comparison parameters are discussed as follows.

Figure 7:

De-noising Performance on MR Image.

(A) Ground truth. (B) Ground truth image corrupted with white Gaussian noise with σ_η=0.3. (C) Image de-noising using Donoho’s soft threshold. (D) Image de-noising using Bayes Shrink. (E) Image de-noising using method in Ref. [8]. (F) Image de-noising using the proposed method.

1. The PSNR is defined as

   (11) PSNR=20 log10255mse,

   where mse is defined as

   (12) mse=1m×n∑i=1m×n(y^(i,j)−y(i,j))2.

2. The SSIM is the quality assessment parameter and shows the similarity between two images. SSIM is given by

   (13) SSIM=(2y(i,j)¯y^(i,j)¯+2.55)(2σyy^+7.65)/(y(i,j)¯2+y^(i,j)¯2+2.55)(σy2+σy^2+7.65),

   where, y(i,j)¯ and y^(i,j)¯ are the expectation of input and recovered image, respectively. σyy^ is the covariance information between input and recovered image. σy2 and σy^2 are the variance of input and recovered image, respectively. For good visual quality, a unit value of SSIM is required.

3. The pratt-FOM is defined as

   (14) FOM=(1/max(IO,ID))∑j=1ID11+αD2,

   where D is the difference between the original image edge points and detected image edge points. I_O is the number of edge points in the original image. I_D is the number of edge points in the detected image. α is the scaling constant and usually taken as a value of 1/9.

4. The BC is defined as

   (15) BC=e[−12Nln(|E||QTCOQ||QTCDQ|)],

   where C_O is the covariance matrix of original image matrix and C_D is the covariance matrix of detected image matrix.

   (16) Q=(CO+CD2).

   E is the diagonal matrix containing the eigenvalues of Q and N is the number of modes required to capture 90% variance of Q.

5. The MSD is defined as

   (17) MSD=(y^(i,j)−y(i,j))2¯,

   where, (−·) denote the expected value.

6. The SNR is defined as

   (18) SNR=10 log10[∑i=1m×ny(i,j)2∑i=1m×n(y^(i,j)−y(i,j))2].

It is observed that the proposed method gains better PSNR, SSIM, pratt-FOM, BC, MSD, and SNR value irrespective of the noise parameter. The proposed approach is based on a multiresolution approach. Compared to other transform techniques, DWT provides a powerful tool for removing noise from the image. Daubechies wavelet defines the scaling signals and wavelets using more values from the signal. This feature provides great improvement in image enhancement. One of the most important properties of Daubechies wavelet transform is that it conserves the energy of signal, thus reducing the losses and improving the de-noising efficiency. All noise and signal parameters and shrinkage factor are estimated using Bayesian estimators. The Bayes approach is based on prior knowledge. NIG PDF and Gaussian PDF are used as priors to find the signal and noise data. HH sub-band contains most of the noise information. Thus, this sub-band is modeled using Gaussian PDF to obtain the correct noise variance information. Knowledge of the suitable prior has a great role in de-noising performance. It is clearly shown in Figures 3, 4, and 5 that the proposed PDFs correctly model the wavelet coefficients in different sub-bands. The obtained results from Tables 1–6 validate the choice of PDFs, and hence the effectiveness of the proposed method. Thus, wavelet transform in the Bayesian environment works effectively to reduce noise from MR image.

Referring to Table 1, it can be observed that the proposed method gains the highest value of PSNR (dB). It is 40.16 for noise standard deviation (σ_η) of 0.1, which is more than the required PSNR for medical images. The improvement of PSNR of the proposed method is 16.85%, 14.79%, and 8.83% over Donoho’s soft thresholding, Bayes Shrink, and method in Ref. [8], respectively. Table 2 shows the SSIM parameters for different de-noising methods. The SSIM of the proposed method is 0.9983, 0.9719, 0.9612, 0.9560, and 0.9431 for noise standard deviation of 0.1, 0.2, 0.3, 0.4, and 0.5, respectively. The SSIM of the proposed method is closer to 1 than that of the other tested methods. The improvement of SSIM of the proposed method is 5.61%, 3.15%, and 2.02% over Donoho’s soft thresholding, Bayes Shrink, and method in Ref. [8], respectively, for noise standard deviation of 0.1. Table 3 shows the comparison of pratt-FOM for different de-noising methods. The proposed method achieves good improvement in edge preservation, which is a desirable property in high-resolution images. For noise standard deviation of 0.1, the pratt-FOM improvement of the proposed method is 19.55%, 8.21%, and 6.61% over Donoho’s soft thresholding, Bayes Shrink, and method in Ref. [8], respectively. Table 4 shows the measured values of BC for different de-noising methods. The improvement of BC of the proposed method is 51.89%, 36.38%, and 30.74% over Donoho’s soft thresholding, Bayes Shrink, and method in Ref. [8], respectively, for noise standard deviation of 0.1. This shows the effectiveness of the proposed method in terms of similarity measures. From Table 5, it can be seen that the lowest value of MSD is 0.000324 (σ_η=0.1) for the proposed method, and it can be concluded that the degree of dissimilarity of the proposed method is less than that of other state-of-the-art methods. Table 6 shows the SNR for the proposed and state-of-the-art methods. The SNR improvement of the proposed method is 5.61%, 17.86%, 11.08%, 12.27%, and 13.79% over the next best method for noise standard deviation of 0.1, 0.2, 0.3, 0.4, and 0.5, respectively. This demonstrates the effectiveness of the proposed method for suppression of white Gaussian noise.

The comparison of de-noising algorithm in terms of computational time is shown in Table 7. The de-noising methods were run using 1.70-GHz Intel Core i3 processor with σ_η=0.3. The run time of the proposed method is better than the method proposed in Ref. [8] with the exception of Donoho’s soft threshold [15] and Bayes Shrink [10] methods. The computation time of the proposed method was reduced by 16% compared with the method proposed in Ref. [8]. However, the proposed method is still recognized as a fast method.

A high value of the SNR parameter proves the effectiveness of the proposed method in reducing Gaussian noise from MR images. The proposed algorithm is simple and computationally less complex.

5 Conclusion

The authors have proposed a de-noising method based on multiresolution and Bayesian approaches. A normal inverse Gaussian probability distribution function was used to model the wavelet coefficients. Modeled coefficients were analyzed by use of MLE to obtain the signal parameters and hence the signal variance. Noise variance was estimated using a median estimator in diagonal sub-bands. The proposed method was able to suppress the additive noise, while keeping the structure of the MR image unaltered. The algorithm simulated on the MR image showed that the combination of the Bayesian approach estimator and multiresolution approach is efficient and easier to implement.

We would like to further evaluate the performance of the method for Rician distributed noise, which will be reported in a future communication.