Not incorporating this knowledge leads inevitably to biased results, in partic- ... Key words: adaptive, anisotropic, diffusion filter, Rice distribution, Maximum.
Adaptive anisotropic noise filtering for magnitude MR data J. Sijbers1 , A. J. den Dekker1 , A. Van der Linden2 , M. Verhoye2 , and D. Van Dyck1 Vision Lab1 and Bio-Imaging Lab2 Department of Physics, University of Antwerp, Groenenborgerlaan 171, 2020 Antwerpen, Belgium http://www.ruca.ua.ac.be/visielab/sijbers
Magnetic Resonance Imaging, Vol. 17, Nr. 10, p. 1533-1539, (1999)
Abstract Conventional noise filtering schemes applied to magnitude magnetic resonance (MR) images tacitly assume Gauss distributed noise. Magnitude MR data, however, are Rice distributed. Not incorporating this knowledge leads inevitably to biased results, in particular when applying such filters in regions with low signal-to-noise ratio. In this work, we show how the Rice data probability distribution can be incorporated so as to construct a noise filter that is far less biased. Key words: adaptive, anisotropic, diffusion filter, Rice distribution, Maximum Likelihood estimation
1
Introduction
In the past, many efforts were devoted to the reduction of noise, or equivalently, the estimation of the underlying, noiseless signal. At present, the most common noise reduction technique is still time averaging. It has the major advantage that the signal-to-noise ratio (SNR) is increased while the spatial resolution remains intact, provided the imaging process is stationary. Since this approach is quite time consuming and hence not feasible in case of conventional 3D MR imaging, time averaging is usually replaced by spatial averaging. However, there are two important (and often taken for granted) assumptions in applying spatial averaging to estimate the underlying, noiseless signal: Stationarity The first assumption is stationarity of the pixels used for the estimation of the underlying signal. Fortunately, Wang and Lei [1] justified the heuristics that an MR image can be regarded as consisting of many regions in which the signal is stationary and ergodic in the mean and variance. The mean ergodic property justifies signal estimation within those regions to estimate the underlying, noiseless signal value. The main problem, however, is to find these stationary regions. The problem of finding the proper area for local signal estimation is partly solved by the so-called sigma filters, where spatial averaging is halted as soon as strong image gradients are detected [2, 3]. A more elaborated approach can be found in the work of Perona et al. , who developed an anisotropic diffusion scheme for image data, where the 1
diffusion and flow functions are guided by local gradient strengths in different directions [4]. The properties of Perona’s filter are: efficient noise removal in homogeneous regions, preservation of object boundaries, and edge sharpening. This filtering technique was successfully applied to 2D and 3D MR data by Gerig et al [5]. Although the performance of the noise filter is excellent, the underlying image model is piecewise constant or slowly varying. As a result, the edge sharpening causes a region with a constant gray value slope to be broken up in constant plateaus. The edge sharpening property is not retained in the approach of Yang et al [6]. In the present paper, Yang’s version of the adaptive, anisotropic diffusion filter is adopted and modified for the filtering of magnitude MR data. Data PDF The second assumption is made on the probability density function (PDF) of the data points at disposal. By computing the spatial average, the image noise is assumed to be zero mean and Gaussian distributed. Indeed, only in that case the expectation value of the spatial average equals the noiseless, underlying signal. Unfortunately, this assumption is not met when processing magnitude MR data as these can be shown to be Rice distributed [7]. As the expectation value of the spatial average of Rice distributed data points within a stationary region does not equal the underlying signal, filtering methods based on this spatial averaging yield biased results. In this paper, it is shown how the bias can be reduced by incorporating the Rice distribution into the noise filtering procedure.
2
Methods
2.1
The probability density function of magnitude MR data
The raw, complex valued MR data are acquired in the Fourier domain. It is generally known that the noise contributions of these data are additive and uncorrelated, characterized by a zero mean Gauss PDF [1, 8]. After the inverse Fourier transform, the complex data are still Gauss distributed due to the linearity and orthogonality of the Fourier transform 1 [9]. The transformation to a magnitude image, however, changes the Gauss PDF of the data into a Rice PDF [10], given by: 2 fM M − M 2 +f 2 pM (M |f ) = 2 e 2σ I0 (1) σ σ2 where I0 is the zeroth order modified Bessel function of the first kind and M denotes the pixel variable of the magnitude image. The underlying (noiseless) signal is represented by f . At high SNR, the Rice distribution approaches a Gaussian distribution, while at low SNR, it approaches a Rayleigh distribution.
2.2
Signal estimation
In this section, the anisotropic, adaptive diffusion filter proposed by Yang is briefly described, after which a scheme is presented based on the incorporation of the Rice distribution into the estimation procedure. 1
It will be assumed that the MR signals are sampled on a uniform grid in K-space. Furthermore, the variance of the noise will be assumed to be equal for each raw data point.
2
2.2.1
Anisotropic adaptive diffusion filter
Suppose g(~r) denotes a 3D image, where ~r = (x1 , x2 , x3 ) is a three-dimensional position vector. The filtering process consists of convolving g(~r) with a Gauss kernel h(~r) of which the shape is pointwise adapted to the local structure within a neighborhood Ω. The resulting filtered function g 0 (~r) can be written as follows: RRR
h(~r0 − ~r)g(~r)d~r
Ω
0
g (r~0 ) =
RRR
,
h(~r0 − ~r)d~r
(2)
Ω
where
3 1X ((~r − ~r0 ) .~ni )2 h(~r − ~r0 ) = exp − 2 i=1 σi2 (~r0 )
"
#
.
(3)
In Eq. (3), the vectors {~ni } correspond to the main axes of the local structure. Thereby, ~n1 represents the dominant spatial orientation. The level contours of an oriented pattern consist of parallel lines. The power spectrum of such a pattern clusters along a line through the origine of the Fourier domain. Hence, by computing the direction of that line, the orientation of the pattern in the spatial domain can be found. Therefore, the 3×3 second moment matrix R of |G(~u)|2 is computed, where |G(~u)|2 denotes the Fourier power spectrum of g(~r). The coefficients of R can as well be calculated in the spatial domain: ZZZ
Rij
ui uj |G(~u)|2 du1 du2 du3
=
(4)
Ω
=
1 (2π)3
ZZZ Ω
∂g ∂xi
∂g ∂xj
!
dx1 dx2 dx3
,
(5)
where i, j = 1, 2, 3. R is positive semi-definite and Hermitian, hence having only positive eigenvalues. The direction of the eigenvector ~n1 , corresponding to the smallest eigenvalue, say λ1 , determines the main direction of the pattern in the neighborhood Ω of the spatial domain. Alternatively, the three eigenvalues λ1 , λ2 , and λ3 , with λ1 ≤ λ2 ≤ λ3 , determine the relative orientation of the pattern in the respective directions. The shape of the kernel h [see Eq. (3)] is controlled by the standard deviations σ1 , σ2 , and σ3 , which are functions of the local gradient strength in the respective directions. From the eigenvalues, anisotropy measures are derived, a12 =
λ 2 − λ1 P λi
and a13 =
i
λ 3 − λ1 P λi
,
(6)
i
which are used to design the kernel’s standard deviations. These standard deviations should be large along the main direction(s) of the pattern, such that the data are only smoothed in homogeneous regions and along instead of across edge surfaces. In addition, corners should be preserved during filtering. A corner is identified as a condition where the pattern is relative isotropic (a12 ≈ 0 ; a13 ≈ 0) while the local gradient strength |∇g(~r)|2 is large. Therefore, a spatial dependent corner strength C is defined as: C(~r) = (1 − a12 − a13 )|∇g(~r)|2
.
(7)
From Eq. (6) and Eq. (7), the standard deviations are designed as follows: σ1 (~r) =
σ 1 + C(~r)
σ2 (~r) =
σ(1 − 2a12 ) 1 + C(~r) 3
σ3 (~r) =
σ(1 − a12 − a13 ) 1 + C(~r)
,
(8)
Notice that, whenever a corner is detected (i.e., when C is large), the standard deviations of the Gaussian kernel become small such that smoothing is halted. In case of anisotropy of the local structure, smoothing is halted in the respective direction. The above filter actually fits a constant plane to the data in Ω where the weights of the data points are governed by the kernel h. This filter was previously shown to be accurate in reducing noise while retaining edge information. 2.2.2
Incorporation of the Rice distribution
Despite of the ability of the above filter to strongly reduce the noise while leaving edges intact, a bias is introduced when applied to magnitude MR data. This bias, which is due to the inaccurate assumption of Gaussian distributed data, was previously shown to be increasing with decreasing SNR [11]. In order to reduce this bias, we propose a modified version of the above filter in which the full knowledge of the data PDF is exploited. Thereby, the underlying signal model is assumed to be piecewise linear. Consider a J-dimensional data set consisting of N Rice distributed data points {Mi }, where i = 1, 2, ..., N . For each data point Mi , a local environment Ω is chosen, centered about Mi . Then, within each region Ω a J − 1 dimensional hypersurface f is fitted to the data points within that region. Thereby, the hypersurface f is parameterized by J + 1 parameters {βj }: f ({βj }) = β0 +
J X
βj xj
(9)
j=1
The parameters of the hypersurface are then estimated through maximization of the likelihood function with respect to the parameter set {βj }: (
)
{βˆj } = arg max (log L)
(10)
{βj }
where L({βj }|{Mi }) =
N Y Mi i=1
!
Mi2 + fi2 ({βj }) Mi fi ({βj }) exp − I0 σ2 2σ 2 σ2
(11)
with σ 2 denoting the noise variance, which can be measured accurately by using data points from the background or homogeneous regions [12] or from two realizations of the same image [13–15]. Maximizing the likelihood function (11) is equivalent to maximizing log L: log L =
N X i=1
"
N X Mi2 Mi fi ({βj })Mi log 2 − 2 + I0 σ 2σ σ2 i=1
#
"
f 2 ({βj }) − i 2 2σ
#
(12)
where the terms dependent on the unknown parameters {βj } are seperated. The Maximum Likelihood estimate is the global maximum of log L. For each data point, the parameter set {βj }, which represent the local intensity gradients, is estimated after which the elements of the matrix R are evaluated [see Eq. (5)]. After computing the eigenvalues {λi }, the weights are determined from Eq. (3) and Eqs. (6-8). Finally, a constant plane is fitted to the data using Maximum Likelihood estimation where use is made of the Rice distribution and the computed weights.
4
3
Results and discussion
Simulation experiments were set up to test the performance of the proposed filter. Thereby, various input functions (one of which is shown in Fig. 1(a)) were polluted with Rice distributed noise (see Fig. 1(b)). Both filters, Yang’s adaptive anisotropic diffusion filter and the presented version, were applied to the noisy signal and their performance was compared in terms of the Mean Squared Error (MSE), defined as: v uP 2 u N u i=1 fˆi − fi MSE = t PN 2 i=1 fi
,
(13)
where fˆi denotes the estimated signal at position i. This performance comparison was done as a function of the SNR, defined as: q
SNR =
1 N
PN
2 i=1 fi
σ
(14)
The performance of both filters was tested by computing the MSE using Eq. (13) for a large number of noise variances. Fig. 3 shows the results of the MSE’s from the two filters as a function of the SNR as defined in Eq. (14). At high SNR, the difference in performance of both filters is generally small. However, with decreasing SNR, the MSE of Yang’s filter becomes significantly larger than the filter incorporating the Rice distribution of the data. Next to the simulation experiments, the filter was tested on real magnitude MR data (Fig. 4). In Fig. 4(c) and 4(b), the result of the Gauss based and the presented Rice based filter applied to a magnitude MR image of a cucumber are shown. One can observe strongly reduced noise while the image details are well retained. As explained above, the difference between Gaussian and Rician based filters is most apparent in regions with low SNR. In particular, the background area (where the SNR approaches zero) is effectively put to zero by the last filter, while Gaussian based filters tend to set this area to a constant, non-zero greylevel, hence losing contrast.
4
Conclusions
Filters applied to magnitude MR data should incorporate the probability density function that characterizes these data. In the case under concern, this is the Rice probability density function. Full use is made of the characteristics of magnitude MR data by incorporating the Rice distribution into the Maximum Likelihood estimation of the filter parameters. In addition to good noise filtering performance, spatial resolution is retained due to the adaptive and anisotropic characteristics of the filter.
Acknowledgements This work was financially supported by the Flemish institute for the promotion of scientific and technological research in the industry (IWT), Brussels.
5
References [1] Wang Y. and Lei T. Statistical analysis of mr imaging and its applications in image modeling. In Proceedings of the IEEE International Conference on Image Processing and Neural Networks, volume I, pages 866–870, 1994. [2] Saint-Marc P., Chen J. S., and Medioni G. Adaptive smoothing: A general tool for early vision. IEEE Transactions on Pattern Analysis and Machine Intelligence, 12:629–639, 1990. [3] Ying K., Clymer B. D., and Schmalbrock P. Adaptive filtering for high resolution magnetic resonance images. Journal of Magnetic Resonance Imaging, 6(2):367–377, 1996. [4] Perona P. and Malik J. Scale space and edge detection using anisotropic diffusion. IEEE Transactions on Pattern Analysis and Machine Intelligence, 12(7):629–639, 1990. [5] Gerig G., Kubler O., Kikinis R., and Jolesz F. A. Nonlinear anisotropic filtering of MRI data. IEEE Transactions on Medical Imaging, 11(2):221–232, 1992. [6] Yang G. Z., Burger P., Firmin D. N., and Underwood S. R. Structure adaptive anisotropic filtering for magnetic resonance images. Lecture Notes in Computer Science, 970:384–391, 1995. [7] Andersen A. H. and Kirsch J. E. Analysis of noise in phase contrast MR imaging. Medical Physics, 23(6):857–869, 1996. [8] Breiman L. Probability. Addison-Wesley, Reading, Massachusetts, 1968. [9] Ross S. A first course in probability. Collier Macmillan Publishers, New York, 1976. [10] Papoulis A. Probability, Random Variables and Stochastic Processes. McGraw-Hill, Tokyo, Japan, 2nd edition, 1984. [11] Sijbers J., den Dekker A. J., Scheunders P., and Van Dyck D. Maximum likelihood estimation of Rician distribution parameters. IEEE Transactions on Medical Imaging, 17(3):357–361, June 1998. [12] Henkelman R. M. Measurement of signal intensities in the presence of noise in MR images. Medical Physics, 12(2):232–233, 1985. [13] Sijbers J., Scheunders P., Bonnet N., Van Dyck D., and Raman E. Quantification and improvement of the signal-to-noise ratio in a magnetic resonance image acquistion procedure. Magnetic Resonance Imaging, 14(10):1157–1163, 1996. [14] Sijbers J., den Dekker A. J., Verhoye M., Van Audekerke J., and Van Dyck D. Estimation of noise from magnitude MR images. Magnetic Resonance Imaging, 16(1):87–90, 1998. [15] Sijbers J., den Dekker A. J., Raman E., and Van Dyck D. Parameter estimation from magnitude MR images. International Journal of Imaging Systems and Technology, 10(2):109– 114, 1999.
6
45
40
35
30
f(x)
25
20
15
10
5
0 0
50
100
150
200
250 x
300
350
400
450
500
350
400
450
500
(a) Input function 60
50
N(x)
40
30
20
10
0 0
50
100
150
200
250
300
x
(b) Noise added
Figure 1: Simulation experiment: Original input function (a) polluted with Rice distributed noise (b).
7
45 original function Gauss expectation values 40
35
30
E[M]
25
20
15
10
5
0 0
50
100
150
200
250 M
300
350
400
450
500
Figure 2: The original function along with the expectation values of the Gaussian based filters. (The expectation values for the Rician based filter equal the values of the original input function.)
Figure 3: The Mean Squared Error (MSE) computed as a function of the signal-to-noise ratio (SNR), defined in Eq. (14).
8
(a) Unprocessed MR image
(b) Gaussian based filtering
(c) Rician based filtering
Figure 4: Result of the adaptive, anisotropic diffusion filter incorporating the Gauss (b) and Rice distribution (c), applied to a magnitude MR image of a cucumber.
9
