Segmentation Based Noise Variance Estimation from Background MRI Data Jeny Rajan, Dirk Poot, Jaber Juntu, and Jan Sijbers Vision Lab University of Antwerp Antwerp 2610, Belgium {jeny.rajan,dirk.poot,jaber.juntu,jan.sijbers}@ua.ac.be

Abstract. Accurate and precise estimation of the noise variance is often of key importance as an input parameter for posterior image processing tasks. In MR images, background data is well suited for noise estimation since (theoretically) it lacks contributions from object signal. However, background data not only suffers from small contributions of object signal but also from quantization of the intensity values. In this paper, we propose a noise variance estimation method that is insensitive to quantization errors and that is robust against low intensity variations such as low contrast tissues and ghost artifacts. The proposed method starts with an automated background segmentation procedure, and proceeds then by correctly modeling the background’s histogram. The model is based on the Rayleigh distribution of the background data and accounts for intensity quantization errors. The noise variance, which is one of the parameters of the model, is then estimated using maximum likelihood estimation. The proposed method is compared with the recently proposed noise estimation methods and is shown to be more accurate. Keywords: Noise Estimation, MRI, Segmentation.

1

Introduction

Noise estimation in Magnetic Resonance Images (MRI) is important as it can play a key role in effective denoising of these images. It also finds application in quality assessment of MR images and various parameter estimation problems. Many techniques have been proposed in the literature to estimate the noise variance from MR images. A survey of these methods is given in [2]. The noise variance in MRI can be estimated from either complex or magnitude images. Usually, the estimation is done in magnitude MR image, since it is the usual output of the scanning process. Most of the methods in this category estimate the image noise variance from the background region of the image. Typical MR images usually include an empty region of air outside the tissue of interest. Especially in multi-slice or 3D images, there is an abundant number of background voxels available for noise estimation. Since the signal in these empty region of air is zero, the noise in these areas will be Rayleigh distributed. A. Campilho and M. Kamel (Eds.): ICIAR 2010, Part I, LNCS 6111, pp. 62–70, 2010. c Springer-Verlag Berlin Heidelberg 2010 

Segmentation Based Noise Variance Estimation from Background MRI Data

63

In previous work the noise variance was estimated by fitting a Rayleigh probability density function to the partial histogram of the MR image [1]. This approach proved to be highly effective as long as the noise variance is not too high. For high noise variance, however, the information in the signal region of the MR data may significantly contribute to the partial histogram and this may lead to a bias in the estimation of the variance. This is particularly applicable in the case of diffusion weighted MRI (DW-MRI) where the signal-to-noise ratio (SNR) is inherently low. The Rayleigh model of the background can also fail when ghosting artifacts are present. In this paper, we propose a method to reduce the influence of low signal-intensity areas (eg: scalp in DW-MRI of brain) and ghost effects in the noise estimation process. The improvement is achieved through background segmentation and there by estimating the noise variance by fitting the Rayleigh PDF to the histogram of the segmented background. The paper is organised as follows. Section 2 discusses the segmentation algorithm proposed for extracting the image background from the MR image. In Section 3, the noise estimation procedure followed for the estimation of noise level from the segmented background is given. Comparative analysis of the proposed method with recently proposed approaches is shown in Section 4. Finally, conclusion and remarks are given in Section 5.

2

Background Segmentation

Segmentation of an image with low SNR is a challenging task. In this work, we combine the wavelet based bivariate shrinkage and morphological operations to achieve this goal. The image is first denoised using the wavelet based method to avoid the segmentation artifacts. Morphological operations are then applied to the denoised image for segmenting the image background. In the following subsection, we will discuss the approaches we followed for segmenting the signal and background from the noisy MR image. 2.1

Noise Reduction

Presence of noise can always lead to wrong segmentation. Considering high noise levels in MR images, especially in diffusion weighted images, an efficient denoising algorithm is a must. Since the algorithm which we are developing for noise estimation is fully automatic, the smoothing operation should be adaptive to varying noise levels. Another requirement for smoothing is anisotropic behavior of filter. i.e, the edges should be preserved while smoothing. Considering these requirements, we choose the bivariate shrinkage with the Dual Tree Complex Wavelet Transform (DTCWT)[9] for denoising the noisy MR image. The DTCWT calculates the complex transform of a signal with two separate Discrete Wavelet Transform (DWT) decompositions. The two main disadvantages of DWT, the lack of shift invariance and poor directional selectivity, can be overcome by using DTCWT. The properties of DTCWT are discussed in detail in [9]. Along with DTCWT, we used the bivariate shrinkage method proposed in [3] to estimate the original wavelet coefficients from the noisy one. For

64

J. Rajan et al.

Fig. 1. Selection of threshold from the histogram of the denoised MR volume. It can be seen from the histograms that the accurate threshold selection for segmentation is not possible from noisy MR image without denosing.

the implementation of DTCWT, the wavelet software from [10] is used. In our work we applied a 4 level wavelet decomposition. The model proposed in [3] is a modification of Bayesian estimation problem where the statistical dependency between adjacent wavelet coefficients is also considered. This model can be written as  √ 2 3σn y1 2 2 y1 + y2 − (1) wˆ1 =  2 σ y1 + y22 +

where (g)+ is defined as  (g)+ =

0, if g < 0 g, otherwise

(2)

and y1 and y2 are noisy observations of adjacent wavelet coefficients w1 and w2 (w2 represents the parent of w1 ). The noise variance σˆn2 is estimated from the finest scale wavelet coefficients [4]. 2.2

Background Extraction

Once the image is denoised, a threshold t is to be estimated for creating a background mask. We computed this threshold from the histogram of the denoised image. The index of the first local minimum occurs after the maximum value (peak) in the histogram is considered as the threshold value t. An example for the selection of t is shown in Fig. 1. The MR volume is converted to a binary image based on this threshold value. Morphological closing is then applied on this binary image to fill the holes in the data area which is defined below C = (A ⊕ B)  B

(3)

Segmentation Based Noise Variance Estimation from Background MRI Data

(a)

(b)

(c)

(d)

(e)

(f)

65

Fig. 2. Automatic background segmentation from noisy MR image (a),(b) and (c) Original simulated images : T1,T2 and PD respectively (d),(e) and (f) Image after background segmentation

where A is the binary image and B is the structuring element. ⊕ denotes the morphological dilation operation and  denotes the morphological erosion operation. We used a disk shape structuring element of size 5 × 5. Morphological closing may still leave some holes inside the signal area and some noisy data in the background area. For an improved segmentation, we applied connected component analysis to select the largest connected component from the binary image C. The resulting mask (with all the holes filled) was then used to extract the background from the foreground areas of the MR image. The result of this operation is shown in Fig. 2. The algorithm was tested for T1, T2, Proton Density (PD), and Diffusion Weighted MR images with various noise levels. One problem earlier reported with the background segmentation of Diffusion Weighted MR images of brain is the improper segmentation of the scalp [2,5]. Most of the conventional algorithms segment scalp as background and this may introduce a bias in the noise estimation. This wrong segmentation is mainly due to the high contrast difference between brain and scalp area. Contrary to conventional MR images, the intensity of scalp is very low in diffusion weighted images which makes the segmentation of scalp difficult here. As the noise level increases, it becomes more and more difficult to differentiate between scalp and the noisy background. Experiments with our proposed background segmentation

66

J. Rajan et al.

(a)

(b)

(c)

Fig. 3. Background segmentation from noisy diffusion weighted MR image (a) Original Image (b) Image segmented with the conventional approach (c) Image segmented with the proposed method

(a)

(b)

Fig. 4. Background segmentation from noisy MR image with ghost artifact (a) Original Image (b) Image segmented with the proposed method

algorithm show that the segmentation results are good for diffusion weighted MR images also. For diffusion weighted MR images, we used b0 images for generating the mask and this mask is used to segment the background from all b value images. This will also help in reducing segmentation artifacts. Fig. 3 shows the segmentation of a DW-MR image with the proposed segmentation algorithm. Another issue reported with the method in [1] is the induction of error, if there is significant ghost artifacts. The proposed method will consider ghost data as part of MR signal, if there is significant ghost effect. The segmentation result of an MR image corrupted with ghost artifacts are shown in Fig. 4.

3

Estimation of the Noise Variance

The raw, complex MR data acquired in the Fourier domain is characterized by a zero mean Gaussian probability density function (PDF) [8]. After the inverse Fourier transform, the noise distribution in the real and imaginary components will be still Gaussian due to the linearity and orthogonality of the Fourier transform. But when it is transformed to a magnitude image, the noise distribution

Segmentation Based Noise Variance Estimation from Background MRI Data

67

will be no longer Gaussian but Rician distributed [7]. If A is the original signal amplitude, then the PDF of the reconstructed magnitude image M will be pM (Mi,j |Ai,j , σn ) =

Mi,j − e σn2

2 +A2 Mi,j i,j 2 2σn

 I0

Ai,j Mi,j σn2

 u(Mi,j )

(4)

where I0 (.) being the 0th order modified Bessel function of the first kind, u(.) the Heaviside step function and σn2 the variance of the noise. Mi,j is the magnitude value of the pixel (i, j) and Ai,j , the original value of the pixel without noise. In the image background, where the SNR is zero due to the lack of waterproton density in the air, the Rician PDF simplifies to a Rayleigh distribution with PDF[5] M2

Mi,j − i,j pM (Mi,j |σn ) = 2 e 2σn2 u(Mi,j ) σn

(5)

Once the segmentation algorithm is applied to the image, there will be only background, and the image variance reduces to the variance of a Rayleigh distribution. One straight forward approach to estimate the noise standard deviation from the segmented background is to use the equation which relates normal distribution and Rayleigh distribution [5] which is given below

 π −1 2 2− (6) σ n = σM 2 2 where σM is the variance of the segmented background. One problem with this approach is the possibility of over-estimation of the noise level, if the segmented area also contains signal contributions. A maximum likelihood (ML) estimation of the noise standard deviation using the joint PDF of the histogram of the Rayleigh data was proposed in [1]

  l2    K l2 l2 l2 i−1 K − 2σ02 − 2σ − 2σ − 2σi2 2 2 −e ni ln e −e − (7) σ n = arg min Nk ln e σ

i=1

where li with i = 0, ..., K denote the set of boundaries of the histogram bins, ni represent the number of observations within the bin [li−1 , li ] which are multinok mially distributed and Nk = i=1 ni . A method to select the optimal number of bins is also described in [1]. In our experiments we used Eq. (7) for the estimation of the noise standard deviation from the segmented background. Usage of the background histogram, instead of the partial histogram, makes the proposed approach more robust to higher noise levels.

4

Results and Discussion

Experiments were conducted on both synthetic and real MR images (2D and 3D) to measure the improvement achieved with the proposed method over the existing ones.

68

J. Rajan et al.

Fig. 5. Comparison of different noise estimators for Rician magnitude MR data : 100 experiments were considered for each σ value. The graph shows the mean of the estimated value divided by the actual value. The proposed method is compared with the estimators suggested in [5].

(a)

(b)

Fig. 6. MR image of a cherry tomato (a) acquired with 1 average (b) acquired with 12 averages

Synthetic data: For simulations, we used the standard MR image phantom (T1, T2 and PD) with a voxel resolution of 1 mm3 (8 bit quantization) from the Brainweb database [6]. Rician distributed data with varying σ were then generated from this noiseless image. The dimensions of the image were 181×217× 181. For computational simplicity (with wavelets), we resized it to 256×256×181 by zero padding. The graph in Fig. 5 shows the mean of the estimated value ( σ) divided by the actual value (σ). The value closer to 1 is the best estimator. The result of the experiments on a 2D slice of the above mentioned MR data is shown in Fig. 5. In the experiment, for every σ, the σ  is estimated from the mean of 100 simulations. The proposed method is compared with the recently proposed estimators mentioned in [5]. In [2] a comparative analysis of these estimators with other popular methods can be seen. The labels in the graph, Aj1, Aj2 and Aj3 refer to the estimators based on local second order moment, local mean, and local variance in the MR data, respectively. These local moments

Segmentation Based Noise Variance Estimation from Background MRI Data

69

Fig. 7. Estimated σ as a function of the number of averages n used during the acquisition. MR image of a cherry tomato was used for this experiment.

were calculated using a 7 × 7 window. It can be seen from the graph that the σ  estimated with the proposed method is closer to the ground truth than other methods. Real data: For the experiments on real data, we used the MR image of a cherry tomato. A set of MR images was reconstructed by averaging 1 to 12 acquired images. Fig. 6 shows the images reconstructed with 1 and 12 averages, respectively. Averaging was done in the complex k-space. The resulting noise variance as a function of the number of averages over n images was then estimated. The theoretical reduction of the noise standard deviation as a function of the number of images n over which the average was taken is known to be √1n . Since the experimental setup for all the acquisitions were the same, except for averaging, √ σ  × n should be constant over n. It can be seen from Fig. 7 that the proposed method exhibits this property.

5

Conclusion

In this paper, we presented a method that needs prior segmentation of MR image background for the estimation of noise level. The proposed method minimizes the artifacts introduced in the segmentation process by conventional approaches. The reliability of the method is proved on both simulated and real MR images at various noise levels. Comparative analysis with the recently proposed methods shows that the proposed approach has an edge over the existing ones. Acknowledgments. This work was financially supported by the IWT (Institute for Science and Technology, Belgium; SBO QUANTIVIAM) and the UIAP.

70

J. Rajan et al.

References 1. Sijbers, J., Poot, D., den Dekker, A.J., Pintjens, W.: Automatic estimation of the noise variance from the histogram of a magnetic resonance image. Physics in Medicine and Biology 52, 1335–1348 (2007) 2. Aja-Fern´ andez, S., Antonio, T.V., Alberola-L´ opez, C.: Noise estimation in single and multiple magnetic resonance data based on statistical models. Magnetic Resonance Imaging 27, 1397–1409 (2009) 3. S ¸ endur, L., Selesnick, I.W.: Bivariate shrinkage with local variance estimation. IEEE Signal Processing Letters 9, 438–441 (2002) 4. Donoho, D.L., Johnstone, I.M.: Ideal spatial adaptation by wavelet shrinkage. Biometrika 81, 425–455 (1994) 5. Aja-Fern´ andez, S., Alberola-L´ opez, C., Westin, C.F.: Noise and signal estimation in magnitude MRI and rician distributed images: A LMMSE approach. IEEE Trans. Medical Imaging 17, 1383–1398 (2008) 6. Cocosco, C.A., Kollokian, V., Kwan, R.K.S., Evans, A.C.: BrainWeb: online interface to a 3D MRI simulated brain database. NeuroImage 5, S425 (1997) 7. Henkelman, R.M.: Measurement of signal intensities in the presence of noise in MR images. Medical Physics 12, 232–233 (1985) 8. Wang, Y., Lei, T.: Statistical analysis of MR imaging and its applications in image modelling. In: Proceedings of the IEEE International Conference on Image Processing and Neural Networks, vol. 1, pp. 866–870 (1994) 9. Kingsbury, N.: The dual tree complex wavelet transform: A new technique for shift invariance and directional filters. In: EUSIPCO, pp. 319–322 (1998) 10. Cai, S., Li, K., Selesnick, I.: Matlab implementation of wavelet transforms, http://taco.poly.edu/WaveletSoftware/index.html