CrossMark

Available online at www.sciencedirect.com

ScienceDirect

Procedía Computer Science 89 (2016) 658 - 665

Twelfth International Multi-Conference on Information Processing-2016 (IMCIP-2016)

Coupled PDE for Ultrasound Despeckling using ENI Classification

R. Soorajkumara, P. Krishnakumara, D. Girishb'* and Jeny Rajana

aNational Institute of Technology, Surathkal Mangalore 575 025, Karnataka, India bDevi Scan Centre Kothamangalam, Cochin, Kerala 686 691, India

Abstract

Speckle is a type of noise which is often present in ultrasound images. Speckle is formed due to constructive or destructive interference of ultrasound waves. Due to the granular pattern of speckle noise, it hides important details in ultrasound images. Many despeckling techniques are proposed in the literature, but most of them fail to reach a balance between the removal of speckle noise and preservation of the fine details in the image. In this work, an improved coupled PDE model is proposed which combines second order selective degenerate diffusion (SDD) model and fourth order PDE model based on the assumption that speckle in ultrasound image follows Gamma distribution. An edge noise interior (ENI) method is used to control the diffusion. With the help of ENI controlling function, the diffusion at edge pixels and noisy pixels are selectively accomplished with varying speed. Thus, the proposed model preserves the edges and fine texture details in the image. The model is tested on simulated images after corrupting the images with various levels of Gamma noise. Further, we have tested it on real ultrasound images also. The performance of the proposed model is compared with other similar techniques and the proposed method outperforms other state-of-the-art methods, both in terms of qualitative and quantitative measures.

© 2016 The Authors.PublishedbyElsevierB.V. This is an open access article under the CC BY-NC-ND license (http://creativecommons.Org/licenses/by-nc-nd/4.0/).

Peer-review underresponsibility of organizing committee of the Organizing Committee of IMCIP-2016 Keywords: Coupled PDE; ENI; Gamma Noise; Speckle; Ultrasound.

1. Introduction

Ultrasound (US) imaging is a popularly known method in medical imaging modalities because of its versatility, non-invasiveness, portability, non-ionizing and low cost nature. However, granular pattern which exists in US image called 'speckle', degrades the performance of post-processing tasks applied to US images such as image segmentation and registration and also influences the visual analysis. Hence, there is a need for designing a filter which suppresses speckle noise without losing relevant image features.

In the past two decades, various methods have been proposed to despeckle the US images. These techniques include linear filtering, non-linear filtering, wavelet filtering, stochastic based filtering and anisotropic diffusion filters. Linear filters like Lee1, Kuan2 and Frost3 filters replaces the center pixel by the weighted average of all the pixels around the neighbourhood pixels within the kernel window. These filters depends on the coefficient of variation (COV), which maintains a balance between the edge and texture features in non-homogeneous and homogeneous regions. Coupe et al.4 introduced a non-local means (NLM) based approach (OBNLM) for speckle reduction, which exploits the

* Corresponding author. Tel.: +91-984-5848743. E-mail address: soorajkumarbhat@gmail.com

1877-0509 © 2016 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY-NC-ND license (http://creativecommons.Org/licenses/by-nc-nd/4.0/).

Peer-review under responsibility of organizing committee of the Organizing Committee of IMCIP-2016 doi:10.1016/j.procs.2016.06.031

data redundancy in the image. Partial differential equation (PDE) based methods have been employed in the past for denoising and enhancing the images. The essential thought of PDE based methods is to warp the pattern of a broken image with a PDE model and get the desired results as the solution of this PDE with the noisy image as input5. Witkin6 proposed a PDE model which follows heat equations. Even though this model diffuses in all directions to remove noise, the edges are not well preserved. To overcome this limitation, later researchers viewed this problem in a different perspective by taking three factor5 into account (i) controlling the speed of diffusion (ii) controlling direction of diffusion (iii) or combination of both5.

There exist many PDE-based methods which were used to remove noise from images. Perona and Malik (PM)7 proposed a second order PDE model, in which the authors introduced the controlling function depending upon image gradient. Catte et al.8 proposed a selective smoothing model by improving the controlling function of PM equation. Further, Alvarez et al.9 improved the controlling function of PM model and added a diffusion direction term to it and named it as selective degenerated diffusion (SDD) model. Here, diffusion takes place with respect to the direction as well as speed. Yu and Acton10 proposed a speckle reducing anisotropic diffusion (SRAD) method to suppress the speckle noise by preserving the edges. Aja-Fernandez et al.11 proposed a detail preserving anisotropic diffusion (DPAD) which uses a diffusion function based on Kuan filter. Both SRAD and DPAD use instantaneous coefficient to preserve the edges.

Most of the second order PDE models proposed so far are trying to achieve a trade off between removing noise and preserving edges. However, second order models usually suffer from blocky staircase effects due to piecewise constant approximation which add an unnatural look to the image. You and Kaveh12 proposed a fourth order PDE to remove speckle noise, where they used a diffusion controlling function which is based on Laplacian of an image. By using piecewise planar approximation, the blocky stair case effects get reduced, resulting in a more natural looking image. However, fourth order PDE methods require more number of iterations to converge.

A coupled PDE model is presented which is a combination of second order SDD9 model and a fourth order12 PDE model. Our diffusion controlling function is based on edge-noise-interior (ENI) computation. In this way, diffusion at edge and noisy pixels are selectively accomplished with various speeds. Since we made the assumption that the speckle in the US image follows the Gamma distribution13, it become essential to remove the bias introduced after denoising. For this purpose, we have estimated shape and scale parameters of the Gamma distribution using maximum likelihood (ML) approach14.

We organised the paper as follows: Speckle noise characteristics are discussed in Section 2. The coupled PDE model is explained in Section 3. In Section 4, results obtained through experiments are discussed and we conclude the paper in Section 5.

2. Speckle Noise Characteristics

Speckle noise exists in US image which degrades its quality. To effectively remove speckle noise, it is essential to know the statistical distribution of noise. Speckle is formed due to overlapping of two or more backscattered US waves. Depending upon the scatter number density (SND) and presence/absence of deterministic components speckle is classified into four types. When SND is large with the absence of the deterministic components, the speckle is categorized as fully developed speckle noise and it follows Rayleigh distribution15. When SND is large with the presence of the deterministic components, the speckle is categorized as fully resolved speckle noise and it follows Rician distribution15. When SND is small with the absence of the deterministic components, the speckle is categorized as a partially developed speckle noise and it follows K-distribution16. When SND is small with the presence of deterministic components, the speckle is categorized as partially resolved speckle and it follows homodyned K-distribution16. Although in practice, the number of scatters are high, the distribution doesn't follow exactly Rayleigh distribution instead it follows the weighted sum of Rayleigh variables which can be approximated with a Gamma distribution13.

3. Proposed Model

3.1 Speckle noise model

In this study, we assume speckle noise in US image follows Gamma distribution which can be modelled as:

I = Io + n

where, I, /0 and n denotes the noisy image, the original image and the Gamma noise respectively. The probability density function (PDF) of Gamma distribution is as follows17:

,a-l ixi-y>

y-^l — y

P{x\y,a,p) =

(x' 1 ]— exp P Vy < x < oo,

0 otherwise.

where a > 0 represents a shape parameter, f > 0 represents a scale parameter and y represents location parameter. The mean ¡ix and variance ax? is given by14:

Mx = Y + af (3)

a? = af2 (4)

and the third parameter y is estimated by14:

3.2 Background of PDE models

y = Hx - aft (5)

Perona and Malik (PM) proposed a PDE model in which the gradient of the image is used to control the speed of diffusion which can be written as7:

dJL = div{c{\Vl\)Vi) (6)

where V denotes the gradient and div denotes the divergence operator. c(.) represents the controlling function which depends on the gradient of an image. I denotes the observed noisy image.

Catt et al.s introduced a selective smoothing PDE model by improving the controlling function c(.) of Eq. (6) of the PM model which is as follows8:

— = div(c(\ VG*/ |)V/) (7)

where G denotes the Gaussian smoothing kernel.

Alvarez etal.9 proposed an improved PDE model which is based on controlling the direction of diffusion and named as degenerate diffusion PDE which is given as9:

d I ( VI \

— = \VI\div (-) (8)

8t 1 1 V IV/I /

In this model diffusion is done in the orthogonal direction of gradient value which preserves the edges. Further, to enhance the edges, they incorporated a diffusion speed term into it and they renamed the model as SDD9:

51 i VI \

— = c(\VG*I\)\VI\div^—j (9)

Controlling the speed and direction of diffusion is introduced in Eq. (9) which enhances the edges by suppressing the noise.

Since second order PDE suffers fromblocky staircase effects, You and Kaveh12 proposed a fourth order PDE model which replaces the controlling function c(.) of PM by Laplacian of an image as follows12:

51 o o

— = — V (c(.)V /) (10)

where c(.) controlling function which depends on Laplacian of an image and k is a constant.

c(|A/|) =-- (11)

The fourth order PDE method reduces the blocky effects, but requires more iterations to converge. To overcome the limitations of fourth order PDE, a new diffusion controlling function is used in this work which is based on the ENI and also proposed a novel coupled PDE method which combines second order SDD9 and fourth order PDE model12. This reduces the blocky effects and requires less number of iterations to converge.

(a) (b) (c)

Fig. 1. (a) Simulated Image Corrupted by Gamma Noise;(b) ENI Image;(c) Controlling Function Output Image.

3.3 ENI based diffusion controlling function

In most of the PDE based methods, the controlling function diffuses the edges along with noise and interior pixels (smooth region). A good controlling function should retain the edges by allowing more diffusion near the noisy pixels than the edge pixels. We have borrowed the concept of edge noisy interior (ENI)5. ENI counts the number of homogeneous pixels in a local neighbourhood of the image. The ENI value of edge pixels is near to zero. The ENI value of noisy pixel is intermediate and that of interior pixels is near to the maximum intensity of the image (255). The ENI calculation for an image is explained below.

Let p(i, j) be the pixel location, then neighbourhood points of the center pixel p of window size (2w + 1) x (2w + 1) for w > 0 is denoted by Nt

Nw = {(m, n) : \m — i | <= w, \n — j | <= w}

and NW is a set of neighbouring points excluding center pixels p. For each pixel qt NW absolute difference is defined

as d(p,q)5:

d(p,q) =\Ip - Iq | (13)

which gives the correlation between neighbouring pixels with the center pixel.

Further, the intensity value of each qtNW is classified into two groups using a global threshold T5:

1 if d(p,q) < T,

0 d( p,q)>T )

Finally, the ENI for pixel p is given as5:

ENI (I, w, T) = ^ Eq

q e N(w)

Hence, the controlling speed function c(.) is given as in5:

1 1 {2n ENI (I, w, T)

c(ENI(I, w, T)) = - + - cos i-N

where ENI (I, w, T) is the ENI image which depends on gray scale intensity I, window size w = 2, and threshold

T = 40. c(ENI(I, w, T)) values lies between 0 and 1.

Figure 1 illustrates the output of controlling function, given a simulated circle image as input which is corrupted

with Gamma noise a = = 6. The corresponding ENI function output and the controlling function output images

are shown in Fig. 1(b) and Fig. 1(c) respectively. Both figures clearly depicts that the ENI technique preserves the

edges and other texture details.

3.4 Proposed PDE

A novel coupled PDE model is proposed which combines fourth order12 PDE and second order SDD9 PDE by

modifying with ENI5 concept. The fourth order PDE reduces the blocky staircase effect and second order SDD PDE

to achieve fast convergence. The proposed coupled PDE is given as:

dt=D\VI\div{—)-V-W-I) (17)

where D = c(ENI(I, w, T)) is a controlling function which depends on ENI of the image which is given in Eq. (16), V denotes gradient operator.

To compute V2(DV21) from Eq. (17), first the Laplacian of the image is calculated by using symmetric boundary conditions as given below:

L = V21 = Ii+1, j + Ii-1, j + Ii, j-1 + Ii,j+1 - 4Ii,j (18)

Then diffusion coefficient D is multiplied by L followed by application of Laplacian operation by using symmetric boundary conditions as follows:

V2(DL) = Di+1,jLi+1,j + Di-1,jLi-hj + Di,j-1 Li,j-1 + D,j+1 L,j+1 - 4D,jL,j (19)

The discrete scheme of the proposed coupled PDE is give as follows:

/»+1 = + a+ _ AfV2(Z)jW m

where, n denotes number of iterations, I^^j denotes the numerical solution in 2-D with i and j as subscripts of the pixel position, At is the time step. Ix, Iy represent the first order partial derivatives along x and y directions respectively and Ixx, Iyy represent the second order partial derivatives of Ii, j.

3.4.1 Maximum likelihood estimation

The ML estimation technique14 is used to estimate the value of a (shape) and ¡3 (scale). After applying Eq. (20), unbiased value for each pixel in the image is computed using the formula14'18:

I = max{I - aft, 0} (21)

INPUT : Image(/) corrupted with Gamma noise, n, w, T OUTPUT : Denoised Ultrasound image (/)

1. INITIALIZATION:

i Window size (w <- 2), N <- (2w + 1) x (2w + 1)

ii Threshold (T <- 40), Number of iterations (Iterate)<- 12, time step (Ai) <- 0.15

l: while lterate>0 do

i Calculate ENI image using Eq. (12), (13), (14) and (15)

ii Calculate the controlling function and output image c(ENI(I, w, T)) using Eq. (16)

iii Compute coupled PDE using equation using the Eq. (17)

iv Update the image using Eq. (20)

v Iterate <— Iterate - 1

2: end while

3: a andp are calculated using ML approach and unbiased values are computed using Eq. (21)

Algorithm 1. Coupled PDE Algorithm

(g) (h) (i)

Fig. 2. (a) Original Phantom Image (ground truth);(b) Noisy Image (a = = 5); (c) Denoised Image using DRAD;(d) using SRAD;(e) using OBNLM;(f) using Proposed Method;(g) Residual Image of SRAD Method;(h) Residual Image of OBNLM Method;(i) Residual Image of Coupled PDE Method (in the scale 0 to 40).

4. Experimental Results and Discussion

The experiments were carried on both simulated images as well real US images. Performance evaluation was measured using standard metrics like peak signal to noise ratio (PSNR)19 and mean structural similarity index metric (SSIM)20. The proposed filter was compared with different filters which includes SRAD10, Lee1, Frost3, Kuan2, DPAD11, SBF21 and OBNLM4 filters. For existing methods, we have set a parameter values which produces high PSNR and SSIM. For SBF filter, we have used an averaging window size of 3 x 3 and total iterations were set to 15. For SRAD filter, the number of iterations was set to 20 and smoothing time step to 0.3. For OBNLM filter, we have used a search window size of 11 x 11, similarity window size of 3 x 3 and the smoothing parameter was set to 0.7.

The experimental results on simulated phantom image of size 256 x 256 are shown in Fig. 2. Figure 2(a) shows the original phantom image. Figure 2(b) shows the image corrupted with Gamma noise (a = = 5) and

Fig. 2(c) to Fig. 2(e) shows the output images of different denoising techniques. Figure 2(f) is the output image of the proposed method which is visually more closer to the ground truth when compared to other denoising techniques. Figure 2(g)-2(i) shows the residual images. From Fig. 2, it can be seen that, the proposed method preserves edges and other fine details better, when compared to other techniques. Further, quantitative analysis was performed on the same phantom image which is corrupted with different levels of Gamma noise. The results are tabulated in Table 1. We have also done experiments on real US images which are shown in Fig. 3(a) to Fig. 3(f). The evaluation of real US images was done by an expert radiologist (with 5 years of experience) by giving a score out of 5. The experts score is shown in Table 2 (0-worst, 5-best). These results further confirm the superior performance of the proposed method over other methods.

Table 1. Illustrating the Change in PSNR and SSIM Values for Different Noise and for Different Methods.

Methods a = 3, ß = 4 a =4, ß = 5 a = 5. ß = 6

PSNR (dB) SSIM PSNR (dB) SSIM PSNR (dB) SSIM

Frost 22.094 0.462 19.92 0.42 17.43 0.39

Lee 22.32 0.45 20.01 0.4012 17.50 0.365

kuan 22.34 0.4546 20.036 0.41 17.7 0.385

SBF 19.72 0.4326 18.6 0.402 17.8535 0.398

DPAD 26.81 0.4725 24.995 0.450 23.36 0.44

SRAD 27.40 0.5053 24.44 0.4605 22.12 0.433

OBNLM 27.43 0.517 24.762 0.4915 22.37 0.4739

Proposed 37.45 0.9855 34.05 0.97 30.58 0.94

(d) (e) (f)

Fig. 3. (a) US Image, (b) Denoised using SBF, (c) using DRAD, (d) using SRAD, (e) OBNLM, (f) Proposed Method.

Table 2. Radiologist Expert Validation (out of 5) of US Images.

Images Original DPAD SBF SRAD OBNLM Proposed Coupled PDE

Image-1 2.5 3.75 3.25 4 4.25 4.75

Image-2 2 3.5 2.75 3.5 3.75 4.5

We have done a performance comparison of the proposed method against other existing methods. We found that local statistics filters like Lee, Frost, and Kuan filters either smooth the edges or preserves the edges depending upon the window size. The performance of the OBNLM filter depends on both searching and similarity window sizes. However, the proposed PDE method outperforms these methods by using a selective diffusion strategy at noisy, edge and interior pixels. This has been made possible with a controlling function. Hence, we showed both quantitative and qualitative improvement compared to previous models.

Second order PDE methods like SRAD and DRAD removes the speckle noise while preserving the edges. But in these models the diffusion control function depends on the coefficient of variation which is calculated based on a user defined homogeneous region. Further, these methods suffer from staircase effects. However, a coupled PDE model can overcome the staircase effects by piecewise planar approximation which results in a less blocky natural looking image. Another advantage is that coupled PDE converges with less number of iterations.

5. Conclusions

A novel coupled PDE model is proposed which combines the SDD model with a fourth order PDE model to remove the speckle in US images. We have assumed that the speckle in US image follows the Gamma distribution. This model is able to reduce the blocky stair case effects produced by second order PDE methods. Further, it converges in less number of iterations. Diffusion was carried out at different speeds at edge, noisy, and interior pixels. This selective diffusion is achieved using controlling speed function which relies on ENI technique. The ML technique was used to estimate the a (shape) and ¡ (scale) parameters of the Gamma distribution. Experiments were carried out on both simulated and real US images to validate the proposed method. Qualitative and quantitative comparative analysis with other similar techniques shows the goodness of the proposed method.

References

[1] J.-S. Lee, Digital Image Enhancement and Noise Filtering by Use of Local Statistics, IEEE Transactions on Pattern Analysis and Machine Intelligence, (2), pp. 165-168, (1980).

[2] D. T. Kuan, A. Sawchuk, T. C. Strand, P. Chavel, et al., Adaptive Noise Smoothing Filter for Images with Signal-Dependent Noise, IEEE Transactions on Pattern Analysis and Machine Intelligence, (2), pp. 165-177, (1985).

[3] V. S. Frost, J. A. Stiles, K. S. Shanmugan and J. C. Holtzman, A Model for Radar Images and its Application to Adaptive Digital Filtering of Multiplicative Noise, IEEE Transactions on Pattern Analysis and Machine Intelligence, (2), pp. 157-166, (1982).

[4] P. Coupé, P. Hellier, C. Kervrann and C. Barillot, Nonlocal Means-Based Speckle Filtering for Ultrasound Images, IEEE Transactions on Image Processing, vol. 18(10), pp. 2221-2229 (2009).

[5] J. Wu and C. Tang, Pde-Based Random-Valued Impulse Noise Removal Based on New Class of Controlling Functions, IEEE Transactions on ¡mage Processing, vol. 20(9), pp. 2428-2438, (2011).

[6] A. P. Witkin, Scale-Space Filtering: A New Approach to Multi-Scale Description, In IEEE International Conference on Acoustics, Speech, and Signal Processing, ICASSP'84, vol. 9, pp. 150-153, (1984).

[7] P. Perona and J. Malik, Scale-Space and Edge Detection using Anisotropic Diffusion, IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 12(7), pp. 629-639, (1990).

[8] F. Catte, P.-L. Lions, J.-M. Morel and T. Coll, Image Selective Smoothing and Edge Detection by Nonlinear Diffusion, SIAM Journal on Numerical Analysis, vol. 29(1), pp. 182-193, (1992).

[9] L. Alvarez, P.-L. Lions and J.-M. Morel, Image Selective Smoothing and Edge Detection by Nonlinear Diffusion, SIAM Journal on Numerical Analysis, vol. 29(3), pp. 845-866, (1992).

[10] Y. Yu and S. T. Acton, Speckle Reducing Anisotropic Diffusion, IEEE Transactions on Image Processing, vol. 11(11), pp. 1260-1270, (2002).

[11] S. Aja-Fernandez and C. Alberola-Lopez, On the Estimation of the Coefficient of Variation for Anisotropic Diffusion Speckle Filtering, IEEE Transactions on Image Processing, vol. 15(9), pp. 2694-2701, (2006).

[12] Y.-L. You and M. Kaveh, Fourth-Order Partial Differential Equations for Noise Removal, IEEE Transactions on Image Processing, vol. 9(10), pp. 1723-1730, (2000).

[13] Z. Tao, H. D. Tagare and J. D. Beaty, Evaluation of Four Probability Distribution Models for Speckle in Clinical Cardiac Ultrasound Images, IEEE Transactions on Medical Imaging, vol. 25(11), pp. 1483-1491, (2006).

[14] P. Sudeep, P. Palanisamy, J. Rajan, H. Baradaran, A. Saba, L Gupta and J. Suri, Speckle Reduction in Medical Ultrasound Images using an Unbiased Non-Local Means Method, Biomedical Signal Processing and Control, vol. 28, pp. 1-8, (2016).

[15] T. Eltoft, Modeling the Amplitude Statistics of Ultrasonic Images, IEEE Transactions on Medical Imaging, vol. 25(2), pp. 229-240, (2006).

[16] R. W. Prager, A. H. Gee, G. M. Treece and L. H. Berman, Decompression and Speckle Detection for Ultrasound Images using the Homodyned k-Distribution, Pattern Recognition Letters, vol. 24(4), pp. 705-713, (2003).

[17] A. Clifford Cohen and B. Jones Whitten, Modified Moment and Maximum Likelihood Estimators for Parameters of the Three-Parameter Gamma Distribution, Communications in Statistics-Simulation and Computation, vol. 11(2), pp. 197-216 (1982).

[18] P. Sudeep, S. I. Niwas, P. Palanisamy, J. Rajan, Y. Xiaojunb, X. Wang, Y. Luo and L. Liu, Enhancement and Bias Removal of Optical Coherence Tomography Images: An Iterative Approach with Adaptive Bilateral Filtering, Computers in Biology and Medicine, vol. 71, pp. 97-107, (2016).

[19] R. C. Gonzalez, Digital Image Processing, Pearson Education India, (2009).

[20] Z. Wang, A. C. Bovik, H. R. Sheikh and E. P. Simoncelli, Image Quality Assessment: From Error Visibility to Structural Similarity, IEEE Transactions on Image Processing, vol. 13(4), pp. 600-612, (2004).

[21] P. C. Tay, S. T. Acton, J. Hossack, et al., A Stochastic Approach to Ultrasound Despeckling, In 3rd IEEE International Symposium on Biomedical Imaging: Nano to Macro, 2006, IEEE, pp. 221-224, (2006).