An Optimal Partial Differential Equations-based Stopping Criterion for Medical Image Denoising.

Improving the quality of medical images at pre- and post-surgery operations are necessary for beginning and speeding up the recovery process. Partial differential equations-based models have become a powerful and well-known tool in different areas of image processing such as denoising, multiscale image analysis, edge detection and other fields of image processing and computer vision. In this paper, an algorithm for medical image denoising using anisotropic diffusion filter with a convenient stopping criterion is presented. In this regard, the current paper introduces two strategies: utilizing the efficient explicit method due to its advantages with presenting impressive software technique to effectively solve the anisotropic diffusion filter which is mathematically unstable, proposing an automatic stopping criterion, that takes into consideration just input image, as opposed to other stopping criteria, besides the quality of denoised image, easiness and time. Various medical images are examined to confirm the claim.

INTRODUCTION

One of the most important fields of medical image analysis involves post acquisition like denoising, because the presence of noise in image is unavoidable. It can be caused by the image formation process, image recording and image transmission. These random distortions make problems to perform any required image processing. Even a small amount of noise is detrimental to the high accuracy analyzing, e.g., in sub-cell image analysis. Thus, removing noise is often the first step in interferograms analysis. In recent years, partial differential equations methods in image processing have been extensively used due to their good performance to improve traditional approaches. These methods cannot only remove the noise but also keep much more details without blurring or changing the location of the edges. Moreover, their scope include not only the significant problem of image denoising but also other restoration tasks such as deblurring,[12] blind deconvolution,[1] inpainting,[345] edge detection[67] and segmentation.[8910] In this paper, we focus on the use of anisotropic diffusion model which has been extensively studied since it is introduced by the Perona and Malik in 1987.[11] Perona-Malik model has been considered as a useful tool for image noise removal and other areas of image processing.[121314] However, if we decide to restore noisy images using some methods starting from an input image, which led to a set of possible filter solutions by gradually removing noise, the crucial question is when to stop filtering in order to get the optimal restoration result. The objective is quite challenging, because the stopping time has a great effect on the output result. Introducing small amount of stopping time gives more confidence to the input data, so more noise in input data is remained; whereas, most parts of the image are destroyed in the large value of stopping time.[15] On the other hand, stopping in the steady state yields an overly smooth image;[16] as a result, the time derivative acts as a regularizing parameter and a fundamental feature of abovementioned reconstruction procedure is the necessity to decide when to stop the iterations. In general, this is done by trial and error. Many studies have been conducted in the area of introducing stopping criterion. Let us briefly review previous works on the stopping criteria. Dolcetta and Ferretti defined a stopping time by finding a minimum of the performance index.[17] They need a constant that is found using an ordinary image with similar details and discontinuities. This is a rather ambiguous requirement; also they need some approximations to determine the constant. Sporring and Weickert proposed the stopping criterion based on the signal to noise ratio, the relative variance at time t and the initial image.[18] Mrazek and Navara extended this idea and introduced the stopping criterion by considering some assumptions to minimize the correlation of the signal and noise.[19] Other methods have been introduced to choose a stopping criterion, but most of them used statistics information[1520] or were applied to one dimensional signal.[21] Some of the previous methods were computationally intensive.[2223] Furthermore, the above techniques did not discuss about advancing the equation in time and the cause of iterative methods on the stopping criterion or frequency approaches. One of the working in these areas is suggested by Ilyevsky and Turkel.[16] In their stopping criterion the Perona-malik model is solved by complex and time consuming multigrid method. Here, the proposed automatic stopping criterion can satisfactory reduce the computational cost besides improving image quality and speed property in which it uses only the corrupted image content itself and does not need any additional information. The performance of proposed algorithm on medical images from the aspects of great improvement in time and denoising trough peak signal to noise ratio (PSNR) is evaluated. An outline of the paper is as follows:

In Section 2, we briefly introduce Perona-Malik model. This helps us to understand the important advantages of this model especially for medical imaging. Section 3 analyzes our new proposed stopping criterion with the important principles behind our approach. Some empirical evidences to show its efficiency based on reduced costs, PSNR analysis, and visual verification for a more systematic comparison is considered in Section 4 and the paper is concluded with a summary in Section 5.

PERONA-MALIK MODEL

The idea behind the use of the diffusion equation in image processing comes from the result of Gaussian filter in multiscale image analysis. Convolving the given image with a Gaussian filter kσ:

with standard deviation σ, is equivalent to the solution of the diffusion equation in two dimensions. Therefore, for an N × N image I the diffusion process defined as:

Where L is the maximal gray level, I(x, y, t) is the image I(x, y) at time t = 0.5 σ2 with initial condition I (x, y, 0) = I0 (x, y) and I0 denotes the original image.

In general, this can be written as:

Where c(x, y, t) is the diffusion coefficient, ∇ represents the gradient operator and div denotes the divergence operator.

Perona and Malik introduced an inhomogeneous diffusivity by considering the following function for c:

Here the parameter k is a threshold to evaluate the gradient of the image. This definition of c provides a good control for the process of denoising and edge preserving. It means that in the area with the high value of ∇I such as edges results in decreasing of c. Therefore, the problem of edge smoothing will be solved by the model of Perona-Malik. Whereas in homogeneous regions where observation of noise is more likely, ∇I is small; thus, c tends to 1; consequently, smoothing is achievable. Interestingly, there exists a relevancy between (3) and the neural dynamics of brightness perception. Cohen and Grossberg[24] presented a model of the primitive visual cortex with similar inhibition effects as in the Perona-Malik model. Furthermore, rapid decay of diffusivity c causes non-monotone flux function. This property of rapidly decreasing diffusivities are explicitly intended in the Perona-Malik method, so they offer the desirable result of blurring small fluctuations and sharpening edges in one single process.[25] The results of Perona and Malik are visually impressive; edges remains stable over a very long time. It is demonstrated[11] that edge detection based on this process clearly outperforms the linear canny edge detector. It is easily seen that many of the preceding results can be generalized to higher dimensions. This can be useful especially for medical image sequences from computed tomography (CT) or magnetic resonance imaging (MRI) images.[25]

Implementation of nonlinear diffusion filters are mostly based on finite difference methods, because they are easy to apply. Furthermore, the pixel structure of digital images provides a natural discretization of a fixed rectangular grid and they are well-suited for parallel architectures.[25] On the other hand, the Perona-Malik process is instable. However, the only observed instability result is the so-called staircasing effect. Weickert and Benhamouda[26] showed that the regularizing effect of a standard finite difference discretization is sufficient to solve instability problem and turn the Perona-Malik filter to a well-posed initial value problem. This is other significant advantage of utilizing finite difference-based explicit scheme.

On the other hand, the problem with the explicit scheme is that due to the stability restrictions the time step must be very small, resulting in very slow convergence. As a result, in this paper a swap technique to tackle this essential difficulty of explicit method is used, so the main objectives in this article are (1) to present a fast and effective stopping criterion with low computational cost, (2) to utilize explicit algorithm due to its capabilities and as a remedy for ill-posedness problem of Perona-Malik model, (3) to provide a software technique with high ability to remove the difficulty of low convergence in the explicit scheme and decrease processing time.

In the following subsection, we introduce an effective time stopping procedure which takes the advantages of using explicit time marching algorithm, very low computationally expenses and improved quality.

NEW STOPPING CRITERION

Definitions

High frequencies are defined as frequencies from

Denote the L2 as a Euclidean norm of the transform of the image in the frequency domain when only high frequencies are considered (indicative of a criterion to evaluate the energy of high frequencies).

We solved Perona-Malik equation by explicit method and computed L2 in every iteration (L2) Then we computed high frequency relative forward difference (HFRFD) to analyze grey value fluctuations within a neighborhood of each image point, for we need information about its derivatives to understand the structure of an image.

It is known that noise is usually represented by high frequencies in the frequency domain. The converse is not true, because all high frequencies are not noise. Interestingly results arise when one probes the diagram of HFRFD. It means that at first, the rate of decreasing in L2 (energy of high frequencies) is fast (during the process of denoising) but after some iterations the rate of decreasing becomes slow and slower and consequently uniform (an asymptotic line) [Figure 1]. There is a region with the maximum cavity in its diagram before the asymptotic line including an important point of undulation point. This point is located close the minimum point of the curve. In this region the race of decreasing in L2 becomes very slow and most parts of the image is denoised; thus, in practice, more denoising after this situation is not useful and even damages the important features of image such as edges and lines.

Figure 1

The L2 and high frequency relative forward difference reduction in every iteration for a noisy abdomen computed tomography contains a Gaussian noise of variance = 0.02 during the denoising process

Furthermore, with more investigation, it will be clear that the thresholds for parameter r (r ≃ 1 to 1.3) which is introduced in[16] as the ratio of the spectral energy of the denoised image and the blurred version of the image with no noise for a stopping criterion are located in this region; consequently, a threshold should be determined to stop denoising in this region. This threshold is selected as 0.09 by various experiments on different images with different details. Figure 1 represents L2, HFRFD diagrams and the situation of r for the abdomen CT which is contaminated by a Gaussian noise of variance = 0.02. It shows the details of the proposed restoration algorithm. Thus, we will have this algorithm:

Solve the Perona-Malik equation and compute HFRFD as

.

Stop when

We use PSNR for quantification comparing and observation manner for qualification comparing.

The PSNR is calculated as:

With,

In the equations (5) and (6) MSE is the mean square error. Where O is an original non-corrupted image, R is a restored image and Max is the maximum possible pixel value of the image. The advantages of proposed method from the aspects of visual and qualification criteria and processing time will be discussed in the next sections.

EXPERIMENTAL RESULTS

In this section, we present some examples to show the effectiveness of the proposed approach. To this aim, we used the set of images shown in Figure 2. including body CT (a), pelvic CT (b) abdomen CT (c), ankle MRI (d), knee MRI (e), and neck MRI (f) images and we added to them Gaussian noise of zero mean with different variances (σ = 0.02, 0.04, 0.08). Then the different stopping times, using proposed and Ilyevsky et al. methods are compared. All images are in size of 512 × 512; however the proposed technique can be used for higher dimensions due to its low processing time. The denoised image, in every case, is evaluated using the PSNR and processing time and optical results. The results of the measurements are shown in Table 1 and Table 2. The optical results can be seen in, Figures 3–20. For each image, methods are executed by a personal computer with these characteristics:

Figure 2

Tested images

Table 1

Results of implementation for proposed method and Ilyevsky algorithm using Perona-Malik model on CT images

Table 2

Results of implementation for proposed method and Ilyevsky algorithm using Perona-Malik model on MRI images

Figure 3

(a) Original image. (b) A noisy image with σ = 0.02. (c) The Ilyevsky filtering result using Perona-Malik model with t = 5639.6 s, peak signal to noise ratio (PSNR) = 27.2691. (d) The proposed method using Perona-Malik model with t = 10.21 s, PSNR = 27.3588

Figure 20

(a) A noisy image with σ = 0.08. (b) The Ilyevsky filtering result using Perona-Malik model with t = 14578.4 s, peak signal to noise ratio (PSNR) = 28.6534. (c) The proposed method using Perona-Malik model with t = 22.07s, PSNR = 28.6575

(a) A noisy image with σ = 0.04. (b) The Ilyevsky filtering result using Perona-Malik model with t = 8676.9 s, peak signal to noise ratio (PSNR) = 27.2095. (c) The proposed method using Perona-Malik model with t = 15.07 s, PSNR = 27.3473

(a) A noisy image with σ = 0.08. (b) The Ilyevsky filtering result using Perona-Malik model with t = 13670.2 s, peak signal to noise ratio (PSNR) = 27.1131. (c) The proposed method using Perona-Malik model with t = 21.98 s, PSNR = 27.3181

(a) Original image. (b) A noisy image with σ = 0.02. (c) The Ilyevsky filtering result using Perona-Malik model with t = 5633 s, peak signal to noise ratio (PSNR) = 28.7256. (d) The proposed method using Perona-Malik model with t = 9.87 s, PSNR = 28.7434

(a) A noisy image with σ = 0.04. (b) The Ilyevsky filtering result using Perona-Malik model with t = 8723 s, peak signal to noise ratio (PSNR) = 28.6074. (c) The proposed method using Perona-Malik model with t = 14.93 s, PSNR = 28.6304

(a) A noisy image with σ = 0.08. (b) The Ilyevsky filtering result using Perona-Malik model with t = 14295 s, peak signal to noise ratio (PSNR) = 28.4347. (c) The proposed method using Perona-Malik model with t = 21.98 s, PSNR = 28.4805

(a) Original image. (b) A noisy image with σ = 0.02. (c) The Ilyevsky filtering result using Perona-Malik model with t = 5674.3s, peak signal to noise ratio (PSNR) = 29.2248. (d) The proposed method using Perona-Malik model with t = 9.99 s, PSNR = 29.2310

(a) A noisy image with σ = 0.04. (b) The Ilyevsky filtering result using Perona-Malik model with t = 11659 s, peak signal to noise ratio (PSNR) = 29.0844. (c) The proposed method using Perona-Malik model with t = 14.49 s, PSNR = 29.0881

(a) A noisy image with σ = 0.08. (b) The Ilyevsky filtering result using Perona-Malik model with t = 14360.4 s, peak signal to noise ratio (PSNR) = 28.8931. (c) The proposed method using Perona-Malik model with t = 21.05 s, PSNR = 28.8892

(a) Original image. (b) A noisy image with σ = 0.02. (c) The Ilyevsky filtering result using Perona-Malik model with t = 6224.5 s, peak signal to noise ratio (PSNR) = 28.5784. (d) The proposed method using Perona-Malik model with t = 8.84 s, PSNR = 28.6108

(a) A noisy image with σ = 0.04. (b) The Ilyevsky filtering result using Perona-Malik model with t = 8674.8 s, peak signal to noise ratio (PSNR) = 28.4865. (c) The proposed method using Perona-Malik model with t = 13.79 s, PSNR = 28.5157

(a) A noisy image with σ = 0.08. (b) The Ilyevsky filtering result using Perona-Malik model with t = 14192.3 s, peak signal to noise ratio (PSNR) = 28.3473. (c) The proposed method using Perona-Malik model with t = 22.47 s, PSNR = 28.3877

(a) Original image. (b) A noisy image with σ = 0.02. (c) The Ilyevsky filtering result using Perona-Malik model with t = 6307.5 s, peak signal to noise ratio (PSNR) = 27.3957. (d) The proposed method using Perona-Malik model with t = 9.75 s, PSNR = 27.4777

(a) A noisy image with σ = 0.04. (b) The Ilyevsky filtering result using Perona-Malik model with t = 10116.1 s, peak signal to noise ratio (PSNR) = 27.3241. (c) The proposed method using Perona-Malik model with t = 14.56 s, PSNR = 27.4482

(a) A noisy image with σ = 0.08. (b) The Ilyevsky filtering result using Perona-Malik model with t = 14213.6 s, peak signal to noise ratio (PSNR) = 27.2263.(c) The proposed method using Perona-Malik model with t = 22.46 s, PSNR = 27.4197

(a) Original image. (b) A noisy image with σ = 0.02. (c) The Ilyevsky filtering result using Perona-Malik model with t = 5956.9 s, peak signal to noise ratio (PSNR) = 28.9074. (d) The proposed method using Perona-Malik model with t = 10.03s, PSNR = 28.9366

(a) A noisy image with σ = 0.04. (b) The Ilyevsky filtering result using Perona-Malik model with t = 9704.6 s, peak signal to noise ratio (PSNR) = 28.8173. (c) The proposed method using Perona-Malik model with t = 15.14 s, PSNR = 28.8316

CPU: AMD Athlon™ 64 × 2 dual core processor 4800 + 2.50 GHz

Memory: 2.00 G.B

Operating system: Windows 7 (32-bit)

Image processing software: Matlab 7.10.0.499 (R2010a).

Numerical Experiments

As it can be seen from the Table 1 that results of proposed method show much less computational cost and better PSNR is achieved for most of the images. This increase of PSNR for images with more details is higher than other cases. For example, 14 iterations of Ilyevsky method need 5633 s, and obtained PSNR is 28.7256 whereas 14 iterations of HFRFD reconstruction procedure need 9.87 s, and obtained PSNR is 28.7434 for a noisy version of Pelvic CT image with σ = 0.02. Moreover, for the body CT image with σ = 0.08, 31 iterations of the multigrid algorithm and 52 iterations of proposed technique produce restored images with PSNR = 27.1131 in 13670.2 s and PSNR = 27.3181 in 21.98 s, respectively and finally 20 iterations of Ilyevsky method is executed in 8676.9 s, and PSNR is 27.2095 but 35 iterations of proposed procedure need 15.07 s, with PSNR = 27.3474 for the Body CT image which Gaussian noise with σ =0.02 is added to it.

It should be noted that HFRFD approach gives visually better results than those of the denoised images carried out by multigrid scheme which they seem blocky. As shown in Table 1, the proposed technique takes in general more iteration to reach its best restored image. However, this observation does not imply that it makes the denoising process slow down because we use swap technique. It is safe to say that the proposed algorithm does speed up the restoration procedure but runs more steps to continuously improve image quality in higher levels.

Table 2 shows PSNRs and processing time of restored images by applying two different mentioned denoising mechanisms on MRI images. It is appeared from the Table 2 that proposed algorithm shows great improvement in computational efficiency and gains better PSNRs compared with another method for all MRI images. To illustrate, Knee MRI image with variance of 0.02 is restored by proposed strategy in t = 9.75 s with PSNR = 27.4777 and using Ilyevsky technique in t = 6307.5 s with PSNR = 27.3957.

Representing the Gradient

Following, the Sobel gradient of images, which is a popular edge detection method, yielding by our method will be presented for better analyzing a tradeoff between edge preserving and denoising. The Sobel operator is used in image processing based on an edge detection algorithm. It is a discrete differentiation operator which computes an approximation of the gradient of the image. It is based on the convolving of the image with the separable and integer valued filter in the horizontal and vertical directions. The filter is an approximation of a Gaussian which is the only filter that does not introduce artifacts. Each direction of Sobel filter is applied to an image. Then two new images are created. One image shows the vertical response and the other shows the horizontal response. Afterwards, two images are combined into a single image. The purpose is to determine the existence and location of edges in an image. The results of Sobel gradient images are shown in Figures 21–23. It can be seen that the difference in denoising for HFRFD = 0.09 and HFRFD = 0.1 is visible, whereas the amount of denoising with compared to degraded edges after HFRFD = 0.09 is not remarkable. This is the reason for very slow decreasing in energy of high frequencies of denoised image after HFRFD = 0.09 as shown in Figure 1.

Figure 21

(a) Sobel gradient of the original image. (b) Sobel gradient of noisy image with σ = 0.02. (c) Sobel gradient of the proposed filtering result with high frequency relative forward difference (HFRFD) = 0.1. (d) Sobel gradient of the proposed filtering result with HFRFD = 0.09. (e) Sobel gradient of the proposed filtering result with HFRFD = 0.06. (f) Sobel gradient of the proposed filtering result with HFRFD = 0.03

Figure 23

(a) Sobel gradient of the original image. (b) Sobel gradient of noisy image with σ = 0.08. (c) Sobel gradient of the proposed filtering result with high frequency relative forward difference (HFRFD) = 0.1. (d) Sobel gradient of the proposed filtering result with HFRFD = 0.09. (e) Sobel gradient of the proposed filtering result with HFRFD = 0.06. (f) Sobel gradient of the proposed filtering result with HFRFD = 0.03

(a) Sobel gradient of the original image. (b) Sobel gradient of noisy image with σ = 0.04. (c) Sobel gradient of the proposed filtering result with high frequency relative forward difference (HFRFD) = 0.1. (d) Sobel gradient of the proposed filtering result with HFRFD = 0.09. (e) Sobel gradient of the proposed filtering result with HFRFD = 0.06. (f) Sobel gradient of the proposed filtering result with HFRFD = 0.03

DISCUSSION AND CONCLUSION

In the field of image processing and analysis, the problem of image denoising forms a significant preliminary step of image restoration. A major concern in image denoising models is to preserve important features such as edges and lines which are easily detected by the human visually system while denoising process. Perona-Malik model as anisotropic diffusion filter aims at degrading noise without removing fundamental parts of the image. This filter has become a valuable tool for many areas of image processing. However, this filter performs successfully as long as some important points are considered. A main feature of anisotropic diffusion procedure is presenting an easy solving method due to its nonlinearity property in addition to removing its ill-posedness problem. The other important point is introducing the stopping time with low computational cost as a regularizing parameter which forms the quality of output image. Finally, using effective and swift software technique to handle processing time when one uses this filter is significant subject that cannot be neglected. Here, we used explicit method for addressing the ill-posedness problem of Perona-Malik model besides its other mentioned advantages.[25] Furthermore, an automatic stopping criterion is proposed that takes into consideration the quality of denoised images and uses only the content of noisy images in terms of its high frequency energies. Moreover, in order to further improvement in denoising procedure, a software technique with high capability to tackle the difficulty of low convergence in explicit scheme and processing time is introduced. We verified effectiveness of the proposed denoising strategy using several images with different levels of noise and Sobel gradient images. The results of the measurements which are shown in the Section 4 indicated that the proposed method achieved better results regarding the quality measure of denoising, visual verification, preserving important features and time.

BIOGRAPHIES

Maryam Khanian received B.Sc. and M.Sc. degrees in applied mathematics in 2006, and 2011 respectively, and a Ph.D. scholarship from Brandenburg University of Technology, Germany, in 2013. Now she works on her thesis in the area of Photometric Stereo based on the Prospective Projection. She has worked in the area of image processing for about four years.

E-mail: mar.khanian@gmail.com

Awat Feizi obtained his B.S degree in Statistic in 2000 from Alameh Tabatabaei University (Tehran-Iran), MS.c in Biostatistics from Tarbiat Modares University (Tehran-Iran) in 2003 and Ph.D in Biostatistics from Tarbiat Modares University (Tehran-Iran) in 2008. He joint to Department of Biostatistics and Epidemiology at school of health in Isfahan University of medical sciences as an assitant professor and now is working as an Associate professor of Biostatistic at Department of Biostatistics and Epidemiology.

E-mail: awat_feiz@hlth.mui.ac.ir

Ali Davari received his Ph.D. degree in 2004 from Kharazmi University, Tehran, Iran. He spent 1.5 years as a Postdoctoral researcher at the University of Isfahan, Isfahan, Iran. He is an assistant professor in the department of mathematics at University of Isfahan, Isfahan, Iran.

E-mail: a_davari@sci.ui.ac.ir