Speckle Noise Reduction of Medical Ultrasound Images in Complex Wavelet Domain Using Mixture Priors Hossein Rabbani, Mansur Vafadust, Purang Abolmaesumi, Saeed Gazor

Abstract—Speckle noise is an inherent nature of ultrasound images, which may have negative effect on image interpretation and diagnostic tasks. In this paper, we propose several multi-scale nonlinear thresholding methods for ultrasound speckle suppression. The wavelet coefficients of the logarithm of image are modeled as the sum of a noise-free component plus an independent noise. Assuming that the noise-free components have some local mixture distribution (MD), and the noise is either Gaussian or Rayleigh, we derive the minimum mean squared error (MMSE) and the averaged maximum a posteriori (MAP) estimators for noise reduction. We use Gaussian and Laplacian MD for each noise-free wavelet coefficient to characterize their heavy-tailed property. Since we estimate the parameters of the MD using the expectation maximization (EM) algorithm and local neighbors, the proposed MD incorporates some information about the intra-scale dependency of the wavelet coefficients. To evaluate our spatially adaptive despeckling methods, we use both real medical ultrasound and synthetically introduced speckle images for speckle suppression. The simulation results show that our method outperforms several recently and the stateof-the-art techniques qualitatively and quantitatively. Index Terms—Ultrasound image, speckle noise, mixture model, statistical modeling, MAP and MMSE estimator, complex wavelet transform.

I. I NTRODUCTION

U

LTRASONOGRAPHY is one of the most powerful techniques for imaging the internal anatomy (e.g., abdomen, breast, liver, kidney and musculoskeletal). It is relatively inexpensive, noninvasive, harmless for human body and portable, but it suffers from a main disadvantage, i.e., contamination by speckle noise. Speckle noise significantly degrades the image quality and complicates diagnostic decisions for discriminating fine details in ultrasound images. Many techniques have been proposed H. Rabbani is with the Department of Physics and Biomedical Engineering, Isfahan University of Medical Sciences, Isfahan, Iran. M. Vafadust is with the Department of Bioelectrical Engineering, Faculty of Biomedical Engineering, Amirkabir University of Technology, Tehran, Iran. Email: [email protected]. S. Gazor and P. Abolmaesumi are with the Department of Electrical and Computer Engineering, Queen’s University, Kingston, ON, K7L 3N6, Canada. Tel: (613) 533-6591, Fax: (613) 533-6615, Email: [email protected], [email protected].

to reduce this noise [1]–[9]. Early methods use various filters in spatial domain such as average, median and Wiener filter [1]–[5], however, they usually do not preserve accurately all the useful information such as anatomical boundaries in the image. Recently, waveletbased despeckling has been considered [6]–[10]. These methods usually include 1) logarithmic transformation, 2) wavelet transformation, 3) modification of noisy coefficients using shrinkage function, 4) inverse wavelet transformation, and 5) exponential transformation. First step is exerted to convert the multiplicative speckle noise into an additive noise [5], [6] and after that, a wavelet shrinkage technique [11]–[19] is employed for noise reduction. Wavelet-based image denoising methods are formulated as an estimation problem in a Bayesian framework. Thus, the employed probability density function (pdf) for noise-free data and noise (for wavelet coefficients of log-transformed data), and the type of estimator have significant impact on the performance of the noise reduction process. Zong et al. [6] use soft thresholding to remove noise energy in the finer scales. Achim et al. [7] propose a minimum mean square error (MMSE) estimator by assuming a symmetric alpha-stable prior distribution. Bhuiyan et al. [10] present a MAP estimator using a symmetric normal inverse Gaussian (SNIG) prior distribution for modeling the noise-free data. In contrast with above methods, Gupta et al. [9] assume a Rayleigh pdf for the log-transformed noise in wavelet domain which characterize the speckle noise better than a Gaussian pdf. In this paper, we propose a MAP/MMSE estimator assuming Laplacian/Gaussian mixture pdfs for each wavelet coefficient of noise-free log-transformed data and Gaussian/Rayleigh pdf for additive noise in complex wavelet domain [22]. The parameters of MD are obtained locally using EM algorithm [19], [23]. The proposed local mixture pdf is able to model some statistical properties of wavelets such as heavy-tailed nature and intra-scale dependency properties [19]. We obtain several closed-form shrinkage functions (for all combined cases of: MAP/MMSE estimator, Laplacian/Gaussian mixture priors and Gaussian/Rayleigh additive noise) in Section

R +∞

II. We apply them for despeckling of synthetically introduced speckle and real ultrasound images in Section III. Finally, in Section IV, we conclude this paper. II. P ROPOSED A PPROACH Assume that we want to estimate the noise-free image (tissue reflectivity) sl from the recorded ultrasound image xl that is mainly corrupted by a multiplicative speckle noise l [6]: xl

≈ s l l ,

In the context of MDs, the MAP estimation approach does not usually lead into simple useful expressions. In [19]–[21], we suggested an averaged MAP (AMAP) estimator as follows: PI bi,k gi (yk ) i=1 ai,k w , (6) AMAP: w bk = P I a i=1 i,k gi (yk )

(1)

where l is the image index. Ultrasound images also have an additive noise (which is neglected in this paper) whose effect is less significant than the multiplicative noise [5], [6]. If xl = log(xl ), sl and εl represent the logarithms of xl , sl and l , the above equation is written as: xl = sl + εl . Since wavelet transform is a linear transform, the wavelet of xl = sl + εl is written as: yk

= w k + nk ,

wk pi (wk )p(yk |wk )dwk

R +∞ Defining w ei,k = −∞ , and gi (yk ) = p (wk )p(yk |wk )dwk −∞ i R +∞ p (wk )p(yk |wk )dwk , the above equation is simpli−∞ i fied to: PI ei,k gi (yk ) i=1 ai,k w MMSE: w ek = P (5) I i=1 ai,k gi (yk )

where w bi,k is the MAP estimator for ith component of the mixture, i.e., w bi,k = arg maxwk [p(yk |wk )pi (wk )], which is obtained by setting the derivative to zero with respect to w bi,k :

(2)

d log p(yk |wk ) d log pi (wk ) + = 0|wk =wbi,k . dwk dwk

where k is the spatial location of the wavelet coefficient, yk denotes the wavelet coefficients of the logarithm of the ultrasound image, wk denote that of the noisefree image and nk represents the noise. Usually, nk is modeled by Gaussian [6]–[8], [10] or Rayleigh [9], [24], [25] pdfs. In [19], it has been shown that local MDs are useful for noise-free wavelet coefficients. PI We assume local MD described by p(wk ) = i=1 ai,k pi (wk ), where pi (wk )√ can be a Laplacian pdf, pi (wk ) = 2 1√ |wk |), or a Gaussian pdf, pi (wk ) = exp(− σi,k σi,k 2 2 PI wk 1√ exp (− 2σ2 ), ai,k ≥ 0, i=1 ai,k = 1 and I σi,k 2π i,k is the number of mixture components. Assuming such distributions, we can use the following MMSE estimator to retrieve noise-free data from noisy observation: R +∞ wk p(wk )p(yk |wk )dwk w ek = E[wk |yk ] = −∞ p(yk ) R +∞ w p(w )p(y |w )dw k k k k k −∞ , (3) = R +∞ p(w )p(y |w )dw k k k k −∞

(7)

A. Proposed Spatially Adaptive Shrinkage Functions Many researchers model nk with a zero-mean Gaussian noise with variance ς 2 [6]–[8], [10]. In this2 case 1 k) exp(− (yk −w ) in by substituting p(yk |wk ) = √2πς 2 2ς 2 the above formulas, we obtain gi (yk ), w bi,k and w ei,k as in Table I for both Gaussian and Laplacian MD (see [19] for details of derivations). Several researchers have modeled nk2 with two sided nk k| Rayleigh pdf, p(nk ) = |n 2α2 exp(− 2α2 ) [9]. In this −wk | −wk )2 case, we have: p(yk |wk ) = |yk2α exp (− (yk2α ). 2 2 In contrast with the case of additive Gaussian noise, for each Local mixture prior, we derive two different shrinkage functions for MMSE and MAP estimators using Rayleigh noise distribution. Gupta et al. [9] obtained the MAP estimator w bi,k for a Gaussian prior. After some simplification, we derived gi (yk ), w bi,k and w ei,k for Gaussian MD and Laplacian MD and summarized the resulting expressions in Table II (see Appendix A for details of derivations). Eight estimators (refereed to as nonlinear local shrinkage functions) are proposed in this paper. These estimators are named by abbreviation in Table III based on the estimation approach (AMAP or MMSE), the employed noise pdf (Gaussian or Rayleigh), and the employed prior distribution (Gaussian mixture or Laplacian mixture). To implement these estimators the corresponding quantities from Tables I and II must be substituted in (5) and (6).

where p(yk |wk ) represents the conditional pdf yk given wk , and p(wk ) and p(yk ) are the pdf PIof wk and yk . For the proposed mixture pdf p(wk ) = i=1 ai,k pi (wk ), (3) is rewritten as: R +∞ PI wk ( i=1 ai,k pi (wk ))p(yk |wk )dwk −∞ w ek = R +∞ PI ( i=1 ai,k pi (wk ))p(yk |wk )dwk −∞ R +∞ PI i=1 ai,k −∞ wk pi (wk )p(yk |wk )dwk = (4) R +∞ PI i=1 ai,k −∞ pi (wk )p(yk |wk )dwk 2

TABLE I gi (yk ), w bi,k AND w ei,k FOR A DDITIVE G AUSSIAN N OISE , AND MAP AND MMSE E STIMATORS . Gaussian Mixture

Laplacian Mixture

2 yk ) 2(ς 2 +σ 2 ) i,k q 2 2 2π(ς +σi,k )

exp(−

gi (yk ) MAP w bi,k

1 2 1+ ς2 σ i,k

MMSE w ei,k

1 2 1+ ς2 σ i,k

2

√1 2σi,k

y2

exp(− 2ςk2 )[erfcx( σ ς

i,k

yk √ ) ς 2

−

2 √2

(yk − ςσ

yk

i,k

)erfcx( σ ς

i,k (σς i,k

where erfcx(x) =

erfcx R∞ 2 √

π

0

2 √2

y ς 2 yk − √ )+ ς 2

k )+(y + ς − √ k σ

erfcx

2

e−t

+

i,k

yk √ )] ς 2

√ ς2 2 , 0) σi,k

sign(yk ) max(|yk | −

yk

+ erfcx( σ ς

i,k (σς i,k

y ς 2

)erfcx( σ ς yk ) + √ ς 2

k ) + √

i,k

−2tx dt.

TABLE II gi (yk ), w bi,k AND w ei,k FOR A DDITIVE R AYLEIGH N OISE , AND MAP AND MMSE E STIMATOR . Gaussian Mixture −

e

gi (yk )

where

Laplacian Mixture y2

√ √ π erfcx(−zk )−zk π erfcx(zk )) 2 q α 2 2(1+ 2 ) 2πσi,k

α,σi,k k|

(2+zk

σ i,k r 2 +α2 |y |− α4 y 2 +4α4 σ 2 +4α2 σ 4 2|yk |σi,k k k i,k i,k max( 2(α2 +σ 2 ) i,k

MAP w bi,k MMSE w ei,k

2 yk 2σ 2 i,k

D|y

√

α,σ Ax i,k

erfcx

√

=

π

3 4σi,k

√

2 + xασ (2α3 + ασi,k i,k 2),

=

√

α,σi,k

Dx

=

σi,k

√

2

(1 −

√ x 2 2σi,k

AMAP estimator MMSE estimator

−

2 yk 2α2

C

A

+ 2A

√ α π erfcx( σα )), σi,k i,k

TABLE III A BBREVIATIONS OF P ROPOSED S HRINKAGE F UNCTIONS . Gaussian Mixture Prior Rayleigh Gaussian Noise Noise

+

A

C

r

α4 2 2σi,k

•

Laplacian Mixture Prior Rayleigh Gaussian Noise Noise

GGP

GRP

LGP

LRP

GGE

GRE

LGE

LGR

+ α2 , 0)

zk =

), C

erfcx(

yk 2 σi,k

s 2 α2

1 + 22

= erfcx(−

C

+

),

.

σ i,k

E-step: In each iteration, the responsibility factors ri,k are updated by, ai,k gi (yk ) ri,k ← PI , i=1 ai,k gi (yk )

•

B

C

erfcx

α,σ Bx i,k

2 x − σi,k e

2

α,−σi,k ) k|

− C|y

√ |y | 2 − k α,σi,k α,σi,k α,−σi,k α,σi,k α,σi,k 2 σi,k α e ( 2 − −|y | |y | − |y | )+e |yk | |yk | σ k k k i,k √ y2 2 − |yk | − k √ √ σi,k α,−σi,k α,σi,k α π α πe 2α2 √ √ )) sign(yk )( e ( |y | − |y | (1− σ ( σ α ))+ σi,k 2 2σ 2 2 i,k i,k k k i,k α,σi,k α,σi,k α2 α α x √ x 2 −x σi,k σi,k σi,k α 2 −

erfcx(−zk )) erfcx(zk ))

α,σi,k k|

(C|y

√α 2σi,k

sign(yk ) max(|yk | −

z2 σ2 σ2 )+ π (1−2 k 2i,k )( (zk )− 2zk 2(2− i,k 2 α2 α r √ √ 1 + 1 (2+z (−zk )−zk π k π α2 σ2 i,k

erfcx

√ − k α πe 2α2 √ 2 2σi,k 2

,0)

sign(yk ) √

+

for i = 1, · · · , I

(8)

M-step: Then, the parameters ai,k are updated by: X 1 ai,k ← ri,j . (9) |N (k)| j∈N (k)

In this step, the update equations for {σi,k }Ii=1 are approximately estimated as: P 2 j∈N (k) ri,j yj 2 σi,k ← max{ P − ς 2 , δ}. (10) j∈N (k) ri,j

B. Parameter Estimation To implement the introduced shrinkage functions GGP, GRP, GRE, LGP, LGE, LRP and LRE, we need to have MD parameters. The EM algorithm is often used to estimate the MD parameters [23]. A simple description of this iterative algorithm is described in Appendix B. Following a local version of the EM algorithm described in [19]–[21], for each local squared window N (k) with size |N (k)| centered at pixel k, the parameters are iteratively estimated as follows:

where δ is a small positive constant number introduced here to avoid numerical errors, since there 2 are several divisions involving σi,k as denominator. The parameter ς (for the Rayleigh noise, we have ς 2 = 2α2 ) can be estimated using a robust median filter for 3

(a) Real image

complex wavelet transform (DCWT) [22] domain that usually outperforms denoising in other sparse transform domains [19]. For this reason we use the proposed 6tap filter in [26].The window size is chosen |N (k)| = N0 + 2(L−s) for sth scale where L is the depth of wavelet and N0 is an initial value. Figure 1 shows a real ultrasound image and denoised images with LGE (that uses local parameters) and LGEg (LGE that uses global parameters). We can see the effect of local modeling of data in this figure. The boundaries of each brachytherapy seed is more pronounced in the image and the speckle effect is much more reduced in Figure 1(c), therefore it should be more straightforward to develop segmentation techniques to extract the seed locations for image-guided therapy applications. Figure 2 shows a comparison between a real ultrasound image from liver, despeckled images with proposed methods in [5] and [8] and GRP for C = 2.97 and C = 2.24. The results demonstrate that our proposed method, GRP, outperforms the ones compared qualitatively, by significantly reducing the speckle pattern while keeping the anatomical features intact. Other real ultrasound despeckling examples with our methods in this paper for a prostate tissue phantom implanted with brachytherapy seeds, carotid artery and liver can be seen in Figures 3-5. Since LRE and GRE have several oscillations, they do not have acceptable performance in comparison with other shrinkage functions, therefore we have not used them in these figures. σ2 Table IV shows the values of PSNR= 10 log( σ2s ) for n a synthetically speckled image that is proposed in [10]. The σs2 is the variance of noise-free image and σn2 is the variance of unit mean complex Gaussian random field used to produce the speckled image [7]. We can see in this table that our proposed methods have higher PSNR in comparison with the ones compared. As indicated by bold-fonts in this table, LGP and LRP outperform other methods respectively for low level noise values σn < 0.5 and for high level noise σn > 0.5. The correlation coefficient (CoC), edge preservation index (EPI) and structural similarity index (SSI) are other quality metrics defined as follows [27], [28]: P (s − s¯)(ˆ s − ¯sˆ) CoC = pP (12) P (s − s¯)2 (ˆ s − ¯sˆ)2 P ¯s) ¯ (∆s − ∆s)(∆ˆ s − ∆ˆ EPI = pP (13) P ¯s)2 ¯ 2 (∆ˆ (∆s − ∆s) s − ∆ˆ (2¯ s¯sˆ + 2.55)(2σsˆs + 7.65) SSI = (14) 2 (¯ s2 + ¯sˆ + 2.55)(σs2 + σsˆ2 + 7.65)

ultrasound (b) Denoising with (c) Denoising with loglobal parameters cal parameters

Fig. 1. Despeckling with LGE and LGEg for an ultrasound image captured from a prostate tissue phantom implanted with brachytherapy seeds.

(a) Real ultrasound (b) Homomorphic (c) Proposed method image from liver Wiener [5] (7 × 7) in [8] for T = 10

(d) Soft Thresholding (e) GRP for C = 2.97 (f) GRP for C = 2.24 [15] Fig. 2. Comparison between despeckling using GRP with other methods for liver.

the finest scale of noisy wavelet coefficients [9], [10]: ς = Cmedian {|yk | |∀k ∈ D1,2 } .

(11)

where C ∈ [2.2, 3] is a smoothing factor and D1,2 is the set of subband HH indexes in first and second scale. We must reduce the value of C where the texture information of image is more important. In contrast, we shall increase the value of C where image discontinuities are more informative. A too large smoothing factor causes blurring of the edges, in contrast a too small smoothing factor prevents effective noise removal. Therefore, the value of C must be tuned according to the clinical preference, during initial trial for each specific medical application. An example is provided in Figure 2 illustrating the effect of C. III. E XPERIMENTAL R ESULTS In this section, we apply the introduced shrinkage functions for reduction of synthetical and real noises. We implement our despeckling algorithms in discrete

where s, sˆ are the original and denoised images respectively, s¯ shows the mean of s, ∆(s) is the highpass 4

[1 0] LG

LR P

G

G RP

G E

LG P

12.67 12.56 12.42 12.26 12.06 11.82 11.60 11.32 11.06

22.58 19.20 16.83 15.07 13.56 12.03 11.32 10.32 9.41

22.18 19.39 17.26 15.59 14.16 12.88 11.91 10.99 10.11

22.41 19.65 17.71 16.29 15.06 14.16 13.29 12.57 11.90

23.33 20.49 18.48 17.07 15.63 14.53 13.49 12.96 12.11

23.29 20.44 18.35 16.91 15.46 14.43 13.39 12.76 11.97

23.16 20.39 18.40 17.11 15.74 14.81 13.83 13.02 12.48

23.35 20.46 18.38 16.92 15.47 14.32 13.30 12.71 11.83

14.25 13.37 12.62 12.16 11.64 11.12 10.75 10.47 10.17

E

Pr op

os ed

ed op os

ho d M et

et h M

k rin Sh ye s

in

in od

5] [1

c ph i om or

Pr

= 0.2 = 0.3 = 0.4 = 0.5 = 0.6 = 0.7 = 0.8 = 0.9 = 1.0

Ba

σn σn σn σn σn σn σn σn σn

H om

N oi

se

Va ria

nc e

W ie ne

r[

[8 ]

5]

TABLE IV C OMPARISON BETWEEN PSNR S OF PROPOSED METHODS AND SIMILAR COMPETITIVE METHODS IN THE LITERATURE .

Fig. 5. Despeckling with the methods proposed in the current paper for an ultrasound image captured from human liver. From top left, clockwise: Real ultrasound image, GRP, LGP, LGE, LRP, GGE. Fig. 3. Despeckling with the methods proposed in the current paper for an ultrasound image captured from a prostate tissue phantom implanted with brachytherapy seeds. From top left, clockwise: Real ultrasound image, LRP, LGP, LGE, GRP, GGE.

filtered s using the discrete Laplacian operator (a 3 × 3 pixel standard approximation) and σsˆs represents the covariance between s and sˆ. The CoC and SSI indicate the measure of similarity between the original and denoised images and EPI shows the measure of edge preservation. The best value of these parameters for good visual quality is 1. We use CoC, EPI and SSI quality metrics to better evaluate our algorithms. Table V shows these metrics for a synthetic image including regions with uniform intensity, strong scatters and sharp edges (as proposed in Figure 4 of [8]) for σn = 0.3 and 0.6. For this reason, we implement our algorithms in redundant wavelet domain [29] as in [28] and compare them with 1) median filter, 2) homomorphic Wiener filter [5], 3) soft thresholding [30], 4) Pizurica’s method which is a robust non-homomorphic method, 5) proposed technique in [31] using a speckle reducing anisotropic diffusion filter, 6) GNDShrink [28] which models the noise-free wavelet coefficients by Gaussian pdf and speckle in

Fig. 4. Despeckling with the methods proposed in the current paper for an ultrasound image captured from carotid artery. From top left, clockwise: Real ultrasound image, GRP, LGP, LGE, LRP, GGE.

5

TABLE V C OMPARISON BETWEEN C O C, SSI AND EPI OF PROPOSED

cases, we prefer to make a tradeoff between the bias and the the variance. For example, the bias using soft thresholding is greater than using hard thresholding; however, soft thresholding is often preferable because of the variance. From Table VI, among our methods, the MMSE estimators (LGE and GGE) have the smallest biases, where as the shrinkage functions obtained using a Rayleigh noise have the largest biases. In addition, comparing N0 = 7 with N0 = 3 in this table, we conclude that the larger the window size, the larger the bias. Figure 6 shows the proposed shrinkage functions in this paper. In this figure, we see that shrinkage functions obtained using MMSE estimator (LGE and GGE) are smoother than other shrinkage functions (LGP, LRP and GRP). As a result, they lead to smaller bias (Table VI). In contrast, according to Table IV, we observe that shrinkage functions in a MAP framework have higher PSNR. We also observe that for Gaussian noise and MMSE estimator, Gaussian mixture prior produces higher PSNR for σn ≥ 0.6, However for σn < 0.6, using Laplacian MD leads to higher PSNR. Using Laplacian mixture prior, Table IV indicates that Rayleigh noise model is preferred for σn ≥ 0.5 in terms of PSNR in contrast, Gaussian noise model is preferred for σn < 0.5. The Gaussian noise assumption results in smaller bias compared with Rayleigh noise assumption.

METHODS IN REDUNDANT DISCRETE WAVELET DOMAIN AND SEVERAL COMPETITIVE METHODS IN THE LITERATURE FOR AN ARTIFICIAL ULTRASONIC TEST IMAGE .

Speckled input Median filter Method in [5] Soft Thresh. [30] Method in [31] Method in [8] GNDShrink [28] GNDThresh [28] LGE GRP GGE LGP LRP

CoC 0.969 0.964 0.971 0.975 0.979 0.985 0.992 0.994 0.988 0.986 0.991 0.989 0.989

σn = 0.6 EPI .785 .613 .678 .838 .853 .926 .951 .960 .950 .950 .974 .961 .975

SSI .596 .752 .785 .642 .865 .921 .953 .958 .986 .985 .989 .989 .987

CoC .992 .973 .978 .993 .983 .991 .998 .998 .997 .997 .998 .998 .998

σn = 0.3 EPI .929 .668 .713 .944 .900 .980 .988 .990 .984 .991 .992 .988 .995

SSI .811 .881 .889 .830 .922 .983 .985 .985 .997 .997 .998 .997 .998

wavelet domain by two-sided Generalized Nakagami distribution and 7) GNDThresh [28] which models the noise-free signal with Laplacian pdf and uses two-sided Generalized Nakagami distribution for speckle noise. From Table V, we understand that LRP has the best EPI results, and GGE outperforms the others in terms of SSI, however, GendShrink usually outperforms in terms of CoC. Table VI compares the proposed methods in this paper in terms of PSNR, CoC, SSI and EPI for a phantom ultrasonic image with σn = 0.7 using window sizes |N (k)| = 3 + 2(5−k) and |N (k)| = 7 + 2(5−k) for k th scale. From this table, we see that the window size has important impact especially on EPI, i.e., the bigger the window size N0 = 7, the smaller EPI. In this table, LGP3 represents the shrinkage function obtained using a Laplacian mixture prior of three components (in contrast to LGP which uses two mixture components). From this table, by comparing the values of LGP, LGP3 (AMAP estimators) and non-mixture LGP (the MAP estimator non-mixture Laplacian model), we understand that AMAP estimators perform better than MAP only, but increasing the mixture components does not lead to major improvement to the results (the performance gain in LGP3 compared with LGP is not significant). Our proposed estimators are biased, since they are implemented in the logarithm domain which involves nonlinear transformations. Table VI also shows the biases and variances of our methods, homomorphic Wiener filter and Pizurica’s method. As we expect, the Pizurica’s method [8] which is a non-homomorphic technique is almost unbiased. However, from the last column of Table VI, we also see that our methods have less variance than Pizurica’s and homomorphic Wiener filter. In many

IV. C ONCLUSION In this paper, we derived several shrinkage functions using local mixture priors for modeling log-transformed speckled images in complex wavelet domain. We used MAP and MMSE estimators, Gaussian mixture and Laplacian mixture priors, and additive Gaussian and Rayleigh noise. The simulation results reveals that our nonlinear local thresholding functions outperform a few of similar competitive methods in the literature, both visually and in terms of quality metrics. The proposed methods in this paper may also be used for processing of other medical imaging modalities corrupted with various types of noise. The results may be improved using more complex prior models which describe more accurately other statistical properties of wavelets, e.g., mixture of local circular symmetric Laplacian pdfs models the inter-scale dependency [32]. A PPENDIX A E STIMATORS FOR TWO - SIDED R AYLEIGH NOISE Here, we drive the expressions of gi (yk ), w bi,k and w ei,k for two-sided Rayleigh noise. 6

TABLE VI C OMPARISON BETWEEN PSNR, C O C, SSI, EPI, B IAS AND VARIANCE OF PROPOSED METHODS IN THIS PAPER FOR A PHANTOM ULTRASONIC IMAGE WITH σn = 0.7 USING WINDOW SIZES |N (k)| = 3 + 2(5−s) AND |N (k)| = 7 + 2(5−s) FOR sth SCALE . PSNR

CoC

Speckled image non-mixture LGP Method in [8] LRP LGE LGP3 LGP GGE GRP

12.53 18.18 16.58 18.86 18.82 18.87 18.85 18.85 18.07

.9732 .9924 .9890 .9935 .9934 .9935 .9935 .9935 .9923

Homomorphic Wiener [5] non-mixture LGP Method in [8] LRP LGE LGP3 LGP GGE GRP

16.14 18.31 16.69 18.87 18.80 18.89 18.86 18.91 18.16

.9879 .9926 .9892 .9935 .9934 .9935 .9935 .9935 .9924

SSI N0 .9731 .9923 .9890 .9927 .9928 .9929 .9928 .9928 .9913 N0 .9872 .9925 .9893 .9927 .9928 .9929 .9929 .9929 .9915

EPI =3 .4227 .8011 .6685 .8271 .8266 .8222 .8218 .8258 .7943 =7 .6447 .8011 .6880 .8223 .8210 .8174 .8165 .8210 .7941

Bias

Variance

.0090 .3686 .0027 .5114 .3802 .3984 .3970 .3753 .5844

32.25 8.82 12.74 7.58 7.61 7.53 7.55 7.56 9.05

.7404 .3744 .0795 .5141 .3815 .4006 .3987 .3808 .5863

14.11 8.56 12.43 7.54 7.60 7.50 7.54 7.46 8.87

Gupta et al. [9] solved the above equation for a Gaussian prior and obtained: max( w bi,k =

2|yk |

2 σi,k α2

r +|yk |− 2(1+

2 +4σ 2 +4 yk i,k σ2 i,k α2

σ4 i,k α2

, 0)

)

sign(yk )

. (16)

For a √Laplacian pdf we have log pi (wk ) √ 2 − log(σi,k 2) − σi,k |wk |. This leads to √ −1 yk − w bi,k 2 − + sign(w bi,k ) = 0, yk − w bi,k α2 σi,k

Fig. 6. paper.

=

(17)

and by solving the above equation, the MAP estimator comes out to be q 4 α2 max(|yk | − √2σ + 2σα2 + α2 , 0) i,k i,k w bi,k = . (18) sign(yk )

A comparison between proposed shrinkage functions in this

B. Distribution of the noisy observations To calculate the proposed shrinkage functions in this paper we need to obtain the distribution of noisy observation: Z +∞ |yk − wk | − (yk −w2k )2 2α gi (yk ) = pi (wk ) e dwk . (19) 2α2 −∞

A. MAP estimator For a two-sided Rayleigh noise, we have p(yk |wk ) = −wk )2 exp(− (yk2α ). Substituting this pdf in (7), we 2 have: |yk −wk | 2α2

−1 yk − w bi,k d log pi (wk ) + + |wk =wbi,k 2 yk − w bi,k α dwk

1) Gaussian prior: P4 In this case, using (19) we have gi (yk ) = Il where Il are given in l=1 r 1 = 0. (15) (20), where ξi,k = yk α22 + σ22 and erfcx(x) = i,k

7

Z I1

yk

= yk −∞

Z

−

e e √ σi,k 2π −

yk

Z

∞

= −yk

I3

yk

Z I4

∞

= yk

R∞ 0

2

e−t

−2tx

−

−

(yk −wk )2 2α2

2 wk 2σ 2 i,k

−

2 wk 2σ 2 i,k

e e √ σi,k 2π −

2 wk 2σ 2 i,k

e e √ wk σi,k 2π

dt. Denoting zk =

ξi,k 2 , σi,k

dwk =

(y −w )2 − k2α2k

−

2α2

2α2

2α2

e 2α2

√ ξi,k q ξi,k πerfcx(− 2 ), σ 2 i,k 2πσi,k 2 yk 2σ 2 i,k

−

e

(−1 +

dwk = − 2(1 +

(yk −wk )2 2α2

(y −w )2 − k2α2k

2 yk 2σ 2 i,k

−

2α2

e e √ = − wk σi,k 2π −∞

I2

√2 π

2 wk 2σ 2 i,k

−

dwk = − −

e

e 2α2

2 yk 2σ 2 i,k

2 yk 2σ 2 i,k

q

2 2πσi,k

(1 +

dwk = 2(1 +

we obtain:

ξ ξi,k √ πerfcx(− σi,k 2 )) α2 i,k α2 2 ) σi,k

q 2 2πσi,k

(20a)

,

√ ξi,k ξi,k πerfcx( 2 ), σi,k

ξi,k √ ξ πerfcx( σi,k 2 )) α2 i,k

q

α2 2 ) σi,k

.

(20b)

(20c)

(20d)

2 2πσi,k

2) Laplacian prior: For a Laplacian prior, the numerator of (25) is an odd function and is equal to D1 + D2 + D3 for positive input values where Dl s are given √ √ α,σ 2 in (28) and Ax i,k = 4σ3π (2α3 + ασi,k + xασi,k 2).

2 yk 2σ 2 i,k

( zk2√π + erfcx(−zk ) − erfcx(zk )) q √i,k . (21) 2 α,σ 2 x 2 i,k 2 − σα2 +2A−x i,k erfcx( σαi,k ), So, defining Bα,σ = 2σ 2(1 + σα2 ) 2σi,k x i,k i,k i,k and using (24) for gi (yk ), we obtain (29). 2) Laplacian prior: In this case, (19) is written as: A PPENDIX B √ Z +∞ − σ 2 |wk | EM A LGORITHM 2 i,k (yk −wk ) |yk − wk |e PI √ gi (yk ) = e− 2α2 dwk . (22) Suppose a MD p(x) = a p (x) where 2 R PI i=1 i,k i 2α σi,k 2 −∞ pi (x; θi )dx = 1, ai ≥ 0, i=1 ai = 1 and I is the The function gi (yk ) is even and is calculated by gi (yk ) = number of mixture components. We want to estimate the E1 + E2 + E3 where Ei is defined in (23). Therefore, parameters {ai , θi }Ii=1 from the observation {xn }N n=1 α,σ α,σi,k x data where N is the number of observed samples. Since = defining Cx i,k = erfcx(− σαi,k + α√ ) and D x 2 √ 2 x the mixture component which has produced each ob− √ π α e σi,k α √ (1 − erfcx( )), we obtain served sample is unknown, our observation set {xn }N σi,k σi,k n=1 σi,k 2 is an incomplete data. In such as case, the Expectation 2 yk √ Maximization (EM) algorithm is often employed for α πe− 2α2 α,σi,k α,σi,k α,−σi,k parameter estimation. This algorithms maximizes the √ gi (yk ) = D|yk | + (C − C ). (24) |yk | |yk | 2 2σi,k 2 expected value of likelihood function by defining e

w ei,k = zk

Q=

C. MMSE estimator

N X I X

ri (n) ln(ai pi (xn ))

(30)

n=1 i=1

The MMSE estimator for two-sided Rayleigh noise is defined as: R +∞ w ei,k =

−∞

−wk | − e wk p(wk ) |yk2α 2

gi (yk )

(yk −wk )2 2α2

dwk

.

where {ri (n)}Ii=1 is an auxiliary variable namely responsibility factor that for each observed data xn represents how likely that observed data was produced by each component {pi (x)}Ii=1 . The EM algorithm is an iterative algorithm and consists of two steps namely E-step and M-step. After the initialization of the parameters, the responsibility factors are calculated in E-step by: ai pi (xn ) , for i = 1, · · · , I (31) ri (n) ← PI i=1 ai pi (xn )

(25)

1) Gaussian prior: The numerator P4 in (25) for a Gaussian pdf can be written as sum l=1 Jl , where Jl s are defined in (26). Thus, using (21) for gi (yk ), the MMSE estimator is obtained as in (27). 8

Z E1

√ 2 σi,k

0

= −∞

(yk − wk )e √ σi,k 2

=

=

=

=

−

yk

wk

√

dwk = −

y2

2 − k2 2α 2 e 4σi,k

√ yk α (−σi,k + α πerfcx( √ + )), (23a) α 2 σi,k

(yk −wk )2

(yk −wk )2

−

e −e √ wk2 σi,k 2π −∞ wk ∞

= yk yk

Z J4

2 wk 2

2 −wk 2σ 2 i,k

yk

Z J3

2

e 2σi,k e− 2α2 √ wk 2α2 σi,k 2π −∞

= yk

Z J2

−

yk

Z J1

(yk −wk )2 2α2

2α2 √

E3

e−

e σi,k e− 2α2 √ dwk (yk − wk ) 2α2 σi,k 2 0 √ √ 2  2y y2 yk √  − 2yk − σ k σi,k − 2αk2 i,k ) (e − e α 2π e σi,k erfcx( σαi,k ) − e− 2α2 erfcx( αy√k2 − σαi,k ) + α√ π − , 2 4σi,k √ (yk −wk )2 √ √ Z ∞ − 2 wk k √ e σi,k e− 2α2 2 − σ2y α i,k (−σ √ dw = − e )) (wk − yk ) k i,k + α πerfcx( 2 2 2α 4σ σ σ 2 i,k yk i,k i,k Z

E2

wk

∞

wk2

= yk

−

(yk −wk )2 2α2

2α2

2 wk 2σ 2 i,k

e −e √ σi,k 2π 2 −wk 2σ 2 i,k

e e √ σi,k 2π

−

−

dwk = e

(yk −wk )2 2α2

2α2

(yk −wk )2 2α2

2α2

α2 2 σi,k

2(1 + r 2 −yk 2σ 2 i,k

3yk 2 σi,k

2 α2

ξi,k ξi,k √ q (−1 + 2 πerfcx(− 2 )), α σi,k 2 ) 2πσi,k

1 +

−

e 2(1 + r

dwk = e

2 σ2 i,k

2 yk α4 ( α12 +

− (1 +

2α2 ( α12 +

dwk = −yk

2 −yk 2σ 2 i,k

3yk 2 σi,k

1 2 σi,k

2 yk 2σ 2 i,k

q

α2 2 ) σi,k

2 α2

(23c)

2 yk 2σ 2 i,k

−

e

dwk = yk

(23b)

1 +

2 σ2 i,k

(1 + 2 2πσi,k

+ (1 +

2α2 ( α12 +

)3/2

q

)

) √2π erfcx(

−ξi,k ) 2 σi,k

,(26b)

2 πσi,k

ξi,k √ ξi,k πerfcx( 2 )), α2 σi,k

2 yk α4 ( α12 +

1 2 σi,k

1 σ2 i,k

√

Z D1

0

= −∞

2

(yk − wk )e σi,k √ σi,k 2 √

)

) √2π erfcx( σi,k 2 ) i,k

q 2 )3/2 πσi,k

D3

2

(27)

σi,k

(yk −wk )2

wk e− 2α2 2α2

wk

. (26d)

2 yk

dwk = e− 2α2 (

(yk −wk )2

α2 yk α α,σi,k √ + )), 2 − A|yk | erfcx( 2σi,k α 2 σi,k

(28a)

(yk − wk )e σi,k wk e− 2α2 √ = dwk (28b) 2α2 σi,k 2 0 √ y y2 − σ k√2 α2 2yk α α yk α,σ α,σ − 2αk2 i,j = e (− 2 + + A−|yi,k erfcx( )) − e (A−|yi,k erfcx( − √ ), | k k| 2σi,k 4σi,k σi,k σi,k α 2 √ − 2 (yk −wk )2 √ Z ∞ w k − −y k 2α2 √ (wk − yk )e σi,k wk e α2 2yk α α,σ √ = dwk = e σi,j 2 (− 2 + + A−|yi,k erfcx( )). (28c) 2 k| 2α 2σ 4σ σ σ 2 i,k i,k yk i,k i,k Z

D2

−

wk

(26c)

ξ

1 σ2 i,k

2 2 √ pπ σ2 zk σi,k 2zk 2(2 − αi,k (1 − 2 J1 + J2 + J3 + J4 2 ) + 2 α2 )(erfcx(zk ) − erfcx(−zk )) q = . w ek = √ √ 1 1 gi (yk ) πerfcx(−zk ) − zk πerfcx(zk )) 2 (2 + zk 2 +

α

(26a)

yk

9

y2

− 2αk2

sign(yk )(e D1 + D2 + D3 w ek = = √ 2 |y | sign(yk )gi (yk ) − k σi,k e

σi,k

√

2

2 ( σα2 i,k

−

(1 −

√ α π α σi,k erfcx( σi,k ))

α,σ α,σ A−|yi,k C|yk |i,k k|

Using these responsibility factors in the M-step, the parameters {ai }Ii=1 are updated using Lagrange multiplier: ai ←

N 1 X ri (n), N n=1

(32)

σi2 ←

i

PN

2 n=1 ri (n)xn . PN n=1 ri (n)

+

√ − k2 2α α πe √ 2 2σi,k 2

+e

−

√ |yk | 2 σi,k

α,σ

B|yk |i,k )

y2

α,σ (C|yk |i,k

−

.

(29)

α,−σ C|yk | i,k )

[10] M. I. H. Bhuiyan, M. O. Ahmad, M. N. S. Swamy, “Waveletbased despeckling of medical ultrasound images with the symmetric normal inverse Gaussian Prior,” in Proc. ICASSP, April 2007, vol. 1, pp. 721-724. [11] D. L. Donoho, I. M. Johnstone, “Ideal spatial adaptation by wavelet shrinkage,” Biometrika, vol. 81, no. 3, 1994. [12] D. L. Donoho, I. M. Johnstone, “Adapting to unknown smoothness via wavelet shrinkage,” Journal of American Statistical Association, vol. 90, no. 432, pp. 1200-1224, 1995. [13] D. L. Donoho, “De-noising by soft-thresholding,” IEEE Trans. on Information Theory, vol. 41, pp. 613-627, May 1995. [14] M. S. Crouse, R. D. Nowak, R. G. Baraniuk, “Wavelet-based statistical signal processing using hidden Markov models,” IEEE Trans. on Signal Processing, vol. 46, no. 4, pp. 886-902, Apr. 1998. [15] S. G. Chang, B. Yu, M. Vetterli, “Adaptive wavelet thresholding for image denoising and compression,” IEEE Trans. on Image Processing, vol. 9, no. 9, pp. 1532-1546, 2000. [16] J. Portilla, V. Strela, M. J. Wainwright, E. P. Simoncelli, “Image denoising using gaussian scale mixtures in the wavelet domain,” IEEE Trans. on Image Processing, vol. 12, pp. 1338-1351, 2003. [17] M. K. Mihcak, I. Kozintsev, K. Ramchandran, P. Moulin, “Low complexity image denoising based on statistical modeling of wavelet coefficients,” IEEE Signal Processing Letters, vol. 6, pp. 300-303, Dec. 1999. [18] L. Sendur, I.W. Selesnick, “Bivariate shrinkage with local variance estimation,” IEEE Signal Processing Letters, vol. 9, no. 12, pp. 438-441, Dec. 2002. [19] H. Rabbani, S. Gazor, M. Vafadoost, “Sparse domain image denoising employing local mixture models,” submitted to IEEE Trans. on Pattern Analysis and Machine Intelligence. [20] H. Rabbani, M. Vafadust, I. Selesnick, “Wavelet based image denoising with a mixture of Gaussian distributions with local parameters,” in Proc. Int. Symp. ELMAR focused on Multimedia Signal Processing and Communications, June 2006, pp. 85-88. [21] H. Rabbani, M. Vafadust, S. Gazor, “Image denoising based on a mixture of Laplace distributions with local parameters in complex wavelet domain,” in Proc. ICIP, Oct. 2006, pp. 2597-2600. [22] I. W. Selesnick, R. G. Baraniuk, N. Kingsbury, “The dualtree complex wavelet transforms - A coherent framework for multiscale signal and image processing,” IEEE Signal Processing Magazine, vol. 22, no. 6, pp.123-151, Nov. 2005. [23] A. P. Dempster, N. M. Laird, D. B. Rubin, “Maximum Likelihood from Incomplete Data via the EM Algorithm,” Journal of the Royal Statistical Society. Series B (Methodological), vol. 39, no. 1, pp. 1-38, 1977. [24] R. F. Wagner, S. W. Smith, J. M. Sandrik, H. Lopez, “Statistics of speckle in ultrasound B-scans,” IEEE Trans. on Sonics and Ultrasonics, vol. 30, no. 3, pp. 156-163, 1983. [25] A. N. Evans, M. S. Nixon, “Mode filtering to reduce ultrasound speckle for feature extraction,” IEE Proc. on Vision, Image and Signal Processing, vol. 142, no. 2, 1995.

The parameters {θi }Ii=1 are obtained by maximization of Q. The update equations for {θi }Ii=1 depend on the pdf functions {pi (x)}Ii=1 . For example for a Gaussian PI x2 MD p(x) = i=1 ai,k σ √1 2π exp(− 2σ 2 ), we have i

−

α,σ α,−σ A|yk |i,k C|yk | i,k )

(33)

The sequence of E-M steps proceeds until the parameters satisfy some convergence conditions. R EFERENCES [1] J. U. Quistgaard, “Signal acquisition and processing in medical diagnostic ultrasound,” IEEE Signal Processing Magazine, pp. 67-74, Jan. 1997. [2] M. Karaman, M. A. Kutay, G. Bozdagi, “An adaptive speckle suppression filter for medical ultrasonic imaging,” IEEE Trans. on Medical Imaging, vol. 14, pp. 283-292, June 1995. [3] J. W. Goodman, “Some fundamental properties of speckle,” Journal of Optical Society of America, vol. 66, pp. 1145-1150, Nov. 1976. [4] T. Loupas, W. N. Mcdicken, P. L. Allan, “An adaptive weighted median filter for speckle suppression in medical ultrasonic images,” IEEE Trans. on Circuits and Systems, vol. 36, pp. 129-135, Jan. 1989. [5] A. K. Jain, Fundamental of Digital Image Processing. Englewood Cliffs, NJ: Prentice-Hall, 1989. [6] X. Zong, A. F. Laine, E. A. Geiser, “Speckle reduction and contrast enhancement of echocardiograms via multi-scale nonlinear processing,” IEEE Trans. on Medical Imaging, vol. 17, pp. 532540, Aug. 1998. [7] A. Achim, P. Tsakalides, A. Bezarianos, “Novel Bayesian multiscale method for speckle removal in medical ultrasound images,” IEEE Trans. on Medical Imaging, vol. 20, pp. 772-783, 2001. [8] A. Pizurica, W. Philips, I. Lemahieu, M. Acheroy, “A versatile wavelet domain noise filtration technique for medical imaging,” IEEE Trans. on Medical Imaging, vol. 22, pp. 323-331, 2003. [9] S. Gupta, R. C. Chauhan, S. C. Saxena, “Locally adaptive wavelet domain Bayesian processor for denoising medical ultrasound images using speckle modelling based on Rayleigh distribution,” IEE Proc. on Vision, Signal and Image Processing, vol. 152, pp. 129-135, 2005.

10

Purang Abolmaesumi received the B.Sc. and M.Sc. degrees from Sharif University of Technology, Tehran, Iran, and the Ph.D. degree from the University of British Columbia, Vancouver, Canada, all in Electrical Engineering. He then joined the School of Computing at Queen’s University, Kingston, Ontario, Canada where he is now an Assistant Professor. Dr. Abolmaesumi’s research interests are in image-guided surgery, medical imaging,

[26] N. G. Kingsbury, “A dual-tree complex wavelet transform with improved orthogonality and symmetry properties,” in Proc. ICIP, Sept. 2000, vol.2, pp. 375–378. [27] Z. Wang, A. C. Bovik, H. R. Sheikh, E. P. Simoncelli, “Image quality assessment: From error visibility to structural similarity,” IEEE Trans. on Image Processing, vol. 13, no.4, pp. 600-612, Apr. 2004. [28] S. Gupta, R. C. Chauhan, S. C. Saxena, “A versatile technique for visual Enhancement of Medical Ultrasound Images,” Digital Signal Processing (Elsevier), vol. 17, no. 3, pp. 542-560, May. 2007. [29] M. Lang, H. Guo, J.E. Odegard, C.S. Burrus, R.O. Wells Jr., “Noise reduction using an undecimated discrete wavelet transform,” IEEE Signal Processing Letters, vol. 3, pp. 10–12, 1996. [30] S. Gupta, R.C. Chauhan, S.C. Saxena, “A wavelet based statistical approach for speckle reduction in medical ultrasound images,” IEE J. Int. Fed. Med. Biol. Eng., vol. 42, pp. 189–192, 2004. [31] Y. Yu, S.T. Acton, “Speckle reducing anisotropic diffusion,” IEEE Trans. on Image Processing, vol. 11, no. 11, pp. 1260–1270, 2002. [32] H. Rabbani, M. Vafadust, I. Selesnick, S. Gazor, “Image denoising employing a mixture of circular symmetric Laplacian models with local parameters in complex wavelet domain,” in Proc. ICASSP, April 2007, vol. 1, pp. 805-808.

robotics and haptics.

Saeed Gazor (S’94-M’95-SM’98) received the B.Sc. degree in electronics and the M.Sc. degree in communication systems from Isfahan University of Technology, Isfahan, Iran, in 1987 and 1989, respectively, both with the highest standing, and the Ph.D. degree (with highest honors) in signal and image processing from D´epartement Signal, Ecole Nationale Sup´erieure des T´el´ecommunications (ENST), Paris (T´el´ecom Paris), France, in 1994. From 1995 to 1998, he was with the Department of Electrical and Computer Engineering, Isfahan University of Technology. From January 1999 to July 1999, he was with the Department of Electrical and Computer Engineering, University of Toronto, Toronto, ON, Canada. Since 1999, he has been on the Faculty at Queen’s University at Kingston, ON, Canada, where he currently holds the position of Associate Professor in the Department of Electrical and Computer Engineering. He is also cross-appointed to the Department of Mathematics and Statistics at Queen’s University. His main research interests are array signal processing, statistical and adaptive signal processing, speech processing, MIMO communication systems, networking, analog adaptive circuits, channel modeling, and information theory. Dr. Gazor received a number of awards, including a Provincial Premier’s Research Excellence Award, a Canadian Foundation of Innovation Award, and an Ontario Innovation Trust Award.

Hossein Rabbani was born in Iran in 1978. He received the B.Sc. degree in communication engineering from Isfahan University of Technology, Isfahan, Iran, in 2000 with the highest honors, and the M.Sc. and Ph.D degrees in bioelectrical engineering in 2002 and 2008, respectively, from Amirkabir University of Technology (Tehran Polytechnic), Tehran, Iran. From January 2007 to July 2007, he was with the Department of Electrical and Computer Engineering, Queen’s University, Kingston, ON, Canada, as a Visiting Researcher. His main research interests are multidimensional signal processing, multiresolution transforms, probability models of sparse domain’s coefficients and image restoration.

Mansur Vafadust received his B.sc. degree in electrical engineering (electronics) from Sharif University of Technology, Tehran, Iran, in 1985, his M.sc. degree in the same field from Amirkabir University of Technology, Tehran, Iran, in 1988, and his Ph.D. in biomedical engineering from the University of New South Wales, Sydney, Australia, in 1994. His current research interests include image enhancement and noise reduction from medical images. He is currently with the Department of Biomedical Engineering, Amirkabir University of Technology, Tehran, IRAN.

11