20092009 13th13th International Irish Machine Machine Vision Vision and and Image Image Processing Processing Conference Conference

Denoising Magnetic Resonance Images using Fourth Order Complex Diffusion Jeny Rajan, Ben Jeurissen, Jan Sijbers IBBT Vision Lab University of Antwerp Wilrijk, 2610 Antwerp, Belgium jeny.rajan, ben.jeurissen, [email protected]

efficients by filtering the observed coefficients drj (k) at any level j, spatial position k, and wavelet orientation r [5]. In PDE based methods, the basic idea is to deform an image, a curve, or a surface with a PDE and obtain the expected results as a solution to this equation [2]. Both wavelet and PDE based methods have their own merits and demerits, even though wavelet based methods seem to be more flexible for MR denoising. A comparative study of both methods is not in the scope of this work. Interested readers can refer to [5],[4]-[8]. In the proposed method, we used a PDE based approach for denoising MR images. PDE based methods, especially anisotropic diffusion (Perona-Malik), have proved to be particularly effective in pre-filtering MR images [5]. A variant of standard anisotropic diffusion method was extended by Yang [15] using both the local intensity orientation and an anisotropic measure of level counters, instead of utilizing local gradients to control the anisotropism of the filters. Even though anisotropic diffusion seems attractive, it assumes image noise to be Gaussian distributed. When processing magnitude MR data, a Gaussian assumption for image noise is not acceptable as it can be shown to be Rice distributed [10]. Using anisotropic diffusion or its variants can generate a bias in the magnitude MR data (which increases with decreasing SNR). Sijbers et al [12] proposed an adaptive anisotropic diffusion method to tackle this problem. However, several papers [3],[16],[14],[7] have noted that 2nd order anisotropic diffusions are ill posed in the sense that images close to each other are likely to diverge during the diffusion process. Solutions like 4th order PDEs and complex PDEs are suggested to solve this problem. In this paper, we propose a 4th order, non linear, complex diffusion for denoising MR images. Complex diffusion is studied as an efficient tool for image denoising, which certainly has an edge over the ordinary PDEs. Simulation and real experiments based on the proposed 4th order complex PDE show that its performance with respect to the Peak signal-to-noise ratio (PSNR) and the Structural Similarity Index Matrix (SSIM) is superior compared to 2nd

Abstract Complex diffusion is a comparatively new Partial Differential Equations (PDE) based method introduced for removing noise from images. The efficiency of 2nd order complex diffusion for image denoising is already proved by many researchers. 2nd order non linear complex diffusion can behave like 3rd and 4th order real PDEs enabling a variety of new options with standard 2nd order numerical schemes. Extending 2nd order non linear complex diffusion to 4th order can produce a much better result. In this paper we present a 4th order non linear complex diffusion. Our experimental results show that this 4th order complex PDE is a good choice for denoising Magnetic Resonance images. The efficacy of the algorithm is demonstrated on both simulated and real Magnetic Resonance images.

1. Introduction Denoising is a crucial step required for correct interpretation of Magnetic Resonance (MR) data. Depending on specific diagnostic tasks, high spatial resolution and high contrast may be required for MR images, whereas for image processing applications, a high SNR is usually necessary because most of the algorithms are very sensitive to noise [5]. It has been shown that noise in Diffusion Weighted (DW) MR Images introduces errors in the estimation and sorting of the diffusion tensor eigenvalues and the derived anisotropy measures [6]. The quality of the estimated tensor field in Diffusion Tensor Imaging (DTI) also depends on the noise level. Hence, it is clear that an efficient denoising algorithm can have a significant impact on MR image processing tasks. Various approaches for removing noise from MRI have been proposed. Most of these methods are based on Wavelets or Partial Differential Equations (PDEs) [2][1],[15]. The goal of wavelet domain filtering is to obtain a better estimate of the noise free image wavelet co-

978-0-7695-3796-2/09 $26.00 © 2009 IEEE DOI 10.1109/IMVIP.2009.29

K. Kannan Medical Imaging Research Group NeST (P) Ltd Trivandrum 81, India [email protected]

123

order complex diffusion. This paper is organized as follows. Section 2 discuss the bias removal in MR images and explains the proposed 4th order complex PDE. In Section 3, we will discuss the results of the proposed method. Finally, conclusions and remarks are drawn in Section 4.

the maximum-minimum principle, it preserves other desirable mathematical and perceptual properties. When an image is processed with complex diffusion, we will get low frequency components (plateaus) of the image in the real plane and high frequency components (edges) in the imaginary plane [9]. The components in the real and imaginary plane are almost equivalent to that of the image convolved with a Gaussian and Laplacian of Gaussian (LOG) at various scales respectively. Here we propose a 4th order non linear complex diffusion method, which is an improvement over its 2nd order counterpart in terms of preserving edges. The method is based on:

2. Proposed Method The proposed 4th order complex PDE is an improvement over the one proposed by Gilboa et al [3]. Since the method assumes the noise to be Gaussian distributed, we first did a bias reduction from the squared magnitude MR image as proposed in [8].

I t = −∇2 [(c(=(I))∇2 I]

where =(.) takes the imaginary part and the diffusivity function c(.)is defined as

2.1. Bias Removal

c(=(I t )) =

Nowak [8] proposed a method to remove the bias, making use of the properties of the squared magnitude image. If M is the Rician distributed magnitude image, S is the unknown, noiseless image of interest, then using the moments of the non central chi square distribution, the mean and variance of the squared magnitude image M 2 are   E M 2 = S 2 + 2σ 2 (1) Var[M 2 ]

=

4σ 2 (S 2 + σ 2 )

(5)

eiθ

(6)

2 1 + ( =(I) kθ )

and with initial conditions: <(I t=0 ) = I

;

=(I t=0 ) = 0

(7)

In Eq. (6) k is the threshold parameter, θ is the phase angle (it can also determine the direction of diffusion) and I is the original image. Eq. (5) can be solved numerically using an iterative approach as given below. First, we calculate the Laplacian of the image I as follows:

(2)

It is important to realize that Eq. (1) has a fixed, signal independent bias [11]. This property can be exploited to reduce the bias by subtracting an estimate of the bias from each pixel in the squared magnitude image as follows:

t ∇2 Ii,j =

1 t t t t t (I + Ii−1,j + Ii,j+1 + Ii,j−1 ) − Ii,j (8) 4 i+1,j

with symmetric boundary conditions c2 = M 2 − 2b M σ2 where σ b2 =

1 2 M 2

,

(3)

t t t t I−1,j = I0,j , IM +1,j = IM,j

,

j = 1, 2, · · · , N (9)

and .

(4)

t t t t Ii,−1 = Ii,0 , Ii,N +1 = Ii,N

Even though this method does not remove the bias in the magnitude image completely, there will be a clear contrast enhancement.

,

i = 1, 2, · · · , M (10)

where M × N is the total image space. Once the Laplacian of the image I is computed, the function c(=(I t ))∇2 I t has to be computed. Let

2.2. Denoising using 4th order complex diffusion

Rt = c(=(I t ))∇2 I t

(11)

Then ∇2 Rt can be computed as

The concept of complex diffusion in image processing was introduced by Gilboa et al. [3] as an alternative to 2nd order anisotropic diffusion, which introduces blocky effects in images while processing. This blocky effect is inherent in the nature of ordinary second order equations; it can be avoided by using complex diffusion. Complex diffusion is derived by combining the standard diffusion equation with the free Schr¨odinger equation [3]. Even though it violates

∇2 R t =

1 t t t t t (R +Ri−1,j +Ri,j+1 +Ri,j−1 )−Ri,j (12) 4 i+1,j

Finally, the numerical approximation to the differential Eq. (5) can be written as t+1 t t Ii,j = Ii,j − 4t∇2 Ri,j

(13)

In the proposed method we used a step size (4t) of 0.1.

124

3. Results and Discussions

Fig. 2 and Fig. 3 shows the plot of PSNR and SSIM, respectively, against various values of the noise standard deviation. For these experiments we used the same image shown in Fig. 1(a). From Fig. 1, it can be observed that the proposed method produces images with best PSNR and MSSIM when compared with 2nd order complex diffusion or 4th order real PDE. Next to the simulation experiments, the proposed filtering method was applied to real data as well. Fig. 4 shows the performance of the proposed method on a T1 weighted MR image (MR Brain, 1.5 T). Fig. 4(b) and Fig. ?? show the results after 25 and 50 iterations, respectively. It can be seen from both figures that there is a clear reduction of noise. In Fig. 5, we tested the algorithm on a diffusion weighted image of a rat brain. These images are of low resolution and the noise level is generally high. The result shows that there is a significant reduction in image noise while the image details are well retained.

Experiments were done on both simulated and real MR images. For all experiments we set k = 1.5, θ = 355◦ , and the step size ∆t = 0.1. The selection of k and θ was based on the parameter settings discussed in [9]. First, an phantom image, shown in Fig. 1(a), was polluted with Rician noise (PSNR -3.87), the result of which is shown in Fig. 1(b). After reducing the bias caused by the Rician distributed noise (using Eq. (3)), the image was processed with the 4th order as well as with the 2nd order complex PDE based method. The results are shown in Fig. 1(c) and Fig. 1(d), respectively. Next, the proposed noise reduction method was quantitatively evaluated in terms of the peak signal-to-noise ratio (PSNR) and the Structural Similarity Index Matrix (SSIM). SSIM [13] is used to evaluate the overall image quality and is in the range 0 to 1. The SSIM works as follows. Let x and y be two non negative images, where one image has a perfect quality (groundtruth). Then, the SSIM can serve as a quantitative measure of the similarity of the second image. The system separates the task of similarity measurement into three comparisons: luminance, contrast and structure. It can be defined as SSIM (x, y) =

(2µx µy + C1 )(2σxy + C2 ) (µ2x + µ2y + C1 )(σx2 + σy2 + C2 )

4. Conclusions In this paper, we proposed a 4th order non linear complex diffusion for denoising MR images. The behavior of the proposed 4th order complex diffusion has to be further explored, but the initial results are promising. Experiments with the proposed method shows that there is a clear improvement in image quality compared to its 2nd order counterpart, and is proved in both real and simulated images. We plan to further improve the proposed method my incorporating the rician behavior of noise in it.

(14)

where µx and µy are the estimated mean intensity along x and y directions and σx and σy are the standard deviations respectively. σxy can be estimated as N

σxy =

1 X (xi − µx )(yi − µy ) N − 1 i−1

(15)

References [1] J. Aelterman, B. Goossens, A. Pizurica, and W. Philips. Removal of correlated rician noise in magnetic resonance imaging. Proceedings of EUSIPCO 2008, 2008. [2] R. Fedkiw, G. Sapiro, and C. W. Shu. Shock capturing, level sets, and pde based methods in computer vision and image processing: A review of osher’s contributions. Journal of Computational Physics, 185(2):309–341, March 2003. [3] G. Gilboa, N. Sochen, and Y. Y. Zeevi. Image enhancement and denoising by complex diffusion processes. IEEE Trans. on PAMI, 26(8):1020–1036, August 2004. [4] L. Jiang and W. Yang. Adaptive magnetic resonance image denoising using mixture model and wavelet shrinkage. Proceedings of 8th Digital Image Computing : Techniques and Applications, pages 831–838, 2003. [5] L. Landini, M. Lombardi, and A. Benassi. Noise filtering methods in mri. In Advanced Image Processing in Magnetic Resonance Imaging. CRC Press, 2005. [6] A. Leemans. Modelling and processing of diffusion tensor magnetic resonance images for improved analysis of brain connectivity. 2006.

C1 and C2 in Eq. (14) are constants and the values are given as C1 = (K1 L)2 (16) and C2 = (K2 L)2

(17)

where K1 , K2  1 is a small constant and L is the dynamic range of the pixel values (255 for 8 bit gray scale images)[13]. The PSNR, given in decibel (dB), is the ratio between the maximum possible power of a signal and the power of the corrupting noise:   M AXx2 P SN R = 10 log10 (18) M SE with M AXx the maximum possible value of x and MSE de mean squared error between x and y.

125

(a)

(b) Figure 3 – SSIM Computed against different noise sigma values.

(c)

(d)

Figure 1 – Comparative study of the proposed method with 2nd order Complex Diffusion (a) Original Image (b) Image corrupted with Rician noise (PSNR -3.87, MSSIM 0.2095) (c) Processed with 2nd order Complex Diffusion (PSNR 3.67, SSIM 0.3471) (d) Processed with (proposed) 4th Order Complex Diffusion (PSNR 6.67, SSIM 0.4890).

(a) Original

(b) 25 iterations

Figure 4 – MR denoising with 4th order non-linear complex diffusion. (a) Original MR brain image (1.5 T), (b) processed with proposed method with 25 iterations.

Figure 2 – PSNR Computed against different noise sigma values.

[7] M. Lysaker, A. Lundervold, and X.-C. Tai. Noise removal using fourth-order partial differential equation with applications to medical magnetic resonance images in space and time. IEEE Trans Imag Proc, 12(12):1579–1590, 2003. [8] R. D. Nowak. Wavelet based Rician noise removal for magnetic resonance images. IEEE Trans. Image Processing, 10(8):1408–1419, 1999. [9] J. Rajan, K. Kannan, and M. R. Kaimal. Smoothening and sharpening effects of theta in complex diffusion for image processing. Proceedings of ICAPR 2009, pages 325–328, February 2009. [10] J. Sijbers, A. J. den Dekker, E. Raman, and D. Van Dyck. Parameter estimation from magnitude MR images. Int J Imag Syst Tech, 10(2):109–114, 1999.

(a) Slice 1

(b) N = 50

(c) N = 200

(d) Slice 2

(e) N = 50

(f) N = 200

Figure 5 – Two diffusion weighted MR slices of the rat brain (a) and (d) are the original (corrupted with noise) slices. (b), (c) and (e), (f) are the processed images with the proposed method with iterations 50 and 200 respectively.

126

[11] J. Sijbers, A. J. den Dekker, P. Scheunders, and D. Van Dyck. Maximum likelihood estimation of Rician distribution parameters. IEEE Trans. Med. Imag., 17(3):357–361, June 1998. [12] J. Sijbers, A. J. den Dekker, M. Verhoye, A. Van der Linden, and D. Van Dyck. Adaptive anisotropic diffusion filter for magnitude MR data. In Proceedings of SPIE’99: Medical Imaging, volume 3661, pages 1418–1425, San Diego, CA, USA, February 1999. [13] Z. Wang, A. Bovik, H. R. Sheik, and E. P. Simoncelli. Image quality assessment: From error visibility to structural similarity. IEEE Trans. on Image Pocessing, 13(4):600–612, April 2004. [14] R. T. Whitaker and S. M. Pizer. A multiscale approach to non uniform diffusion. CVGIP: Image Understanding, 57(1):99–110, 1993. [15] G. Z. Yang, P. Burger, D. N. Firmin, and S. R. Underwood. Structure adaptive anisotropic filtering for magnetic resonance image enhancement. In Proceedings of CAIP: Computer Analysis of Images and Patterns, pages 384–391, 1995. [16] Y. L. You and M. Kaveh. Fourth order partial differential equations for noise removal. IEEE Trans. on Image Pocessing, 9(10):1723–1730, October 2000.

127