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
where the speckle image y_{ij} is the product of the original image X_{ij,} and the non-Gaussian noise n_{ij.} 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 y_{ij} [7]. The noise component n_{ij} is then approximated as an additive zero mean gaussian noise.
The Discrete Wavelet Transform (DWT) is then applied to ln y_{ij} 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
where
Mean Square Error: The MSE measures the quality change between the original image (
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
where,
Peak Signal-to-Noise Ratio: The PSNR is computed as
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
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 (T_{c}) 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, T_{c} 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 T_{c} is better than a filter having higher T_{c} 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)
where Y_{ij} 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
and then summed with
where MxM is the size of the kernel and K_{uv} is the pixel value within the kernel at indices u and v,
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):
The weighting function is computed from estimated noise variation coefficient of the image, C_{u} given by the Eq.(12):
and the variation coefficient C_{i} of the image given by the Eq. (13):
where ENL is given by the Eq.(14) :
In the Eq.(14), the
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
where Y_{ij} denotes the despeckled image,
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.
Despeckling filters | Computational 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 |
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
for some suitable sequence of coefficients a_{L}(k). Once
Wavelet analysis leads to perfect reconstruction filter banks using the coefficient sequences a_{L}(k) and a_{H}(k). The input sequence X is convolved with high pass (HPF) and low pass (LPF) filters
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, X_{w} denotes the input wavelet coefficients and Y_{t} denotes the output wavelet coefficients after thresholding, we define the following thresholding functions:
Hard thresholding:
Soft thresholding:
Semi soft thresholding:
(23) |
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 N_{k} is the size of the subband in the wavelet domain, then λ_{i}
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
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 W_{k} of the image X, is given by
where
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
and
in which N_{k} is the number of the wavelet coefficients W_{k} 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 (σ_{n} /σ_{x}) << 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 (σ_{n} /σ_{x})>> 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
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.
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 d_{j,} j=1, 2,...,J and a coarse approximation image
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.
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.
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 (d_{j} [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 (c_{j} [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.
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.
The frequency bands obtained by using optimal level L2-11 of contourlet decomposition are as follows: the 2^{nd} 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, S_{i,j} is the 2D circular shift in i^{th} row and j^{th} column directions,
where B is the series of bit shifts in the i^{th} row and j^{th} 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 i^{th} row and j^{th} 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.
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.
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].
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)
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]..
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).
where X_{ij} represents the noisy pixel in the image X, f_{ij} represents the noise free pixel, n_{ij} and g_{ij} 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)
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.
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.
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
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.
The data points (PSNR_{i,}m_{i}) and (PSNR_{i},σ_{i}), 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 X_{ref} 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 X_{ref}. 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) :
where Y_{i}=actual mean value of Gaussian noise,
n_{1}= total number of ultrasound images in the dataset.
The RMSE is given by Eq. (40) :
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.
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.
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.
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 methods | PSNR (in dB) | Computational Time (in Secs.) |
Contourlet transform using cycle spinning | 35.91 | 24.09 |
Proposed method based on linear regression model | 36.98 | 0.34 |
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.
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.