IOP PUBLISHING

PHYSICS IN MEDICINE AND BIOLOGY

Phys. Med. Biol. 55 (2010) N441–N449

doi:10.1088/0031-9155/55/16/N02

NOTE

Noise measurement from magnitude MRI using local estimates of variance and skewness Jeny Rajan, Dirk Poot, Jaber Juntu and Jan Sijbers Vision Lab, Department of Physics, University of Antwerp, Belgium E-mail: [email protected], [email protected], [email protected] and [email protected]

Received 10 May 2010, in final form 2 July 2010 Published 3 August 2010 Online at stacks.iop.org/PMB/55/N441 Abstract In this note, we address the estimation of the noise level in magnitude magnetic resonance (MR) images in the absence of background data. Most of the methods proposed earlier exploit the Rayleigh distributed background region in MR images to estimate the noise level. These methods, however, cannot be used for images where no background information is available. In this note, we propose two different approaches for noise level estimation in the absence of the image background. The first method is based on the local estimation of the noise variance using maximum likelihood estimation and the second method is based on the local estimation of the skewness of the magnitude data distribution. Experimental results on synthetic and real MR image datasets show that the proposed estimators accurately estimate the noise level in a magnitude MR image, even without background data. (Some figures in this article are in colour only in the electronic version)

1. Introduction Estimation of the noise variance from magnetic resonance images (MRI) is often of key importance as an input parameter for image post-processing tasks. The estimated noise variance gives a measure of the quality of the MR data. Moreover, the noise variance is often a crucial parameter in image processing algorithms such as noise reduction, segmentation, parameter estimation or clustering (Nowak 1999, Basu et al 2006, Koay et al 2009). Many methods have been proposed in the literature for the estimation of the noise level from MRI (Brummer et al 1993, Nowak 1999, Sijbers and den Dekker 2004, Sijbers et al 2007, AjaFern´andez et al 2008). A survey of these methods is given in Aja-Fern´andez et al (2009). Most of the methods proposed earlier estimate the noise level from the background area of the magnitude MR image, which is known to be Rayleigh distributed (Sijbers et al 2007, Aja-Fern´andez et al 2008). Unfortunately, these methods cannot be used for images where no background information is available. For MR images other than the brain, like cardiac 0031-9155/10/160441+09$30.00

© 2010 Institute of Physics and Engineering in Medicine

Printed in the UK

N441

N442

J Rajan et al

or lung images, background data may not be available, for example in the case that the field of view (FOV) is small, such that noise assumptions based on the Rayleigh distribution fail (Aja-Fern´andez et al 2009). Also, the new scanning techniques and software eliminate most of the noisy background, which in fact affect the methods based on the Rayleigh model that need a certain amount of background pixels to perform proper estimation (Aja-Fern´andez et al 2010). Recently, in the work of Aja-Fern´andez et al (2008), a method was proposed to estimate the noise in MR images without the background region. Their method was based on the assumption that the variance of the Rician distribution approaches a Gaussian distribution at high SNR. However, at low SNR, the Rician distribution is not approximated well by a Gaussian distribution, causing a bias in the estimated noise level. In this note, we propose two different approaches to address the aforementioned issues related to the Rician noise estimation. The first method is based on maximum likelihood (ML) estimation of the local variance for each pixel of the image using a local neighborhood. The noise variance can then be computed from the mode of distribution of the estimated local variance. The ML method simultaneously estimates the underlying true signal and noise variance from the Rician distributed data if the underlying signal is constant (Sijbers et al 1998). It will be shown that this method is highly accurate in estimating the noise level, but has a rather high computational complexity. Therefore, a computationally more efficient method is also proposed. This method is based on the local skewness. A correction factor for the variance estimated with the Gaussian assumption is introduced based on the estimated skewness for the actual computation of the Rician noise variance. 2. Noise distribution in MRI The raw, complex MR data acquired in the Fourier domain is characterized by a zero mean Gaussian probability density function (PDF). After the inverse Fourier transform, the noise distribution in the real and imaginary components will still be Gaussian due to the linearity and the orthogonality of the Fourier transform. However, due to the subsequent transform to a magnitude image, the noise distribution will be no longer Gaussian but Rician distributed. If A is the original signal amplitude, then the PDF of the reconstructed magnitude image M will be (Sijbers and den Dekker 2004)   AM M − M 2 +A2 2 (M), (1) pM (M|A, σ ) = 2 e 2σ I0 σ σ2 where I0 (·) denotes the zeroth-order modified Bessel function of the first kind, (·) the Heaviside step function and σ 2 denotes the variance of the Gaussian noise in the complex MR data. The νth raw moment of the Rice PDF can be analytically expressed as a function of the confluent hypergeometric function of the first kind, denoted by 1 F1 :    ν A2 ν − , (2) ; 1; − F E[M ν ] = (2σ 2 )ν/2  1 + 1 1 2 2 2σ 2 where  represents the Gamma function. When the underlying signal intensity A equals zero, the Rice PDF simplifies to a Rayleigh distribution (Henkelman 1985): M2 M pM (M|σ ) = 2 e(− 2σ 2 ) (M). (3) σ The moments of this distribution depend only on σ and can be written as  ν E[M ν ] = (2σ 2 )ν/2  1 + . (4) 2

Noise measurement from magnitude MRI using estimates of variance and skewness

N443

If MB is a non-signal background area in the MRI, then σ can be directly estimated from the moments of the Rayleigh distribution with or without explicit segmentation of the background. If an explicit segmentation algorithm is applied to extract the background, then the noise standard deviation can be estimated as 2 ¯ MB , σn = (5) π ¯ B is the mean of the segmented background. It is also possible to estimate the noise where M standard deviation without an explicit segmentation by using the local statistics (Aja-Fern´andez et al 2008): 2 ¯ i,j }), mode ({M (6) σn = π ¯ i,j corresponds to the local mean computed for each pixel at (i, j ). where M At high SNR, the Rician PDF approaches a Gaussian PDF with a mean A and variance σ 2: (M−A)2 1 pM (M|σ ) = √ e− 2σ 2 (M) (7) 2π σ 2 and the moments will be  ν−1   A2 A2 d 2σ 2 E[M ν ] = σ 2(ν−1) (8) A e e− 2σ 2 , ν−1 dA i.e. at high SNR, the noise variance can be estimated by simply computing the variance of a homogeneous area in the image. This approach is normally followed in the literature for the estimation of the noise in the absence of the background area. A homogeneous area in an image can be selected without manual intervention using the method suggested in Aja-Fern´andez et al (2008). The noise can be estimated from such a homogeneous area by

 σ2 = mode σ 2 , (9) n

Mi,j

2 where is the variance of the estimated noise and σM corresponds to the local variance i ,j computed around each pixel at (i, j ). However, at low SNR, equation (9) will lead to an underestimation of the noise level. In section 3, we discuss how to measure the noise variance from an MR image without the need for high SNR regions or a background region.

σn2

3. Noise estimation using local maximum likelihood estimation Let m = (m1 , . . . , mN ) be the N Rician distributed magnitude data points for which the underlying, noiseless intensity value is A. Then the joint PDF pm can be written as (Sijbers and den Dekker 2004)   N  AMn Mn − Mn2 +A2 2 e 2σ I0 pm = , (10) σ2 σ2 n=1 where {Mn } are the magnitude variables corresponding to the magnitude observations {mn }. The ML estimate of A and σ 2 can then be constructed by substituting the available observations {mn } in equation (10) and maximizing the resulting likelihood function L(A, σ 2 ) or equivalently log L(A, σ 2 ), with respect to A and σ 2 :   N N N m    Amn m2n + A2  n . (11) ln L = − ln + ln I0 σ2 2σ 2 σ2 n=1 n=1 n=1

N444

J Rajan et al

The ML estimate is then found from the global maximum of ln L w.r.t. A and σ 2 (Sijbers and den Dekker 2004):   (12) AML , σ2 ML = arg{max(ln L)}. A,σ 2

The above procedure assumes the underlying signal to be constant within the local neighborhood from which the signal and noise variance is estimated. If this assumption is valid for most of the local pixel neighborhoods, a robust estimator of the noise variance is given by the mode of all ML estimated local noise levels:

 2 2 = mode σ , (13) σ m

MLi,j

2 2 is the estimated noise variance and σ where σ m MLi,j is the ML estimate of the noise variance at each pixel (i, j ). In fact, most of the MR images are piecewise constant with a reasonably small number of classes (e.g. MR image of brain) (Zhang et al 2001). Since the noise is estimated from the available piecewise constant regions in the image, this method neither depends on the image background nor on the SNR of the image.

4. Noise estimation using measurement of local skewness Even though the method based on local ML estimation is accurate, a drawback is its time complexity. In this section, we propose a method based on the local computation of the skewness of the magnitude data distribution to estimate the noise variance. From equations (2), (4) and (8), it can be inferred that the moments of the Rician distribution are always in between the moments of Rayleigh and Gaussian distributions. The relationship between σ 2 and the variance of a Rician distribution σM2 at low and high SNR can be written as  π −1 (14) σ 2 = σM2 2 − 2 and σ 2 = σM2 ,

(15)

respectively. In general, σ 2 in terms of σM2 can be written as σ 2 = σM2 × ϕ,

(16) −1

where ϕ is a correction factor in the range [1; (2 − π/2) ], i.e. when the Rician distribution approaches a Rayleigh distribution (at low SNR), the correction factor tends to (2 − π/2)−1 and when the Rician distribution approaches a Gaussian (at high SNR), the correction factor tends to 1. The proximity of the Rician distribution toward Rayleigh or Gaussian can be measured using its skewness, which is defined by γ =

2E[M]3 − 3E[M]E[M 2 ] + E[M 3 ]

, (17) 3 (−E[M]2 + E[M 2 ]) 2 which can be analytically computed using equation (2). The skewness of a Rician distribution is a monotonically decreasing function of the SNR, defined as A/σ with values ranging from 0.631 (for SNR = 0) to 0 (for SNR = ∞). If we fix A, the skewness depends only on the noise variance σ 2 of the underlying Gaussian distribution. On the other hand, the variance of the Rician distributed data, σM2 , is given by σM2 = E[M 2 ] − E[M]2 .

(18)

Noise measurement from magnitude MRI using estimates of variance and skewness

N445

Correction Factor

2.5

2

1.5

1

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

Skewness

Figure 1. Relation between the skewness and variance correction factor for Rician distributed data.

If A is fixed, σM2 also solely depends on σ 2 . To obtain the true value of σ 2 , σM2 needs to be multiplied by a correction factor: ϕ=

σ2 . 2 σM

(19)

Hence, the correction factor ϕ at various SNR can be computed by exploiting the skewness of the Rician distribution. Figure 1 shows the relationship between the correction factor ϕ and the skewness γ of a Rician distribution. A lookup table can be created to represent the relationship between ϕ and γ for a Rician distribution from low to high SNR (A/σ ) (i.e. from Rayleigh to Gaussian) by varying the value of A from zero to a high value and by keeping σ constant. The corresponding values for ϕ and γ can then be computed by substituting the value of A and σ in equations (17)–(19). Now, for every γ , a corresponding value for ϕ can be found from the lookup table. A polynomial model that computes an approximation of the lookup table is given in the appendix. Now, to compute the noise variance in an MR image, the sample skewness and sample variance are computed for each pixel (i, j ) within a neighborhood . The correction factor for each pixel can then be found from the lookup table and the noise variance for each pixel is estimated using equation (16). Following a similar reasoning as in the previous section, the noise variance can now be estimated as the mode of all local estimates of noise variance around each pixel, which can be written as

 , (20) σ2 = mode σ2 S

L i,j

where σS2 is the estimate of the noise variance in the image and σL2 i,j is the estimated noise variance for a neighborhood around pixel (i, j ). 5. Experimental results Experiments were conducted on both real and synthetic MR images with and without the background region. The window size used in our experiments for local estimates was 9 × 9. For simulations we used the standard MR image phantom of the brain obtained from the Brainweb database (Cocosco et al 1997) and a cardiac MR image both with intensity values in the range 0–255. The cardiac image, in contrast to the brain image, contained no background areas. Both images were artificially corrupted with Rician noise with the noise level σ as 50. Figure 2 shows the distribution of the local estimates of σ computed using ML and skewness-based methods for both the brain and the cardiac images. It can be observed from the distributions that their mode is almost equal to the true standard deviation of the noise.

N446

J Rajan et al 2500 2000 1500 1000 500 0

0

10

20

30

40

50

60

70

80

90

100

110

120

130

140

70

80

90

100

110

120

130

140

Distribution of local estimates of σ

(a)

(b) 1200 1000 800 600 400 200 0

(c)

0

10

20

30

40

50

60

Distribution of local estimates of σ

(d) Figure 2. Noisy MR images of the brain (a) and heart (c), both corrupted with Rician noise with σ = 50. On the right, the distribution of local estimates of σ computed using both ML (solid red line) and skewness-based methods (dashed blue line) are shown for both the brain (b) and the cardiac (d) images.

This experiment demonstrates the independence of the proposed estimators to the image background. Now, to show the reliability of the proposed methods on both low and high SNR, we conducted the experiments on the cardiac MR image after corrupting the image with Rician noise with σ ranging from 10 to 100. Results of this experiment are shown in figure 3. The mean of 100 experiments divided by the actual value of σ is depicted. The value closer to 1 is the best estimate. The results are compared with the estimator given in equation (9). Since, at high SNR the Rician PDF approaches a Gaussian PDF, the estimation based on equation (9) will be closer to the true σ . However, as SNR drops, as expected, it can be seen from the figure that the Gaussian assumption introduces a bias. It can also be observed from the graph that the proposed estimators are significantly less biased especially for low values of the SNR. Even though both the proposed estimators are independent of the image background or SNR, the proposed estimator based on the local ML estimation of σ is closer to the true σ than the one based on the skewness. However, the time complexity of the skewness-based method is less than that of the ML-based method. The estimation method based on the skewness is almost ten times faster than the ML-based method. For the experiments on real data we used the MR image of a cherry tomato (with a mean intensity of around 10 000). A set of MR images was reconstructed by averaging 1 to 12 images. Averaging was done in the complex k-space. The resulting noise variance

Estimated σ / Actual σ

Noise measurement from magnitude MRI using estimates of variance and skewness

N447 Eq. (9) Eq. (13) Eq. (20)

1.6 1.4 1.2 1 0.8 0.6 0.4 0.2 10

20

30

40

50

σ

60

70

80

90

100

Estimated σ * sqrt(n averages)

Figure 3. The estimated σ /actual σ for various values of σ ranging from 10 to 100. This experiment was done on the cardiac MR image. The red and blue lines correspond to the proposed methods using ML and skewness, respectively, and the green line corresponds to the estimator given in equation (9). The shaded area shows the standard deviation for 100 experiments. 1500 Eq. (13) Eq. (20) 1000

500

0

0

2

3

4

5

6

7

8

9

10

11

12

number of averages

Figure 4. 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.

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 1/ n. Since the experimental setup for all acquisitions √ was identical except for averaging, the estimated noise standard deviation  σ multiplied by n is expected to be constant as a function of n. It can be seen from figure 4 that the proposed method exhibits this property. 6. Conclusion Two different approaches based on the local estimates of the variance and skewness were proposed to address the estimation of the noise from MR images when the Rayleigh background data or high SNR image regions are not available. The estimation of σ based on local ML is slightly more accurate than the method based on the skewness. However, the time complexity of the ML-based method is significantly greater than that of the skewness-based method. Experimental results on synthetic and real MR images show the reliability of the proposed methods. Even though the proposed methods are based on the noise characteristics of magnitude image data acquired with single coil MRI, they can, under certain conditions, be extended for multiple coil imaging methods like SENSE.

N448

J Rajan et al Table A1. Values of the coefficients in equation (A.1).

Coefficient p0 p1 p2 p3 p4 p5 p6 p7 p8 p9

Value 1 2.89 −7.29 × 10 1.16 × 103 −9.84 × 103 4.78 × 104 −1.38 × 105 2.31 × 105 −2.09 × 105 7.86 × 104

Acknowledgments This work was financially supported by the Inter-University Attraction Poles Program 6-38 of the Belgian Science Policy, and by the SBO-project QUANTIVIAM (060819) of the Institute for the Promotion of Innovation through Science and Technology in Flanders (IWTVlaanderen). Appendix. Expression for the lookup table The lookup table mentioned in section 4 can be approximated by the following polynomial equation: ϕ(γ ) =

9 

pi γ i ,

(A.1)

i=0

where ϕ is the correction factor and γ is the skewness. The values for the coefficients (with 95% confidence interval) are given in table A1. References Aja-Fern´andez S et al 2008 Noise and signal estimation in magnitude MRI and Rician distributed images : a LMMSE approach IEEE Trans. Imag. Proc. 17 1383–98 Aja-Fern´andez S et al 2009 Noise estimation in single and multiple coil magnetic resonance data based on statistical models Magn. Reson. Imaging 27 1397–409 Aja-Fern´andez S et al 2010 About the background distribution in MR data: a local variance study Magn. Reson. Imaging 28 739–52 Basu S et al 2006 Rician noise removal in diffusion tensor MRI Medical Image Computing and Computer-Assisted Intervention—MICCAI 2006 LNCS 4190 pp 117–25 Brummer M E et al 1993 Automatic detection of brain contours in MRI data sets IEEE Trans. Med. Imaging 12 153–68 Cocosco C A et al 1997 Brainweb: online interface to a 3D MRI simulated brain database NeuroImage 5 S425 http://www.bic.mni.mcgill.ca/brainweb/ Henkelman R M 1985 Measurement of signal intensities in the presence of noise in MR images Med. Phys. 12 232–3 Koay C et al 2009 Simultaneous identification of noise and estimation of noise standard deviation in MRI Proc. Int. Soc. Mag. Reson. Med. (Honolulu, HI) vol 17 p 4691 Nowak R D 1999 Wavelet based Rician noise removal for magnetic resonance images IEEE Trans. Image Processing 10 1408–19

Noise measurement from magnitude MRI using estimates of variance and skewness

N449

Sijbers J and Dekker A J den 2004 Maximum likelihood estimation of signal amplitude and noise variance from MR data Magn. Reson. Med. 51 586–94 Sijbers J et al 1998 Maximum likelihood estimation of Rician distribution parameters IEEE Trans. Med. Imag. 17 357–61 Sijbers J et al 2007 Automatic estimation of the noise variance from the histogram of a magnetic resonance image Phys. Med. Biol. 52 1335–48 Zhang Y et al 2001 Segmentation of brain MR images through a hidden Markov random field model and the expectation maximization algorithm IEEE Trans. Med. Imaging 20 45–57