Guided trilateral filter and its application to ultrasound image despeckling

Biomedical Signal Processing and Control 55 (2020) 101625

Contents lists available at ScienceDirect

Biomedical Signal Processing and Control journal homepage: www.elsevier.com/locate/bspc

Guided trilateral filter and its application to ultrasound image despeckling Wenchao Cui ∗ , Mengmeng Li, Guoqiang Gong, Ke Lu, Shuifa Sun, Fangmin Dong College of Computer and Information Technology, China Three Gorges University, Yichang 443002, China

a r t i c l e

i n f o

Article history: Received 19 January 2019 Received in revised form 17 June 2019 Accepted 21 July 2019 Keywords: Speckle reduction Ultrasound image Bilateral filter Guided filter Fisher-Tippett distribution

a b s t r a c t Speckle reduction has been a hot research topic in ultrasound imaging community. However, designing a method with better despeckling performance and lower algorithm complexity is still an active pursuit. In this paper, we propose a novel local spatial filtering framework named as the guided trilateral filter (GTF) to implement restoration of a noisy image with a well-fitting statistical distribution model. The GTF is derived as the local maximum likelihood estimation from the probabilities of the residuals between the filtered/guided image and the noisy image. The resulting iteration algorithm is essentially a local weighted filtering whose weights come from the guided image, including spatial distance, range discrepancy and statistical distribution. Choosing specific functions for these trilateral weights can direct to some classical local spatial filters, such as classical bilateral filter, robust bilateral filter, joint/cross bilateral filter, speckle reduction bilateral filter, etc. As a despeckling application of the GTF, we embed the Fisher-Tippett distribution model of ultrasound image into the GTF and thereby present the guided trilateral despeckling filter (GTDF). Experimental results on synthetic and real ultrasound images have demonstrated that the GTDF not only can achieve superior performance against state of the art methods, but also has fast and robust convergence and parameter setting insensitivity. © 2019 Elsevier Ltd. All rights reserved.

1. Introduction Ultrasound imaging has been one of the most widely adopted imaging modalities in medical imaging community owing to its predominant characteristics such as non-invasiveness, harmlessness, portability, instantaneity and low cost [1]. Unfortunately, the coherent nature of the acquisition system leads to a phenomenon that there exists a random granular pattern across the whole image. This pattern, called as speckle and modeled as multiplicative noise, more or less influences the interpretability of the imaging tissues [1,2]. Therefore, speckle in ultrasound imaging should be attenuated as much as possible. In the past decades, several algorithms have been proposed in the literature for speckle reduction [1–5], including local spatial filters [6–9], anisotropic diffusion (AD) filters [10–13], nonlocal means (NLM) filters [14–17], total variation (TV) methods [18–22], multiscale methods [23–26] and homomorphic approaches [27], etc. Local spatial filters are usually designed based on local statistics. For example, the adaptive weighted median filter [6] tailored clas-

∗ Corresponding author. E-mail address: [email protected] (W. Cui). https://doi.org/10.1016/j.bspc.2019.101625 1746-8094/© 2019 Elsevier Ltd. All rights reserved.

sical median filter through the introduction of weight coefficients according to the local statistics around each pixel of the ultrasound image. Guo et al. [7] used local intensity texture information to determine the nonhomogeneous pixels and then adopted the directional average filters to reduce speckle. To be different, Thakur et al. [8] used gray level co-occurrence matrix to compute four texture features and then determined the appropriate mean filter to smooth the nonhomogeneous regions. An improved algorithm known as speckle reduction bilateral filter (SRBF) [9] redesigned the range weight function of the classical bilateral filter (BF) [28] based on local Rayleigh distribution, in order to resist the interference of speckle in ultrasound imaging. Being easy to implement and control, local spatial filters are broadly used for despeckling task. However, the inherent speckle usually influences the reliability of the local image statistics used for these filters and hence leads to unsatisfactory restoration results. AD filters mostly improved the P-M model [29] by reformulating a partial differential equation (PDE), aiming to achieve a balance between details preservation and speckle suppressing. A well-known method named as speckle reduction AD (SRAD) filter [10] exploited the instantaneous coefficient of variation to implement the edge-sensitive diffusion for speckled images. AjaFernández and Alberola-López [11] discussed the estimation of the coefficient of variation and proposed details preserving AD (DPAD)

2

W. Cui, M. Li, G. Gong et al. / Biomedical Signal Processing and Control 55 (2020) 101625

filter. Krissian et al. [12] extended SRAD by allowing different levels of filtering along the image contours and their principal curvature directions. Recently, Ramos-Llordén et al. [13] proposed an AD filter with a probabilistic-driven memory mechanism to overcome over-filtering problem by following a tissue selective philosophy. To get the numerical solutions of the PDEs, AD filters are likely to eliminate meaningful structural details during iterations and tend to produce staircase effect in the filtered results. NLM filters compute the weighted average for each pixel not in its neighborhood but throughout the whole image. Meanwhile, the corresponding weights are calculated based on blockwise similarity, not pointwise similarity. For instance, Coupe et al. [14] presented an optimized Bayesian NLM filter for speckle reduction using the Pearson distance for block comparison. Guo et al. [15] used maximum likelihood (ML) estimation to obtain an initial noise-free image, and then applied NLM filter to refine it. Zhan et al. [16] updated the weights iteratively in a lower dimensional subspace using principal components analysis to improve accuracy of the weights and reduce its computational complexity. An unbiased NLM filter [17] further enhanced filtering performance by using a three-parameter Gamma distribution function to fit the ultrasound images. Since NLM filters usually involve excessive computation in terms of measurements on blockwise similarity, they require massive computational burden that may obstruct their real-time running. TV methods formulate the despeckling problem as minimization of an energy function, which usually includes a data fitting term and a TV regularization term. For example, the AA model proposed by Aubert and Aujol [18] assumed the multiplicative speckle to obey Gamma distribution, and then used maximum a posteriori (MAP) estimation to derive a minimization function. Subsequently, two different convex optimization extensions to the AA model had been proposed to fast solve the TV regularization function, based on an inverse space scale algorithm [19] and a variable splitting method [20], respectively. Converting the multiplicative model into an additive one by taking logarithm, a TV method called as multiplicative image denoising by augmented Lagrangian (MIDAL) [21] had also presented a MAP based minimization function, which could be solved using the augmented Lagrangian framework. Huang and Yang [22] used a generalized Kullback-Leibler distance to build an objective function and its solution was reached using variable splitting and the Bregman iterative method. TV methods control the balance between speckle smoothing and edge preserving by the key coefficient of the TV regularization term, such that an unsuitable setting value will lead to either over-smoothing or under-smoothing. Multiscale methods firstly adopt wavelet transformation to decompose the original image or its logarithmic form into subband coefficients, which are then processed by means of some well-designed scheme that can eliminate the ingredient of speckle. Finally, the filtered image is synthesized from the processed wavelet coefficients by inverse wavelet transformation. In [23], Achim et al. used the heavy-tailed alpha-stable distribution to model subband decompositions, and then designed a Bayesian processor to remove speckle. Based on Daubechies complex wavelet transformation, Khare et al. [24] presented a despeckling method that detected strong edges using the imaginary components of complex scaling coefficients and then applied shrinkage on the magnitudes of complex wavelet coefficients at non-edge points. Zhang et al. [25] used a fast bilateral filter (FBF) to process the lowpass coefficients while a Bayesian shrinkage algorithm trimmed the others. This method will be called as wavelet shrinkage and FBF (WSFBF) algorithm in the rest of this paper. Subsequently, a new trial that combined the wavelet shrinkage and the guided filter [30] could be found in [26]. Owing to the multiresolution and directionality characteristics of wavelet transformation, multiscale

methods can effectively preserve texture details, though they tend to produce ringing artifacts. Homomorphic approaches usually convert the multiplicative speckle into the additive noise by logarithmic transformation such that an additive noise removal algorithm can be employed. Recently, two additive Gaussian denoisers, denominated as dualdomain image denoising (DDID) [31] and block matching and 3-D filtering (BM3D) [32], were used in a homomorphism fashion for speckle reduction [27], hereinafter named as HomoDDID and HomoBM3D, respectively. HomoDDID was a hybrid method based on the joint bilateral filtering and the shrinkage of Fourier coefficients. HomoBM3D searched the similar image patches and then performed the three-dimensional wavelet shrinkage. Homomorphic approaches need to perform logarithmic transformation of the original image, such a nonlinear operation may lead to the biased characteristics computed from the image and thereby affect the despeckling performance to some extent. Although the above-mentioned methods are all designed to reduce speckle and enhance the visual effect, the optimum balance between smoothing speckle and preserving image details remains an intractable problem. Additionally, it may be a difficult problem how to reach the best compromise between despeckling performance and algorithm complexity, especially when processing ultrasound 3D data or videos in the real-time environment. Therefore, researchers devote their time to find a better method with lower algorithm complexity to deal with speckle for ultrasound images. To this end, we propose a novel local spatial filtering framework that is derived as the local ML estimation from the probabilities of the residuals between the filtered image and the noisy image. The resulting formulation contains three aspects of information, that is, spatial distance, range discrepancy and statistical distribution. The corresponding weights are computed from the filtered image like the guided filter [30]. Hence, we call it as the guided trilateral filter (GTF). It will be explained in the following sections that classical BF, robust BF, joint/cross BF [33,34] and SRBF can be regarded as special cases of the GTF. For its application to ultrasound image despeckling, the Fisher-Tippett distribution is introduced into the GTF such that the guided trilateral despeckling filter (GTDF) is generated. We give extensive experimental results on synthetic and real ultrasound images to demonstrate the superior performance of GTDF by comparing it with several despeckling methods. The rest of this paper is organized as follows: In Section 2, we introduce the proposed GTF framework and its four special cases of filters. Based on the GTF framework and the statistical distribution of ultrasound image, we give a detailed description of the proposed GTDF method in Section 3. Section 4 introduces the quantitative indices of despeckling performance estimation. We report the experimental results on synthetic and real ultrasound images, comparing the GTDF method with some state-of-the-art methods in Section 5. Section 6 concludes this paper.

2. Guided trilateral filter (GTF) framework As shown in Fig. 1, let  be the whole image domain, Nr (x) be the local spatial neighborhood centered at pixel location x (x ∈ ) with size [-r, r]×[-r, r]. I and ˆI denote a noisy image and its expected filtered image, respectively. Then the probability of the residual, bt = ˆI (x)— I(x+t), between ˆI (x) and I(x+t) for each neighborhood pixel t [t ∈ Nr (x)], can be described as [35]:

 2  ⎞qt bt s 1 ⎠ p(bt ) ∝ ⎝ e ⎛

−

s

(1)

W. Cui, M. Li, G. Gong et al. / Biomedical Signal Processing and Control 55 (2020) 101625

3

of the overall residuals b in Nr (x) can be written as

p (b) ∝



p (bt ) =

t ∈ Nr (x)



⎛  ⎜1 ⎝se

=



−

 2  ⎞qt bt s ⎠ ⎝1e ⎛

t ∈ Nr (x) ˆI (x)−I(x+t) s

−

s

2  ⎞gt t ⎟ ⎠

(5)

t ∈ Nr (x)

which can be also called as the likelihood function with a unknown variable ˆI (x). Following the maximum likelihood (ML) estimation for ˆI (x), one can perform a minus natural logarithm operation on (5), resulting in the following minimization function





F=

gt t 

t ∈ Nr (x)

Fig. 1. An illustration for the different neighborhoods in the GTF framework. The asterisk marker represents the location of the pixel x to be processed, whose filtered value will be generated by its neighborhood pixels in Nr (x) (the square [-r, r]×[-r, r] centered at x with the solid line). Therein two neighborhood pixels, x+t and x+t , t, t ∈ Nr (x), are indicated by the dot markers. Their neighborhood regions Nr (x+t) and Nr (x+t ) are drawn by squares with the dotted line and the dashed line, respectively. Note that the PDF of the distribution weight t in Eq. (3) should be estimated by fitting the pixel values of the neighborhood Nr (x+t) with some distribution function.

where s is the scale parameter, ϕ (·) is a differentiable kernel function with monotonically increasing property, qt is an indicator of the similarity of the noiseless intensities at pixels location x and x+t [35]. In this paper, the similarity indicator qt records both the spatial similarity and the distribution similarity. The former is measured by the usual Gaussian weight defined as



gt = exp

−

||t||2

ˆI (x) − I(x + t) s

−



gt t ln

1

t ∈ Nr (x)

s

(6)

Obviously, the optimum estimation ˆI (x) can be obtained by setting the derivative of F with respect to ˆI (x) to zero, that is





gt t 



t ∈ Nr (x)

ˆI (x) − I(x + t) s

2 





ˆI (x) − I(x + t) = 0

(7)

where ϕ (·) is the differential of the kernel function ϕ (·). Finally, by solving the Eq. (7), we can derive the optimum estimation ˆI (x), which is also the filtered result of our proposed guided trilateral filter (GTF) framework, as the following iteration form:





(2)

2s2

2 

gt t 



t ∈ Nr (x)

ˆI k (x)−I(x+t) s



2 

I(x + t)

2 

where  s is the standard deviation. It can be observed that the nearer the spatial distance is, the bigger the similarity is. The latter is measured by the probability that ˆI (x) belongs to the distribution that generated from the local neighborhood Nr (x+t) centered at x+t, as shown in Fig. 1. The corresponding weight is computed as

ˆI k+1 (x) = GTF

P ˆI (x) ∈ t (y t ) = t (ˆI (x) t )

where k is the iteration number, ˆI 0 (x) = I(x) at the initialization (k = 0). Actually, Eq. (8) represents a local spatial filtering process which includes three aspects of weights: spatial distance weight gt , statistical distribution weight t and range discrepancy weight ϕ . Inspired by the bilateral filter (BF) having two weights, we call it the trilateral filter. Additionally, the weights t and ϕ need to be computed through the intermediate filtered image. This behavior is similar to the guided filter [30] whose weights are calculated not from the original noisy image but from the guided image. When considering the intermediate filtered image as the guided image, the proposed framework (8) can be named as the guided trilateral filter (GTF). Remark. It should be clarified that our proposed GTF is different from the guided bilateral filter (GBF) proposed in [35]. The GBF has also three weights: one is the spatial weight and the other two are range discrepancy weights that computed from the guided image and the noisy image, respectively. However, our proposed GTF introduces a distinct statistical distribution weight, which has been verified in [9] that it is suited to perform filtering task on the premise of a known distribution model concerning the noisy image.





(3)

where t (y| t ) is the probability density function (PDF) fitted with pixel values in the neighborhood Nr (x+t) of the filtered image ˆI .  t denotes the parameter set of the PDF. When ˆI (x) and ˆI (x+t) belong to a similar distribution, the weight t yields high values and promotes smoothing. On the contrary, t has a small value when a contour is present, thus preserving the edges. Obviously, the distribution weight t has flexible alternatives depending on the known appropriate statistical distribution models with respect to different imaging modes. It is noted that the same idea can be found in the SRBF algorithm [9] in which the PDF matched with the Rayleigh distribution. Therefore, the similarity indicator qt in Eq. (1) is defined as qt = (gt )˛ (t )ˇ

(4)

where the parameters ˛ and ˇ are used to adjust the relative importance of the spatial similarity and the distribution similarity. For simplicity, we set ˛ =ˇ = 1 in this paper. Assuming the residuals bt in the whole neighborhood Nr (x) are independent and identically distributed, then the joint probability



gt t 

ˆI k (x)−I(x+t) s

(8)

t ∈ Nr (x)

4

W. Cui, M. Li, G. Gong et al. / Biomedical Signal Processing and Control 55 (2020) 101625

2.1. Special case 1: robust bilateral filter

where  t is the scale parameter of the Rayleigh distribution and it can be computed using the ML estimation (see Appendix A):

If setting t = 1and ϕ (u)=exp(—u/2), then Eq. (8) can be simplified as

  k

 ˆI k+1 (x) RBF

=

ˆI (x)−I(x+t)

gt exp −

t ∈ Nr (x)

2 

I(x + t)

2s2

 



gt exp −

2 

(9)





1



2

Nr (x + t)

ˆ t =

I 2 (j)

(14)

j ∈ Nr (x+t)

where || Nr (x + t) || represents the total number of the pixels in neighborhood Nr (x + t).

ˆI k (x)−I(x+t)

3. Guided trilateral despeckling filter (GTDF)

2s2

t ∈ Nr (x)

This iterative algorithm was named in [35] as the robust BF. 2.2. Special case 2: classical bilateral filter When doing only one iteration (k = 0), Eq. (9) will become a local weighted average same to the classical BF. Therefore, substituting ˆI 0 (x) = I(x) into Eq. (9) can derive the formulation of the classical BF as





2

gt exp − (I(x)−I(x+t)) 2



2s

ˆIBF (x) =

t ∈ Nr (x)



I(x + t)



2

gt exp − (I(x)−I(x+t)) 2



(10)

2s

t ∈ Nr (x)

2.3. Special case 3: joint/cross bilateral filter In some situations, the information about the structure of the filtered image is available as the form of a similar image, usually called as the guided image [30]. Computing the weights from the guided image is equivalent to exploring effectively a prior knowledge on the specific processing task [33,34]. On the other hand, if the weights are computed directly from the noisy image, then large estimation bias will arise due to the seriously corrupted correlations of pixels under strong noise [36]. Therefore, some researchers [31,37,38] attempted to alleviate the estimation bias by calculating the weights from the intermediate filtered image. Such a variant of the classical BF is called as the joint or cross BF with the expression:



 ˆIJBF (x) =

gt exp

(x)−IG (x+t)) − (IG 2s2

t ∈ Nr (x)





gt exp − (IG

2



I(x + t)

(x)−IG (x+t))2 2s2



(11)

t ∈ Nr (x)

where IG represents the guided image or the intermediate filtered image.

In this section, we will give a specific application of the GTF framework on speckle reduction for ultrasound image. We name this instantiation as the guided trilateral despeckling filter (GTDF). As previously mentioned, the weight t is related to the statistical distribution of the noisy image to be processed. Therefore, we firstly discuss the distribution model for ultrasound image. 3.1. Statistical distribution for ultrasound image There are several PDFs to model the intensities of a ultrasound image, such as Rayleigh distribution [39,40], Lognormal distribution [41], Gamma distribution [42,43], Weibull distribution [43] and Fisher-Tippett distribution [1,44,45], etc. In this paper, starting from the Rayleigh distribution of the ultrasound envelope, we derive the log-compression ultrasound image to satisfy with the Fisher-Tippett distribution. It is well-known that for the fully developed speckle, the amplitude of the RF ultrasound signal envelope in the local neighborhood Nr (x + t) may be well fitted by a Rayleigh random variable Y such that p(Y ) =

Y 2t



exp

−

where t is the tissue reflectivity for the specific resolution cell [45]. In order to adjust the large dynamic range of the RF image to the normal range of the display equipment, one may use the logarithm transformation of Y and define a new random variable Z for the ultrasound image, which can be expressed as Z = log(Y + 1)

(16)

Note that by adding one in Eq. (16), we conveniently keep the variable Z in the same support interval [0, +∞) with the variable Y since there should be no negative intensities in the log-transformed image [45]. Thus, the distribution of the variable Z in the local neighborhood Nr (x + t) has a double exponential or Fisher-Tippett shape, whose PDF can be given as

  Z 2 

e −1

dY eZ − 1 p(Z) = p(Y ) exp Z −

= 2 2 dZ



(12)

where the statistical distribution weight t becomes as



 I(x) I 2 (x)

t I(x)  = 2 exp − 2 t 

2t

(17)

ˆ eI (x)

−1

2t



⎜ exp ⎝ˆI (x) −

2 ⎞

ˆ

eI (x) − 1 22t

⎟ ⎠

(18)

where the parameter t can be estimated in the local neighborhood Nr (x + t) based on the ML method. One can derive the corresponding estimator as (see Appendix B):

t ∈ Nr (x)

t

⎛

gt t

2t

Considering the fitting PDF in Eq. (3) as the Fisher-Tippett distribution (17), the statistical distribution weight t in the GTF framework (8) can be written as t (ˆI (x) t ) =

gt t I(x + t)

t ∈ Nr (x)

(15)

t

If setting ϕ (u) = 1 and defining t with the Rayleigh distribution, then the GTF in Eq. (8) at only one iteration (k = 0, ˆI 0 (x) = I(x)) reduces to SRBF:

ˆISRBF (x) =

Y ≥0

,

22t

2.4. Special case 4: speckle reduction bilateral filter (SRBF)





Y2

(13)

 ˆ 2t =

1



2

Nr (x + t)

  j ∈ Nr (x+t)

ˆ

2

eI (j) − 1

(19)

W. Cui, M. Li, G. Gong et al. / Biomedical Signal Processing and Control 55 (2020) 101625

5

Fig. 2. Comparison of fitting (b) the histogram of a homogeneous region (the white rectangular box) manually selected from (a) a B-mode CCA ultrasound image with the Fisher-Tippett distribution (the red solid line) and the Rayleigh distribution (the blue dashed line).

To demonstrate the validity of the Fisher-Tippett distribution model, a histogram fitting experiment is conducted using a Bmode common carotid artery (CCA) ultrasound image. As shown in Fig. 2(a), the white rectangular box represents the selected homogeneous region whose histogram is drawn in Fig. 2 (b). The Rayleigh distribution is joined in the comparison as a widely used fitting model for ultrasound image. It can be viewed that the PDF of the Fisher-Tippett distribution traces the profile of the histogram better than that of the Rayleigh distribution, especially at the peak and the right slope of the fitting curve.

Besides the statistical distribution weight t defined in Eq. (18), we choose the spatial distance weight gt same to Eq. (2) and the range discrepancy weight ϕ (u)=exp(—u/2). In this case, the GTF in Eq. (8) will evolve into the formulation of the GTDF as follows

 ˆI k+1 (x) GTDF

=

gt t exp −

t ∈ Nr (x)



2 

ˆI (x)−I(x+t)

I(x + t)

2s2

  gt t exp −

2 

(20)

ˆI k (x)−I(x+t) 2s2

t ∈ Nr (x)

where the scale parameter s is set in proportion to the parameter t of the distribution weight t (18), with the simple form: s2 =  ˆ 2t

4. Quantitative evaluation of despeckling performance To evaluate the performance of a despeckling algorithm quantitatively, two types of indices: full-reference indices and non-

reference indices are usually adopted depending on whether the standard despeckled image exists.

3.2. Implementation of the GTDF method

  k

where is the proportional coefficient. For clarity, the block diagram of the proposed GTDF method is drawn in Fig. 3 and the detailed implementation procedures are given in Algorithm 1.

(21)

4.1. Full-reference indices In an experiment on a synthetic image, the filtered image can be evaluated based on the available noise-free image. In this paper, we select three full-reference indices, namely peak signal-to-noise ratio (PSNR), feature similarity index (FSIM) [46] and edge preservation index (EPI) [45,47], to investigate the objective performance of the methods. All three indices reflect better despeckling performance with a higher value. Let X, Xˆ be the noise-free image and the filtered image, respectively, then PSNR is computed as follows:

ˆ = 10 log PSNR(X, X)

2552 MSE

 (22)

6

W. Cui, M. Li, G. Gong et al. / Biomedical Signal Processing and Control 55 (2020) 101625

Fig. 3. The block diagram of the proposed GTDF method.

where log(·) denotes the common logarithm, MSE represents the mean squared error and is defined as

2 1  MSE = Xi,j − Xˆ i,j MN M

EPI is a measure of the edge preservation capability of the despeckling algorithm, which is defined as [45]



N

(23)

ˆ = EPI(X, X)

i=1 j=1

where M×N is the image size. In view of the fact that human visual system understands an image mainly according to its low-level features, FSIM measures the structure similarity by two low-level features: phase congruency (PC) and gradient magnitude (GM). Let {PCi (x), Gi (x)|i = 1, 2, x ∈ } records all the computed PCs and GMs of the noise-free image and the filtered image, respectively, then the similarity measure for PC1 (x) and PC2 (x) is defined as SPC (x) =

2PC1 (x)PC2 (x) + T1 PC12 (x) + PC22 (x) + T1

(24)

where T1 is a positive constant to increase the stability of SPC , depending on the dynamic range of PC values [46]. Likewise, the similarity measure for G1 (x) and G2 (x) is defined as SG (x) =

2G1 (x)G2 (x) + T2 G12 (x) + G22 (x) + T2

ˆ = FSIM(X, X)

x∈˝

(26) PCm (x)

x∈˝

where PCm (x)=max[PC1 (x), PC2 (x)].

x∈˝

DX2 (x)



(27) D2ˆ (x) X

x∈˝

where D(·) = (·) − E [ (·)] and herein  means the high-pass filtered version of the image that can be obtained with a 3 × 3 pixel standard approximation of the Laplacian operator, E(·) indicates the mean operator. 4.2. Non-reference indices For an experiment on a real ultrasound image, the objective assessment can only be done between the original ultrasound image (noisy image) and the filtered result, due to the lack of the ground truth. Let In denote the noisy image, If denote the filtered image and VAR(·) indicate the variance operator. The speckle strength in ultrasound image can be defined as speckle index (SI) [48], which is the ratio of the standard deviation and mean in a homogeneous region: SI =

VAR(In ) E(In )

(28)

Meanwhile for a specific homogeneous region, the speckle reduction in the filtered image can be quantified by inherent signalto-noise ratio (ISNR) [49]:



ISNR =

SPC (x)SG (x)PCm (x)







(25)

where T2 is a positive constant depending on the dynamic range of GM values. As the overall similarity, FSIM is defined as



DX (x)DXˆ (x)

x∈˝

2

E(If )

VAR(If )

(29)

It is clear that the higher the ISNR value is, the better the effect of improving the signal-to-speckle ratio over homogeneous areas will be. In addition, speckle suppression index (SSI), which normalizes SI of the filtered image by that of the original noisy image, has been

W. Cui, M. Li, G. Gong et al. / Biomedical Signal Processing and Control 55 (2020) 101625 Table 1 The key parameters for the six rival algorithms.

Table 2 Performance comparison for seven algorithms on the synthetic ultrasound image with moderate speckle.

Algorithm

Experiment parameters

SRAD DPAD SRBF

Window size 7 × 7, time step t =0.2,iterations k = 100 Window size 7 × 7, time step t =0.2,iterations k = 100 Window size 9 × 9, standard deviation  c = 1.9 of the spatial kernel Regularization parameter = 4.0, penalty parameter = 4.0, iterations k = 20 Window size 15 × 15, standard deviation  s = 7 of the spatial kernel, bilateral range parameters r ={100, 8.7, 0.7}, wavelet shrinkage parameters f ={4.0, 4.0, 0.7} Window size 9 × 9, wavelet=’db8’, decomposition level J = 3

MIDAL HomoDDID

WSFBF

widely used to evaluate the despeckling performance [45,49,50]. SSI has the expression as



SSI =

VAR(If ) E(If )

·



E(In )

(30)

VAR(In )

Obviously, reducing the speckle will lead to SSI < 1 and the smaller the SSI value is, the stronger despeckling ability the filtering algorithm has. To make the despeckling performance be in direct proportion to the SSI value, we define a variant of SSI as SSI-dB = −20 log(SSI)

(31)

Algorithm

PSNR

FSIM

EPI

Noisy image SRAD DPAD SRBF MIDAL HomoDDID WSFBF GTDF

24.0077 27.7493 27.7478 28.2504 27.9198 28.5653 28.0862 28.6386

0.7873 0.8770 0.8776 0.8743 0.7998 0.8424 0.8804 0.8814

0.3942 0.5493 0.5508 0.4397 0.5328 0.5008 0.3974 0.5533

The bold values indicate the best value in the corresponding index.

Table 3 Performance comparison for seven algorithms on the synthetic ultrasound image with heavy speckle. Algorithm

PSNR

FSIM

EPI

Noisy image SRAD DPAD SRBF MIDAL HomoDDID WSFBF GTDF

18.4672 23.9896 24.1092 25.4501 25.3917 24.8993 23.7720 25.7349

0.6253 0.7960 0.7980 0.7938 0.7469 0.8235 0.8184 0.8300

0.2119 0.3472 0.3492 0.2804 0.3524 0.3933 0.2356 0.3535

The bold values indicate the best value in the corresponding index.

which will give a conclusion that a higher value represents better despeckling performance. If the filtering algorithm overestimates the mean of the filtered image, then SSI may fail to value the despeckling performance. To avoid this adverse effect, speckle suppression and mean preservation index (SMPI) has been adopted and its definition [49] is



7



SMPI = 1 + E(In ) − E(If )





VAR(If )

(32)

VAR(In )

Note that like SSI, a smaller SMPI represents better filtering performance. Similar to SSI-dB, the variant of SMPI, SMPI-dB = −20 log(SMPI)

(33)

will be used to indicate better despeckling performance with a higher value. Therefore in this paper, the despeckling performance for real ultrasound images will be contrasted through ISNR, SSI-dB and SMPI-dB. These indices will be calculated based on the homogeneous regions manually selected from the original noisy image, as shown later in Figs. 6(a) and 8(a). 5. Experimental results In this section, we will evaluate the performance of the proposed GTDF against several despeckling algorithms, including SRAD [10], DPAD [11], SRBF [9], MIDAL [21], HomoDDID [27] and WSFBF [25]. All the algorithms will be comprehensively assessed through the aforementioned indices based on synthetic and real ultrasound images. All the experiments have been conducted on Matlab 2014a by a HP Z440 workstation with a 3.5 GHz Xeon E5-1620 CPU and a 16 G RAM. The corresponding parameters of the rival algorithms have been determined following the setting rules recommended by the authors or set by trial and error. Table 1 gives out the key parameters for the six algorithms. 5.1. Experiments for synthetic images According to [23,47,51], ultrasound image can be simulated by a reference noise-free image multiplying with a unit mean random

field representing the variations in speckle amplitude. The reference image used in this paper is a kidney image downloaded from http://field-ii.dk, as shown in Figs. 4(a) and 5 (a). The spatially correlated speckle can be generated by lowpass filtering a complex Gaussian random field and taking the magnitude of the filtered output [47]. By changing the variance of the underlying complex Gaussian random field, one can generate images with different levels of speckle [51]. In our experiments, we perform the lowpass filtering by averaging the complex values in a 3 × 3 sliding window. It has been verified that such an operation is sufficient to model the real ultrasound image well [47]. By setting the variance of the underlying complex Gaussian random field as 1.0 and 4.0, we generate two synthetic kidney images with moderate and heavy speckle, as shown in Figs. 4(b) and 5 (b), respectively. For the proposed GTDF algorithm, the corresponding four parameters: the neighborhood radius r, the standard deviation  s of the spatial weight, the total iteration number k and the proportional coefficient are set as follows: r = 7,  s = 1.3, k = 2 and = 0.5 for Fig. 4(b); r = 9,  s = 2.5, k = 2 and = 1.5 for Fig. 5(b). Fig. 4 shows the filtered results obtained by seven algorithms in the case of moderate corruption. It can be seen that SRAD and DPAD remain little speckle while MIDAL fades some details due to over smoothing. This can be verified from Table 2 that the corresponding PSNR values are all less than 28.0 dB. WSFBF gets the second best FSIM that means this algorithm preserves the visual features well, whereas the other two indices are unsatisfactory, especially that the EPI is only slightly more than the noisy image. SRBF has comparable despeckling performance on PSNR and FSIM against HomoDDID, whereas its lower EPI reveals poor edge preservation property. As a contrast, the proposed GTDF achieves the highest PSNR and the best performance on preserving the visual feature (FSIM) and edge feature (EPI). In the presence of strong speckle, SRAD, DPAD and WSFBF reveal considerable speckle that can be observed from Figs. 5(c), 5(d) and 5(h), respectively. Although HomoDDID almost entirely remove the speckle, its filtered result brings out serious wavy artifacts, as shown in Fig. 5(g). These adverse effects cause their corresponding PSNRs in Table 3 to be less than 25.0 dB. However, HomoDDID is still an acceptable algorithm due to the second best FSIM and the best

8

W. Cui, M. Li, G. Gong et al. / Biomedical Signal Processing and Control 55 (2020) 101625

Fig. 4. Comparison of despeckling results for seven algorithms namely (c) SRAD, (d) DPAD, (e) SRBF, (f) MIDAL, (g) HomoDDID, (h) WSFBF and (i) GTDF on (b) a synthetic kidney ultrasound image with moderate speckle contaminating (a) the original image.

EPI. Comparing with MIDAL, SRBF has subtle advantages on PSNR and FSIM while it has worse edge preservation property. Finally, it is reasonable to conclude from Table 3 that the GTDF method has the best despeckling performance by virtue of the highest PSNR, the best FSIM and the second best EPI.

5.2. Experiments for real ultrasound images An open ultrasound RF data set provided by the Ultrasonic Imaging Laboratory in University of Illinois at Urbana-Champaign can be downloaded from http://ultrasonics.bioengineering.illinois. edu/data phantom.asp. From the RF data, the B-mode breast tumor ultrasound images can be generated. We select a benign tumor image shown in Fig. 6(a) as the first experiment of real ultrasound image. The corresponding parameters of the proposed GTDF algorithm are set as follows: r = 9,  s = 2, k = 2 and =15. Observing the filtered results in Fig. 6, one can find that SRAD, DPAD and MIDAL fail to suppress the influence of the speckle effectively. Besides, HomoDDID degrades the readability of the image due to the advent of the wavy artifacts. Their poor results can also be

Table 4 Performance comparison for seven algorithms on the tumor ultrasound image in Fig. 6(a). Algorithm

ISNR

SSI-dB

SMPI-dB

SRAD DPAD SRBF MIDAL HomoDDID WSFBF GTDF

126.6587 129.3145 191.0469 160.4573 173.6574 196.0820 204.7093

3.9760 4.0657 5.7605 5.0022 5.3457 5.8742 6.0606

2.2894 2.3256 3.9678 0.5810 −2.5551 2.9602 4.1634

The bold values indicate the best value in the corresponding index.

reflected from the quantitative estimation values listed in Table 4. To be specific, the ISNRs are less than 174.0, the SSI-dBs are less than 5.4 and the SMPI-dBs are less than 2.4. Furthermore, it can be judged from the recorded results that SRBF and WSFBF demonstrate comparable despeckling capabilities, though they are both inferior to the GTDF method. Further experiments have been done with the first 10 frames of images generated from the tumor ultrasound RF data set. The cor-

W. Cui, M. Li, G. Gong et al. / Biomedical Signal Processing and Control 55 (2020) 101625

9

Fig. 5. Comparison of despeckling results for seven algorithms namely (c) SRAD, (d) DPAD, (e) SRBF, (f) MIDAL, (g) HomoDDID, (h) WSFBF and (i) GTDF on (b) a synthetic kidney ultrasound image with heavy speckle contaminating (a) the original image.

responding evaluation indices with regard to each image for seven algorithms are shown in Fig. 7. It is clear to be observed from Figs. 7(a) and 7(b) that the ISNRs and the SSI-dBs of the GTDF method have dominant advantages among the seven algorithms except that MIDAL has subtle superiority on the Frame008 image. Besides, the SMPI-dBs shown in Fig. 7(c) also confirm that the GTDF algorithm has the best despeckling resuts for these tumor ultrasound images. Another open ultrasound image database is provided by the Signal Processing Laboratory in Brno University of Technology. This database contains 84 B-mode ultrasound images of common carotid artery (CCA) in longitudinal section, which can be downloaded from http://splab.cz/en/download/databaze/

ultrasound. The selected CCA ultrasound image is displayed in Fig. 8(a). The settings of the corresponding parameters are same to those for the benign tumor image in Fig. 6(a). After conducting the experiments on the real CCA image, we record the despeckled images in Fig. 8 and the corresponding objective indices in Table 5. It can be found that SRAD, DPAD, MIDAL and HomoDDID get poor ranking with the scope of ISNRs less than 3.0, SSI-dBs less than 0.7 and SMPI-dBs less than 0.1. Besides, SRBF and WSFBF can be classified into the second rank depending on their ISNR values, both of which lie in the interval [3.0, 4.0]. Finally, the GTDF algorithm possesses the dominant advantage in contrast to the other six despeckling algorithms.

10

W. Cui, M. Li, G. Gong et al. / Biomedical Signal Processing and Control 55 (2020) 101625

Fig. 6. Comparison of despeckling results for seven algorithms namely (b) SRAD, (c) DPAD, (d) SRBF, (e) MIDAL, (f) HomoDDID, (g) WSFBF and (h) GTDF on (a) a real B-mode tumor ultrasound image. The white boxes in (a) indicate the manually selected homogeneous regions in which the quantitative indices will be computed.

Fig. 7. Quantitative evaluation on (a) ISNR, (b) SSI-dB and (c) SMPI-dB for seven algorithms with 10 real B-mode tumor ultrasound images.

We also choose 10 images from the CCA ultrasound image database to verify the despeckling abilities of the seven algorithms further. The x-coordinates of Fig. 9 are marked by the corresponding names of these selected images. According to Fig. 9(a), we can observe that seven algorithms have approximate ISNR values for the first six images while the GTDF has an obvious enhancement over the other six algorithms for the last four images. As

to the SSI-dB indices, the GTDF obtains almost the highest values ¨ except for the image 09-54-44 ¨. In Fig. 9(c), the SMPI-dB values with respect to different algorithms for each image clearly demonstrate that the GTDF surpasses all the rival algorithms. Therefore, the experimental results on the CCA ultrasound images validate again that the proposed GTDF has superior despeckling effect.

W. Cui, M. Li, G. Gong et al. / Biomedical Signal Processing and Control 55 (2020) 101625

11

Fig. 8. Comparison of despeckling results for seven algorithms namely (b) SRAD, (c) DPAD, (d) SRBF, (e) MIDAL, (f) HomoDDID, (g WSFBF and (h) GTDF on (a) a real B-mode CCA ultrasound image. The white boxes in (a) indicate the manually selected homogeneous regions in which the quantitative indices will be computed.

Fig. 9. Quantitative evaluation on (a) ISNR, (b) SSI-dB and (c) SMPI-dB for seven algorithms with 10 real B-mode CCA ultrasound images.

12

W. Cui, M. Li, G. Gong et al. / Biomedical Signal Processing and Control 55 (2020) 101625

Fig. 10. PSNR values at different (a) neighborhood size r, (b) iteration number k, (c) standard deviation  s and (d) proportional coefficients by running GTDF on Figs. 4(b) and 5 (b).

Table 5 Performance comparison for seven algorithms on the CCA ultrasound image in Fig.8(a). Algorithm

ISNR

SSI-dB

SMPI-dB

SRAD DPAD SRBF MIDAL HomoDDID WSFBF GTDF

2.2783 2.3028 3.1258 2.4584 2.4703 3.4629 4.1761

0.3050 0.3510 1.6783 0.6350 0.6565 2.1237 2.9370

−0.5690 −0.5527 1.0946 0.0777 0.0961 1.9545 2.1759

The bold values indicate the best value in the corresponding index.

5.3. Parameter setting insensitivity for the GTDF In this section, we will give a detailed discussion on the parameters of the proposed GTDF, namely the local neighborhood size r, the total iteration number k, the standard deviation  s and the proportional coefficient . The test images are the synthetic images previously used in Figs. 4(b) and 5(b). PSNR is used as the objective index. Fig. 10(a) gives the PSNR values at different neighborhood radius r from 3 to 19 with step size 2. It can be observed that the two curves reach the peaks at r = 7 and r =9, respectively. That is to say, the optimum neighborhood radius should be larger if the corrupting speckle becomes more serious. Additionally, the PSNR values rapidly approximate to the stable values for the two images even

though the size r continues to increase. This demonstrates that the GTDF is insensitive to the parameter r. As to the total iteration number k, both test images have the peak values at the same second iteration and almost no variation with the increase of the iteration when k is larger than 5, as can be seen in Fig. 10(b). Therefore, the proposed GTDF algorithm requires only few iterations to converge. As the important parameter of the spatial weight gt , the standard deviation  s has a key influence on despeckling performance. Such a conclusion can be drawn from the variation curves in Fig. 10(c). It can be observed that for the image with moderate speckle (Fig. 4(b)), the optimum setting is  s = 1.3 while  s = 2.5 for the image with heavy speckle (Fig. 5(b)). Thus, this parameter also needs to be adjusted proportionally depending on the speckle strength. Besides, the flat tails of the two curves demonstrate again that GTDF is insensitive to the parameter  s . The proportional coefficient directly influences the scale parameter s in Eq. (20). It can be concluded from Fig. 10(d) that an unsuitable small value is apt to output a poor PSNR because the PSNR values have large variations when is less than 0.5. According to its insensitiveness to a larger , we suggest the setting value increases gradually from an initialization value = 0.5. 5.4. Comparison for statistical distribution weight In order to test the influence of the distribution model of ultrasound image, we substitute the Rayleigh distribution weight

W. Cui, M. Li, G. Gong et al. / Biomedical Signal Processing and Control 55 (2020) 101625

Fig. 11. Quantitative evaluation on (a) ISNR, (b) SSI-dB and (c) SMPI-dB for GTF-Ra and GTDF with 10 real B-mode tumor ultrasound images.

Fig. 12. Quantitative evaluation on (a) ISNR, (b) SSI-dB and (c) SMPI-dB for GTF-Ra and GTDF with 10 real B-mode CCA ultrasound images.

13

14

W. Cui, M. Li, G. Gong et al. / Biomedical Signal Processing and Control 55 (2020) 101625

Table 6 Running time (s) for seven algorithms on four test images. Algorithm

Fig.4(b) 256 × 256

Fig.5(b) 256 × 256

Fig.6(a) 356 × 284

Fig.8(a) 384 × 399

Average time

SRAD DPAD SRBF MIDAL HomoDDID WSFBF GTDF

1.4445 1.5376 0.4488 9.4599 35.0343 0.3141 1.1236

1.3872 1.3502 0.5430 9.4965 36.5364 0.3249 1.4884

2.0141 2.0808 0.7548 17.5483 57.3789 0.3671 2.3412

3.2800 3.2299 1.1489 22.4081 85.9975 0.5559 3.2854

2.0315 2.0496 0.7239 14.7282 53.7368 0.3905 2.0597

adopted by SRBF (see Section 2.4) for the Fisher-Tippett distribution weight of GTDF to generate another instantiation of the GTF framework, hereinafter called as GTF-Ra. Consequently, GTF-Ra is different from GTDF only in terms of the statistical Rayleigh distribution weight t , which should be calculated using Eqs. (13) and (14) instead. The experiments have been performed on both the tumor images and the CCA images and the quantitative indices obtained by GTDF and GTF-Ra have been plotted in Figs. 11 and 12. By contrast, the proposed GTDF surpasses GTF-Ra over five of the total six sub-figures, except for Fig. 11(c). Therefore, we can make a conclusion that for despeckling ultrasound images, the Fisher-Tippett distribution weight has better performance than the Rayleigh distribution weight. 5.5. Comparison for running time To compare the running speed for seven algorithms, we list the average CPU time in Table 6 by individually executing each algorithm three times for each test image. HomoDDID processes the logarithmic image using the DDID algorithm, which includes three steps and in each step, there are M×N (total image pixels) complex neighborhood operations in both the spatial domain and the frequency domain. As a result, HomoDDID requires the longest elapsed time, as shown in Table 6. SRBF and WSFBF belong to noniterative algorithms and just do some uncomplicated operations, thus they have faster running speed. Among the remaining four iterative algorithms, SRAD and DPAD require less time due to their simple computations within each iteration while MIDAL is the case that more sophisticated operations result in more running time. The proposed GTDF has a rapid implementation in less than four seconds because its convergence requires only few iterations. Furthermore, the proposed GTDF can strengthen its applicability by replacing the MATLAB programing with a C program.

such as CT imaging and MRI imaging by introducing their respective statistical distribution model. Declaration of Competing Interest None. Acknowledgments This work was supported in part by the National Key Research Program (2016YFB0800403), the National Natural Science Foundations of China (61871258, 61272237, 61102155, U1401252) and the Research Fund for Excellent Dissertation of China Three Gorges University (2019SSPY073). The authors would like to thank anonymous reviewers for constructive comments, which improved quality of the paper. Appendix A Assuming that the intensities I(j) j ∈ Nr (x+t) in the local neighborhood Nr (x+t) are samples of independent and identically distributed observations coming from the Rayleigh random variable M with the PDF p(M) =

t2



−

exp

M2



(A.1)

2t2

where  t is the scale parameter of the Rayleigh distribution. Then the likelihood function is given by the following joint density function: L1 =

 I(j) j ∈ Nr (x+t)

6. Conclusions In this paper, we have presented a local spatial filtering framework called as the guided trilateral filter (GTF) to deal with the restoration problem when a noisy image has a known statistical distribution model. The derivation of the GTF ascribes to the local ML estimation from the probabilities of the residuals between the filtered image and the noisy image. Comparing with classical bilateral filter, the GTF supplements the statistical distribution information and computes the corresponding weights in the analogous fashion with the guided filter. As a result, classical bilateral filter (BF), robust BF, joint/cross BF and speckle reduction BF can be considered as special cases of the GTF. For despeckling of ultrasound image, we introduce the Fisher-Tippett distribution model into the GTF and thus propose the guided trilateral despeckling filter (GTDF). We have conducted some experiments on synthetic and real images to contrast quantitatively with several classical algorithms. Experimental results have validated that the GTDF has superior despeckling performance, fast and robust convergence as well as insensitiveness to parameter setting. In future, we will continue further work about the GTF to serve new applicable fields

M

t2

−

exp

I 2 (j)

 (A.2)

2t2

Applying the log transformation to Eq. (A.2), we can derive the log-likelihood function as Ln(L1 ) =





ln I(j) −

j ∈ Nr (x+t)

 

ln t2 −

j ∈ Nr (x+t)

 I 2 (j) j ∈ Nr (x+t)

2t2

(A.3)

Taking the derivative of the above Ln(L1 ) with respect to t2 and equating to zero, we obtain



1

−

Nr (x + t)

t2

+

1



2t4 j ∈ Nr (x+t)

I 2 (j) = 0

(A.4)

where || Nr (x + t) || represents the total number of the pixels in the neighborhood Nr (x + t). Finally the ML estimator of t2 can be obtained as the solution of the above Eq. (A.4): ˆ t2 =

1



2

Nr (x + t)

 j ∈ Nr (x+t)

I 2 (j)

(A.5)

W. Cui, M. Li, G. Gong et al. / Biomedical Signal Processing and Control 55 (2020) 101625

or equivalently,



ˆ t =



1



2

Nr (x + t)

I 2 (j)

(A.6)

j ∈ Nr (x+t)

Appendix B In the local neighborhood Nr (x+t), the intensities ˆI (j) j ∈ Nr (x+t) can be deemed as samples of independent and identically distributed observations coming from the Fisher-Tippett random variable Z in Eq. (17). The likelihood function is given by the following joint density function:

⎛

L2 =



ˆ eI (j)

j ∈ Nr (x+t)

−1

2t



⎜

exp ⎝ˆI (j) −

2 ⎞

ˆ

eI (j) − 1

⎟ ⎠

22t

(B.1)

Applying the log transformation to Eq. (B.1), we can derive the log-likelihood function as





Ln(L2 ) =



ˆ

j ∈ Nr (x+t)



+



ln eI (j) − 1 −

j ∈ Nr (x+t)



ˆI (j) −

j ∈ Nr (x+t)

 

ln 2t



2

ˆ

eI (j) − 1

(B.2)

22t

j ∈ Nr (x+t)

Taking the derivative of the above Ln(L2 ) with respect to 2t and equating to zero, we obtain



1

−

Nr (x + t)

2t

+

1

 

24t j ∈ Nr (x+t)

2

ˆ

eI (j) − 1

=0

(B.3)

where || Nr (x + t) || represents the total number of the pixels in the neighborhood Nr (x + t). Finally the ML estimator of 2t can be obtained as the solution of the above Eq. (B.3):  ˆ 2t =

1



2

Nr (x + t)

 

ˆ

2

eI (j) − 1

(B.4)

j ∈ Nr (x+t)

References [1] O.V. Michailovich, A. Tannenbaum, Despeckling of medical ultrasound images, IEEE Trans. Ultrason. Ferroelectr. Freq. Control 53 (1) (2006) 64–78, http://dx.doi.org/10.1109/TUFFC.2006.1588392. [2] S.H.C. Ortiz, T. Chiu, M.D. Fox, Ultrasound image enhancement: a review, Biomed. Signal Process. Control 7 (5) (2012) 419–428, http://dx.doi.org/10. 1016/j.bspc.2012.02.002. [3] J. Park, J. Kang, J. Chang, Y. Yoo, Speckle reduction techniques in medical ultrasound imaging, Biomed. Eng. Lett. 4 (1) (2014) 32–40, http://dx.doi.org/ 10.1007/s13534-014-0122-6. [4] A. Perperidis, Postprocessing approaches for the improvement of cardiac ultrasound b-mode images: a review, IEEE Trans. Ultrason. Ferroelectr. Freq. Control 63 (3) (2016) 470–485, http://dx.doi.org/10.1109/TUFFC.2016. 2526670. [5] X. Feng, X. Guo, Q. Huang, Systematic evaluation on speckle suppression methods in examination of ultrasound breast images, Appl. Sci. 7 (1) (2017) 37, http://dx.doi.org/10.3390/app7010037. [6] T. Loupas, W.N. McDicken, P.L. Allan, An adaptive weighted median filter for speckle suppression in medical ultrasonic images, IEEE Trans. Circuits Syst. 36 (1) (1989) 129–135, http://dx.doi.org/10.1109/31.16577. [7] Y. Guo, H. Cheng, J. Tian, Y. Zhang, A novel approach to speckle reduction in ultrasound imaging, Ultrasound Med. Biol. 35 (4) (2009) 628–640, http://dx. doi.org/10.1016/j.ultrasmedbio.2008.09.007. [8] A. Thakur, R. Anand, Speckle reduction in ultrasound medical images using adaptive filter based on second order statistics, J. Med. Eng. Technol. 31 (4) (2007) 263–279, http://dx.doi.org/10.1080/03091900600718402. [9] S. Balocco, C. Gatta, O. Pujol, J. Mauri, P. Radeva, SRBF: speckle reducing bilateral filtering, Ultrasound Med. Biol. 36 (8) (2010) 1353–1363, http://dx. doi.org/10.1016/j.ultrasmedbio.2010.05.007.

15

[10] Y. Yu, S. Acton, Speckle reducing anisotropic diffusion, IEEE Trans. Image Process. 11 (11) (2002) 1260–1270, http://dx.doi.org/10.1109/TIP.2002. 804276. [11] S. Aja-Fernández, C. Alberola-López, On the estimation of the coefficient of variation for anisotropic diffusion speckle filtering, IEEE Trans. Image Process. 15 (9) (2006) 2694–2701, http://dx.doi.org/10.1109/TIP. 2006.877360. [12] K. Krissian, C.F. Westin, R. Kikinis, K. Vosburgh, Oriented speckle reducing anisotropic diffusion, IEEE Trans. Image Process. 16 (5) (2007) 1412–1424, http://dx.doi.org/10.1109/TIP.2007.891803. [13] G. Ramos-Llordén, G. Vegas-Sánchez-Ferrero, M. Martin-Fernandez, C. Alberola-López, S. Aja-Fernández, Anisotropic diffusion filter with memory based on speckle statistics for ultrasound images, IEEE Trans. Image Process. 24 (1) (2015) 345–358, http://dx.doi.org/10.1109/TIP.2014.2371244. [14] P. Coupe, P. Hellier, C. Kervrann, C. Barillot, Nonlocal means-based speckle filtering for ultrasound images, IEEE Trans. Image Process. 18 (10) (2009) 2221–2229, http://dx.doi.org/10.1109/TIP.2009.2024064. [15] Y. Guo, Y. Wang, T. Hou, Speckle filtering of ultrasonic images using a modified non local-based algorithm, Biomed. Signal Process. Control 6 (2) (2011) 129–138, http://dx.doi.org/10.1016/j.bspc.2010.10.004. [16] Y. Zhan, M. Ding, L. Wu, X. Zhang, Nonlocal means method using weight refining for despeckling of ultrasound images, Signal Process. 103 (2014) 201–213, http://dx.doi.org/10.1016/j.sigpro.2013.12.019. [17] P. Sudeep, P. Palanisamy, J. Rajan, H. Baradaran, L. Saba, A. Gupta, J. Suri, Speckle reduction in medical ultrasound images using an unbiased non-local means method, Biomed. Signal Process. Control 28 (2016) 1–8, http://dx.doi. org/10.1016/j.bspc.2016.03.001. [18] G. Aubert, J. Aujol, A variational approach to removing multiplicative noise, SIAM J. Appl. Math. 68 (4) (2008) 925–946, http://dx.doi.org/10.1137/ 060671814. [19] J. Shi, S. Osher, A nonlinear inverse scale space method for a convex multiplicative noise model, SIAM J. Imaging Sci. 1 (3) (2008) 294–321, http:// dx.doi.org/10.1137/070689954. [20] Y. Huang, M. Ng, Y. Wen, A new total variation method for multiplicative noise removal, SIAM J. Imaging Sci. 2 (1) (2009) 20–40, http://dx.doi.org/10. 1137/080712593. [21] J.M. Bioucas-Dias, M.A. Figueiredo, Multiplicative noise removal using variable splitting and constrained optimization, IEEE Trans. Image Process. 19 (7) (2010) 1720–1730, http://dx.doi.org/10.1109/TIP.2010.2045029. [22] J. Huang, X. Yang, Fast reduction of speckle noise in real ultrasound images, Signal Process. 93 (4) (2013) 684–694, http://dx.doi.org/10.1016/j.sigpro. 2012.09.005. [23] A. Achim, A. Bezerianos, P. Tsakalides, Novel Bayesian multiscale method for speckle removal in medical ultrasound images, IEEE Trans. Med. Imaging 20 (8) (2001) 772–783, http://dx.doi.org/10.1109/42.938245. [24] A. Khare, M. Khare, Y. Jeong, H. Kim, M. Jeon, Despeckling of medical ultrasound images using Daubechies complex wavelet transform, Signal Process. 90 (2) (2010) 428–439, http://dx.doi.org/10.1016/j.sigpro. 2009.07. 008. [25] J. Zhang, G. Lin, L. Wu, C. Wang, Y. Cheng, Wavelet and fast bilateral filter based de-speckling method for medical ultrasound images, Biomed. Signal Process. Control 18 (2015) 1–10, http://dx.doi.org/10.1016/j.bspc.2014.11. 010. [26] J. Zhang, G. Lin, L. Wu, Y. Cheng, Speckle filtering of medical ultrasonic images using wavelet and guided filter, Ultrasonics 65 (2016) 177–193, http://dx.doi. org/10.1016/j.ultras.2015.10.005. [27] C.A. Deledalle, L. Denis, S. Tabti, F. Tupin, MuLoG, or how to apply Gaussian denoisers to multi-channel SAR speckle reduction? IEEE Trans. Image Process. 26 (9) (2017) 4389–4403, http://dx.doi.org/10.1109/TIP.2017.2713946. [28] C. Tomasi, R. Manduchi, Bilateral filtering for gray and color images, Proceedings of the Sixth IEEE International Conference on Computer Vision (ICCV) (1998) 839–846, http://dx.doi.org/10.1109/ICCV.1998.710815. [29] P. Perona, J. Malik, Scale-space and edge detection using anisotropic diffusion, IEEE Trans. Image Process. 12 (7) (1990) 629–639, http://dx.doi.org/10.1109/ 34.56205. [30] K. He, J. Sun, X. Tang, Guided image filtering, IEEE Trans. Pattern Anal. Mach. Intell. 35 (6) (2013) 1397–1409, http://dx.doi.org/10.1109/TPAMI.2012.213. [31] C. Knaus, M. Zwicker, Dual-domain image denoising, Proceedings of the 20th IEEE International Conference on Image Processing (ICIP) (2013) 440–444, http://dx.doi.org/10.1109/ICIP.2013.6738091. [32] K. Dabov, A. Foi, V. Katkovnik, K. Egiazarian, Image denoising by sparse 3-D transform-domain collaborative filtering, IEEE Trans. Image Process. 16 (8) (2007) 2080–2095, http://dx.doi.org/10.1109/TIP.2007.901238. [33] G. Petschnigg, R. Szeliski, M. Agrawala, M. Cohen, H. Hoppe, K. Toyama, Digital photography with flash and no-flash image pairs, ACM Trans. Graph. 23 (3) (2004) 664–672, http://dx.doi.org/10.1145/1015706.1015777. [34] E. Eisemann, F. Durand, Flash photography enhancement via intrinsic relighting, ACM Trans. Graph. 23 (3) (2004) 673–678, http://dx.doi.org/10. 1145/1015706.1015778. [35] L. Caraffa, J.P. Tarel, P. Charbonnier, The guided bilateral filter: when the joint/cross bilateral filter becomes robust, IEEE Trans. Image Process. 24 (4) (2015) 1199–1208, http://dx.doi.org/10.1109/TIP.2015.2389617. [36] T. Dai, W. Lu, W. Wang, J. Wang, S. Xia, Entropy-based bilateral filtering with a new range kernel, Signal Process. 137 (2017) 223–234, http://dx.doi.org/10. 1016/j.sigpro.2017.02.005. [37] H. Yu, L. Zhao, H. Wang, Image denoising using trivariate shrinkage filter in the wavelet domain and joint bilateral filter in the spatial domain, IEEE Trans.

16

[38]

[39]

[40]

[41]

[42]

[43]

[44]

W. Cui, M. Li, G. Gong et al. / Biomedical Signal Processing and Control 55 (2020) 101625 Image Process. 18 (10) (2009) 2364–2369, http://dx.doi.org/10.1109/TIP.2009. 2026685. K. Zhang, G. Lafruit, R. Lauwereins, L. Gool, Constant time joint bilateral filtering using joint integral histograms, IEEE Trans. Image Process. 21 (9) (2012) 4309–4314, http://dx.doi.org/10.1109/TIP.2012.2198220. T.C. Aysal, K.E. Barner, Rayleigh-maximum-likelihood filtering for speckle reduction of ultrasound images, IEEE Trans. Med. Imaging 26 (5) (2007) 712–727, http://dx.doi.org/10.1109/10.1109/TMI.2007.895484. A. Sarti, C. Corsi, E. Mazzini, C. Lamberti, Maximum likelihood segmentation of ultrasound images with Rayleigh distribution, IEEE Trans. Ultrason. Ferroelectr. Freq. Control 52 (6) (2005) 947–960, http://dx.doi.org/10.1109/ TUFFC.2005.1504017. Y. Zimmer, R. Tepper, S. Akselrod, A lognormal approximation for the gray level statistics in ultrasound images, Proceedings of the 22nd Annual International Conference of the IEEE Engineering in Medicine and Biology Society (2000) 2656–2661, http://dx.doi.org/10.1109/IEMBS.2000.901405. Z. Tao, H.D. Tagare, Tunneling descent level set segmentation of ultrasound images, Proceedings of the 19th International Conference on Information Processing in Medical Imaging (IPMI) (2005) 750–761, http://dx.doi.org/10. 1007/11505730 62. Z. Tao, H.D. Tagare, J.D. Beaty, Evaluation of four probability distribution models for speckle in clinical cardiac ultrasound images, IEEE Trans. Med. Imaging 25 (11) (2006) 1483–1491, http://dx.doi.org/10.1109/TMI.2006. 881376. G. Slabaugh, G. Unal, M. Wels, T. Fang, B. Rao, Statistical region-based segmentation of ultrasound images, Ultrasound Med. Biol. 35 (5) (2009) 781–795, http://dx.doi.org/10.1016/j.ultrasmedbio.2008.10.014.

[45] C.A.N. Santos, D.L.N. Martins, N.D.A. Mascarenhas, Ultrasound image despeckling using stochastic distance-based BM3D, IEEE Trans. Image Process. 26 (6) (2017) 2632–2643, http://dx.doi.org/10.1109/TIP.2017.2685339. [46] L. Zhang, L. Zhang, X. Mou, D. Zhang, FSIM: a feature similarity index for image quality assessment, IEEE Trans. Image Process. 20 (8) (2011) 2378–2386, http://dx.doi.org/10.1109/TIP.2011.2109730. [47] F. Sattar, L. Floreby, G. Salomonsson, B. Lovstrom, Image enhancement based on a nonlinear multiscale method, IEEE Trans. Image Process. 6 (6) (1997) 888–895, http://dx.doi.org/10.1109/83.585239. [48] C.P. Loizou, C.S. Pattichis, Despeckle filtering for ultrasound imaging and video Algorithms and Software, Volume I, second ed., Morgan & Claypool, California, 2015, http://dx.doi.org/10.2200/S00641ED1V01Y201504ASE014. [49] S.G. Dellepiane, E. Angiati, Quality assessment of despeckled SAR images, IEEE J. Sel. Top. Appl. Earth Observ. Remote Sens. 7 (2) (2014) 691–707, http://dx. doi.org/10.1109/JSTARS.2013.2279501. [50] Y. Sheng, Z.G. Xia, A comprehensive evaluation of filters for radar speckle suppression, Proceedings of International Geoscience and Remote Sensing Symposium (IGARSS) (1996) 1559–1561, http://dx.doi.org/10.1109/IGARSS. 1996.516730. [51] A. Pizurica, W. Philips, I. Lemahieu, M. Acheroy, A versatile wavelet domain noise filtration technique for medical imaging, IEEE Trans. Med. Imaging 22 (3) (2003) 323–331, http://dx.doi.org/10.1109/TMI.2003.809588.