A New Adaptive Diffusive Function for Magnetic Resonance Imaging Denoising Based on Pixel Similarity.

Although there are many methods for image denoising, but partial differential equation (PDE) based denoising attracted much attention in the field of medical image processing such as magnetic resonance imaging (MRI). The main advantage of PDE-based denoising approach is laid in its ability to smooth image in a nonlinear way, which effectively removes the noise, as well as preserving edge through anisotropic diffusion controlled by the diffusive function. This function was first introduced by Perona and Malik (P-M) in their model. They proposed two functions that are most frequently used in PDE-based methods. Since these functions consider only the gradient information of a diffused pixel, they cannot remove noise in noisy images with low signal-to-noise (SNR). In this paper we propose a modified diffusive function with fractional power that is based on pixel similarity to improve P-M model for low SNR. We also will show that our proposed function will stabilize the P-M method. As experimental results show, our proposed function that is modified version of P-M function effectively improves the SNR and preserves edges more than P-M functions in low SNR.

INTRODUCTION

Magnetic resonance imaging (MRI) has very important role on current medical and research procedures. However, these images are normally corrupted by random noise from the measurement process. Therefore, removing noise in MRI is one of the major issues in the study of medical imaging. The traditional methods of image denoising use smoothing (such as average filters and Gaussian filtering) for achieving good results, but the edge information and texture may be lost.[1] However, medical images consist of many details and subtle features such as angular and edge, and the lack of edge information of image brings difficulties for subsequent analysis. From above, in the image restoration field, we extremely need an image restoration process technology which can remove the noise and keep edge information also. So, partial differential equation (PDE)-based methods are one of the best selections to MRI denoising.

In addition to PDE-based methods, there are many methods that are based on nonlocal means (NLM) filter[23] such as multi component NLM and its variants,[4] unbiased NLM,[5] adaptive NLM,[67] nonparametric principal component analysis NLM,[8] and MRI denoising using nonlocal PCA.[9]

The main drawback of the NLM-based algorithms is computational burden due to calculate weighted average for all pixels.[10] Although this complexity value is reduced in some modified version (by blocking image[6] or using PCA), but still the calculation volume is high. On the other hand, they cannot remove the noise and simultaneously preserve edges, as well as PDE-based methods, also. In these methods, oversmoothing is observed in some regions which results in a loss of edges and fine structures in the image.[11] In additional, PDE-based methods are one of the best selections to denoising smooth area and preserving edges simultaneously. These filtering are still widely used in the preprocessing of MRIs because of its efficiency and simplicity in implementation.

Several PDE-based denoising methods were proposed such as Perona and Malik (P-M) method.[12] They proposed a nonlinear equation which replaced isotropic diffusion expressed through a linear heat equation with an anisotropic diffusion (AD). On the other hand, AD is controlled by the diffusion coefficient (also called conductance or diffusive function). The diffusive function g tends to “0” while the absolute of image gradient (|∇u|) is large (edges) and “1” while |∇u| is small (flat regions). Therefore the P-M model can effectively remove the noise, as well as containing edge information.

Many revised AD models have been proposed done as stated by P-M diffusion model. Alvarez et al.[13] introduced a mean curvature motion approach for AD. This method performs the diffusion process only in the direction of minimum rate of gradient change. Therefore, the image will be smoothed in a direction parallel to edges, not across it. On the other hand, Catte et al.[14] showed that the P-M AD model is ill-posed, that is, very similar images could produce divergent solutions. To solve the problem correctly, they showed that the diffusion coefficient should be calculated from isotropically smoothed versions of the image.

For noisy images with a low signal-to-noise (SNR) ratio, however, noise variations may be comparable to or even greater than edge variations. In this case, the discrimination between edge and noise can be difficult and then the conventional AD usually fails to produce satisfactory results for low SNR images.[1516] Therefore, the kernel AD (KAD) is proposed to improve the performance of AD for low SNR images. In this method, the input space is mapped to a higher order feature space, and then discriminates between edges and noise in the feature space.[171819]

The classical P-M's AD model only considers the local gradient information of each pixel in the image. However, a high gradient magnitude is generally a good indication of edges, whereas a low gradient magnitude may not always point to nonedge regions or noise. For example, the important local details along with edges in the image may have low gradient. In order to retain fine details while removing noise, the local gray-level variance along with the gradient is added to a modified AD model.[20] According to above discussion, it is obvious that the gradient magnitude alone cannot be a good criterion to discriminate edge, smooth area, and noise. Because of this, Wang et al.[21] introduced a pixel similarity measurement into the PDE-based on mean curvature flow of curve. Using this measurement, each pixel can be classified in five categories: Edge, strong noise, median noise, weak noise, and smooth area. In this paper we utilize similarity measurement to determine edge, noise, and smooth area. In this work, we improve P-M diffusive function by pixel similarity. Our novel diffusive function can stable P-M model. Experimental results show that the proposed function can effectively improves the SNR and preserves both edges and fine details of objects in the restored image even for low SNR.

In Section 2, the AD is introduced. In this section the P-M model and KAD model that is appropriate for low SNR is studied. Section 3 is dedicated to the pixel similarity measurement. In Section 4, we present an improved diffusive function based on pixel similarity with fractional and adaptive power and then the stability of our function is analyzed. Section 5 introduce two stabled AD. The experimental and comparison results are described in Section 6.

ANISOTROPIC DIFFUSION MODEL

Perona and Malik Model

Diffusion algorithms remove noise from an image by modifying the image via aPDE. The diffusion driven by heat equation is isotropic and image features will be blurred by noise removing. One of the advanced properties of AD is to integrate prior knowledge of the image to the diffusion coefficient. Perona and Malik[12] first introduced the idea of nonlinear diffusion that is preferred within a smooth region to diffusion near an edge. They proposed the enhancement of noisy image I (x, y) by the solution of the following PDE:

in Ω×(0,T)      (1)

u(x, y; 0)=I(x, y) in Ω

on ∂Ω

where u (x, y;t) is the enhanced version of the image at time t, Δu denotes the local gradient of u, Ω ⊂ R2 is typically a square domain, ∂Ω means the edges of Ω, represents outward normal direction to ∂Ω and g: (0, +∞) → (0, 1) so called conductance or diffusive function is a nonnegative function of the magnitude of local gradient. This function is chosen to be close to zero where gradient is large, so that the diffusion is stopped across edges. g should be large in pixels with low gradient variation, so that the diffusion is maximal within smooth regions. Functions satisfying these assumptions are commonly called edge stopping functions, and examples of such a function were firstly proposed by Perona and Malik,[12] e.g.,

g1(s)=exp(–(s / k)2)      (2)

where k is the gradient magnitude threshold parameter that controls the rate of the diffusion and serves as a soft threshold between the image gradients that are attributed to noise and those attributed to edges. According to Perona and Malik, Tsiotsios and Petrou,[1222] the g1 favors high-contrast edges over low-contrast ones, while the g2 function favors wide regions over smaller ones.

In fact, above AD equation (Eq. 1) can be obtained by minimizing following functional.

where

Kernel Anisotropic Diffusion

In low SNR image that the noise magnitude is large, it may be that the diffusive function cannot distinguish between noise and edge. Therefore, the conventional AD usually fails to produce satisfactory results for low SNR images.[1516] For low SNR images, the discrimination between edges and noise can be viewed as a nonlinearly separable classification problem. According to Cover's theorem on the separability of patterns, nonlinearly separable patterns in an input space are linearly separable with high probability if the input space is transformed nonlinearly to a feature space with high dimensionality.[192425] While the gradient operation in the original AD is ineffective for low SNR images, the input space can be mapped to a higher order feature space, and then discriminate between edges and noise in the feature space.

If the input data is represented by U (U ⊆ R), and feature space is denoted by F (F ⊆ R), the nonlinear mapping function Φ is defined as:

Φ: U → F, u → Φ(u)      (5)

where u is an input in U, which is mapped into a feature space with dimensionality of n. By substituting the gradient magnitude |∇u| in the P-M's AD with ‖∇Φ(u)‖, the KAD can be obtained as following:[19]

where ‖∇Φ(u)‖ is the gradient magnitude in feature space that is defined by the mapping function Φ. It is called as a kernelized gradient operator or a kernelized edge detector, and plays a major role in KAD.

Suppose ||∇u|| denotes the norm of gradient for pixel p, we can write:

where ∂p represents the spatial neighborhood of pixel p (eight neighboring pixels around p), |∂p| refers to cardinality of ∂p (|∂p| = 8) and up denotes the intensity of p. Similarly, ‖∇Φ(u)‖ can be represented as Yu et al.:[19]

where K is the kernel function and thus the mapping function Φ is defined by kernel function K. Commonly used kernel functions are polynomial kernels, radial-basis function kernels, and two-layer perceptron kernels:[25]

K(x, y) = (x + 1)      (9)

K(x, y) = tanh(β0x + β1)      (11)

where d, σ, β0 and β1 are specified a priori by the user.

PIXEL SIMILARITY AND NEIGHBORING INCONSISTENCY

The pixel similarity measures the similarity between two pixels. By calculating the pixel similarity of center pixel with its neighbors, we can obtain the neighboring consistency of the center pixel. Informally, if a pixel is not consistent enough with its surrounding pixels, it can be considered as a noise, but an extremely strong inconsistency suggests an edge. The pixel similarity is a probability measure; its value is between 0 and 1. Accordingly, the pixel similarity can be used to determine neighboring inconsistency. It indicates the inconsistency of a pixel with its neighbors. Thus, not only we distinguishing noise and nonnoise pixels, but also we divide pixels into five categories namely smooth area, weak noise, medium noise, strong noise, and edge.[26]

If pi (i = 1, …, 8) denotes the neighboring pixels of pixel p, the pixel similarity between p and pi is defined as:

f(p, p) = exp(–λ)|u – u)      (12)

where up, upi are the intensity of pixels p, pi and λ, a are constants (typically λ = 1 and a = 2). The pixel similarity measures the similarity between two pixels. Then the degree of consistency of the center pixel with its surrounding pixels can also be calculated. The idea is, if a pixel is not consistent enough with its surrounding pixels; it can be considered as a noise. Therefore, by using a threshold we can obtain the consistency or inconsistency of center pixel with its neighbors. This threshold that is named confidence level is a predetermined constant (0 ≤ β ≤ 1 and typically β = 0.5). Using this threshold we can define the consistency discrimination function d:

where d (p, pi) = 1 means p and pi are inconsistent. By defining the degree of inconsistency of center pixel with its neighbors is determined such that if M > L (L is a predetermined constant and typically L = 5), then p is low consistent (high inconsistent) with its neighboring pixels and is considered as a noise. In this case the center pixel p is called a low consistency pixel under the confidence level β.[26]

Finally, in MRIs, if a pixel is not consistent enough with its surrounding pixels, it can be considered as a noise, but an extremely strong inconsistency suggests an edge. By detecting all kinds of low consistency pixels under different levels of confidence, the following thresholding rule that is shown in Table 1 is found:

Table 1

Category division for low consistency pixels under different confidence levels

DIFFUSIVE FUNCTION BASED ON PIXEL SIMILARITY

In Figure 1 diffusive functions g1 and g2 are plotted in terms of |∇u|. As can be deduced from this figure, the value of g1 and g2, even for not very large values of |∇u|, sharply (strictly) decrease. Thus, the diffusive process will stop and noise will remain. This result is due to the exponential operator in g1 and second ordering of |∇u|/K in g2. According to[1222] function g1 is appropriate for high contrast images, while noisy MRIs especially in high detail regions have not high contrast.

Figure 1

Diffusive functions g1 (a) and g2 (b) in terms of |∇u|

The power of |∇u|/K in g2 is a constant value, whereas this value can be more flexible and can be changed. According to this fact, we propose a new diffusive function with fractional power. In our function the real number υ is used for this power. Thus, function g2 will change as follows:

If υ is negative, the diffusion process will stop in smooth regions and will fully perform in edges. The value of υ is obtained by stability condition for P-M model that is discussed following.

When we use PDEs to carry through image processing, by supposing that ∇u≠0, =∇u/∣∇u∣ and ⊥, where and is unit tangent and normal vector respectively, we can get the following expressions:

Then we can get:[27]

where uηη, uξξ and gη, gξ are the second order derivatives and diffusive functions in and direction, respectively. The Eq. 17 is called as the global scheme for PDE-based restoration approaches.[28]

Obviously, g(|∇u|) = g(|∇u|) in Eq. 17 is always positive which means diffusion forward and making image smooth. However, g(|∇u|) = g′(|∇u|)|∇u|+g(|∇u|) can be either positive or negative leading to unstable results, this causes instability of the diffusion process.[272930] Considering P-M functions of Eq. 2 and Eq. 3, it appears that their corresponding g functions, in the global scheme of Eq. 17, can sometimes takes negative values [Figure 2]. This leads to local instabilities of the P-M's PDE which degrades the processed image instead of denoising it. As shown in Figure 2, g2 is more stable than g1 and for this reason; diffusive function g2 is utilized in this paper.

Figure 2

Investigation of the instability of Perona and Malik functions for different value of K: (a) gη for Eq. 2, (b) gη for Eq. 3

As mentioned in above, for the stability of P-M model, gη and gξ must be nonnegative. gξ is always positive and by doing some mathematical calculations, it can be found that gη will be nonnegative if:

(1 - v) (s/K)v ≤ 1      (18)

This condition is valid for |∇u| < K and practical value of μ, but for |∇u| ≥ K that higher than the previous case occurs, the following relationship must be established:

1 ≤ v ≤ 2      (19)

As describe in Section 2, the main weakness of P-M model is laid in discriminating between edges and noise such that the noise may be assigned to the edge for low SNR images and then the noise remains. In order to solve this problem, υ can be a function of neighboring inconsistency. So that, if pixel p assigned to edge, υ should be tend to 2. Thus the value of diffusion function more rapidly tends to zero. On the other hand, if p assigned to noise υ should be tending to 1 that prevents from reducing diffusive function value. In this case, since the value of υ decreases, the diffusion process cannot be stopped completely when |∇u| is large but the pixel is assigned to noise (in low SNR images).

According to above discussion, we propose the following decreasing function for υ:

and then our diffusive function is obtained as:

We name this function as the pixel similarity-based diffusive function. As can be deduced from Eq. 21, the proposed function is not only a function of the gradient magnitude and its value is related to neighboring inconsistency. Therefore, the slope of diffusive function is controlled with neighboring inconsistency. It will prevent from remaining noise in relatively low SNR images.

OTHER STABLED ANISOTROPIC DIFFUSION

As can be deduced from subsection 2.2, the main weakness of P-M model is occurred in high noise level images where the diffusive function cannot distinguish between noise and edge. Therefore, we use the consistency discrimination function to classify pixels into five groups: Smooth area, weak noise, medium noise, strong noise, and edge. This causes noisy pixels and edges to be easily distinguished from each other, even in noisy images with high noise level. On the other hand, we use consistency discrimination function in diffusive function to stable diffusion process. In other words, our proposed method will causes the conventional AD is stable. For this reason, it is essential that the proposed method is compared with the some stabled method, in addition to the classical P-M and KAD models.

Nonlocal Orientation Diffusion Partial Differential Equation

In Qiao and Zou[28] In order to achieve more stability and yet preserve edge information of the image, the diffusion process is only performed in the direction that is parallel to the image edges and can be limited in gradient direction. In this method, g is supposed zero and therefore does not take the negative values. Then the Eq. 1 may be changed to:

Where α(r) that is called nonnegative monotony digression function is related with NLM of gradient gη. This function can control diffusion velocity. When the nonlocal area of image is flat, there may be less nonlocal gradient operator, and consequently the value of α(r) is increased. This will make to a faster diffusion. When the image contains more edge features in the nonlocal areas, there may be larger nonlocal gradient operator, and it may decrease diffusion velocity to preserve the edge information. By selecting an appropriate value for α(r) to control diffusion velocity in Eq. 22, we can wipe off the image noise better and preserve the edge information effectively.

Diffusive Function Based on Double Well Potential Function

In Histace and Menard,[31] authors propose a new diffusive function to stable the P-M model:

g(|∇u|) = 1 – ψ(|∇u|)      (23)

where |∇u| denotes the normalized gradient magnitude and ψ(.) is double well potential. In Figure 3, some graphical representations of ψ(.) for different values of α are shown.

Figure 3

Double well potential function in terms of normalized |∇u| for different value of α

Now, if the mathematical expression of g is calculated, one can obtain that:

g(|∇u|) = g’(|∇u|)|∇u| + g(|∇u|)      (24)

Since |∇u| ∈ [0,1] and g’(|∇u|) is a one-order-less polynomial function than g(|∇u|) it can be concluded that:[31]

g(|∇u|) ≈ g(|∇u|) = g(|∇u|)      (25)

By taking Eq. 23, it is clear that corresponding functions in the tangential and orthogonal directions never takes negative values. Thus, diffusive process stays stable for all gradient values of processed image which is of primary importance.

According to Eq. 23, g and g are almost the same. Therefore, g tends to zero when g is close to zero. As a result, by using a right choice for α, it becomes possible to completely stop the diffusion process in the tangential and orthogonal directions of the contours, simultaneously.

RESULTS AND ANALYSIS OF EXPERIMENT

In this section, we compare our proposed diffusive model with P-M, KAD, nonlocal orientation diffusion (NLOD) and double well potential (DWP) models. We perform experiments on T1-weighted MRIs with 8-bit quantization (1000 images) to show the performance of our proposed model. Our samples are real MRIs from McConnell Brain Imaging Centre of the Montreal Neurological Institute. These images are destroyed by additive Gaussian noise with zero mean and various standard deviations.

Noise in MRI can be Gaussian or Rician distributed depending on the image type (real or magnitude data).[432] Since all images in our experiments are real data, they are corrupted with Gaussian-distributed random noise. Although MRIs are usually magnitude data and their noise obeys a Rician distribution[203334] but at high SNR (let us define SNR as I/σ), the Rician noise will be described by the Gaussian distribution. In other word, If the SNR is sufficiently high (I/σ ≥ 2,[29] 3[2033] or 4[35] or 5[30] which are usually sufficient for any practical purposes), the Rician distribution can be approximated by the Gaussian distribution.

To measure the quality of the recovered image, the criterions of SNR and improved SNR (ISNR) are employed here. For any given image u, SNR and ISNR are characterized by:

where A, I and u are the original, noisy and restored image, respectively. Generally, the larger the SNR value, the better the image quality. Another criterion that we use in this paper is mean square error (MSE) that is defined as:

where N is the number of pixels in the image. In order to evaluate the edge-preserving ability, we propose new measurement which calculates the MSE of edge pixels (edge MSE):

In which, B denotes the edges of original image A and NB is the number of pixels in B. We also adopt the structure similarity index measure (SSIM) to measure the similarity between two images,[30] which is defined as:

where μA and μu are the average of A and u, σA and σu are the variance of A and u, σA, u is the covariance of A and u, finally c1, c2 are two constants to avoid instability.

In Figure 4, we compare the proposed model with P-M, KAD, NLOD, and DWP methods. In this comparison, the gradient threshold and iteration number are selected as K = 6 and T = 10, respectively. Figure 4a is the original T1-weighted MRIs and Figure 4b presents the noisy Image by random Gaussian white noise with SNR = 10 db. Figure 4c shows the restoration result from the P-M model, in which the noise cannot be effectively eliminated. The result of KAD model with kernel function from Eq. 10 is shown in Figure 4d. The result of NLOD model with α(r)=1/(1+r2) is shown in Figure 4e. Figure 4f shows the denoising results with DWP function with α = 0.7. Figure 4g shows the denoising results with proposed model by setting a = 2, λ = 0.1 and β = 0.2. It can be observed that our model performs better than the other methods to eliminate noise and preserve edges. As can be seen in the Figure 4c–g, the ability of the proposed method in edge preserving and noise reduction is much more than the others.

Figure 4

Comparison of proposed diffusive function with Perona and Malik function and kernel anisotropic diffusion, nonlocal orientation diffusion, double well potential methods for T = 10 and K = 6: Original image (a), noisy image by Gaussian white noise with signal-to-noise = 10 db (b), output of Perona and Malik model (c), restored image by kernel anisotropic diffusion (d), output of nonlocal orientation diffusion method (e), restored image by double well potential function with α = 0.7 (f), result of proposed diffusive function (g)

Table 2 summarizes the SNR, ISNR, MSE, EMSE, and SSIM performance of the P-M, KAD, NLOD, DWP, and proposed methods for 1000 MRIs from McConnell Brain Imaging Centre of the Montreal Neurological Institute (MNI, McGill University). These images are composed of 20 different clean (without noise) images with 50 slices that are captured from the same machine. The experimental results show that the proposed model can significantly improve the measurements. From these data, the absolute superiority of the proposed method is obvious. In our method, measurements ISNR, MSE, EMSE, and SSIM are improved more than 1.6 db, 1.6 db, 1 db, and 2%, respectively. However, these improvements are given to about 5.2 db, 5.2 db, 3.8 db and 21% for ISNR, MSE, EMSE, and SSIM, respectively, in comparison with P-M model. In addition, these data indicate that the proposed method can simultaneously reduce noise and preserve edge, much more than the other methods. These advantages are due to the selection of a proper method for stabilizing the diffusion process. In Table 3, we employ P-M, KAD, NLOD, DWP, and proposed models for 1000 MRIs, and then average the values of ISNR, EMSE, and SSIM value for all images. As expected, this table confirms the previous results. We have also performed additional experiments in Table 4 for different σ values and 1000 MRIs to show the efficiency of the proposed models in higher noise levels. In this table, the additive noise is supposed white Gaussian noise. Experimental results show that the proposed method can effectively remove noise while well preserving the sharp edges and fine details of a noisy image even if the noise level is high.

Table 2

Comparison of proposed function with P-M, KAD, NLOD and DWP methods

Table 3

Comparison of proposed function with P-M function and KAD, NLOD, DWP methods for 1000 MR images

Table 4

Comparison of proposed function with P-M function and KAD, NLOD, DWP methods for high noise levels (these results are the averaged values for 1000 MR images)

CONCLUSIONS

In this paper, the instability of classical P–M model is analyzed and an improved diffusive function is proposed. This new function is based on pixel similarity and it is an improved version of P-M function with fractional and adaptive power. In the proposed function, the decreasing rate of diffusive function is controlled with the power of |∇u|/K. Such that this power (υ) increases in edge and decreases for smooth and noisy areas.

Since the P-M functions are only the function of gradient magnitude and |∇u| is large for low SNR images, the value of g tends to zero and then the noise remains. On the one hand, the noisy image is mapped to feature space and therefore, the edges of the image are not available. As can be seen from Figure 4d, the restored image in KAD is blurred.

The proposed method uses two criteria for the classification of pixels: Gradient magnitude and pixel similarity. Unlike the P-M functions that a pixel is classified in two categories, a pixel is classified into five categories in our diffusive function. Therefore, if the gradient magnitude of pixel is large but it is not noise, the diffusive function does not decrease and then the diffusion process is not stopped. Since proposed function is directly performed in input space (nor the feature space), it has much better edge-preserving than KAD. Our function allows the classical P-M model to be stable, also. Experimental results show that the proposed function can effectively improve the ISNR, MSE, and SSIM and preserve both edges and fine details of objects in the restored image. Advantages of the proposed function are apparent even for higher noise levels.

BIOGRAPHIES

Mostafa Heydari has received his Bs.c & Ms.c in electrical and electronic engineering 2004 & 2006 respectively. Now he is PhD candidate in Babol Noshirvani Institute of Technology, department of electrical and electronic engineering. He published 9 articles in journals and conferences. His research interests are digital signal processing, biomedical signal processing and image processing.

E-mail: m.heydari@stu.nit.ac.ir

Mohammad Reza Karami has received the Bs.c in electrical and electronic engineering in 1992, Ms.c of signal processing in 1994, and PhD in 1998 in biomedical engineering from I.N.P.L d’Nancy of France. He is now the Associate professor with the Department of Electrical and Computer Engineering, Babol University of Technology. Since 1998 his research is in signal and speech processing. He published more than 100 articles in journals and conferences. He teaches Digital signal processing, biomedical signal processing and speech processing in university. His research interests include Speech, Image and signal processing.

E-mail: mkarami@nit.ac.ir