成子Biomedical Signal Processing and Control 6 (2011) 129–138
Contents lists available at ScienceDirect
Biomedical Signal Processing and
Control
j o u r n a l h o m e p a g e :w w w.e l s e v i e r.c o m /l o c a t e /b s p
c
Speckle filtering of ultrasonic images using a modified non local-bad algorithm
Y.Guo,Y.Wang ∗,T.Hou
Department of Electronic Engineering,Fudan University,No.220,Handan Rd.,Shanghai 200433,China
a r t i c l e i n f o Article history:
Received 26February 2010
Received in revid form 9October 2010Accepted 21October 2010
载不动许多愁Available online 18 November 2010Keywords:
Speckle reduction
Maximum likelihood estimator Non-local means algorithm Ultrasonic image
厦门岛内a b s t r a c t
The speckle reduces the contrast and obscures important details in ultrasound images,which limits the accuracy in clinical diagno.Many methods have been developed to reduce the speckle.How-ever,they suffer from veral limitations,difficult to reach a balance between the noi attenuation and edge prerving.Recently,a newly developed filter,namely Non-Local Means (NL-means),works well in the image denoising,jointly enhancing the edge information.Since the NL-means filter is origi-nally designed for Gaussian noi removal,it cannot be directly applied to ultrasonic image despeckling with the Rayleigh distribution noi.In this paper,we propod a modified non local-bad (MNL)filter to adapt for the speckle reduction.Experiments were carried out on synthetic images to find optimum parameters.Besides,this filter has been compared with other six state-of-the-art denoising methods on synthetic data and real ultrasonic images.Results show that the MNL method is capable of effectively reducing the speckle noi while well prerving tissue boundaries for ultrasonic images.
© 2010 Elvier Ltd. All rights rerved.
1.Introduction
The ultrasound imaging has become widely ud in medical diagno becau of its noninvasive nature,low cost,and contin-uing improvements in the image quality [1].The speckle,a form of locally correlated multiplicative noi,is caud by the interfer-ence of the energy from randomly distributed structure scatters,who resolution is beyond image system capabilities.The speckle degrades the contrast resolution in ultrasound images,limiting the detectability of small,low-contrast lesions,thus making ultrasonic images generally difficult for the non-specialist to interpret [2].
Prior to the further image processing,the speckle needs to be removed.The goal of the enhancement is to attenuate the speckle without affecting important image features.As a multi-plicative noi,the speckle is not truly a noi becau its texture often carries uful information of the image.In some conditions,despeckling may be retroactive.Many methods have been devel-oped to reduce the speckle in ultrasonic images,including local statistics methods such as the Lee filter [3],the Frost filter [4]and the Kuan filter [5],the median filtering [6],the speckle reducing anisotropic diffusion (SRAD)[7],and the wavelet filtering [8].
The Lee filter [3]forms an output by a weighted average cal-culation using sub-region statistics over
different pixel windows.The Kuan filter [5]has the same form as the Lee filter but with a different weighting function,for it makes no approximation to
∗Corresponding author.Tel.:+862165642756;fax:+862165643526.
E-mail address:guoyi@ (Y.Guo),yywang@ ,yywang@ (Y.Wang),082021028@ (T.Hou).the original model.The Frost filter [6]is an exponentially damped convolution kernel which adapts to regions containing edges by exploiting local statistics.The locally bad filters compromi between the averaging (in homogeneous regions)and prerving (at edges and features).Such coefficients depend on local statistics in the moving windows.In ca of the low coefficients,the filter tends to be average.Given the high coefficients,it will rerve more sharp features.
The median filter [6]is a nonlinear operator,often ud in medi-cal image denoising.It replaces the middle pixel in the window with the median value of its neighbors.A single and unreprentative pixel in a neighborhood will not affect the median value signifi-cantly.The main problem is that the median filter would blur edges and tiny details.
The SRAD [7]is a non-linear and space-variant filter.It produces a family of resulting images bad o
n an iterative diffusion process.The diffusion coefficient rves as the edge detector,exhibiting high values at edges or features and low values in homogenous regions.With the diffusion coefficient and the number of iterative times,the SRAD achieves a balance between the despeckling and edge rerving.
The wavelet transform is widely ud in the image denoising.The discrete wavelet transform captures the non-stationary nature of the signal and noi,by parating the image into a t of scales.It applies soft thresholds to wavelet coefficients of different scales to eliminate the noi.Gupta et al.[8]propod a wavelet-bad statistical approach for speckle reduction in ultrasonic images.It is bad on the generalized Gaussian distributed modeling of sub-band coefficients,using the BayesShrink method.In the wavelet filters,the threshold plays an important role in the denoising.A
1746-8094/$–e front matter © 2010 Elvier Ltd. All rights rerved.doi:10.1016/j.bspc.2010.10.004
130Y.Guo et al./Biomedical Signal Processing and Control6 (2011) 129–138
small threshold will retain the noi whereas a large one leads to the loss of signal details.
Though the image denoising algorithms operates well under some situations,they have veral limitations.First,as the locally bad image denoising methods,different window sizes and shapes greatly affect thefiltering quality.A large window leads to over-smoothing,making edges blurred.Contrastively,a small window reduces the smoothing capability.Second,the algorithms require thresholds during thefiltering process.An inappropriate choice of a threshold may lead to averagefiltering,thus leaving sharp features unfiltered.Third,they ldom enhance edges and small structures.
A tradeoff always exists between the denoising and prerving.
Buades et al.propod a new denoising method,named Non-Local Means(NL-means)algorithms[9].It is a weighted Gaussian filter,taking advantage of the high degree of redundancy in any nature image.Compared with above denoising methods,the NL-means algorithm us the region comparison instead of pixel comparison and the pattern redundancy is not restricted to be local. Hence,it performs well in Gaussian noi reduction and sharp edges prervation.However,the speckle model in ultrasonic images is subject to the Rayleigh distribution,the NL-meansfilter cannot be applied to despeckling directly.
In this paper,we prent a modified motivation for the NL-meansfilter with maximum likelihood estimation as thefirst step. The main contribution of this modified NL-means(MNL)algorithm is the adaptation for ultrasonic images.The paper is organized as follows.Section2introduces the speckle model and a brief descrip-tion of the NL-means method.In Section3,we will explain our new method in detail.In Section4,parameters optimization and a com-parison of our method with other previously propodfilters will be prented.Finally,we will draw the conclusion and prent future rearch directions in Section5.
2.Speckle model and NL-means algorithm
2.1.Speckle model
The speckle is a coherent interference artifact who verity depends on the relative pha between two overlapping returning ultrasonic echoes.It produces small structures,masking true inter-faces.It is generally assumed that the speckle is fully developed in ultrasonic images.The speckle follows a Rayleigh density probabil-ity distribution.Then,the multiplicatively corrupted backscattered signal is modeled as[10],
g(i)=f(i)∗Á(i),i∈I(1) where g(·),f(·)andÁ(·)denote the obrved noisy image,noi-free image and fadin
g variable.I is a t of pixels in the image.Here Á(·)is a stationary random process independent of f(·),following
a
Fig.2.Similar windows in a nature image.
Rayleigh distribution,
p(Á, Á)=
Á
Á
·exp
−Á2
2 Á
(2) (Á)= Á
2
(3) (Á)2=
Á(4− )
2
(4) where Áis the shape parameter, (Á)is the mean and (Á)2is the variance.The real and imaginary parts ofÁare independent, orthogonal,identically distributed Gaussian random variable with zero means.Here,we restrict the mean ofÁto1.
With this model,a synthetic image with the speckle noi can be simulated.For an example,given Fig.1(a)as an original image, a“noisy image”shown as Fig.1(b)can be simulated by the original image multiplied with a Rayleigh distribution fading variable.
2.2.NL-means algorithm
For a certain small window p in an image,it has many similar windows in the same image.Sometimes,the“similar windows”are not clo to p,spreading anywhere in the whole image.An example is shown in Fig.2.According to the intensity and struc-ture,the window p is similar with the window q1and the window q3,but unlike with the window q2.
In the locally bad image denoising methods,we only u the“neighborhood window”forfii
n Fig.2,q1and q3 are far away from p in the location,thus they will be omitted in the denoising.It is obviously far from enough.If the conventional “neighborhood”is extended to the“whole image”,the window
p Fig.1.Speckle simulation(a)the original image,(b)the noisy image.
Y.Guo et al./Biomedical Signal Processing and Control 6 (2011) 129–138131
may be predicted by all similar windows in the whole image,which will result in the improvement of the denoising performance.That is the basic idea of the NL-means filter.
In the NL-means filter,for a given pixel,the restored gray value of each pixel is obtained by the weighted average of gray values of all pixels in the image.Each weight is proportional to the simi-larity between the local neighborhood of the pixel being procesd and the neighborhood corresponding to the other pixels [11].Given a discrete noisy image g ={g (i )|i ∈I },the filtered value NL(g (i ))is calculated as a weighted average of all pixels in the image,NL (g (i ))=
批评之后j ∈I
w (i,j )g (j )(5)
where i is the pixel being filtered,j reprents each one of pixels in
the noisy image,w (i ,j )is the similarity between the pixel i and j ,sat-isfying the conditions 0≤w (i ,j )≤1and
j ∈I
w (i,j )=1.The weight w (i ,j )depends on the similarity between the neighborhoods N i and N j of pixel i and j ,where N i denotes as a square neighborhood of size R sim *R sim ,centered at the pixel i .Pixels with a similar neigh-borhood will have larger weights.For an example,in Fig.2,w (p ,q 1)and w (p ,q 3)are higher than w (p ,q 2).The neighborhood similar-ity is measured as a decreasing function of the weighted Euclidian distance,
d (i,j )=G ||g (N i )−g (N j )||2
(6)
where G is a normalized Gaussian weighted function with zero
mean and standard deviation.Then,w (i ,j )is calculated as,w (i,j )= 1Z (i )
·exp
−d (i,j )h
(7)Z (i )=
j ∈I
exp
−d (i,j )h 2
(8)
Here Z (i )is the normalized constant.The parameter h acts as
a degree of filtering,controlling the decay of the exponential function.When i =j ,the lf similarity is high enough to yield over-weighting effect.In this ca,we define,w (i,i )=max(w (i,j ),∀i /=j )
(9)
Ideally,the NL-means filter may arch similar windows in the whole image.However,it will produce the high quality denoising,but computationally expensive and inefficient.Hence,we restrict this “arch area”to the window ˝with the size of R arch *R arch pixels to reduce the operation time.
3.The modified non local-bad speckle filter
The kernel of the NL-means algorithm is normalized Gaus-sian weighted function as shown in (6).So,
this filter is originally designed for Gaussian noi reduction,which is improper in remov-ing the Rayleigh distributed speckle.Therefore,modification should be made in the despeckling of ultrasonic images.
Our MNL method consists of two steps.First,the maximum likelihood (ML)estimation is employed to calculate the initial noi-free intensity.Then,we u the NL-means algorithm to restore details.
3.1.Maximum likelihood estimator
The maximum likelihood (ML)estimation is a robust filtering method to provide estimation for the statistical model’s parameters [12].Given the noi-free pixel f (i ),we define (i )is the neigh-borhood of pixel i .For a small window size,f (i )can be assumed
constant.The noisy pixel g (i )is a Rayleigh distribution with the parameter g , g = Á*f (i ).
p g (g (i ))=
g (i )
2g
·exp
−g 2(i )
2 2g
(10)
For a t of n statistically independent obrvations in i ,where n
is the number of pixels in i ,the likelihood function is given by,L (g (i 1),g (i 2),...,g (i n )| g )=
n k =1
p g (g (i k ))
(11)
In practice,it is always more convenient to work with the log-likelihood function,
ln(L (g (i 1),g (i 2),...,g (i n )| g ))=
n k =1
ln
g (i k )
2g
−
1
2 2g
n
k =1
g 2(i k )
(12)
Then,the ML estimator of g is calculated as,ˆ g =arg max[ln(L (g (i 1),g (i 2),...,g (i n )| g ))](13)
According to (10)–(13),ˆ g is given by,求婚歌
ˆ g =
12n
·
n
k =1
g 2(i k )
1/2
(14)
and f (i )is,
f (i )=ˆ
g ( Á)−1=
12n ( Á)
2
·
n
k =1
g 2(i k )
1/2
(15)
3.2.Spatially adaptive pixels
In the ML estimator,all pixels in i are ud to compute f (i ),
regardless of their similarity with the being filtered pixel i .It is quite ineffective.As shown in Fig.3,the intensity and neighborhood structure of the pixel k are contrary to tho of the pixel i .If the pixel k is ud to estimate the value of i ,the prediction would be wrong.
The mean intensity is the first order statistics,showing the char-acter of the window,which indicates the similarity between two pixels.Thus,we apply it as the criterion to discard tho irrelevant pixels in
a preliminary step.We define M k as a square neighbor-hood of pixel k .The mean intensity of M k is:¯M k =1/m m k =1
g (k ),where m is the number of pixels in M k .If ¯M
k is clo to ¯M i ,that is 1<¯M
红鼻头k /¯M i < 2,the pixel k is considered,where 1<1and 2>1are two constants clo to one.Otherwi,the pixel k will be eliminated.Therefore,n can be adaptive according to local spatial
contexts.
Fig.3.One blood cell in Fig.2.
132Y.Guo et al./Biomedical Signal Processing and Control6 (2011) 129–138 3.3.MNL specklefilter
The ML estimator is a kind of locally bad image denoising
method.It does not retainfine structure details and usually makes
edge blurred.Inspired by the NL-means algorithm,we combine
the two methods together.First,ML estimator with spatially
adaptive method is applied to obtain the initial noi-free pixels.
Then,we u NL-means algorithm tofilter and reconstruct edges
or textures.
MNL specklefilter scheme
(1)ML estimator:
For each pixel i in the noi image g,
(a)take the window i and M i;
(b)compare the average intensity in i to discard unwanted ones;
唐玄宗年号(c)compute the initial noi-free value f ML(i)using(15).
(2)NL-meansfilter:
拍摄方案
For each pixel i in the MLfiltered image f ML,
(a)take the arch window˝i and the neighborhood window N i;
(b)for each pixel j in the arch window,compute d(i,j),Z(i)and
w(i,j)using(6)–(8);
(c)u(5)to compute thefinalfiltered value NL(g(i)).
4.Experiments and results
We implemented the propod algorithm using Matlab
7.4(R2007a)and performed experiments on a PC with2.66GHz
Intel Core2processor.To evaluate the performance of the MNL,
we optimized three parameters of the MNL,tested it on synthetic
images and clinical ultrasonic images.For all test images,the per-
formance of the MNL is compared with tho of other six previously
propodfilters,including the conventional NL-meansfilter,the
ML estimator,the Leefilter,the Medianfilter,the SRAD and the
Med-waveletfilter.
4.1.Synthetic images
Fig.1(a)was utilized as an original image in the simulation.
Realistic spatially correlated speckle noi in ultrasonic images
can be simulated by low-passfiltering a complex Gaussian ran-
domfield and taking the magnitude of thefiltered output.We
performed the low-passfiltering by averaging the complex val-
ues in a3*3sliding window[13].The noisy images with different
SNRs were generated by the original image multiplied with differ-
ent Rayleigh-distributed noi levels.Here SNRs of the generated
noisy images are2,4,6,8,10,12,14,16dB respectively.SNR is
defined as,
SNR(dB)=10log10
N
i=1
g2(i)
N
i=1
(g(i)−f(i))2
(16)
where g(·)is the noisy image,f(·)is thefiltered image,and N is the number of pixels in g(·).The higher the SNR,the less noi is.
4.1.1.Parameter optimization
The MNL algorithm has three parameters:h(the decay of the exponential function),R sim(the radius of the similar neighborhood) and R arch(the radium of the arch window),controlling thefil-tering results.To measure the quality of thefilter with different parameters,the Root Mean Squared Error(RMSE)was
ud.Fig.4.The optimum h values test,(a)the relationship between the RMSE and h at a certain R sim,(b)optimum h at different SNRs and R sim.
Given the noisy image g(·)and thefiltered image f(·),the RMSE is,
RMSE=
N
i
(f(i)−g(i))2
N
(17)
The less RMSE is,the better the imagefiltering quality is.
First,we tested the optimum parameter h.Due to the fast decay of the exponential kernel in w(i,j),large Euclidean distances lead to nearly zero weights.h controls the decay of the weights,theref
ore, the degree of the smoothing.If h is too small,little noi will be filtered,or h is too large,the edge will be blurred.
We t afixed R arch to obrve the variation of R sim and h.At R arch=10,for each SNR of noisy images and R sim(R sim∈[2,6]),an exhaustive arch for optimum h was performed.Fig.4is one of the experiments with SNR=10dB,where 2is the variance of the speckle.In Fig.4(a),as h grew,the RMSE was from descending to ascending.We marked the turning point,corresponding to the least RMSE,as our optimum h under this condition.
Fig.4(b)is the optimum h at different SNRs.As obrved,there is an inver linear relation between the SNR and the optimum h. The higher the SNR,the lower the h.Experimentally,h can be taken values between9 2and14 2,obtaining a high visual quality.Here, we t h to10 2.
Y.Guo et al./Biomedical Signal Processing and Control6 (2011) 129–138
133
Fig.5.Filtered images at different R sim(SNR=8dB,R arch=10,h=10 2),(a)–(f)corresponding to R sim=2to R sim=6.
The cond parameter is R sim,the radius of the similar neighbor-hood N i.A small R sim is not robust to the noi,while a large one will lead to fewer neighborhoods to be found.We t the optimum h (h=10 2)and thefixed R arch(R arch=10)to test the performance of different values of R sim.From the references on image denoising [3–6],R sim usually varies from1to8.Thus,we cho R sim from2to 6in the experiment.
Fig.5illustrates thefiltered output at different R sim.At R sim=2or R sim=3,nois are still visible,especially in the background.When R sim increas,nois arefiltered,with the edge dimmed.We cho R sim=4as our optimum value.This window size is large enough to be robust to the noi and at the same time to be able to prerve fine details.
The third parameter is R arch,the radium of the arch window to reduce the computational complexity of the MNL specklefilter. We t the optimum h(h=10 2)and the optimum R sim(R sim=4) to experiment on different values of R arch.Since R sim=4,we t R arch in the range of[7,12].
In Fig.6,as R arch grows,there are slightly differences among the RMSE,but sharply increasing in t
he operation period.The run-time of R arch=12is nearly three times of that of R arch=7.A good compromi between the accuracy and computational load was found to be R arch=10.
From the above results over the synthetic data,we suggested the arch window˝i, i of21*21pixels,a similarity square neigh-borhood N i,M i of9*9pixels and h was manually tuned to be10 2. Spatially adaptive parameters were t as 1=0.8, 2=1.2.All the parameters gave satisfying denoising results.4.1.2.Methods comparison
For the visual and quantitative comparisons,five criteria,includ-ing the visual display,the signal-to-noi(SNR),the mean structure similarity(MSSIM)[14],thefigure of merit(FOM)[7]and the method noi[15]were taken into account.Each criterion reflects one aspect of the denoising method.The performance of the MNL is compared with tho of other sixfilters,including the conventional NL-meansfilter,the ML estimator,the Leefilter,the Medianfilter, the SRAD and the Med-waveletfilter.The window size of the Lee and the Median is both5*5,the iterative times of the SRAD is350.
The visual comparison will directly show how well thefiltered image is.Figs.7and8are results of venfilters at SNR=12dB and SNR=6dB respectively.Apparently,our propod method has the best despeckling effect,ideally sharp boundaries and no artifacts. The SRAD can also attenuate the speck
le in the background,how-ever,it blurs edges.The Lee,the Median,the Med-wavelet,the ML estimatorfilters hardly remove the noi.As to the conventional NL-means,it works well with the Gaussian noi,but fails in the Rayleigh noi.The speckle is clearly visible in the background.
The SNR is employed to estimate the effectiveness of the noi reduction.Table1shows the SNR of venfilters under different noi conditions.This table ems to corroborate the visual obr-vation mentioned before.The SRAD,the ML,the NL-means and our propod MNL denoi speckle well,all of them greatly increas-ing original SNR of noisy images.Whereas,at low SNRs(SNR=6dB, 4dB and2dB),the SNRs of the ML and the NL-means decrea sig-nificantly,hardly eliminating the speckle.But the MNL still holds high SNR values,surpassing other three methods.
Table1
The SNR comparison under different noi conditions.
Noisy image SNR=16dB SNR=14dB SNR=12dB SNR=10dB SNR=8dB SNR=6dB SNR=4dB SNR=2dB
MNL18.60918.53419.08219.81320.58120.88418.67015.124 Lee18.17417.26415.59514.20512.5131
0.6138.822 6.973 Median22.30220.61418.77216.82414.78812.86310.8478.911 Med wavelet18.40616.77515.00313.24711.4289.307 5.362 3.889 SRAD18.31818.22018.08017.83517.53417.06916.27215.642 ML16.24716.19415.90715.41514.63813.52111.87410.092 NL means18.90818.09917.18316.05914.98312.87310.8838.892
