Denoising Medical Images using Calculus of Variations.

We propose a method for medical image denoising using calculus of variations and local variance estimation by shaped windows. This method reduces any additive noise and preserves small patterns and edges of images. A pyramid structure-texture decomposition of images is used to separate noise and texture components based on local variance measures. The experimental results show that the proposed method has visual improvement as well as a better SNR, RMSE and PSNR than common medical image denoising methods. Experimental results in denoising a sample Magnetic Resonance image show that SNR, PSNR and RMSE have been improved by 19, 9 and 21 percents respectively.

1. INTRODUCTION

Most of images captured by imaging instruments have some extra information or noises that are not concerned to the main image. These noises appear in different form in images, which are usually random information added or multiplied by the main image.

The noise in medical images makes interpretation of images more difficult. Denoising is often necessary before analyzing (Segmentation, Classification and Detection of diseases or injury) medical images. Denoising methods can be used for reducing multiplicative or additive noise. For reducing the noise of images, there are some approaches like Spatial domain filtering and Transform domain filtering. Spatial filtering can be further classified into non-linear and linear filters. Transform domain can be classified into spatial-frequency filtering and wavelet filtering. All of these methods have some limitations. For example most filters could not preserve edge and textures of images as good as possible. Although most of filters use different quality evaluation metrics like RMSE (Root-Mean-Square Error), SNR (Signal-to-Noise Ratio) and PSNR (Peak Signal-to-Noise Ratio) for evaluating the performance of filters, there are additional assessments like visual assessment by experts and texture analysis that could be investigated.

In[1] there is a survey of different denoising techniques and their limitations. For medical image denoising, there are different methods. In[2] there is a comparative study of ultrasound image denoising methods.

Variational denoising known as “Inverse problem” was initially considered by using an energy optimization approach in[3] and further developed by Rudin et al.[4] to introduce the total variation (TV) method. This method was generalized in[5] and was used in different denoising methods like[6].

Magnetic Resonance (MR) images usually corrupted by Rician noise. Considering[7], we can assume Rician noise as additive Guassian or Raily noise in special cases. Ultrasound images also are corrupted by one type of multiplicative noise named speckle. By finding logarithm of these images, multiplicative noise is converted to additive noise[8]. We are going to introduce one variational approach for denoising medical images that can preserve fine scale details and textures of medical images and also can preserve edges of images better. So we introduce a method that is based on[6] which can preserve textures and edges of images based on local variance estimations by shaped windows and we use it for denoising both MR and ultrasound images.

The structure of the paper is as follows. In section 2, the basic concepts are described. The implementation method is outlined in section 3, and experimental results are presented in section 4. Finally, conclusions are drawn in section 5.

2. BASIC CONCEPTS

A. Calculus of Variations and Denoising

Calculus of variations is one of the mathematical branches[9] that is applied in functional analysis. Functional is a kind of function like (1):

S(f) = ∫ΩF(x, y, f, f)dxdy,      (1)

Where

The arguments of functional, also are functions. we can assume f is an image like I, The most calculus of functionals is finding minimum and maximum of them. For minimizing functional (1) the Euler-Lagrange equation is used as following:

Where

Total variational process is introduced by rudin-osher-fatemi in[4]. This process minimizes one energy functional like (3). This functional has two terms, a) smoothing term E and b) term of preserving similarity of denoised image with original image E:

Where J is noisy image, I is denoised image, Ω is image domain and λ is Lagrange coefficient (or lagrange multiplier).

The rudin-osher-fatemi process generalized in[5] to

Where Φ is a function that its properties are shown in[10], Ω is the image domain. The purpose of this process is reducing noise of image I, this can be done by minimizing above functional.

(4) can be modified as:

Where P = σ2 is noise variance. This problem could be solved by using Lagrange multiplier λ and Euler-Lagrange Equation. This problem is known as scaler Φ problem (because λ is unique for all image pixels and does not depend on each pixel of image separately), after minimizing (5) λ and I are as following:

Where I is used in steepest (gradient) descent method in (25).

B. Cartoon pyramid model

The cartoon pyramid model has been defined in[1112] and is used in many image denoising techniques. This model makes a pyramid of images in different scales. The cartoon of scale S according to[6] is defined as:

Where I is response of (5) by steepest descent method. Residue is defined as difference between two scales of cartoons:

R = C – C(p

Noncartoon part of scale S is defined as residue from level zero:

NC = R0, = C0–C      (10)

In Fig. 1 the concept of these definitions is shown.

Figure 1

The concept of cartoon, noncartoon and residue for a sample image[6].

Some basic properties of this model have been proved in[6]. For denoising, one decomposition level can be used which should contain the noise and the texture at a similar or below scale of the noise. For selecting proper scale, an estimation of noise variance can be used. Therefore, a constrained problem like (11) can be selected:

Where α controls the selected scale. The proposed model consists of three components: I where I is the cartoon, I is the noncartoon part and I is additive noise. Note that noncartoon part can consist of textures (small scale details). Thus the residue of the image is I = J – I = Ĩ + Ĩ where Ĩ shows approximation of I.

C. Local variance estimation

Estimation of variance for I (Introduced in subsection 2.B) is obtained by two steps:

In the first step, by applying scaler Φ process, strongly denoised image I is obtained from initial noisy image J, then estimation of primary variance of each image pixel is obtained from:

Where Ω’ is a subset of image pixels in a rectangle around the pixel (x,y), I are pixels in this subset. w(x̄,ȳ)=w(∣x̄–x∣,∣ȳ–y∣) is a normalized (∫w(x̄,ȳ)dx̄dȳ=1) and radially symmetric window. η[.] is the mean value with respect to pdf w(x̄,ȳ)/∣Ω’∣ on the set Ω’*Ω’ of quadruples (x,y,x̄,ȳ).

In the second step, variance of each image pixel is estimated from primary variance. For this purpose, one method like[13] is used. In this method for each pixel a locally adaptive window is found on a region. By assuming a matrix of variances with size of image, then a window with size m*m is taken around of each element of this matrix. The set of variances in each window is called A (13). Then by selecting proper variances from this set by condition (14) and finding weighted averages of them, final variance of the pixel in position (p,q) is obtained from (17). So for each pixel of I, a separate set of variances named A exists which is defined as:

where U is window with size m*m around (p,q), Qz is obtained from (12).

If |a – a0|> τ then r=0

Else r=1      (14)

where a's are members of A is variance of pixel at position (p,q), r is coefficient for selecting appropriate variances that is used in (17) and total average threshold τ is obtained from:

Where Sx,Syare image dimensions and average threshold in each window is calculated by:

So the final variance of each pixel (p,q) is obtained by:

where β is the weighting parameter. In this paper, we assume that β = 1.

3. IMPLEMENTATION

The constraint in (5) can be generalized by imposing a spatially varying variance constraint. Scalar Φ method (5) can be reformulated as adaptive Φ method as following:

Where By Lagrange coefficients, above relation can be reformulated as following:

After using Euler-Lagrange equation and applying steepest descent method we have:

where local average is defined as:

(20) is used in (27). Like[6] λ(x, y) is obtained from following equation:

where Q(x, y) is:

In this method like[6] we assign:

We assume that noise variance (σ2) in whole image is determined before denoising. For more studying, you can see methods like[2] for noise variance estimation of ultrasound images. And[14] in which an automatic histogram-based noise variance estimation technique for MR image is described.

Then the algorithm can be written as:

Set initial parameters and noise variance manually.

Separate the noise and relevant textures by minimizing ∫Ω Φ(|∇I|) subject to (11) and set I = J – I and I=I, where I is strongly denoised image that is obtained from noisy image J.

Compute local variance of I by (17) and then compute the local constraints S(x, y) by (24).

Solve (18) by iteratively evolving (20) and update λ(x, y) according to (22) and then set I=I, where I is denoised image.

Calculate I =(1–ξ)I1 + (ξ)I2 where ξ is coefficient that controls SNR value.

For implementing this algorithm in MATLAB 2009, we had some assumptions for step 1 of algorithm:

Function Φ was selected as and P for relation (11) was P = 1.5σ2, noise variance σ2 was specified before, Average window w(x, y) was selected to be a Gaussian of standard deviation σ = 5, and β = 1 for relation (17). We demonstrate psudo codes for some of the important relations:

For step 2 we used steepest descent method like this:

For k=1 to itrnno

I = I+lambda*I_t      (25)

End For

Where I initially is noisy image J, I_t is like (7), itrno is number of iterations, and lambda is computed from (6) with following code:

Lambda = mean(mean(Div.*(I-J)))/var_n      (26)

where Div is , var_n is variance of noise.

For step 4 we use steepest descent method as:

For k = 1 to itrno

I = I+lambda_xy*I_t      (27)

End For

where I initially is noisy image J, I_t is like (20), itrno is number of iterations, and lambda_xy is computed from (22) with following code:

Lambda_xy = Div.*(I-J)/Sxy      (28)

where Div is

The convergence criteria for step 2 and step 4 is :

Mean(mean(I-I_old))

Where ccval is a small value, and I_old is I before running steepest descent methods.

different shaped windows obtained by window size (WS) m*m=5*5.

The flowchart of this algorithm is shown in Fig. 3.

Figure 3

Flowchart of the proposed algorithm.

4. EXPERIMENTAL RESULTS

To compare the performance of the proposed method with various denoising methods, first, it must be shown that this method can improve image quality better than older methods. Thus, we apply this method to some simple images with specified noise before applying it on medical images.

We compare this method with TV method and adaptive method proposed in[6]. The results are shown in Fig. 4 and table (1). The results show that if we set parameters in this method properly, it can denoise image and can also preserve edges and textures of image better than the other methods.

Figure 4

Top, Original image. Middle-left, Noisy image with guassian noise and σ = 25, Middle-Right, TV method. Down-left, Adaptive Method[6], Down-right, Proposed method.

Table 1

Comparing proposed method with TV and adaptive Method proposed in[6]

The proposed method in table (1) show greater SNR and PSNR and lower RMSE than the other methods.

The drawback of our proposed method is selecting proper parameters.

We compared effect of window size on denoising top image in Fig. 4 that contaminated by Gaussian noise with σ = 35 by proposed algorithm. The effects of changing window size on quality measures are shown in table (2). It seems that window size 3*3 is better than the other sizes, but further investigation is required to distinguish relationship of window size and quality measures.

Table 2

effect of window size (WS) on quality measures

The result of denoising radial image like Fig. 5 and other images showed that when ξ is increasing from 0 to 1, the image denoised better consequently, but the edge and textures of it is destroyed. Because, when ξ is 1 then the effectiveness of image I that obtained from step 2 of process (strongly denoised method) in final I (step 5) is decreased.

Figure 5

Top-left, original image. Top-right: Noisy image with Gaussian noise (σ = 20), Middle-left, Proposed algorithm with ξ = 1, Middle-right, Proposed algorithm with ξ = 0.7, Down-left, Proposed algorithm with ξ = 0.3, Down-right, Proposed algorithm with ξ = 0

The SNR, RMSE and PSNR are shown In table (3) for comparison.

Table 3

Effect of parameter ξ on denoing.

In table (4), the results of comparing proposed algorithm and other despeckling methods[2] (lsmv: Mean and variance local statistics, lsminsc: Minimum speckle index homogeneous mask, wiener, median, homog: Most homogeneous neighborhood, gf4d: geometric, homo: homomorphic, fca: Linear scaling of the gray-level, fad: Perona and Malik anisotropic diffusion ,wavelet, Proposed) for Fig. 6 are shown. In this table SNR, RMSE, PSNR and physician visual evaluation are shown. The SNR of the proposed algorithm is greater than the others. RMSE of the proposed algorithm is lower than the others. Visual evaluation by Physician of the proposed method is after the improved median method.

Table 4

Comparing quality measures of different despeckling methods for Fig. 6.

Figure 6

(Row,Column) Method, (1,1) Breast ultrasound noisy image. (1,2) Lsmv. (1,3) Lsminsc. (2,1) Wiener. (2,2) Median. (2,3) Homog. (3,1) Gf4d. (3,2) Homo. (3,3) Fca. (4,1) Fad. (4,2) Waveltc. (4,3)Proposed.

For results of Fig. 6, the following parameters were used:

Window size m*m=3*3, σ = 5 and ξ = 0.5

In Fig. 7, one frame of heart video film and denoised version of it by the proposed algorithm are shown.

Figure 7

Left: One frame of echocardiography, Denoised frame by proposed method.

For results of Fig. 7: Window size m*m=1*1, σ = 8 and ξ = 0.5

This figure shows that the noisy image improved properly by proposed method. In Fig. 8the results of comparing the proposed method for MR image are shown. The proposed method had better performance than the other methods.

Figure 8

Top row from right to left: First image, Noisy image, TV, Bottom row from right to left: Adaptive method[6], scaler Φ, Proposed method.

For results of Fig. 8: Window size m*m=3*3, σ = 15 and ξ = 0.

Table (5) shows quality evaluation for noisy image and four variational methods by SNR, RMSE and PSNR. The first row is evaluation of noisy image, second row is result of running TV method with λ = 0 and third row is result of scaler Φ method. The fourth row is result of method[6] and last row is the result of the proposed method.

Table 5

quality evaluation of different variational methods in denoising MR image.

Further investment of running the proposed algorithm showed that this algorithm in addition to denoising image can preserve edges and textures of images better than many other methods.

5. CONCLUSIONS

The result of executing the proposed algorithm on some MR and ultrasound images showed that this algorithm can improved the method[6], specially for some figures like MR image, the size of shaped window introduced in section (2.C) can affect on SNR values, the performance of this algorithm can be further improved by using some other statistics that can be used for adjusting β parameter in (17). For separating noise and textures of image in the proposed algorithm we can use other variational method like[15] instead of minimizing of ∫ΩΦ(|∇I|) in second step of algorithm.

The complexity of this method depends on running steepest method that is an iterative numerical method for finding minimum of functional and the running time of it depends on convergence speed of steepest descent method. That is why the time of running this method is high in comparing with other methods, but if we choose parameters for algorithm properly, the convergence speed can be increased effectively and the running time decreases.

This algorithm can be used in other methods like[16] that use Lagrange coefficient for preserving textures of images.

In future works, we are going to find better energy functional for denoising medical images