Open access peer-reviewed chapter

Speckle Noise Reduction in Medical Ultrasound Images

By P.S. Hiremath, Prema T. Akkasaligar and Sharan Badiger

Submitted: May 2nd 2012Reviewed: April 8th 2013Published: June 5th 2013

DOI: 10.5772/56519

Downloaded: 3476

1. Introduction

The use of ultrasound imaging in medical diagnosis is well established because of its non-invasive nature, low cost, capability of forming real time imaging and continuing improvement in image quality. However, it suffers from a number of shortcomings and these include: acquisition noise from the equipment, ambient noise from the environment, the presence of background tissue, other organs and anatomical influences such as body fat, and breathing motion. Therefore, noise reduction is very important, as various types of noise generated limits the effectiveness of medical image diagnosis.

Ultrasound is a sound wave with a frequency that exceeds 20 kHz. It transports energy and propagates through several means as a pulsating pressure wave [1]. It is described by a number of wave parameters such as pressure density, propagation direction, and particle displacement. If the particle displacement is parallel to the propagation direction, then the wave is called a longitudinal or compression wave. If the particle displacement is perpendicular to the propagation direction, it is a shear or transverse wave. The interaction of ultrasound waves with tissue is subject to the laws of geometrical optics. It includes reflection, refraction, scattering, diffraction, interference, and absorption. Except from interference, all other interactions reduce the intensity of the ultrasound beam.

Ultrasound technique is mainly based on measuring the echoes transmitted back from a medium when sending an ultrasound wave to it. In the echo impulse ultrasound technique, the ultrasound wave interacts with tissue and some of the transmitted energy returns to the transducer to be detected by the instrument [2]. Further, the reflected waves are picked up by the transducer probe and relayed to the machine. The machine calculates the distance from the transducer probe to the tissue or organ (boundaries) using the speed of sound in tissue (1,540 m/s) and the time of the each echo's return ( millionths of a second). The machine displays the distances and intensities of the echoes on the screen, forming a two dimensional image. Superficial structures such as muscles, tendons, testes, breast and the neonatal brain are imaged at a higher frequency (7- 18 MHz), which provides better axial and lateral resolution. Deeper structures such as liver and kidney are imaged at a lower frequency 1-6 MHz with lower axial and lateral resolution but greater penetration.

The usefulness of ultrasound imaging is degraded by the presence of signal dependent noise known as speckle. Speckle noise is multiplicative in nature. This type of noise is an inherent property of medical ultrasound imaging and because of this noise the image resolution and contrast become reduced, which effects the diagnostic value of this imaging modality [3]. So, speckle noise reduction is an essential pre processing step, whenever ultrasound imaging is used for medical imaging. Therefore, image despeckling is a very important task, and should be filtered out [4-6], without affecting important features of the image.

In ultrasound images, the noise content is multiplicative and non Gaussian. Such noise is generally more difficult to remove than additive noise, because the intensity of the noise varies with the image intensity. A model of multiplicative noise is given by

yij= XijniE1

where the speckle image yij is the product of the original image Xij, and the non-Gaussian noise nij. The indices i, j represent the spatial position over the image. In most applications involving multiplicative noise, the noise content is assumed to be stationary with unitary mean and unknown noise variance σ2. To convert multiplicative noise into an additive noise, as given in the Eq.(2), a logarithmic transformation is applied to the speckle image yij [7]. The noise component nij is then approximated as an additive zero mean gaussian noise.

ln yij= ln Xij+ ln nijE2

The Discrete Wavelet Transform (DWT) is then applied to ln yij and the wavelet transformed image is subjected to thresholding. After applying the inverse DWT, the processed image is subjected to an exponential transformation, which is the inverse logarithmic operation,that yields a denoised image.

1.1. Image quality assessment

Image quality is important when evaluating or segmenting ultrasound images, where speckle obscures subtle details in the image [8]. In a recent study [9] it is shown that speckle reduction improves the visual perception of the expert in the assessment of ultrasound imaging of the human organs. The statistical parameters like Signal to Noise Ratio (SNR), Correlation Coefficient (CC), variance, Mean Square Error (MSE) and Peak Signal to Noise Ratio (PSNR) for image quality assessment are described below. The following metrics are calculated using the original image X and the despeckled image Y.

Variance: It determines the average dispersion of the speckle in the image. A lower variance gives a cleaner image as more speckles are reduced. The formula for calculating the variance is

σ2=1N2i,j=0N1(XijX¯)2E3

where X-is the mean intensity value of the image X of size N×N.

Mean Square Error: The MSE measures the quality change between the original image (X) and denoised image (Y) of size N×N. It is given by

MSE=1N2i,j=0N1(XijYij)2E4

The MSE has been widely used to quantify image quality and, when used alone, it does not correlate strongly enough with perceptual quality. It should be used, therefore, together with other quality metrics and visual perception.

Signal to Noise Ratio: The SNR compares the level of desired signal to the level of background noise. The higher the ratio, the less obtrusive the background noise is. It is expressed in decibels (dB) as

SNR=10log10(σ2σe2)E5

where, σ2is the variance of the original image and σe2is the variance of error (Difference between the original and denoised image i.e. X-Y).

Peak Signal-to-Noise Ratio: The PSNR is computed as

PSNR=10 log10(S2MSE)E6

where S is the maximum intensity in the original image. The PSNR is higher for a good quality image and lower for a poor quality image. It measures image fidelity, that is, how closely the transformed image resembles the original image.

Correlation Coefficient: It represents the strength and direction of a linear relationship between two variates. The best known is the Pearson product moment correlation coefficient, which is obtained by dividing the covariance of the two variables by the product of their standard deviations, and it is given by

CC=N2XijYijXijYijN2Xij2-(Xij)2N2Yij2(Yij)2E7

where the summations are done over both the indices i and j from 0 to N-1. If the correlation coefficient is near to +1, then there exists stronger positive correlation between the original image and despeckled image.

Computational time: The computational time (Tc) of a filter is defined as the time taken by a digital computing platform to execute the filtering algorithm when no other software, except the operating system, runs on it. Normally, Tc depends on the computing system’s clock time period. But, in addition to the clock period, it depends on the memory size, the input data size, and the memory access time, etc. The computational time taken by a filter should be low for online and real time image processing applications. Hence, a filter with lower Tc is better than a filter having higher Tc value while all other performance measures remain identical.

1.2. Image data set

For the purpose of experimentation, a medical ultrasound image database is prepared in consultation with a medical expert for the present study. The images are acquired using the instrument GE LOGIQ 3 Expert system with 5 MHz transducer frequency, in JPEG format. The data set consists of 70 ultrasound images of size 512x512, of kidney and liver.

In the present Chapter, the aim of the study is to address, novel issues related to despeckling medical ultrasound images. It is envisaged that the results of this investigation would be used as a pre processing step for effective image segmentation or image registration techniques in other applications also.

The remaining part of this Chapter is organised into five sections. The section 2 deals with common speckle filters. The section 3 examines wavelet transform methods, while the section 4 investigates contourlet transform methods. The section 5 describes the Gaussian model of speckle noise. Finally, the section 6 gives the conclusions.

2. Common speckle filters

Speckle reducing filters are originated from the synthetic aperture radar community [10]. Later these filters are applied to ultrasound imaging since the early 1980s [11]. There are two major classifications of speckle reduction filters, viz. single scale spatial filters and transform domain multiscale filters. The spatial filter acts on an image by smoothing it; that is, it reduces the intensity variation between adjacent pixels. The simple sliding window spatial filter replaces the center value in the window with the average of all the neighboring pixel values including itself. By doing this, it replaces pixels, that are unrepresentative of their surroundings. It is implemented with a convolution mask, which provides a result that is a weighted sum of the values of a pixel and its neighbors. It is also called a linear filter. The mask or kernel is a square. Often a 3× 3 square kernel is used. If the coefficients of the mask sum up to one, then the average brightness of the image is not changed. If the coefficients sum to zero, the average brightness is lost, and it returns a dark image.

The common speckle filters such as Lee, Kuan and Wiener filters are considered for the study. The brief definition and mathematical description of the standard spatial filters are discussed below:

Lee filter:

The Lee filter [12] is based on the approach that the smoothing is performed on the area having low variance. However, smoothing will not be performed on area of high variance, which is near edges. The Lee filter assumes that the image can be approximated by a linear model represented by Eq.(8)

Yij=K-+W*(C-K-)E8

where Yij is the gray scale value of the pixel at ( i, j) after filtering. If there is no smoothing, the filter will output, only the mean intensity value K-of the kernel K, otherwise, the difference between the centre pixel C and K-is calculated and multiplied with a weighting function W given in Eq.(9) :

W=σk2(σk2+σ2)E9

and then summed with K-, where σk2 is the variance of the pixel values within the kernel given by the Eq. (10):

σk2=1M2u,v=0M-1 Kuv-K- 2E10

where MxM is the size of the kernel and Kuv is the pixel value within the kernel at indices u and v, K-is the mean intensity value of kernel. The parameter σ2is the variance of the image X, which is given by the Eq.(3). The main disadvantage of Lee filter is that it tends to ignore speckle noise in the areas closest to edges and lines.

Kuan filter:

The Kuan filter [13] is a generalization of the Lee filter. The Kuan filter converts the multiplicative model of speckle into an additive linear form. However, it is based on the equivalent number of looks (ENL), which is computed from an ultrasound image to determine a different weighting function W given by the Eq.(11):

W=1-Cu/Ci1+Cu.E11

The weighting function is computed from estimated noise variation coefficient of the image, Cu given by the Eq.(12):

Cu=ENL-12E12

and the variation coefficient Ci of the image given by the Eq. (13):

Ci=σk/K-E13

where ENL is given by the Eq.(14) :

ENL=(K¯σk)2E14

In the Eq.(14), the σkis the standard deviation of the kernel and K- is the mean intensity value of the kernel. The only disadvantage of the Kuan filter is that the ENL parameter needs to be computed.

Wiener filter:

The Wiener filter [14] is a linear spatial domain filter. There are two alternatives : (i) Fourier transform method ( frequency domain) (ii) mean squared method (spatial domain), for implementing the Wiener filter. The first alternative is used for denoising and deblurring, whereas the second alternative is used for denoising only. The frequency domain alternative of Wiener filtering requires a prior knowledge of noise power spectra and the original image. But, in the spatial domain alternative, no such prior knowledge is required. It is based on statistical least squared principle and minimizes the mean squared error between actual signal sequence and desired signal sequence.

In an image, the statistical properties differ too much from one region to another region. Thus, both global statistics (mean, variance, and higher order moments of entire image) and local statistics (mean, variance, and higher order moments of kernel) are important. Wiener filtering is based on both, global and local statistics and is given by

Yij=K-+σk2σk 2+σ2 Kuv-K-E15

where Yij denotes the despeckled image,  K- is the local mean, σk2 is the local variance,  Kuvis u,vthpixel in the kernel K and σ2is the global variance. Let us consider kernel of size MxM, then local variance  σk 2is defined by Eq.(10). From the Eq.(15), it is observed that the filter output is equal to local mean if the centre pixel value equals local mean, or else it outputs the modified value different from local mean. Thus, filter output varies from the local mean depending upon the local variance and thus tries to hold the true original value as far as possible.

The Lee filter and Wiener filter are implemented using kernel size 3x3, 5x5, 7x7 and Kuan filter using kernel size 3x3 and 5x5.The classical Wiener filter, is not adequate for removing speckle, since it is designed mainly for additive noise suppression. To address the multiplicative nature of speckle noise, a homomorphic approach is developed in [15], which converts the multiplicative noise into additive noise, by taking the logarithm of image and then applies the Wiener filter. The PSNR, SNR, CC, variance and MSE are considered as filter performance measures. The Figures 1.-4. show the average results obtained for 70 ultrasound images, which are despeckled using Kuan, Lee and Wiener filter. The optimality is determined by the criteria, namely (i) higher SNR and PSNR values, (ii) lower variance, MSE values and (iii) Correlation Coefficient is nearly equal to one. From the Figures 1.-4., it is observed that Wiener filter with kernel size 3x3 gives better results than other despeckling filters. The computational time of different filters are given in the Table 1. The filter having less computational time is usually required for online and real time applications. The least value of computation is highlighted. From the Table 1, it is observed that Wiener filter with kernel size 3x3 is better among all the filters compared here, for despeckling medical ultrasound images.

For proper judgement of performance of filters, the subjective evaluation should be taken into consideration. For subjective evaluation, the despeckled images of various filters are shown in the Figure 5. From the Figure 5, it is observed from visual inspection that all the three methods achieved good speckle suppression performance. However, Lee and Kuan filters lost many of the signal details and the resulting images are blurred. Further, Wiener filter with kernel size 3x3 yielded better visual enhancement of medical ultrasound images. However, the Lee filter smoothes away noise in flat regions, but leaves the fine details such as lines and texture unchanged.

Figure 1.

Performance of various despeckling filters, in terms of PSNR, SNR.

Figure 2.

Performance of various despeckling filters, in terms of variance

Figure 3.

Performance of various despeckling filters, in terms of MSE.

Figure 4.

Performance of various despeckling filters, in terms of Correlation Coefficient

Figure 5.

Performance comparison of various despeckling filters by visual inspection of an ultrasound image of kidney.

Despeckling filtersComputational Time (in Secs.)
Lee (3x3)25.60
Lee (5x5)25.63
Lee (7x7)25.64
Kuan(3x3)40.89
Kuan(5x5)41.57
Wiener(3x3)4.14
Wiener(5x5)5.14
Wiener(7x7)6.15

Table 1.

Performance comparison of various despeckling filters based on computational time.

Thus the main disadvantage of Lee filter is that, it tends to ignore speckle noise in the area closest to edges and lines. The Kuan filter is considered to be more superior to the Lee filter. It does not make an approximation on the noise variance within the filter window. The only limitation of Kuan filter is the high computational time due to estimation of ENL parameter. The Wiener filter with kernel size 3×3 is effective in preserving the edges and other detailed information upto some extent. Further, when the various spatial domain filters are compared by visual inspection, it is observed that Wiener filter with kernel size 3×3 yielded better visual enhancement of medical ultrasound images. Further, for the complete removal of speckle without losing any data is not possible at the moment. This is because all of these filters rely on local statistical data related to the filtered pixel. An alternative approach is to use wavelet transform.

3. Wavelet transform method

The primary goal of speckle reduction is to remove the speckle without losing much detail contained in an image. To achieve this goal, we make use of the wavelet transform and apply multiresolution analysis to localize an image into different frequency components or useful subbands and then effectively reduce the speckle in the subbands according to the local statistics within the bands. The main advantage of the wavelet transform is that the image fidelity after reconstruction is visually lossless.

A wavelet is a mathematical function used to decompose a given function or continuous-time signal into different frequency components and study each component with a resolution that matches its scale. A wavelet transform is the representation of a function by wavelets. The wavelets are scaled and translated copies (known as daughter wavelets) of a finite length or fast decaying oscillating waveform (known as mother wavelet). Wavelet transforms are classified into continuous wavelet transform (CWT) and discrete wavelet transform. The CWT analyzes the signal through the continuous shifts of a scalable function over a time plane. Because of computers discrete nature, computer programs use the discrete wavelet transform. The discrete transform is very efficient from computational point of view.

Image denoising using wavelet techniques is effective because of its ability to capture most of the energy of a signal in a few significant transform coefficients. Another reason of using wavelet transform is due to development of efficient algorithms for signal decomposition and reconstruction [16] for image processing applications such as denoising and compression. A survey of despeckling techniques is discussed in [17, 18] and many wavelet domain techniques are already available in the literature. In [19], the authors have presented a novel speckle suppression method for medical ultrasound images, in which it is shown that the subband decompositions of ultrasound images have significantly non Gaussian statistics that are best described by families of heavy tailed distributions such as the alpha stable. Then, a Bayesian estimator is designed to exploit these statistics. Alpha stable model is used to develop a blind noise removal processor that performs a nonlinear operation on the data. In [20], the authors have proposed a novel technique for despeckling the medical ultrasound images using lossy compression. In [21], authors have proposed a new wavelet based image denoising technique, in which the different threshold functions, namely universal threshold, Visu shrink, sure shrink, Bayes shrink and normal shrink are considered for the study. The threshold value is calculated using circular kernel, mean max threshold, nearest neighbour and new threshold function.

Any decomposition of an image into wavelets involves a pair of waveforms, one to represent the high frequencies corresponding to the detailed parts of an image (wavelet function ψ ) and one for high frequencies are transformed with short functions (low scale). The result of WT is a set of wavelet the low frequencies or smooth parts of an image (scaling function ϕ) coefficients, which measure the contribution of the wavelets at different locations and scales. The WT performs multiresolution image analysis [22]. The scaling function for multiresolution approximation can be obtained as the solution to a two scale dilatational Eq.(16):

ϕ(x)=kaL(k)ϕ(2x-k)E16

for some suitable sequence of coefficients aL(k). Once ϕhas been found, an associated mother wavelet is given by a similar looking Eq.(17):

ψ(x)=kaH(k)ϕ(2x-k).E17

Wavelet analysis leads to perfect reconstruction filter banks using the coefficient sequences aL(k) and aH(k). The input sequence X is convolved with high pass (HPF) and low pass (LPF) filters aHk, aLkrespectively. Further, each result is down sampled by two, yielding the transform signals xHand  xL. The signal is reconstructed through upsampling and convolution with high and low synthesis filters sH (k) and sL (k). By cascading the analysis filter bank with itself a number of times, digital signal decomposition with dyadic frequency scaling known as DWT can be formed. The DWT for an image as a 2D signal can be derived from 1D DWT. The easiest way for obtaining scaling and wavelet functions for two dimensions is by multiplying two 1D functions. The scaling function for 2D DWT can be obtained by multiplying two 1D scaling functions: ϕ(x,y)=ϕ(x)ϕ(y)representing the approximation subband image (LL). The analysis filter bank for a single level 2D DWT structure produces three detail subband images (HL, LH, HH) corresponding to three different directional orientations (Horizontal, Vertical and Diagonal) and a lower resolution subband image LL. The filter bank structure can be iterated in a similar manner on the LL channel to provide multilevel decomposition. The separable wavelets are also viewed as tensor products of one dimensional wavelets and scaling functions. If ψ(x) is the one dimensional wavelet associated with one dimensional scaling function ϕ(x), then three 2D wavelets associated with three subband images, called as vertical, horizontal and diagonal details, are given by

ψV(x,y)=φ(x)ψ(y)E18
ψH(x,y)=ψ(x)φ(y)E19
ψD(x,y)=ψ(x)ψ(y)E20

which correspond to the three subbands LH, HL and HH, respectively [23]. The wavelet equation produces different types of wavelet families like Daubenchies, Haar, Symlets, Coiflets and Biorthogonal wavelets [24].

3.1. Thresholding techniques

There are two approaches to perform the thresholding after computation of the wavelet coefficients, namely, subband thresholding and global thresholding [25]. In subband thresholding, we compute the noise variance of the horizontal, vertical and diagonal sub bands of each decomposition level, starting from the outer spectral bands and moving towards inner spectral bands (decomposition from higher levels towards lower levels) and calculate threshold value using Bayes shrinkage or Visu shrinkage rule. In global thresholding, we determine the threshold value from only the diagonal band but we apply this threshold to the horizontal, vertical and diagonal sub bands. This approach assumes that the diagonal band contains most of the high frequencies components; hence the noise content in diagonal band should be higher than the other bands. Thresholding at the coarsest level is not done, because it contains the approximation coefficients that represent the translated and scaled down version of the original image. Thresholding at this level will cause the reconstruction image to be distorted.

Shrinkage scheme:

The thresholding approach is to shrink the detail coefficients (high frequency components) whose amplitudes are smaller than a certain statistical threshold value to zero while retaining the smoother detail coefficients to reconstruct the ideal image without much loss in its details. This process is sometimes called wavelet shrinkage, since the detail coefficients are shrunk towards zero. There are three schemes to shrink the wavelet coefficients, namely, the keep-or-kill hard thresholding, shrink-or-kill soft thresholding introduced by [26] and the recent semi soft or firm thresholding. Shrinking of the wavelet coefficient is most efficient if the coefficients are sparse, that is, the majority of the coefficients are zero and a minority of coefficients with greater magnitude that can represent the image [27]. The criterion of each scheme is described as follows. Given that λ denotes the threshold limit, Xw denotes the input wavelet coefficients and Yt denotes the output wavelet coefficients after thresholding, we define the following thresholding functions:

Hard thresholding:

Yt=Thard(Xw)     ={Xw     , for |Xw|λ0         , for |Xw|<λ  E21

Soft thresholding:

Yt=Tsoft(Xw)    ={sign{Xw}(|Xw|λ)   , for  |Xw|λ0                                  , for  |Xw|<λ  E22

Semi soft thresholding:

Yt=Tsemisoft(Xw) ={0                                               , for         |Xw|λ sign{Xw}λ1( |Xw|λ)λ1λ           , for   λ< |Xw|λ1Xw                                          , for         |Xw|>λ1           where   λ1=2λ.E23

The hard thresholding procedure removes the noise by thresholding only the wavelet coefficients of the detail sub bands, while keeping the low resolution coefficients unaltered. The soft thresholding scheme shown in Eq. (22) is an extension of the hard thresholding. It avoids discontinuities and is, therefore, more stable than hard thresholding. In practice, soft thresholding is more popular than hard thresholding, because it reduces the abrupt sharp changes that occurs in hard thresholding and provides more visually pleasant recovered images. The aim of semi soft threshold is to offer a compromise between hard and soft thresholding by changing the gradient of the slope. This scheme requires two thresholds, a lower threshold λ and an upper threshold λ1 where λ1 is estimated to be twice the value of lower threshold λ.

3.2. Shrinkage rule

A very large threshold λ will shrink almost all the coefficients to zero and may result in over smoothing the image, while a small value of λ will lead to the sharp edges with details being retained but may fail to suppress the speckle. We use the shrinkage rules, namely, the Visu shrinkage rule and Bayes shrinkage rule for thresholding which are explained in the following:

3.2.1. Visu shrinkage rule

Visu shrinkage rule [28] is thresholding by applying universal threshold. The idea is to find

each threshold λi to be proportional to the square root of the local noise variance σ2 in each subband of the ultrasound image after decomposition. If Nk is the size of the subband in the wavelet domain, then λi

λi=σ2log(Nk)E24

The estimated local noise variance, σ2, in each subband is obtained by averaging the squares of the empirical wavelet coefficients at the highest resolution scale as

σ2=1Nk2i,j=0Nk1Wij2E25

The threshold of Eq.(24) is based on the fact that, for a zero mean independent identically distributed (i.i.d.) Gaussian process with variance σ2, there is a high probability that a sample value of this process will not exceed λ. Thus, the Visu shrink is suitable for applications with white Gaussian noise and in which most of the coefficients are zero. In such cases, there is a high probability that the combination of (zero) coefficients plus noise will not exceed the threshold level λ.

3.2.2. Bayes shrinkage rule

In Bayes shrink [29], an adaptive data driven threshold is used for image denoising. The threshold on a given subband Wk of the image X, is given by

λK=σn2σxE26

where  σn, the estimated noise variance found as the median of the absolute deviation of the diagonal detail coefficients on the finest level (sub band HH1), is given by

σn=median({WijHH1})0.67452E27

This estimator is used when there is no a priori knowledge about the noise variance. The σx, which is the estimated signal variance in the wavelet domain, is given by

σx=max(σy2n2,0)E28

and  σy2, an estimate of the variance of the Wk, Since Wk is modeled as zero mean, σy2 can be found empirically by,

σy2=1Nk2i,j=0Nk1Wij2E29

in which Nk is the number of the wavelet coefficients Wk on the subband considered. In the Eq. (27), the value 0.67452 is the median absolute deviation of normal distribution with zero mean and unit variance. In the Eq.(26), if (σnx) << 1, the signal is much stronger than the noise. The normalized threshold is chosen to be small in order to preserve most of the signal and remove the noise. If (σnx)>> 1, the noise dominates the signal. The normalized threshold is chosen to be large to remove the noise more aggressively. In the Eq.(28), if σn2σy2, σx will become zero i.e. λ becomes infinity. Hence, for this case, λ =max ({| Wij |}).

The experimentation is carried out for each filter order, for each family and the decompositions are performed upto 4 levels. The statistical features: variance, PSNR, MSE, SNR and correlation coefficient, computed for different wavelets: DW (db1, db5), CW (coif1, coif2, and coif5), SW (sym1, sym3, sym5) and BW (bior2.2, bior6.8). The BW family, where filters are of order 6 in decomposition and order 8 in reconstruction (BW6.8), gives the better results among all the filter types. The PSNR is calculated for bior 6.8 filter up to 4 decompositions. The PSNR value increases up to 3 decompositions and thereafter reduces. Hence, the optimal level of decomposition is 3. The results obtained for the different thresholding schemes are given in the Table 2 [27]. From Table 2, it is found that, there is significant improvement of the subband threshold approach in terms of image quality assessment parameters over the global threshold approach. This is because subband thresholding approach employs an adaptive thresholding approach to respond to the changes in the noise content of the different subbands. In contrast, the global thresholding relies on a Visu threshold to threshold all subbands. The label (written in the parenthesis) in the Table 2. indicates the global(I) and subband thresholding (II), Bayes’ (1) or Visu (2) shrinkage rule and the hard (i), soft (ii) or semi soft (iii) thresholding function employed to generate that image. The optimal thresholding scheme with shrinkage rule and shrinkage function are determined by the criteria, namely, higher SNR and PSNR values, lower Variance, MSE values and Correlation Coefficient is nearly equal to one. From Table 2, it is observed that subband decomposition (II) with soft thresholding (ii) using Bayes shrinkage rule (1) gives better results than other techniques. Bayes shrinkage is better than the universal threshold because the universal threshold tends to be high, killing many signal coefficients along with the noise. The Figure 6. shows the visual quality of proposed method based on wavelet filter with Wiener filter. It is observed that the wavelet filter (bior6.8 with level 3(L3) ) yields better visualization effect and denoised image than the Wiener filter. The Table 3. shows the performance comparison of proposed method based on wavelet filter with despeckle filter, namely Wiener filter. From Table 3, it is observed that the wavelet filter (bior6.8 with level 3) yields better visualization effect and denoised image than the wiener filter in terms of variance and computational time. But the method based on despeckling using wavelet transforms needs to be improved interms of image quality assessment parameters.

Bayes’ Shrinkage rule(1)Global thresholding (I)Subband thresholding (II)
Hard (i)Soft (ii)Semi soft (iii)VarianceMSESNRCCPSNRVarianceMSESNRCCPSNR
0.03960.002316.890.99126.380.04820.001518.940.99428.48
0.03960.002514.830.99226.020.03680.001219.060.99529.94
0.03700.002416.040.99226.170.04560.001617.940.99328.56
Visu shrinkagerule (2)Global thresholding (I)Subband thresholding (II)
Hard (i)Soft (ii)Semi soft (iii)0.04310.002815.940.99325.520.03920.002916.150.99227.33
0.04610.003214.120.99224.920.03720.002317.790.99428.78
0.04810.003114.150.99225.000.03850.003216.170.99327.51

Table 2.

Performance evaluation of different thresholding methods interms of variance, MSE,PSNR,SNR and CC values.

Denoising methodsPSNRSNRMSEVarianceCCComputational time (in Secs.)
WT-L3-ST (bior 6.8)29.9419.060.001210.03680.99542.80
Wiener (3x3)32.3421.530.000580.04890.99774.14

Table 3.

Performance comparison of wavelet based despeckling method and Wiener filter.

Figure 6.

Despeckled images using Wiener filter and Wavelet filter

3.3. Laplacian pyramid transform

Several speckle reduction techniques based on multi scale methods (e.g. Wavelet transform, Laplacian pyramid (LP) transform) have been proposed [30-33]. The LP has the distinguishing feature that each pyramid level generates only one bandpass image (even for multidimensional cases), which does not have “scrambled” frequencies. This frequency scrambling happens in the wavelet filter bank when a high pass channel, after down sampling, is folded back into the low frequency band, and thus its spectrum is reflected. In the LP, this effect is avoided by down sampling the low pass channel only.

A speckle reduction method based on non linear diffusion filtering of band pass ultrasound images in the Laplacian pyramid domain has been proposed in [34], which effectively suppresses the speckle while preserving edges and detailed features. In [31], the authors have implemented a nonlinear multiscale pyramidal transform, based on non overlapping block decompositions using the median operation and a polynomial approximation. It is shown that this structure can be useful for denoising of one and two dimensional (1-D and 2-D) signals. It can be used for the selection of thresholds for denoising applications.

In [33] the comparison of two multiresolution methods: Wavelet transform and Laplacian pyramid transform, for simultaneous speckle reduction and contrast enhancement for ultrasound images is given. As a lot of variability exists in ultrasound images, the wavelet method proves to be a much better method than the Laplacian one for an overall improvement. However, the Laplacian pyramid scheme need to be explored for achieving better despeckling results.

3.3.1. Laplacian pyramid scheme

One way of achieving a multiscale decomposition is to use a Laplacian pyramid (LP) transform [35]. In the first stage of the decomposition, the original image is transformed into a coarse signal and a difference signal. The coarse signal has fewer samples than the original image but the difference signal has the same number of samples as the original image. The coarse signal is a filtered and down sampled version of the original image. It is then up sampled and filtered to predict the original image. The prediction residual constitutes the detail signal. The coarse signal can be decomposed further and this process can be repeated a few times iteratively. Assuming the filters in LP are orthogonal filters, an image X is decomposed into J detail images dj, j=1, 2,...,J and a coarse approximation image cJ. Then, we have

X2=j=1Jdj2+cJ2E30

The Laplacian is then computed as the difference between the original image and the low pass filtered image. This process is continued to obtain a set of detail filtered images (since each one is the difference between two levels of the Gaussian pyramid). Thus the Laplacian pyramid is a set of detail filters. By repeating these steps several times, a sequence of images are obtained. If these images are stacked one above another, the result is a tapering pyramid data structure and, hence the name the Laplacian pyramid.

A speckle reduction method based on Laplacian pyramid transform for medical ultrasound image is illustrated using the block diagram shown in the Figure 7. In the Figure 7, a homomorphic approach such as the log transformation of the speckle corrupted image, converts the multiplicative noise of the original image into additive noise. Homomorphic operation simultaneously normalizes the brightness across an image and increases contrast. For every difference signal of N-level of Laplacian pyramidal decomposition a threshold value is calculated using Bayes’ shrinkage rule. Further, thresholding is performed to reduce speckle. The exponential operation is performed on the filtered output to obtain the despeckled image.

Figure 7.

Block diagram of speckle noise suppression using Laplacian pyramid transform

The Laplacian pyramid transform is performed on the log transformed image. The Laplacian pyramidal decompositions up to six levels are obtained using biorthogonal filters with sufficient accuracy numbers such as the “9/7” and “5/3”. Further, thresholding schemes such as hard thresholding, soft thresholding and semi soft thresholding is performed to reduce speckle. The threshold value is calculated using Bayes’ shrinkage rule. The experimentation is carried out on 70 ultrasound images of liver and kidney. The performance evaluation of the proposed method is done in terms of variance, MSE, SNR, PSNR, CC values that are computed from despeckled image. The Laplacian pyramid transform with 1 level decomposition and hard thresholding is observed to be better than other thresholding methods.

The Table 4. shows the performance comparison of the proposed LP transform based despeckling method with the wavelet transform based despeckling method [36]. It is noticed that, in comparison with the despeckling medical ultrasound images based on WT, the despeckling based on LP method yields poor results. Because multiplicative noise is a particular type of signal dependent noise, in which the amplitude of the noise term is proportional to the value of the noise free signal having nonzero mean. Therefore, a band pass representation like LP is not suitable for multiplicative noise model, So the method needs to be improved.

In order to capture smooth contours in the images, the contourlet transforms, which allow directional decompositions, are employed for despeckling medical ultrasound images in the next section.

Despeckling methodsPSNRSNRMSEVarianceCCComputational time (in Secs.)
LP transform based despeckling method (LP-HT-L1)28.4418.230.001410.04160.99722.24
Wavelet transform based despeckling method (WT-L3-ST)29.9419.060.001210.03680.99542.80

Table 4.

Performance comparison of the LP transform method and the wavelet transform based despeckling method

4. Contourlet transform method

The contourlet transform (CT) is a multiscale and multidirectional framework of discrete image. It is the simple directional extension for wavelet that fixes its subband mixing problem and improves its directionality. Among the “beyond wavelet” techniques, contourlet allows for different and flexible number of directions at each scale, while achieving nearly critical sampling. The desirable properties of CT for image representation includes multiresolution, allowing images to be approximated in a coarse to fine fashion; localization of the basis vectors in both space and frequency; low redundancy, so as not to increase the amount of data to be stored; directionality, allowing representation with basis elements oriented in a variety of directions; and anisotropy, the ability to capture smooth contours in images, using basis elements that are a variety of elongated shapes with different aspect ratios [37].

The contourlet transform has been developed to overcome the limitations of the wavelets, and hence, the new algorithms based on the contourlet transform are more efficient than wavelet methods. In [38], the authors have presented a contourlet based speckle reduction method for denoising ultrasound images of breast. The double iterated filter bank structure and a small redundancy at most 4/3 using two thresholding methods shows a great promise for speckle reduction. In [39], the despeckling medical ultrasound images using contourlet transform using Bayes’ shrinkage rule is investigated. The algorithm is also tested on ovarian ultrasound images to demonstrate improvements in the segmentation that yields good classification for follicle detection in an ovarian image [40].

In [41], speckle reduction based on contourlet transform using scale adaptive threshold for medical ultrasound image has been examined, where in the subband contourlet coefficients of the ultrasound images after logarithmic transform are modelled as generalized Gaussian distribution. The scale adaptive threshold in Bayesian framework is applied. The method is tested on both synthetic and clinical ultrasound images interms of S/MSE and edge preservation parameter. The proposed method exhibits better performance on speckle suppression than the wavelet based method, while it does well preserve the feature details of the image.

The contourlet transform can be divided into two main steps: Laplacian pyramid decomposition and directional filter banks. Contourlet transform is a multi scale and directional image representation that uses first a wavelet like structure for edge detection, and then a local directional transform for contour segment detection. A double filter bank structure of the contourlet obtains sparse expansions for typical images having smooth contours. In the double filter bank structure, the Laplacian pyramid is used to decompose an image into a number of radial subbands, and the directional filter banks decompose each LP detail subband into a number of directional subbands. The band pass images (dj [n]) from the LP are fed into a DFB so that directional information can be captured. The scheme can be iterated on the coarse image (cj [n]). The combined result is a double iterated filter bank structure, named pyramidal directional filter bank (PDFB), which decomposes images into directional subbands at multiple scales. The general model for despeckling an image using contourlet transform is shown in the Figure 8.

Figure 8.

The general model for speckle reduction using contourlet transform.

In the Figure 8, the CT based despeckling method consists of the log transformed original ultrasound image being subjected to contourlet transform, to obtain contourlet coefficients. The transformed image is denoised by applying thresholding techniques on individual band pass sub bands using a Bayes shrinkage rule, derived from the local statistics of the signal in the transform domain. Bayes’ shrink was proposed by [29]. The goal of Bayes’ shrinkage method is to minimize the Bayesian risk, and hence its name, Bayes’ shrink. Further, thresholding schemes such as hard thresholding, soft thresholding or semi soft thresholding is performed to reduce speckle. The exponential operation is performed on the filtered output to obtain the despeckled image.

The experimentation is carried out on 70 ultrasound images of liver and kidney. The six levels of Laplacian pyramidal decompositions are performed using biorthogonal filters with sufficient accuracy numbers such as the “9/7”. The directional decompositions up to eight are performed in all the pyramidal levels, using two dimensional ladder filters. The contourlet transform uses the “9/7” filters in LP stage because, in the multiscale decomposition stage, it significantly reduces all inter scale, inter location and inter direction mutual information of contourlet coefficients. Similarly, in directional decomposition stage, the ladder structure PKVA filters [42] are more effective in localizing edge direction as these filters reduce the inter direction mutual information. Further, thresholding schemes such as hard thresholding, soft thresholding or semi soft thresholding is performed to reduce speckle. The threshold value is calculated using Bayes’ shrinkage rule. The PSNR is calculated up to 6 LP decompositions. The PSNR value increases up to 2 decompositions using HT, ST and SST, and thereafter reduces. Hence, the optimal level of LP decomposition is 2. Further, it is observed from Table 5, that the 2-level Laplacian pyramidal decomposition and 4 directional bandpass subbands (2 at level 1, 2 at level 2) using hard thresholding yield better results than soft thresholding and semi soft thresholding techniques.

Thresholding methodsLevelsPSNRVarianceMSECCSNRComputational time(in Secs.)
Hard thresholdingL2-1134.210.02030.000370.999027.741.64
Soft thresholdingL2-3132.490.03820.000560.998324.141.68
Semi-soft thresholdingL2-1133.310.02920.000460.998725.911.74

Table 5.

The results obtained for the optimal optimal decomposition of LP levels and directional decompositions in terms of image quality assessment parameters using contourlet method based on different thresholding techniques with Bayes’ shrinkage rule.

Figure 9.

a) Original ultrasound image. (b) Despeckled image using wavelet transform using subband Bayes soft thresholding (level 3). (c) Despeckled image using contourlet transform using soft thresholding. (d) Despeckled image using contourlet transform using hard thresholding. (e) Despeckled image using contourlet transform using semi soft thresholding.

The frequency bands obtained by using optimal level L2-11 of contourlet decomposition are as follows: the 2nd level has 1 approximation band of size 128 x128 and 4 detail components (2 of 128 x 256, 2 of 256 x 512). The reconstructed image is the despeckled image. The hard thresholding is better than other thresholding methods, because small coefficients are removed while others are left untouched in HT, while in ST or SST coefficients above the threshold are shrunk by absolute value of threshold. Further, it is found that the despeckling based on contourlet transform gives better results than the speckle reduction method based on wavelet transform in particular. The wavelet based Bayes’ shrink thresholding method is based on separable 2D wavelet transform that has limited directions (Horizontal, Vertical and Directional). Speckle noise in medical ultrasound images will generate significant coefficients in wavelet domain just like true detail features, such as edges. However, the speckle noise is less likely to generate significant coefficients in the contourlet method, and thus, it directly leads to better performance in suppressing noise than the wavelet based Bayes’ shrink thresholding scheme.

Another way to analyze the effects of despeckling techniques is to study the despeckled images. In the Figure 9., the resultant images of a sample medical ultrasound image are presented to compare the results of different despeckling techniques by visual inspection. In the Figure 9(b) and (c), the speckle is reduced considerably, but the structures are blurred and some visible artifacts are introduced. However, in Figure 9(e), the speckle is reduced well and structures are enhanced. But some details are lost and some are over enhanced. It is encouraging to note that in the Figure 9(d), the speckle is effectively reduced and also structures are enhanced with almost no loss or noticeable artifact.

4.1. Cycle spinning based contourlet transform

The CT is not translation invariant. This means that the errors after denoising will be sensitive to the positions of discontinuities in the data. In order to avoid such effects, it is necessary to build translation invariant version of the transform. Translation invariance is achieved through several ways. For example, in [43], the time invariant schemes of wavelet based decompositions have been proposed and have been often referred to as cycle spinning. Unfortunately, due to the downsamplers and upsamplers present in the directional filter banks of CT, the CT is not shift invariant, which is important in image denoising by thresholding and normally causes pseudo-Gibbs phenomenon. In [44], the cycle spinning algorithm is utilized in developing a translation invariant contourlet based denoising technique. The experimental results clearly demonstrate the capability of the proposed scheme in image denoising, especially for detailed texture images. It is shown that most of the visual artifacts resulting from the contourlet transform denoising process are eliminated. In [45], a cycle spinning method is used to compensate for the lack of translation invariance property of sharp frequency localized contourlet. Experimental results demonstrate that cycle spinning is a simple and efficient way to average out the pseudo-Gibbs phenomena, which are around singularities and produced by the down sampling and up sampling of directional filter banks, and improve the denoising performance interms of visual quality and PSNR.

To compensate for the lack of translation invariance property of the contourlet transform, we apply the principle of cycle spinning to contourlets. Suppose X and Y are original and despeckled images, F and F−1 are forward and inverse contourlet transform, Si,j is the 2D circular shift in ith row and jth column directions, λis the threshold operator in contourlet transform domain. The cycle spinning based contourlet transform for image denoising could be described as

Y=1B2i,j=1BCS-i,-j(F-1(λ(F(CSi,j(X)))))   E31

where B is the series of bit shifts in the ith row and jth column directions. If one decomposes an image of size (N, N) using the contourlet transform, the maximum number of decomposition levels in the LP stage will be B, and therefore, the maximum number of shifts are (B, B) in the row and column directions. After a B number of bit shifts, which depends on the level of decomposition, the transform output degrades. Hence, the cycle spinning has to be stopped after a certain number of bit shifts. The Figure 10. shows the block diagram for speckle reduction method based on contourlet transform with cycle spinning. The cycle spinning is applied to the log transformed image. It performs two dimensional circular shift in ith row and jth column directions. The circular shifting is performed up to B number of bit shifts, where B depends on the level of decomposition. The transform output degrades as B increases. Hence, the cycle spinning has to be stopped after a certain number of bit shifts. Then contourlet transform is performed using double filter bank structure. The six levels of Laplacian pyramidal decompositions are performed using biorthogonal filters with sufficient accuracy numbers such as the ”9-7”.The directional decompositions up to six is performed in the lowest pyramidal level, using two dimensional ladder filters designed in [42]. Further, a thresholding scheme either hard thresholding, soft thresholding or semi soft thresholding, is performed to reduce speckle. The threshold value is calculated using Bayes’ shrinkage rule. The results obtained for the optimal bit shifts of cycle spinning using contourlet method based on hard thresholding, soft thresholding and semi soft thresholding using Bayes’ rule are presented. The results obtained for different reconstruction methods are shown in Figures 11-15, which exhibit graphs of statistical features PSNR, SNR, variance, MSE and CC, respectively, for different levels of Laplacian pyramid decompositions and directional decompositions corresponding to optimal results of contourlet transform with cycle spinning based despeckling method [46]. The optimal reconstruction method is determined by the criteria, namely, lower variance and MSE, higher SNR and PSNR values, Correlation Coefficient is nearly equal to one. The contourlet transform with 2 level of pyramidal decomposition and two directional decompositions in the finest scale and hard thresholding technique with Bayes’ shrinkage rule has yielded better results in comparison contourlet transform based methods [47]. In the Figures 11-15, the horizontal axis label CYC-HT-Bn indicates cycle spin (CYC), thresholding (HT,ST,SST), n number of bit shifts in cycle spinning. From Figures 11-15, it is observed that, 4 bit cycle spinning, having the 2-level of Laplacian pyramidal decomposition with 4 directional bandpass sub bands (2 at level 1, 2 at level 2) subject to soft thresholding, yields optimal results for speckle reduction. The computational time (in Secs.) of the cycle spinning based CT method is shown in Figure 15. The CT based despeckling method takes less computational time as compared to cycle spinning based CT method.

Figure 10.

The block diagram of despeckling method based on contourlet transform with cycle spinning.

The Figure 16. illustrates the resultant despeckled images of an ultrasound image obtained by the cycle spinning based CT method using hard, soft and semi soft thresholding with Bayes’ shrinkage rule, and also that obtained by the CT method [47], for comparison by visual inspection The despeckling method based on cycle spinning using contourlet coefficient shrinkage (Figure16.(b)) performs better and appears to be an improvement over direct contourlet transform based despeckling method.

Among the transform domain filters developed in the previous sections, the contourlet transform with cycle spinning yields better visual quality enhancement of the despeckled images. However, there is still a need to remove Gaussian noise inherent in the medical ultrasound images, which is addressed in the next section.

Figure 11.

The values of PSNR and SNR for various thresholding methods obtained by cycle spinning based CT method and that by direct CT based method.

Figure 12.

The values of MSE for various thresholding methods obtained by cycle spinning based CT method and that by direct CT based method.

Figure 13.

The values of Correlation Coefficient for various thresholding methods obtained by cycle spinning based CT method and that by direct CT based method.

Figure 14.

The values of variance for various thresholding methods obtained by cycle spinning based CT method and that by CT based method.

Figure 15.

The values of computational time for various thresholding methods obtained by cycle spinning based CT method and that by CT based method.

5. Gaussian model for speckle noise

Gaussian probability density is often used in image filtering. The Gaussian has the unique ability not to create new edges as its scale (standard deviation) is increased. This property enables the extraction of edges that represent different levels of detail in an image. As the scale is increased, the number of extracted weak and false edges reduces. However, at the same time edges shift from their true positions. The amount of shift an edge makes not only depends upon the scale of the Gaussian filter, but also on the intensity distributions of the underlying image [48]. Mathematically, a Gaussian filter modifies the input signal by convolution with a Gaussian function. Smoothing is commonly undertaken using linear filters such as the Gaussian function (the kernel is based on the normal distribution curve), which tends to produce good results in reducing the influence of noise with respect to the image. The 2D Gaussian distributions, with standard deviation σ for image X, is given by Eq. (32) [49].

GX,m,σ=1σ2πe-12X-mσ2E32

Figure 16.

a) Original ultrasound image (b) Despeckled image using cycle spinning based CT using ST method. (c)Despeckled image using cycle spinning based CT using HT method. (d)Despeckled image using cycle spinning based CT using SST method. (e) Despeckled image using direct CT method.

5.1. Pre- or post –processing

Usually, medical ultrasound images are affected by the mixed noise, which is the combination of speckle noise and Gaussian noise. There are two factors that influence the usefulness of a smoothing filter. The first reduces the range of resolutions over which variations in the output appear by the filter variation ∆w, in the frequency domain be small and second factor is increase in spatial localization by a small spatial variance ∆x. These localization requirements in the spatial and frequency domain are conflicting and related by the uncertainty principle given in Eq. (33)

w x  14πE33

It has been shown that Gaussian functions are the only ones that provide the optimal trade off between these conflicting requirements constrained by the Eq. (34) [50]..

wx= 14πE34

So Gaussian filters are widely used in image filtering. In [51], the removal of mixed noise using order statistic filter and wavelet domain Wiener filter is proposed. Authors have evaluated two methods. The first method comprises, order statistic prefilter and empirical Wiener filter, which is used to reduce the Gaussian noise. The disadvantage of this method is the higher time consumption. The second method is, order statistic filter for each decomposition level, where decomposition is carried out by the wavelet transform, followed by thresholding. The drawback of this method is that its efficiency is less than that of the first method (about 1dB) in removing the mixed noise. In [52], denoising of mixed noise in ultrasound images is presented. Combined Bayesian maximum a posterior (MAP) estimator and ST-PCNN (Soft threshold pulse coupled neural networks) method has been used for mixed noise reduction. The method removes the speckle noise considerably than the Gaussian noise that degrades the ultrasound images. The drawback of the method is either Gaussian noise or speckle noise is removed. Hence we present a method to remove residual Gaussian noise from despeckled image.

Two alternative algorithms are developed for reducing mixed noise in medical ultrasound images [53]. In the first alternative, the denoising method reduces the Gaussian noise by applying Gaussian filter in pre processing stage, then despeckling is performed using either wavelet transform, Laplacian pyramid transform or contourlet transform. The noise model for the first alternative (i.e. Gaussian noise removal in pre processing followed by despeckling) is given by the Eq.(35).

 Xij=fijnij+gijE35

where Xij represents the noisy pixel in the image X, fij represents the noise free pixel, nij and gij represent the multiplicative speckle noise and additive Gaussian noise, respectively. The indices i, j represent the spatial position over the image. We use transform domain filtering techniques [36,54,47] for despeckling along with Gaussian filter in pre processing stage for removal of Gaussian noise.

In the second alternative, the despeckling of medical ultrasound images is performed either using wavelet transform, Laplacian pyramid transform or contourlet transform and, then, it is followed by postprocessing stage in which Gaussian filter removes Gaussian noise from the despeckled image. However, the second alternative assumes the noise model given by the Eq.(36)

Xij=fij+gij nij E36

The second alternative is investigated for image quality enhancement, due to noise removal, in an ultrasound image.

The experimentation is carried out using various kernel sizes and different values of σ. Larger values of σ produce a wider peak influencing the greater blurring. Kernel size is increased with increasing σ to maintain the Gaussian nature of the filter. Gaussian kernel coefficients depend on the value of σ. The Figure 17. shows different convolution kernels that approximate a Gaussian with σ. Gaussians are locally sensitive and can be made more spatially localized by decreasing parameter σ. It is observed that the kernel size 3×3 with σ = 0.5 yields better results than other kernels. It is found that larger kernels of size 5×5 or 7×7 produce better denoising effect but make the image more blurred. Thus, the empirically determined kernel size 3×3 and σ=0.5 are used in two alternative methods (Gaussian filter in Pre or Post processing).The two alternative methods are evaluated in terms of filter assessment parameters, namely, PSNR, SNR, MSE, variance and CC. The comparisons of the performance of the both alternatives with the despeckling methods discussed in [36,54, 47] are given in the Table 6. From the Table 6, it is observed that the Gaussian filter in pre processing stage is found to be more effective than that in despeckling based on Laplacian pyramid transform and contourlet transform. However, the Gaussian filter in postprocessing stage is found to be more effective in despeckling based on wavelet transform. Thus, the Gaussian filter improves the performance of despeckling methods, because Gaussian noise is characterized by adding to each image pixel a value from a zero mean Gaussian distribution. The zero mean property of the distribution allows such noise to be removed by locally averaging pixel values [55]. Further, it is observed that, the Gaussian filter in pre processing stage followed by contourlet transform based despeckling method yields better visual enhancement than the other denoising methods, which is illustrated in the Figure 18. The denoising and visual enhancement techniques developed in this study lead to improvement in the accuracy and reliability of automatic methods for medical ultrasound imaging systems.

Figure 17.

The 3×3 kernel with (a) σ = 0.5, (b) σ = 1.

Denoising methodsSNRPSNRCCMSEVarianceComputational time (in Secs.)
First alternative with wavelet transform25.2230.230.99920.000500.01280.91
First alternative with LP transform20.0731.660.99650.001060.04260.38
First alternative with CT transform29.9935.150.99980.000310.00711.62
Second alternative with WT transform26.3232.540.99950.000500.01200.89
Second alternative with LP transform19.7630.120.99590.001330.04400.76
Second alternative with CT transform28.8934.350.99960.000360.00882.25
Despeckling method (WT transform ) [36]19.0629.940.99540.001210.03682.80
Despeckling method (LP transform) [54]18.2328.440.99720.001410.04162.24
Despeckling method (CT transform) [47]27.7434.210.99900.000370.02031.64

Table 6.

Comparison of performance of denoising methods based on Gaussian filtering with despeckling methods.

5.2. Linear regression model

We present a linear regression based approach for clinical ultrasound image despeckling in the spatial domain. We propose a linear regression model for Gaussian noise representation of speckle noise for medical ultrasound images. This approach introduces an adaptive filter, well preserving edges and structures in the image. The parameters in the model are estimated through an efficient iterative scheme.

Figure 18.

a) Original ultrasound image,(b) Denoised image using Gaussian filter, (c) Denoised image using 1st alternative with WT method (d) Denoised image using 1st alternative with LP method. (e) Denoised image using 1st alternative with CT method. (f) Denoised image using 2nd alternative with WT method. (g) Denoised image using 2nd alternative with LP method. (h) Denoised image using 2nd alternative with CT method.

In [56], the authors have developed the adaptive weights smoothing algorithm, which is an iterative procedure in which the size of a neighbourhood is adaptive to the surface smoothness. In [57], the estimation of jump surfaces by local piecewise linear kernel smoothing is examined. In [58, 59], an anisotropic diffusion algorithm has been proposed for Gaussian noise removal. In [60], a bilateral filter to remove Gaussian noise is developed. In [61], a window based linear regression filter for echo cardiographic image denoising is proposed. The main draw backs of the above algorithms are that they need more computational time and complex circuit to implement them.

We consider a medical ultrasound image X and the corresponding despeckled image Y obtained by using the contourlet transform with cycle spinning [46]. The subtracted image Z=X-Y is the error image containing speckle noise. We find the mean m and standard deviation σ of Z and then simulate Gaussian noise G with these values of m and σ. The removal of this Gaussian noise G from despeckled image Y yields the new despeckled image Y^, i.e. Y^=Y-G, which is further subtracted from the original image X to obtain the new error image Z containing the residual speckle noise. This procedure is repeated until the percentage of black pixels in error image Z reaches 99.9. We determine the maximum value of PSNR and the corresponding values of mean m and standard deviation σ using the iterated despeckled images Y^. This procedure is applied for all the medical ultrasound images Xi, i=1,..., 63, in the dataset, yielding the two sets of data points (PSNRi,mi) and (PSNRii), i=1,...,n1 exhibits linear correlation. Using the method of least square errors, we obtain the lines of best fit for these data, namely:

m=a * PSNR + bE37
σ =c * PSNR + dE38

where a, b, c and d are the constants, which are the linear regression model parameters for Gaussian representation G(m, σ) of speckle noise in the medical ultrasound images. The Eqs.(37) and (38) represent the linear regression model for Gaussian representation G(m, σ) of speckle noise in a medical ultrasound image with its PSNR value determined already [62]. The steps involved in this procedure are given in the Algorithm 1.

Algorithm 1. Linear regression model for Gaussian representation G (m, σ) of speckle noise.

Input : Medical ultrasound image.
Output:Linear regression model parameters for Gaussian representation G (m, σ) of speckle noise.
Start
Step 1:Input medical ultrasound image X.
Step 2:Input despeckled ultrasound image Y obtained by using contourlet transform with cycle spinning.
Step 3:Find the error image Z as difference between X and Y
(i.e. Z = X – Y).
Step 4:Find the percentage of black pixels in Z.
Step 5:Find the mean (m) and standard deviation (σ) of Z.
Step 6:Simulate the Gaussian noise G with mean m and standard deviation σ
obtained in the Step 5.
Step 7:Subtract the simulated Gaussian noise G from despeckled image Y to obtain the resultant image (Y^), i. e.Y^=Y-G.
Step 8:Find the PSNR, mean and standard deviation of the resultant image
(Y^) obtained in the Step 7. Set Y =Y^.
Step 9:Repeat the Steps 3-8 until the percentage of black pixels in Z reach 99.9%.
Step 10:Find the maximum of the PSNR values obtained for the iterated despeckled images
(Y^) and the corresponding values of m and σ. Store the values of maximum PSNR and corresponding values of m and σ.
Step 11:The Gaussian noise with mean m and standard deviation σ determined in the Step 10 is the modelled Gaussian noise for the image X.
Step 12:Repeat the Steps 1-11 for all the ultrasound images X and their corresponding despeckled images Y in the data set.
Step 13:Obtain the lines of best fit for the both PSNR vs. m and PSNR vs. σ obtained for all the ultrasound images (n1) in the Step 12 using the method of least square errors. i. e. determine the constants a, b, c and d of the Eqs.(37) and (38), which are the lines of best fit that form the linear regression model for representing the speckle noise in ultrasound image as Gaussian noise. The lines of best fit are found using goodness of fit static.
Step 14:Output linear regression model parameters a, b, c and d of Eqs.(37) and (38).
Stop

The data points (PSNRi,mi) and (PSNRii), i=1,...,63 obtained for 63 typical ultrasound images of the image data set are stored and then used for regression of PSNR on m and also PSNR on σ. In the ultrasound image data set used for building regression model for Gaussian representation of speckle noise, we select a reference image Xref for which the PSNR value is minimum. Given an arbitrary input medical ultrasound image X, we compute the PSNR of X with respect to the reference image Xref. Using linear regression model (Eqs. (37) and (38)), we estimate the values of mean m and standard deviation σ of the Gaussian noise, which is removed from the input original image. The resultant image is the despeckled image (Y). The ‘goodness of fit’ statistic for the lines of regression is given by two quantities, namely, the sum of squares due to error of the fit (SSE) and root mean squared error (RMSE).The SSE is given by Eq. (39) :

SSE=i=1n(Yi-Y^i)2E39

where Yi=actual mean value of Gaussian noise, Y^i= estimated mean value of Gaussian noise, and

n1= total number of ultrasound images in the dataset.

The RMSE is given by Eq. (40) :

RMSE=SSEn1E40

If the SSE and RMSE values are closer to zero, they indicate better fit. The general model for Gaussian noise estimation and removal in despeckling ultrasound image is shown in the Figure 19.

Figure 19.

The general model for Gaussian representation of speckle noise.

The steps involved in the denoising process, shown in Figure 19 are given in the Algorithm 2.

Algorithm 2. Despeckling based on linear regression model for Gaussian noise estimation.

Input: Medical ultrasound image.
Output: Despeckled image.
Start
Step 1:Input medical ultrasound image X.
Step 2:Compute PSNR of image X with respect to the reference image Xref prescribed by
the regression model.
Step 3:Obtain Gaussian noise estimation G(m, σ) corresponding to the PSNR computed in
the Step 2 using the regression model (Eqs.(37) and (38))
Step 4:Apply Gaussian filter with estimated mean (m) and standard deviation (σ), obtained
n the Step 3, to the input image X, which yields the despeckled image Y.
Step 4:Apply Gaussian filter with estimated mean (m) and standard deviation (σ), obtained in the Step 3, to the input image X, which yields the despeckled image Y.
Step 5:Output the despeckled image Y.
Stop.

The linear regression model parameters a, b, c and d for Gaussian representation of speckle noise are computed for the dataset of 63 ultrasound images and the linear regression model equations are (Eqs.(37) and (38)),where a=-6.129e-007, b=2.742e-005, c=-0.0002192, d=0.01004, with the measures of ’best fit’ are SSE=4.682e-009, RMSE=8.833e-006 for mean vs. PSNR and SSE=0.0006471, RMSE=0.003284 for standard deviation vs. PSNR. The Figures(20 and 21) show the lines of best fit for mean vs. PSNR and standard deviation vs. PSNR, respectively, which are used for Gaussian noise estimation and removal.

Figure 20.

Linear regression of mean on PSNR

Figure 21.

Linear regression of standard deviation on PSNR

The comparison of the results of the proposed method with the contourlet transform method (with cycle spinning) is given in the Table 7. It is observed that the image quality enhancement obtained by the despeckling method based on linear regression model is better than that obtained by the contourlet transform method in terms of PSNR and computational time required for denoising.

Denoising methodsPSNR(in dB)Computational Time(in Secs.)
Contourlet transform using cycle spinning35.9124.09
Proposed method based on linear regression model36.980.34

Table 7.

Comparison of performance of despeckling based on contourlet transform and proposed method based on linear regression model.

The Figure 22. shows a sample medical ultrasound image, its despeckled image using contourlet transform with cycle spinning and the denoised image using the linear regression model respectively. The visual quality of image enhancement can also be observed from the sample image and its denoised image. The anatomical structures are more clearly visible in the Figure 22.(c) than that in Figure 22.(b). The box indicates the region of image in (b) and (c) showing prominent visual enhancement due despeckling methods.

Figure 22.

a) Original ultrasound image (b) Despeckled image using the contourlet transform with cycle spinning (c) Denoised image using the proposed linear regression model. The box Indicates the region of image in (b) and (c) showing prominent visual enhancement due to the despeckling.

Figure 23.

a) A sub image of original ultrasound image (b) sub image of despeckling using Wiener filter (c) sub image of despeckling using WT method (d) sub image of despeckling using CT method (e) sub image of despeckling using cycle spinning based CT method (f) sub image of despeckling using proposed linear regression model.

The Figure 23 shows the visual enhancement due to various despeckling methods for comparison. (a) shows the sub image of original image, The Figure 23 (b)-(f) indicates the sub image showing visual enhancement due to different despeckling methods namely, Wiener filter with (3X3),wavelet transform method, contourlet transform method, cycle spinning based contourlet transform method and proposed linear regression model. The prominent visual enhancement is observed using the proposed linear regression model.

The proposed method estimates the Gaussian noise content in the input medical ultrasound image for denoising the image efficiently. Hence, it is easily amenable for building embedded system software for ultrasound imaging equipments in order to display the high quality images, which helps the medical expert in the diagnosis with greater accuracy.

6. Conclusion

In this chapter, a despeckling method, based on a 2D directional non separable transform known as contourlet transform is presented. Conventional 2D wavelet transform is separable and thus cannot sparsely represent non separable structures of the image, such as directional curves. It is found that pyramidal directional filter bank feature of contourlet transform makes it a good choice for representation of curves and edges in the image. But, the contourlet transform, one of the recent geometrical image transforms, lacks the feature of translation invariance due to sub sampling in its filter bank structure. In cycle spinning, CT is improved by averaging the estimation of all translations of the degraded image. The Gibbs effect is considerably reduced by the contourlet transform with cycle spinning, because the average of different estimations of the image reduces the oscillations. In the literature, the authors [33,41,45,54,55,61] have considered ultrasound images (natural/synthetic) with artificially added speckle noise content and have proposed methods for despeckling such images. However, in the present study, we considered ultrasound images captured by the ultrasound equipment which contain inherent speckle noise and have proposed methods for removing the speckle noise more effectively.

When the noise characteristics of the images are unknown, it is proposed to denoise by a linear regression model, which is cost effective compared to the other methods. We have proposed a novel linear regression model for Gaussian noise estimation and removal in despeckling medical ultrasound images. The experimental results demonstrate its efficacy both in terms of speckle reduction and computational time required for denoising. Further, the proposed regression model is simple, generic and computationally inexpensive. Hence, it is easily amenable for building embedded system software for ultrasound imaging equipments in order to display the high quality images, which help the medical experts for speedy accurate image analysis and diagnosis. Further, the proposed regression model is simple, generic and computationally inexpensive.

Acknowledgments

Authors are grateful to the reviewers for their helpful comments which improved the quality of the paper. Further, authors are thankful to Dr. Ramesh Mankare, Radiologist, Sangameshwar Scanning Centre, Bijapur, Karnataka, India, for providing the ultrasound images of kidney, liver and also for helpful discussions.

© 2013 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution 3.0 License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

How to cite and reference

Link to this chapter Copy to clipboard

Cite this chapter Copy to clipboard

P.S. Hiremath, Prema T. Akkasaligar and Sharan Badiger (June 5th 2013). Speckle Noise Reduction in Medical Ultrasound Images, Advancements and Breakthroughs in Ultrasound Imaging, Gunti Gunarathne, IntechOpen, DOI: 10.5772/56519. Available from:

chapter statistics

3476total chapter downloads

14Crossref citations

More statistics for editors and authors

Login to your personal dashboard for more detailed statistics on your publications.

Access personal reporting

Related Content

This Book

Next chapter

Strategies for Hardware Reduction on the Design of Portable Ultrasound Imaging Systems

By D. Romero-Laorden, J. Villazón-Terrazas, O. Martínez-Graullera and A. Ibáñez

Related Book

First chapter

Introduction to Infrared Spectroscopy in Life and Biomedical Sciences

By Theophile Theophanides

We are IntechOpen, the world's leading publisher of Open Access books. Built by scientists, for scientists. Our readership spans scientists, professors, researchers, librarians, and students, as well as business professionals. We share our knowledge and peer-reveiwed research papers with libraries, scientific and engineering societies, and also work with corporate R&D departments and government entities.

More About Us