Despeckling of ultrasound medical images using nonlinear adaptive anisotropic diffusion in nonsubsampled shearlet domain

Biomedical Signal Processing and Control 14 (2014) 55–65

Contents lists available at ScienceDirect

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

Despeckling of ultrasound medical images using nonlinear adaptive anisotropic diffusion in nonsubsampled shearlet domain Deep Gupta ∗ , R.S. Anand, Barjeev Tyagi Department of Electrical Engineering, Indian Institute of Technology Roorkee, Roorkee 247667, India

a r t i c l e

i n f o

Article history: Received 5 March 2014 Received in revised form 8 June 2014 Accepted 26 June 2014 Keywords: Denoising Diffusion Edge preservation Nonsubsampled shearlet transform (NSST) Ultrasound Speckle

a b s t r a c t Despeckling is of great interest in ultrasound medical images. The inherent limitations of acquisition techniques and systems introduce the speckles in ultrasound images. These speckles are the main factors that degrade the quality and most importantly texture information present in ultrasound images. Due to these speckles, experts may not be able to extract correct and useful information from the images. This paper presents an edge preserved despeckling approach that combines the nonsubsampled shearlet transform (NSST) with improved nonlinear diffusion equations. As a new image representation method with the different features of localization, directionality and multiscale, the NSST is utilized to provide the effective representation of the image coefficients. The anisotropic diffusion approach is applied to the noisy coarser NSST coefficients to improve the noise reduction efficiency and effectively preserves the edge features. In the diffusion process, an adaptive gray variance is also incorporated with the gradient information of eight connected neighboring pixels to preserve the edges, effectively. The performance of the proposed method is evaluated by conducting extensive simulations using both the standard test images and several ultrasound medical images. Experiments show that the proposed method provides an improvement not only in noise reduction but also in the preservation of more edges as compared to several existing methods. © 2014 Elsevier Ltd. All rights reserved.

1. Introduction Currently, the research in medical imaging has produced many different imaging modalities for the clinical purpose. Among the different imaging modalities, ultrasound imaging is of a particular interest in the medical diagnosis of neck, chest, liver, abdominal cavity, gallbladder, pancreas, spleen, adrenal glands, kidney, prostate and scrotum due to its cost effectiveness, portability, acceptability and safety [1]. However, the ultrasound images are of relatively poor quality due to speckles (considered as a multiplicative noise) present in them. In addition to multiplicative noise, sometimes ultrasound (US) images also suffer from random additive Gaussian noise. The presence of speckles affects the human interpretation as well as accuracy of computer-assisted methods. Therefore, detection and enhancement of the boundaries between different cavities and organs are of great need in ultrasound images and is considered as a challenging problem. Sometimes, this process may suppress the important details of the ultrasound images.

∗ Corresponding author. Tel.: +91 9358190782. E-mail addresses: [email protected] (D. Gupta), [email protected] (R.S. Anand), [email protected] (B. Tyagi). http://dx.doi.org/10.1016/j.bspc.2014.06.008 1746-8094/© 2014 Elsevier Ltd. All rights reserved.

Thus, noise reduction algorithms should be designed in such a manner that they suppress the noise as much as possible without any significant loss of information presented in the US images. Speckle reduction can be done in two categories viz. image averaging and image filtering methods [2]. Image averaging methods usually lead to the loss of spatial resolution. Filtering methods are practical alternative for most of the clinical applications. Furthermore, it can be classified as single scale spatial filtering such as linear [3,4], nonlinear, adaptive methods [4,5], multiscale spatial filtering such as diffusion based methods [2,6–8] and others multiscale methods in different transform domains such as pyramid [9], wavelet [10], ridgelet [11] and curvelet [12]. Simple mathematical linear filters such as mean filter degrade the sharp transitions (line and edges) of the image [4]. Most popular nonlinear filters such as median filter are employed to all the pixels whether they are corrupted or not. The weighted median filter is also used for the noise reduction. It is able to retain more edges than the classical median filter [5]. However, there is a loss of resolution by suppressing the fine details. In the category of multiscale spatial filtering, partial differential equations (PDE) are used to implement some filters for denoising purpose. The most popular methods include Lee filter [13], Kuan method [14] and anisotropic diffusion equations proposed by Perona and Malik [6]. So it is called P–M equations that

56

D. Gupta et al. / Biomedical Signal Processing and Control 14 (2014) 55–65

provide a method for image smoothing. A lot of work has been done with anisotropic diffusion equations in such a way that the important structural information can be retained in the denoised images [8,15–17]. Besides the diffusion methods, various multiscale methods have been employed for noise reduction [10,12,18–21]. Currently, lots of research work on image processing is concentrated in the transform domain. In that series, wavelet thresholding has been presented as a true signal estimation technique that exploits the capabilities of wavelet transform (WT) for signal denoising [18,19]. The main strength of wavelet thresholding is its capability to process the different frequency components of an image, separately [22]. Thus, many efforts have been made to suppress the different types of noise and to overcome the limitations of spatial domain filtering by wavelet thresholding methods [19–21]. However, it may lead to the formation of some visual artifacts around sharp discontinuities. The WT is able to efficiently represent a function with one dimensional singularity. However, it is less efficient in reflecting the sharp transitions such as line and curve singularities due to its limitation of direction [23]. To overcome this limitation, ridgelet transform has been proposed to provide the information about the orientation of the linear edges [11]. However, it does not represent the two dimensional singularities. Donoho et al. [12] has presented curvelet transform (CT) used to represent two dimensional singularities with the smooth curve and provides better denoising and edge preservation results. Contourlet transform proposed by Vetterli et al. [24] performs well in noise reduction due to application of multiscale Laplacian pyramid (LP) followed by directional filter banks. However, it has less directional features than curvelets. To represent the edges more efficiently, Labate et al. introduced a new multiscale analysis tool called shearlet that has all properties like other MGA tools as multiscale, localization, anisotropy and directionality [25,26]. The decomposition of the shearlet transforms (ST) consists of multiscale and multidirectional decomposition that are similar to contourlets except that there is no limitation on the number of directions. Shearlets can also be constructed in discrete domain realized by combination of the Laplacian pyramid (LP) and directional filters, but still the lack of shift invariance problem cannot be overcome. Easley et al. [27] proposed nonsubsampled shearlet transform (NSST) that is realized by nonsubsampled Laplacian pyramid (NSLP) and several shearing filters. The NSST also provides the flexible directional selectivity and shift invariance [27,28]. Based on the above concept, the present work combines the NSST with improved nonlinear diffusion equations and thresholding scheme for despeckling of the US medical images. The paper is structured as follows. Section 2 presents the methodologies used for the proposed method. Section 3 illustrates the implementation steps of the proposed despeckling algorithm. In Section 4, various experimental results are discussed and compared with the several other methods in terms of the different performance measures. Conclusions are drawn in Section 5. 2. Methodology 2.1. Anisotropic diffusion Anisotropic diffusion is modeled by Perona and Malik for defining a scale space image [6]. This model is an extension of the heat equation that is based on the partial differential equation (PDE). Let s(x, y; t) is an image with coordinates (x, y) at time t, and then the continuous anisotropic diffusion is defined as ∂s(x, y; t) = div[g(x, y; t)∇ s(x, y; t)] ∂t

(1)

where div is the divergence operator, g is the diffusion coefficient and  is a gradient operator with respect to space variables. The

diffusion model becomes isotropic, if g is a constant parameter. When g is a function of directional parameters, the diffusion model becomes anisotropic. Perona–Malik (PM) suggested the two well known diffusion coefficients defined as, g(f ) =

1

1+

 f 2

(2a)



  2  f

g(f ) = exp −

(2b)



where f = | ∇ s| and the parameter  serves as a threshold of gradient size. A smaller gradient is diffused, while positioning a larger gradient is treated as edges. Instead of having many computational and theoretical properties, there is one serious problem with the diffusion method. It is very sensitive to the noise which may introduce large oscillations of the gradient. Furthermore, the PM method cannot differentiate between true edges and noises. Another problem is the stair casing effects that arise around the smooth edges [29]. To provide the solution of this problem, Catte et al. [30] proposed that a Gaussian kernel G is convolved with the images to reduce the effect of noise and provide better estimation of the local gradient. However, it is very sensitive to the number of diffusion iteration by considering only the gradient information of the pixel. Normally, large gradient values are treated as edges, but sometimes the important details along with edges available in the image may have low gradient magnitude [17]. Therefore, the gray level variance is incorporated along with the gradient of the pixels to evaluate the diffusion coefficients. In the PM method, the derivative term (s) is calculated using a template of four closest neighbors of a pixel (x, y). This term can be evaluated more accurately by considering the large number of the neighboring pixels within a template and the corresponding algorithm improves the quality of the images. Moreover in this work, eight nearest neighboring pixels are used within 3 × 3 template to evaluate the gradient term and the adaptive gray variance is also included along with the gradient to estimate the diffusion coefficients. 2.2. Nonsubsampled shearlet transform The NSST is an extension of the WT in multidimensional and multidirectional case which combines the multiscale and direction analysis, separately. Firstly, the NSLP is used to decompose an image into low and high-frequency components, and then direction filtering is employed to get the different subbands and different direction shearlet coefficients. Direction filtering is achieved using the shear matrix which provides many more directions. The introduction of the NSST is given as follows: Consider a two-dimensional affine system with composite dilations as [27,28] ADS =



j,k,m (x)

= | det D|j/2 (S k Dj x − m) : j, k ∈ Z, m ∈ Z2



(3)

where D refers to the anisotropic matrix, S denotes the shear matrix and j, k and m are scale, direction and shift parameter, respectively. The D and S are both 2 × 2 invertible matrices and |det S| = 1, the elements of the system are called composite wavelet, if it forms a Parseval frame for L2 (R2 ), which is also an affine like system. For any f ∈ L2 (R2 )



j,k,m

  f,

j,k,m

2  = ||f ||2

(4)









d 0 d1/2 0 or where 1/2 0 d 0 d d > 0 controls the scale of shearlets, which ensures that the frequency support of shearlets becomes increasingly stretched at finer The anisotropic dilation matrix

D. Gupta et al. / Biomedical Signal Processing and Control 14 (2014) 55–65

57

Fig. 1. An illustration of three level multiscale and multidirectional NSST decomposition.

 scales. The shear matrix S =

1 0

s 1



 or

1 0 s 1

 controls only the



direction of shearlets. Generally, assuming d = 4 and s = 1, where D = D0 =



D1 =



S1

2 0 0 4

1 0 1 1

(0) (x) j,k,m





is anisotropic dilation matrix and S = S0



4 0 0 2 1 1 0 1

 or

 or

is shear matrix [28]. The shearlet transform function is,

3

= 2j 2

(0)

j

(S0k D0 x − m) and

(1) (x) j,k,m

3

= 2j 2

(1)

j

(S1k D1 x − m) (6)

2.3. NSST thresholding In the recent years, various thresholding (shrinkage) schemes are considered for image denoising. The thresholding scheme provides the threshold coefficients by comparing the transformed coefficients against a threshold to remove the noise from a signal while preserving the most important information of the original signal. In the present study, we perform the hard thresholding using the NSST with cycle spinning approach, independently in each finer NSST noisy subbands. The fundamental principle of hard thresholding using the NSST is same as the thresholding using wavelet and curvelet transform. To achieve the threshold coefficients from the noisy NSST coefficients (sNSST (x, y)), the thresholding expression is given as, sˆNSST (x, y) = thr (sNSST (x, y))

ˆ (0) () = ˆ (0) (1 , 2 ) = where j≥0, −2j ≤ k ≤ 2j − 1, m ∈ Z2 , ˆ 1 (1 ) ˆ 2 (2 /1 ) and ˆ (1) () = ˆ (1) (1 , 2 ) = ˆ 1 (2 ) ˆ 2 (1 /2 )  2 : |1 |≥ 18 , |2 /1 | ≤ 1 For any ( 1 ,  2 ) ∈ D0 , D0 = (1 , 2 ) ∈ R and D1 = (0) (x) j,k,m





2 : |2 |≥ 18 , |1 /2 | ≤ 1 consist of supports of (1 , 2 ) ∈ R (1) (x). j,k,m

and As mentioned above, the NSLP and shearing filters (ShF) are utilized to provide the multiscale and multidirectional decomposition. At each level of the NSLP decomposition, one high frequency and one low frequency sub images are produced and further the low frequency subband is decomposed, iteratively. At the decomposition level m = 3, an image is decomposed into (m + 1) = 4 subbands with the same size of the source image in which one (first) is the low frequency subband and others m = 3 are the high frequency subband images. Shearing filter is also used to decompose the high frequency sub images using the NSST decomposition without sub-sampling that satisfy the shift invariance property. Using the Shf at level K = [2, 3, 3], the high frequency subband images obtained from the NSLP at each decomposition level, 2K = [4, 8, 8] directional subband image coefficients are produced with the same size as the source images. Three level decomposition is shown in Fig. 1 which illustrates the NLSP and its corresponding directional decompositions. In this study, the numbers of shearing directions are taken to be 8, 8 and 4 from finer to coarser scale using three level decomposition of an image.

(7)

where sˆNSST (x, y) is the estimator of the unknown noiseless coefficients using hard thresholding and thresholding function is considered as,



sˆNSST (x, y)

=

sNSST (x, y),

|sNSST (x, y)|≥ˇ

0,

otherwise

(8)

where  is the variance of the noisy image and 2 is the variance of noisy subband decomposed using the NSST. The noise variance of each subband is estimated from the noisy NSST coefficients using Monte Carlo techniques as,

 



1  2 = sNSST (x, y)sNSST∗ (x, y) L L

L

(9)

x=1 y=1

where sNSST∗ (x, y) is the complex conjugate of sNSST (x, y) and L is the length of the subband at the  = 1, 2, 3, . . ., K − 1th scale. For each high frequency NSST subbands, the scale dependent parameter ˇ is computed by Eq. (10). ˇ=



log(L )

(10)

In this study, the hard thresholding combined with the concept of cycle spinning approach [31,32] is applied on every detailed

58

D. Gupta et al. / Biomedical Signal Processing and Control 14 (2014) 55–65

Fig. 2. Process flow of the proposed method based on NSST.

coefficients of the NSST decomposition. A summary of the algorithm is given as follows: 1. Apply the circular shifts on the noisy image. s (x, y) = circular shift(s(x, y), [xshift , yshift ]) The number of shifts depends on the length of the input vector scale. 2. Perform the multiscale decomposition of the shifted copies of an image using the NSST to obtain the noisy NSST coefficients. sNSST (x, y) = NSST(s (x, y)) 3. Apply the thresholding scheme on the noisy NSST coefficients to get the threshold NSST coefficients.

Let s(x, y) be an observed image that is a combination of the original image with some noise. The implementation steps of the proposed model are as follows: Step 1: Perform the nonsubsampled shearlet transform to an image s(x, y) at K different scales. Let sKNSST (x, y; t) denotes the approximation coefficients at Kth coarser scale and (sNSST (x, y)) denotes the detail subband coefficients of the image, where  = 1, 2, 3, . . ., K − 1. Step 2: Apply the thresholding function thr (·) to each detail coefficients (sNSST (x, y)) of the NSST decomposition and evaluate the threshold coefficients (ˆsNSST (x, y) = thr (sNSST (x, y))) for each scale from 1 to K − 1. Step 3: Now the approximation NSST coefficients are processed with modified nonlinear anisotropic diffusion. Let sKNSST (x, y; t) be the approximation NSST coefficients at coordinates (x, y) and iteration t. The discrete implementation of the anisotropic diffusion in Eq. (1) using the four nearest neighbor is given as,

⎡

sˆNSST (x, y) = thr (sNSST (x, y)) 4. Invert the multiscale decomposition to reconstruct the denoised image.

⎢ ⎥ ⎢ +gW (x, y; t)∇ W sKNSST (x, y; t) ⎥ ⎥ ⎥ ⎣ +gN (x, y; t)∇ N sKNSST (x, y; t) ⎦

sKNSST (x, y; t + 1) = sKNSST (x, y; t) +  ⎢ ⎢

+gS (x, y; t)∇ S sKNSST (x, y; t)

sˆ(x, y) = NSST−1 (ˆsNSST (x, y)) 5. After taking the NSST−1 , perform the inverse shift and resulting denoised image are shifted back to the original position and average the translated result to get the approximated image. ˆ S(x, y) = circular shift(sˆ (x, y), [−xshift , −yshift ])

3. Proposed method In the proposed method, the NSST based thresholding and modified nonlinear anisotropic diffusion equations are considered to increase the performance of the existing methods. As mentioned in the NSST based approach, the high-frequency detail coefficients are thresholded, but low-frequency approximation coefficients remain unchanged. It preserves more edges than the WT based approach, but it adds some unwanted fuzzy edges in a homogenous region like wavelets. In order to solve this problem, coarser approximation NSST coefficients are processed with the modified diffusion model. The process flow of the proposed method is shown in Fig. 2.

⎤

gE (x, y; t)∇ E sKNSST (x, y; t)

(11)

where 0 ≤  ≤ 1/4 and E, W, N and S refer to east, west, north and south, respectively. The gradient (∇ j ) is evaluated as nearest neighbor differences as follows:

⎫ ⎪ ⎪ ⎪ ⎪ d ⎪ ⎪ NSST (x − 1, y; t) − sNSST (x, y; t) ⎪ ⎪ s ⎪ K K ⎪ NSST ⎬ ∇ W sK (x, y; t) = ∇ E sKNSST (x, y; t) =

sKNSST (x + 1, y; t) − sKNSST (x, y; t)

d

⎪ ⎪ ⎪ ⎪ ⎪ d ⎪ ⎪ ⎪ NSST (x, y − 1; t) − sNSST (x, y; t) ⎪ ⎪ s ⎭ K K NSST ∇ s sK (x, y; t) =

∇ N sKNSST (x, y; t) =

sKNSST (x, y

+ 1; t) − sKNSST (x, y; t)

(12)

d

Here, d represents a spatial step size between two successive pixels in horizontal and vertical direction on the digital image. In the first stage of modified diffusion, all the eight nearest neighbors are utilized to evaluate the derivative term in 3 × 3 template. Fig. 3 shows the neighborhood pixel representation in 3 × 3 template. The derivation term is calculated in two groups. In the

D. Gupta et al. / Biomedical Signal Processing and Control 14 (2014) 55–65

59

4. Experimental results and discussion

Fig. 3. The discrete computational approach to evaluate the diffusion coefficient.

first group, we consider only those neighboring pixels which have a distance d from the center pixel sKNSST (x, y; t) as shown in the Eq. (12), √ while in the second group, the pixels which have the distance d 2 from the center pixel are as follows:

⎫ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ NSST NSST ⎪ (x + 1, y + 1; t) − s (x, y; t) s ⎪ K K ⎪ ∇ NW sKNSST (x, y; t) = √ ⎬ sNSST (x + 1, y + 1; t) − sKNSST (x, y; t) ∇ NE sKNSST (x, y; t) = K √ d 2 d 2

⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ NSST NSST sK (x − 1, y − 1; t) − sK (x, y; t) ⎪ ⎪ ⎭

sNSST (x + 1, y − 1; t) − sKNSST (x, y; t) ∇ SE sKNSST (x, y; t) = K √ d 2

∇ SW sKNSST (x, y; t) =

(13)

In the second stage, gray level variance is incorporated with the gradient term to evaluate the diffusivity term. The gray level variance is estimated at each pixel of the approximation NSST coefficients sKNSST (x, y; t) in block size 3 × 3 and normalized it. After including the variance (Varb ) and gradient in the diffusion equations, Eqs. (2a) and (2b) become as, g(f, Varb ) = 1+



f ·Varb 

2

(16)

 variance(noise free image) 

(17)

variance(remaining noise)

(2 S × Sˆ + C1 )(2S,Sˆ + C2 )

(18)

( 2S × 2ˆ + C1 )(S2 +  ˆ2 + C2 ) S

S

L

(14b)

where f = |∇ j sKNSST (x, y; t)| and j = E, W, N, S, NE, NW, SE, SW which indicate the east, west, north, south, north-east, north-west, southeast and south-west direction, respectively from the center pixels. By incorporating the above formulations, the modified nonlinear anisotropic diffusion model is given as,

1 where C1 = (K1 × L)2 and C2 = (K2 × L)2 , S,Sˆ = L−1 (S − k=1 k S )(Sˆ k − ˆ ) is the covariance in moving window, s is the original S

image and sˆ is the reconstructed image with mean S , Sˆ and variance  S and Sˆ , respectively. The EKI is used as another objective criterion to evaluate the edge preservation capability of despeckling methods, and is given by, EKI =

sKNSST (x, y; t + 1) = sKNSST (x, y; t)

⎡ gE (x, y; t)∇ E sNSST (x, y; t) + gW (x, y; t)∇ W sNSST (x, y; t) K

⎤

⎢ +gN (x, y; t)∇ N sKNSST (x, y; t) + gs (x, y; t)∇ s sKNSST (x, y; t) ⎥ ⎥ ⎣ +gNE (x, y; t)∇ NE sNSST (x, y; t) + gNW (x, y; t)∇ NW sNSST (x, y; t) ⎦

+ ⎢

K

√ MSE

A higher SNR value indicates a better denoising capacity. The SSIM index is used to study the structural and perceptual closeness between the original and filtered image. For the default value of K1 = 0.01 and K2 = 0.03 and the dynamic range (L) of the pixel is 255, the SSIM index is evaluated as,

(14a)



K

 MAX 

where MAX is the maximum pixel value of an image and the MSE is the mean squared error, which is equal to mean difference between the original and the restored image. Signal to noise ratio is a commonly used measure to quantify the noise suppression quality of despeckled image. The SNR is the ratio of the variance of noise free image to noise variance. It is defined as,

SSIM =

  2  f · Varb

g(f, Varb ) = exp −

PSNR = 20 log10

SNR = 10 log10

√ d 2

1

In this section, the proposed algorithm was applied to both the simulated images synthesized using MATLAB and real ultrasound images. The denoising abilities of the proposed method are presented qualitatively and quantitatively and also compared with some related noise reduction techniques. To quantitatively evaluate the noise suppression and edge preservation capability, different performance measures are used in our experimentation that measures the performance of proposed algorithm and others. The peak signal to noise ratio (PSNR) [4] and signal to noise ratio (SNR) [33] are the parameters to measure the noise suppression capability of the despeckling methods. However, they do not reflect the performance of the structure and edge preservation capability of the denoising methods. For the quantitative evaluations of the structural similarity and edge preservation in the denoised image, another indices such as structural similarity index (SSIM) [34], and edge keeping index (EKI) [2] are used. The PSNR analysis uses a standard mathematical model to measure an objective difference between two images. It estimates the noise suppression performance in the reconstructed image with respect to an original image. The PSNR is defined as,

(15)

K

+gSE (x, y; t)∇ SE sKNSST (x, y; t) + gSW (x, y; t)∇ SW sKNSST (x, y; t)

Step 4: Finally, an enhanced denoised image is obtained after taking the inverse NSST on the processed subband coefficients from steps 2 and 3.



I( s − s, ˆs − ˆs)

(19)

I( s − s, s − s) · I( ˆs − ˆs, ˆs − ˆs)

where s and ˆs represent the high pass filtered versions of s and sˆ, respectively, calculated by 3 × 3 standard approximation of Laplacian operator with the mean intensity value as s and ˆs. The function I is defined as I(s1 , s2 ) =

N M



s1 (x, y)s2 (x, y), where

x=1 y=1

M × N is the size of the image. This quantitative performance measure should be close to unity for better preservation of edges.

60

D. Gupta et al. / Biomedical Signal Processing and Control 14 (2014) 55–65

Fig. 4. Visual comparisons of the various noise suppressing methods on simulated kidney image. (a) Noisy image degraded with speckle noise. (b) Method 1 [35]. (c) Method 2 [36,37]. (d) Method 3 [37,38]. (e) Method 4 [6]. (f) Method 5 [39]. (g) Method 6 [40]. (h) Method 7 [11]. (i) Method 8 [12,31]. (j) Method 9 [41]. (k) Method 10 [7]. (l) Method 11 [42]. (m) Method 12 [43]. (n) Method 13 [44]. (o) Method 14 [45]. (p) Method 15 [46]. (q) Method 16 [47]. (r) Method 17.

4.1. Results and discussions To assess the performance of the proposed method, different experiments were conducted on the standard test images and several ultrasound images degraded with the different noises such as additive Gaussian noise and multiplicative speckle noise. The results obtained from the proposed method were compared with the following noise suppression techniques. Method 1: The kuan filter as described in [35]. Method 2: The linear homogeneous mask area filter (LHMAF) as described in [36,37] with 5 × 5 square neighborhood window and 15 iterations. Method 3: The maximum homogeneity over pixel neighborhood filter (MHOPNF) as described in [37,38] with a square neighborhood of 3 × 3. Method 4: Anisotropic diffusion filter as described in [6] with 50 iterations and  = 0.25. Method 5: Fourth order PDE filter as described in [39] with 250 iterations, k = 0.5 and time step t = 0.25.

Method 6: Translational invariant wavelet transform based thresholding approach as described in [40]. Method 7: The ridgelet transform based thresholding approach as described in [11]. Method 8: The curvelet based thresholding with cycle spinning approach as discussed in [12,31]. To perform this approach, scale dependent constant value (k) is chosen as k = 4 for the first scale and k = 3 for the others. Method 9: The shearlet based thresholding with cycle spinning algorithm as discussed in [41]. The thresholding approach in shearlet domain is same as in wavelet or curvelet domain. Method 10: The SRAD (speckle reducing anisotropic diffusion) approach as discussed in [7] with second diffusivity equation, 35 averaged iterations and  = 0.2. Method 11: Nonlinear complex diffusion filter (NCDF) as described in [42] with = ␲/30 and = 20. Method 12: Improved adaptive complex diffusion as described in [43] with  max = 10 s, min = 2, max = 28, and = ␲/30.

D. Gupta et al. / Biomedical Signal Processing and Control 14 (2014) 55–65

61

Fig. 5. Visual comparisons of the various noise suppressing methods on simulated fetus image. (a) Noisy image degraded with speckle noise. (b) Method 1 [35]. (c) Method 2 [36,37]. (d) Method 3 [37,38]. (e) Method 4 [6]. (f) Method 5 [39]. (g) Method 6 [40]. (h) Method 7 [11]. (i) Method 8 [12,31]. (j) Method 9 [41]. (k) Method 10 [7]. (l) Method 11 [42]. (m) Method 12 [43]. (n) Method 13 [44]. (o) Method 14 [45]. (p) Method 15 [46]. (q) Method 16 [47]. (r) Method 17.

Table 1 Performance measures obtained by various noise reduction techniques for the simulated kidney image degraded with speckle noise illustrated in Fig. 4a. Metrics

 = 0.1

Methods

PSNR

SNR

SSIM

EKI

PSNR

SNR

SSIM

EKI

PSNR

SNR

SSIM

EKI

Method 1 [35] Method 2 [36,37] Method 3 [37,38] Method 4 [6] Method 5 [39] Method 6 [40] Method 7 [11] Method 8 [12,31] Method 9 [41] Method 10 [7] Method 11 [42] Method 12 [43] Method 13 [44] Method 14 [45] Method 15 [46] Method 16 [47] Method 17

21.84 19.29 20.49 21.74 23.95 27.48 24.54 26.77 27.54 25.86 25.22 26.84 26.23 26.05 28.12 28.37 28.55

11.22 8.67 9.87 11.11 13.33 16.86 13.93 16.15 16.92 15.24 14.61 16.22 15.61 15.43 17.49 17.75 17.93

0.7596 0.7406 0.7435 0.7610 0.7323 0.7412 0.7375 0.7666 0.7603 0.7197 0.6869 0.7286 0.7512 0.7186 0.7653 0.7829 0.8259

0.3333 0.3235 0.3301 0.3181 0.3734 0.4319 0.4136 0.4340 0.4567 0.4598 0.4277 0.4962 0.3719 0.4685 0.4801 0.4827 0.5276

20.03 18.55 20.11 20.79 23.15 27.31 23.75 26.51 27.28 25.04 24.27 25.13 25.42 25.11 26.99 27.37 27.64

9.41 7.94 9.48 10.17 12.54 16.69 13.13 15.95 16.66 14.42 13.65 14.51 14.81 14.49 16.37 16.75 17.02

0.6770 0.7151 0.7161 0.7190 0.7249 0.7391 0.7325 0.7633 0.7548 0.6944 0.6368 0.7531 0.7593 0.7114 0.7458 0.7797 0.7892

0.3150 0.3177 0.3235 0.3041 0.3293 0.4307 0.3782 0.4318 0.4349 0.3332 0.3139 0.3234 0.3397 0.3284 0.4402 0.4519 0.4619

19.37 18.07 19.19 20.46 23.13 26.91 22.09 26.27 27.05 24.69 23.54 24.99 24.77 24.83 25.74 26.20 27.12

8.75 7.44 8.57 9.84 12.51 16.29 11.47 15.65 16.43 14.07 12.92 14.37 14.14 14.21 15.13 15.58 16.50

0.6103 0.6822 0.6847 0.6974 0.7208 0.7353 0.7146 0.7407 0.7560 0.6873 0.6256 0.7197 0.7264 0.7045 0.7229 0.7612 0.7625

0.3046 0.3177 0.3226 0.3011 0.3291 0.4218 0.3672 0.4193 0.4231 0.3211 0.3024 0.4029 0.3242 0.3387 0.4312 0.4236 0.4538

 = 0.2

Method 13: The squeez box filter (SBF) as described in [44] with 7 × 7 square window and 500 iterations. Method 14: The speckle reducing bilateral filter as discussed in [45] with  d = 1.8 and  r = 2. Method 15: A nonlinear total variation approach (TV) as described in [46].

 = 0.3

Method 16: A nonlocal mean based speckle (NLMS) filter as described in [47]. Method 17: The proposed method based on the modified anisotropic diffusion equations for improving the denoising performance of the multiscale NSST based thresholding.

62

D. Gupta et al. / Biomedical Signal Processing and Control 14 (2014) 55–65

Table 2 Performance measures obtained by various noise reduction techniques for the simulated fetus image with speckle noise illustrated in Fig. 5a. Metrics

 = 0.1

Methods

PSNR

SNR

SSIM

EKI

 = 0.2 PSNR

SNR

SSIM

EKI

 = 0.3 PSNR

SNR

SSIM

EKI

Method 1 [35] Method 2 [35,37] Method 3 [37,38] Method 4 [6] Method 5 [39] Method 6 [40] Method 7 [11] Method 8 [12,31] Method 9 [41] Method 10 [7] Method 11 [42] Method 12 [43] Method 13 [44] Method 14 [45] Method 15 [46] Method 16 [47] Method 17

17.12 15.21 15.73 17.06 17.63 18.31 18.71 18.03 18.80 18.33 17.63 18.69 17.95 18.58 20.08 19.56 20.70

12.13 10.22 10.75 12.08 12.64 13.32 13.69 13.04 13.81 13.34 12.64 13.71 12.97 13.60 15.09 14.57 15.71

0.4722 0.5249 0.5635 0.5925 0.6573 0.6077 0.6001 0.6381 0.6490 0.6286 0.5683 0.6381 0.6267 0.6157 0.6660 0.6800 0.6820

0.3580 0.3334 0.3451 0.3572 0.4603 0.3785 0.5242 0.4277 0.4387 0.3828 0.3337 0.4714 0.4072 0.3758 0.6574 0.6336 0.6792

16.22 12.89 13.29 15.72 16.39 17.40 16.71 17.15 17.75 16.66 16.11 16.85 16.79 16.72 18.96 18.13 19.22

11.23 7.91 8.31 10.73 11.41 12.41 11.72 12.16 12.76 11.67 11.12 11.86 11.80 11.75 13.96 13.14 14.23

0.3945 0.4345 0.5110 0.5375 0.6200 0.5758 0.5678 0.6164 0.5764 0.5512 0.4789 0.5893 0.5948 0.5689 0.6279 0.6369 0.6383

0.2944 0.3120 0.2975 0.3246 0.3540 0.3579 0.4824 0.4089 0.4192 0.3866 0.3498 0.4066 0.3646 0.3951 0.5753 0.5454 0.5897

15.17 10.79 10.82 15.58 15.82 16.50 16.21 16.38 16.82 15.39 15.56 15.94 16.34 16.54 17.59 17.91 18.43

10.18 5.81 5.83 10.59 10.81 11.51 11.22 11.39 11.83 10.40 10.57 10.95 11.36 11.55 12.61 12.92 13.45

0.3536 0.3804 0.4732 0.4978 0.5344 0.5445 0.5412 0.5816 0.5095 0.4445 0.4334 0.5079 0.5296 0.4763 0.5678 0.5752 0.6050

0.2754 0.2952 0.2773 0.3092 0.3188 0.3167 0.3951 0.3730 0.3879 0.3625 0.3131 0.3863 0.3154 0.3609 0.4947 0.4541 0.5474

Fig. 6. Visual comparisons of the various noise suppressing methods on prostate ultrasound image. (a) Original image. (b) Method 1 [35]. (c) Method 2 [36,37]. (d) Method 3 [37,38]. (e) Method 4 [6]. (f) Method 5 [39]. (g) Method 6 [40]. (h) Method 7 [11]. (i) Method 8 [12,31]. (j) Method 9 [41]. (k) Method 10 [7]. (l) Method 11 [42]. (m) Method 12 [43]. (n) Method 13 [44]. (o) Method 14 [45]. (p) Method 15 [46]. (q) Method 16 [47]. (r) Method 17.

D. Gupta et al. / Biomedical Signal Processing and Control 14 (2014) 55–65

63

Fig. 7. Visual comparisons of the various noise suppressing methods on splenic cyst ultrasound image. (a) Original image. (b) Method 1 [35]. (c) Method 2 [36,37]. (d) Method 3 [37,38]. (e) Method 4 [6]. (f) Method 5 [39]. (g) Method 6 [40]. (h) Method 7 [11]. (i) Method 8 [12,31]. (j) Method 9 [41]. (k) Method 10 [7]. (l) Method 11 [42]. (m) Method 12 [43]. (n) Method 13 [44]. (o) Method 14 [45]. (p) Method 15 [46]. (q) Method 16 [47]. (r) Method 17.

Fig. 8. Real ultrasound images with two selected image regions as Regions 1 and 2 for quantitative.

4.1.1. Experiment 1 To perform the experiment, two different sets of simulated images from [48] and a speckle noise model available in MATLAB, are considered. The testing gray scale images of kidney and fetus are corrupted with simulated speckle noise model [37]. Three different noise levels were used by setting the different variance of noise  = (0.1, 0.2, 0.3) to analyze and compare the performance of

the different denoising methods. Noisy speckled kidney and fetus images are shown in Figs. 4a and 5a, respectively. The outcomes of kidney and fetus images processed with various speckle reduction techniques corresponding to  = 0.1 and 0.2 are presented in Figs. 4 and 5, respectively. From these figures, it is observed that the proposed method provides better visual results compared to others. Apart from the visual assessment, Tables 1 and 2

64

D. Gupta et al. / Biomedical Signal Processing and Control 14 (2014) 55–65

20

25

15

20 15

10

10

5

(b)

MVR1 MVR2

MVR1 MVR2

30

25

25

20

20

15

15

10

10

5

0

0 Noisy Method 1 Method 2 Method 3 Method 4 Method 5 Method 6 Method 7 Method 8 Method 9 Method 10 Method 11 Method 12 Method 13 Method 14 Method 15 Method 16 Method 17

5

(c)

Noisy Method 1 Method 2 Method 3 Method 4 Method 5 Method 6 Method 7 Method 8 Method 9 Method 10 Method 11 Method 12 Method 13 Method 14 Method 15 Method 16 Method 17

0

MVR1 MVR2

(d)

Noisy Method 1 Method 2 Method 3 Method 4 Method 5 Method 6 Method 7 Method 8 Method 9 Method 10 Method 11 Method 12 Method 13 Method 14 Method 15 Method 16 Method 17

(a)

Noisy Method 1 Method 2 Method 3 Method 4 Method 5 Method 6 Method 7 Method 8 Method 9 Method 10 Method 11 Method 12 Method 13 Method 14 Method 15 Method 16 Method 17

0

5

MVR1 MVR2

Fig. 9. Plot of MVR values obtained by the different denoising methods for two image regions shown in Fig. 8a–d.

summarize all the performance indices evaluated by all the 17 different techniques at three different noise levels for both the images, respectively. From the objective results listed in Tables 1 and 2, it is observed for low level of noise that the proposed method produce the competitive speckle reduction performance as compared to other transform based methods, TV and NLMS filter by achieving the good SNR and PSNR values with better structure and edge preservation indicated by higher SSIM and EKI values, respectively. The proposed method also outperforms the other spatial domain filters. In case of high levels noise, the proposed method provides the competitive edge preservation results compared to methods 6–17 with better noise suppression performance. Moreover, the proposed method suppresses the speckle noise of different levels while preserving the more edges and fine details that is also supported with the quantitative results. 4.1.2. Experiment 2 In this section, different real ultrasound images acquired from open source medical image database1 2 are used to evaluate the performance of the proposed method. The denoising methods were tested on 50 different ultrasound images. The outcomes of the experiments on different ultrasound images such as prostate and splenic cyst are shown in Figs. 6 and 7, respectively. It is very difficult to analyze the denoising results of real ultrasound images processed with the different methods in terms of different performance measures since there are not any noise free reference images. Mean to variance ratio (MVR) is an index for estimating the speckle noise level in ultrasound images over different image regions [37,49]. A larger value of MVR represents a better quantitative performance of different despeckling methods on real

1 2

http://www.ultrasoundcases.info/. http://rad.usuhs.edu/medpix/parent.php3?mode=home page#top.

Table 3 Performance comparisons of the various denoising methods using the averaged values of MVR over 100 different regions obtained on 50 different ultrasound images. Methods

MVR

Method 1 [35] Method 2 [36,37] Method 3 [37,38] Method 4 [6] Method 5 [39] Method 6 [40] Method 7 [11] Method 8 [12,31] Method 9 [41] Method 10 [7] Method 11 [42] Method 12 [43] Method 13 [44] Method 14 [45] Method 15 [46] Method 16 [47] Method 17

16.56 15.19 16.70 18.26 19.27 18.16 17.47 18.22 20.94 19.45 19.32 21.05 20.12 19.69 21.45 21.51 21.95

± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ±

4.98 3.46 3.49 4.23 4.46 3.14 3.69 4.57 3.88 3.97 3.69 3.15 4.87 3.23 4.95 4.13 4.26

ultrasound medical images. Fig. 8 shows the four different ultrasound images in which two different image regions are marked with red and blue rectangles for the quantitative analysis. Red and blue rectangles represent as region 1 and region 2, respectively. The MVR values obtained for both the regions 1 and 2 as MVR1 and MVR2, respectively, are shown in Fig. 9a–d, which show that the proposed method outperforms the others exhibited in terms of larger values of MVR. In addition to further assess the performance of the proposed method, 100 measurements are taken on 50 different ultrasound images (two measurements at different location for each image) to evaluate the MVR values. Table 3 shows the averaged MVR values obtained for all the seventeen different noise reduction techniques.

D. Gupta et al. / Biomedical Signal Processing and Control 14 (2014) 55–65

The results of Table 3 show the superiority of the proposed method in effective noise suppression as compared to others. 5. Conclusions In this paper, a despeckling method for the US images is presented by combining the NSST thresholding approach with nonlinear modified anisotropic diffusion equations. The proposed method is adapted to the multiplicative speckle noise. In the proposed method, the NSST has several advantages over the other platforms. The NSST also provides the multiscale and direction analysis of the noisy images. The large-amplitude noise components are suppressed by applying diffusion process on lowfrequency approximation coefficients and thresholding provides the modified coefficients, which improves the denoising efficiency with better edge/texture preservation. We have also compared the despeckling capabilities of the proposed method with other techniques using a number of test images corrupted with simulated speckle noise and demonstrated better denoising performance of the proposed method with more edge preservation. Besides, the proposed and other methods are also tested on a large database of the real ultrasound images. Finally, it is concluded that the proposed method ensures an improvement in effective noise suppression and preservation of more edges, thus providing the despeckled images with good visual appearance. References [1] A.P. Dhawan, Medical Image Analysis, John Wiley and Sons, Hoboken, NJ, 2003. [2] D. Mittal, V. Kumar, S. Saxena, N. Khandelwal, N. Kalra, Enhancement of the ultrasound images by modified anisotropic diffusion method, Med. Biol. Eng. Comput. 48 (2010) 1281–1291. [3] W.K. Pratt, Digital Image Processing, John Wiley and Sons, New York, 2006. [4] R.C. Gonzalez, R.E. Woods, Digital Image Processing, Prentice-Hall, India, 2001. [5] 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 (1989) 129–135. [6] P. Perona, J. Malik, Scale-space and edge detection using anisotropic diffusion, IEEE Trans. Pattern Anal. Mach. Intell. 12 (1990) 629–639. [7] Y. Yongjian, S.T. Acton, Speckle reducing anisotropic diffusion, IEEE Trans. Image Process. 11 (2002) 1260–1270. [8] X. Liu, J. Liu, X. Xu, L. Chun, J. Tang, Y. Deng, A robust detail preserving anisotropic diffusion for speckle reduction in ultrasound images, BMC Genomics 12 (2011) S14. [9] Z. Fan, Y. Yang Mo, M. Koh Liang, K. Yongmin, Nonlinear diffusion in Laplacian pyramid domain for ultrasonic speckle reduction, IEEE Trans. Med. Imag. 26 (2007) 200–211. [10] D.L. Donoho, I.M. Johnstone, Adapting to unknown smoothness via wavelet shrinkage, J. Am. Stat. Assoc. 90 (1995) 1200–1224. [11] M.N. Do, M. Vetterli, The finite ridgelet transform for image representation, IEEE Trans. Image Process. 12 (2003) 16–28. [12] J.L. Starck, E.J. Candes, D.L. Donoho, The curvelet transform for image denoising, IEEE Trans. Image Process. 11 (2002) 670–684. [13] J.-S. Lee, Digital image enhancement and noise filtering by use of local statistics, IEEE Trans. Pattern Anal. Mach. Intell. 2 (1980) 165–168. [14] D.T. Kuan, A.A. Sawchuk, T.C. Strand, P. Chavel, Adaptive noise smoothing filter for images with signal-dependent noise, IEEE Trans. Pattern Anal. Mach. Intell. PAMI-7 (1985) 165–177. [15] M.J. Black, G. Sapiro, D.H. Marimont, D. Heeger, Robust anisotropic diffusion, IEEE Trans. Image Process. 7 (1998) 421–432. [16] C. Tauber, H. Batatia, A. Ayache, A robust speckle reducing anisotropic diffusion, in: Proceedings of International Conference on Image Processing, 2004, pp. 247–250. [17] S.-M. Chao, D.-M. Tsai, An improved anisotropic diffusion model for detail- and edge-preserving smoothing, Pattern Recogn. Lett. 31 (2010) 2012–2023. [18] D.L. Donoho, De-noising by soft-thresholding, IEEE Trans. Inf. Theory 41 (1995) 613–627. [19] S. Gupta, R. Chauhan, S. Sexana, Wavelet-based statistical approach for speckle reduction in medical ultrasound images, Med. Biol. Eng. Comput. 42 (2004) 189–192.

65

[20] A. Thakur, R.S. Anand, Image quality based comparative evaluation of wavelet filters in ultrasound speckle reduction, Digital Signal Process. 15 (5) (2005) 455–465. [21] C.S. Anand, J.S. Sahambi, Wavelet domain non-linear filtering for MRI denoising, Magn. Reson. Imaging 28 (2010) 842–861. [22] S. Roy, N. Sinha, A. Sen, Fuzzy soft thresholding based hybrid denoising model, in: D. Nagamalai, E. Renault, M. Dhanuskodi (Eds.), Advances in Digital Image Processing and Information Technology, Springer, Berlin Heidelberg, 2011, pp. 1–10. [23] M. Nikpour, H. Hassanpour, Using diffusion equations for improving performance of wavelet-based image denoising techniques, IET Image Process. 4 (2010) 452–462. [24] M.N. Do, M. Vetterli, The contourlet transform: an efficient directional multiresolution image representation, IEEE, Trans. Image Process 14 (2005) 2091–2106. [25] D. Labate, W.-Q. Lim, G. Kutyniok, G. Weiss, Sparse multidimensional representation using shearlets, in: Proceedings of SPIE 5914, 2005, pp. 254–262. [26] K. Guo, D. Labate, Optimally sparse multidimensional representation using shearlets, SIAM J. Math. Anal. 39 (2007) 298–318. [27] G. Easley, D. Labate, W.-Q. Lim, Sparse directional image representations using the discrete shearlet transform, Appl. Comput. Harmon. Anal. 25 (2008) 25–46. [28] B. Hou, X. Zhang, X. Bu, H. Feng, SAR image despeckling based on nonsubsampled shearlet transform, IEEE J. Sel. Topics Appl. Earth Observ. 5 (2012) 809–823. [29] M. Jianwei, G. Plonka, Combined curvelet shrinkage and nonlinear anisotropic diffusion, IEEE Trans. Image Process. 16 (2007) 2198–2206. [30] F. Catté, P. Lions, J. Morel, T. Coll, Image selective smoothing and edge detection by nonlinear diffusion, SIAM J. Num. Anal. 29 (1992) 182–193. [31] W. Aili, Z. Ye, M. Shaoliang, Y. Mingji, Image denoising method based on curvelet transform, in: Proceedings of IEEE Conference on Industrial Electronics and Applications, 2008, pp. 571–574. [32] D. Gupta, R.S. Anand, B. Tyagi, Edge preserved enhancement of medical images using adaptive fusion-based denoising by shearlet transform and total variation algorithm, J. Electron. Imaging 22 (2013) 1–18. [33] M.R. Hajiaboli, An anisotropic fourth-order diffusion filter for image noise removal, Int. J. Comput. Vis. 92 (2011) 177–191. [34] Z. Wang, A.C. Bovik, H.R. Sheikh, E.P. Simoncelli, Image quality assessment: from error visibility to structural similarity, IEEE Trans. Image Process. 13 (2004) 600–612. [35] D.T. Kuan, A. Sawchuk, T.C. Strand, P. Chavel, Adaptive restoration of images with speckle, IEEE Trans. Acoust. Speech Signal Process. 35 (1987) 373–383. [36] M. Nagao, T. Matsuyama, Edge preserving smoothing, Comput. Graph. Image Process. 9 (1979) 394–407. [37] J. Seabra, J. Sanches, Ultrasound speckle/despeckle image decomposition for tissue analysis, in: J.M. Sanches, A.F. Laine, J.S. Suri (Eds.), Ultrasound Imaging, Springer, US, 2012, pp. 73–95. [38] S.M. Ali, R.E. Burge, New automatic techniques for smoothing and segmenting SAR images, Signal Process. 14 (1988) 335–346. [39] Y.L. You, M. Kaveh, Fourth-order partial differential equations for noise removal, IEEE Trans. Image Process. 9 (2000) 1723–1730. [40] S. Qin, C. Yang, B. Tang, S. Tan, The denoise based on translation invariance wavelet transform and its applications, in: Proceedings of Conference on Structural Dynamics, 2002, pp. 783–787. [41] D. Gupta, R.S. Anand, B. Tyagi, Enhancement of medical ultrasound images using multiscale discrete shearlet transform based thresholding, in: Proceedings of IEEE International Symposium on Electronic System Design (ISED), 2012, pp. 286–290. [42] G. Gilboa, N. Sochen, Y.Y. Zeevi, Image enhancement and denoising by complex diffusion processes, IEEE Trans. Pattern Anal. Mach. Intell. 26 (2004) 1020–1036. [43] R. Bernardes, C. Maduro, P. Serranho, A. Araújo, S. Barbeiro, J. Cunha-Vaz, Improved adaptive complex diffusion despeckling filter, Opt. Express 18 (2010) 24048–24059. [44] P.C. Tay, S.T. Acton, J.A. Hossack, Ultrasound despeckling using an adaptive window stochastic approach, in: Proceedings of IEEE International Conference on Image Processing, 2006, pp. 2549–2552. [45] J. Tang, S. Guo, Q. Sun, Y. Deng, D. Zhou, Speckle reducing bilateral filter for cattle follicle segmentation, BMC Genomics 11 (2010) S9. [46] A. Chambolle, An algorithm for total variation minimization and applications, J. Math. Imaging Vis. 20 (2004) 89–97. [47] P. Coupe, P. Hellier, C. Kervrann, C. Barillot, Nonlocal means-based speckle filtering for ultrasound images, IEEE Trans. Image Process. 18 (2009) 2221–2229. [48] Field II ultrasound simulation program, 2013, Available from: http://field-ii.dk/examples (accessed 10.03.13). [49] Q. Huang, M. Lu, Y. Zheng, Z. Chi, Speckle suppression and contrast enhancement in reconstruction of freehand 3D ultrasound images using an adaptive distance-weighted method, Appl. Acoust. 70 (2009) 21–30.