**3. Wavelet transform method**

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

**Figure 4.** Performance of various despeckling filters, in terms of Correlation Coefficient

**Figure 2.** Performance of various despeckling filters, in terms of variance

208 Advancements and Breakthroughs in Ultrasound Imaging

**Figure 3.** Performance of various despeckling filters, in terms of MSE.

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


calculated using circular kernel, mean max threshold, nearest neighbour and new threshold

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 *ϕ*) coeffi‐ cients, 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):

f

y

wavelet is given by a similar looking Eq.(17):

horizontal and diagonal details, are given by

<sup>L</sup> ( ) a (k) (2x-k) *<sup>k</sup>*

<sup>H</sup> ( ) a (k) (2x-k) *<sup>x</sup>* . *<sup>k</sup>*

(,) () () *<sup>V</sup>*

(,) ()() *<sup>H</sup>*

 yj

 jy

y

y

 f

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 aH(k), aL(k) respectively. Further, each result is down sampled by two, yielding the transform signals xH and 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 multi‐ plying 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,

 f

for some suitable sequence of coefficients aL(k). Once *ϕ* has been found, an associated mother

*<sup>x</sup>* <sup>=</sup> <sup>å</sup> (16)

Speckle Noise Reduction in Medical Ultrasound Images

http://dx.doi.org/10.5772/56519

211

<sup>=</sup> <sup>å</sup> (17)

*xy x y* <sup>=</sup> (18)

*xy x y* <sup>=</sup> (19)

function.

**Table 1.** Performance comparison of various despeckling filters based on computational time.

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 *ϕ*) coeffi‐ cients, 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):

$$\boldsymbol{\phi}(\mathbf{x}) = \sum\_{\mathbf{k}} \mathbf{a}\_{\mathbf{k}}(\mathbf{k}) \boldsymbol{\phi}(2\mathbf{x} \cdot \mathbf{k}) \tag{16}$$

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

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

**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

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.

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

The discrete transform is very efficient from computational point of view.

**Table 1.** Performance comparison of various despeckling filters based on computational time.

fidelity after reconstruction is visually lossless.

210 Advancements and Breakthroughs in Ultrasound Imaging

$$
\psi(x) = \sum\_{\mathbf{k}} \mathbf{a}\_{\mathbf{H}}(\mathbf{k}) \phi(2\mathbf{x} \cdot \mathbf{k}).\tag{17}
$$

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 aH(k), aL(k) respectively. Further, each result is down sampled by two, yielding the transform signals xH and 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 multi‐ plying 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

$$\text{hyp}^V(x, y) = \wp(x)\wp(y) \tag{18}$$

$$\text{If } \boldsymbol{\psi}^H(\mathbf{x}, \boldsymbol{y}) = \boldsymbol{\psi}(\mathbf{x})\boldsymbol{\phi}(\boldsymbol{y})\tag{19}$$

$$\boldsymbol{\Psi}^{\text{D}}(\mathbf{x}, \boldsymbol{y}) = \boldsymbol{\Psi}(\mathbf{x})\boldsymbol{\Psi}(\mathbf{y})\tag{20}$$

Yt=*Tsoft*(*X* w)

sign{Xw}(|Xw | −*λ*) , for | *X* <sup>w</sup> | ≥λ 0 , for | *X* <sup>w</sup> | <λ

0 , for |Xw | ≤λ

*X* <sup>w</sup> , for |Xw | >λ<sup>1</sup>

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

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:

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

subband of the ultrasound image after decomposition. If Nk is the size of the subband in the

to be proportional to the square root of the local noise variance σ<sup>2</sup>

λ σ 2log N <sup>i</sup> ( ) <sup>k</sup> <sup>=</sup> (24)

<sup>λ</sup><sup>1</sup> −λ , for λ<sup>&</sup>lt; |Xw <sup>|</sup> ≤λ<sup>1</sup>

λ1( |Xw | −λ)

where λ<sup>1</sup> =2λ.

(22)

213

Speckle Noise Reduction in Medical Ultrasound Images

http://dx.doi.org/10.5772/56519

(23)

in each

={

Yt=*Tsemisoft*(Xw)

sign{Xw}

={

Semi soft thresholding:

lower threshold λ.

**3.2. Shrinkage rule**

*3.2.1. Visu shrinkage rule*

wavelet domain, then λ<sup>i</sup>

each threshold λ<sup>i</sup>

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 thresh‐ olding, 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 thresh‐ old 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-orkill 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:

$$\begin{aligned} \boldsymbol{Y}\_{\text{t}} &= \boldsymbol{T}\_{hard} \{ \mathbf{X}\_{\text{w}} \} \\ &= \begin{cases} \boldsymbol{X}\_{\text{w}} & \text{, for } \mid \boldsymbol{X}\_{\text{w}} \mid \ge \lambda \\ 0 & \text{, for } \mid \boldsymbol{X}\_{\text{w}} \mid < \lambda \end{cases} \end{aligned} \tag{21}$$

Soft thresholding:

#### Speckle Noise Reduction in Medical Ultrasound Images http://dx.doi.org/10.5772/56519 213

$$\begin{aligned} \mathbf{Y\_t =\_{soft}} \{ \mathbf{X\_w} \} \\ = \begin{cases} \text{sign}\{ \mathbf{X\_w} \} (\mid \mathbf{X\_w} \mid -\lambda \rangle & \text{, for } \mid \mathbf{X\_w} \mid \geq \lambda \\ 0 & \text{, for } \mid \mathbf{X\_w} \mid < \lambda \end{cases} \end{aligned} \tag{22}$$

Semi soft thresholding:

(,) () () *<sup>D</sup>*

*xy x y* <sup>=</sup> (20)

 yy

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

There are two approaches to perform the thresholding after computation of the wavelet coefficients, namely, subband thresholding and global thresholding [25]. In subband thresh‐ olding, 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 thresh‐ old 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.

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-orkill 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

denotes the output wavelet coefficients after thresholding, we define the following

(21)

*Y*t=*Thard* (*X* w)

={*<sup>X</sup>* <sup>w</sup> , for <sup>|</sup> *<sup>X</sup>* <sup>w</sup> <sup>|</sup> ≥λ 0 , for | *X* <sup>w</sup> | <λ

y

and Biorthogonal wavelets [24].

212 Advancements and Breakthroughs in Ultrasound Imaging

**3.1. Thresholding techniques**

*Shrinkage scheme:*

and Yt

thresholding functions:

Hard thresholding:

Soft thresholding:

$$\begin{aligned} \mathbf{Y}\_{\text{t}} &= \mathbf{T}\_{\text{semisoff}}(\mathbf{X}\_{\text{w}})\\ &= \begin{cases} 0 & \text{, for } \qquad | \, \mathbf{X}\_{\text{w}} \mid \leq \lambda \\ \text{sign}\{\mathbf{X}\_{\text{w}}\} \frac{\lambda\_{1} \{ | \, \mathbf{X}\_{\text{w}} \mid -\lambda \}}{\lambda\_{1} - \lambda} & \text{, for } \quad \lambda < | \, \mathbf{X}\_{\text{w}} \mid \leq \lambda\_{1} \\ \text{,} & \text{, for } \qquad | \, \mathbf{X}\_{\text{w}} \mid > \lambda\_{1} \end{cases} \\ &\text{where } \quad \lambda\_{1} = 2\lambda. \end{aligned} \tag{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 λ<sup>i</sup> to be proportional to the square root of the local noise variance σ<sup>2</sup> in each subband of the ultrasound image after decomposition. If Nk is the size of the subband in the wavelet domain, then λ<sup>i</sup>

$$
\lambda\_{\mathbf{i}} = \sigma \sqrt{2 \log \left( \aleph\_{\mathbf{k}} \right)} \tag{24}
$$

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

$$\sigma^2 = \frac{1}{\text{N}\_{\text{k}}^2} \sum\_{\text{i}, \text{j}=0}^{\text{N}\_{\text{k}}-1} \text{W}\_{\text{i}\text{j}}^2 \tag{25}$$

mean and unit variance. In the Eq.(26), if (σ<sup>n</sup> /σ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 (σ<sup>n</sup> /σx)>> 1, the noise dominates the signal. The normalized threshold

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 decom‐ position 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

2 ≥ σ<sup>y</sup> 2 , σx will

http://dx.doi.org/10.5772/56519

215

Speckle Noise Reduction in Medical Ultrasound Images

is chosen to be large to remove the noise more aggressively. In the Eq.(28), if σ<sup>n</sup>

become zero i.e. λ becomes infinity. Hence, for this case, λ =max ({| Wij |}).

to be improved interms of image quality assessment parameters.

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

**3.3. Laplacian pyramid transform**

sampling the low pass channel only.

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

$$
\lambda\_{\mathbf{k}} = \frac{\sigma\_{\mathbf{n}}^2}{\sigma\_{\mathbf{x}}} \tag{26}
$$

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

$$\sigma\_{\rm n} = \frac{\text{median}\left(\left<\mathcal{W}\_{\vec{\rm q}} \text{erHH}\_1\right>\right)}{0.67452} \tag{27}$$

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

$$\sigma\_{\mathbf{x}} = \sqrt{\max\left(\sigma\_{\mathbf{y}}^2 \cdot \sigma\_{\mathbf{n}}^2, 0\right)}\tag{28}$$

and *σ <sup>y</sup>* <sup>2</sup> , an estimate of the variance of the Wk, Since Wk is modeled as zero mean, *σ<sup>y</sup>* <sup>2</sup> can be found empirically by,

$$\sigma\_{\mathbf{y}}^{2} = \frac{1}{\mathbf{N}\_{\mathbf{k}}^{2}} \sum\_{i,j=0}^{N\_{\mathbf{k}}-1} \mathbf{W}\_{ij}^{2} \tag{29}$$

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 (σ<sup>n</sup> /σ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 (σ<sup>n</sup> /σ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 σ<sup>n</sup> 2 ≥ σ<sup>y</sup> 2 , σ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 decom‐ position 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**

The estimated local noise variance, σ<sup>2</sup>

214 Advancements and Breakthroughs in Ultrasound Imaging

threshold level λ.

and *σ <sup>y</sup>*

found empirically by,

*3.2.2. Bayes shrinkage rule*

of the empirical wavelet coefficients at the highest resolution scale as

threshold on a given subband Wk of the image X, is given by

N 1 <sup>k</sup> <sup>2</sup> <sup>2</sup> i,j 0 <sup>k</sup>

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

In Bayes shrink [29], an adaptive data driven threshold is used for image denoising. The

2 n <sup>K</sup> <sup>x</sup>

where *σn*, the estimated noise variance found as the median of the absolute deviation of the

({ ij <sup>1</sup>})

Î

This estimator is used when there is no a priori knowledge about the noise variance. The σx,

<sup>2</sup> , an estimate of the variance of the Wk, Since Wk is modeled as zero mean, *σ<sup>y</sup>*

N 1 <sup>k</sup> <sup>2</sup> <sup>2</sup> <sup>y</sup> <sup>2</sup> ij i,j=0 <sup>k</sup> 1 σ= W N


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

median W HH <sup>σ</sup> 0.67452

2 2 σ max σ -σ ,0 <sup>x</sup> y n

diagonal detail coefficients on the finest level (sub band HH1), is given by

n

which is the estimated signal variance in the wavelet domain, is given by

2 1 σ W N ij - =

, in each subband is obtained by averaging the squares

<sup>=</sup> <sup>å</sup> (25)

<sup>σ</sup> <sup>λ</sup> <sup>σ</sup> <sup>=</sup> (26)

<sup>=</sup> (27)

<sup>æ</sup> <sup>ö</sup> <sup>=</sup> <sup>ç</sup> <sup>÷</sup> <sup>è</sup> <sup>ø</sup> (28)

<sup>å</sup> (29)

<sup>2</sup> can be

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.


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


**Table 3.** Performance comparison of wavelet based despeckling method and Wiener filter.

(a) Original ultrasound image (b) Wavelet filter(Bior 6.8,L3) (c) Wiener(3x3)

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.

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

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

> <sup>2</sup> <sup>2</sup> <sup>2</sup> j J <sup>1</sup> X dc *<sup>J</sup> j*=

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

A speckle reduction method based on Laplacian pyramid transform for medical ultra‐ sound 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 con‐ trast. For every difference signal of N-level of Laplacian pyramidal decomposition a threshold value is calculated using Bayes' shrinkage rule. Further, thresholding is per‐ formed to reduce speckle. The exponential operation is performed on the filtered output

. Then, we have

= + <sup>å</sup> (30)

Speckle Noise Reduction in Medical Ultrasound Images

http://dx.doi.org/10.5772/56519

217

It can be used for the selection of thresholds for denoising applications.

results.

*3.3.1. Laplacian pyramid scheme*

dj, j=1, 2,...,J and a coarse approximation image cJ

and, hence the name the Laplacian pyramid.

to obtain the despeckled image.

**Figure 6.** Despeckled images using Wiener filter and Wavelet filter

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*

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

(a) Original ultrasound image (b) Wavelet filter(Bior 6.8,L3) (c) Wiener(3x3)

**Global thresholding (I) Subband thresholding (II)**

**Global thresholding (I) Subband thresholding (II)**

0.0431 0.0028 15.94 0.993 25.52 0.0392 0.0029 16.15 0.992 27.33 0.0461 0.0032 14.12 0.992 24.92 0.0372 0.0023 17.79 0.994 28.78 0.0481 0.0031 14.15 0.992 25.00 0.0385 0.0032 16.17 0.993 27.51

(in Secs.)

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

Denoising methods PSNR SNR MSE Variance CC Computational time

WT-L3-ST (bior 6.8) 29.94 19.06 0.00121 0.0368 0.9954 2.80 Wiener (3x3) 32.34 21.53 0.00058 0.0489 0.9977 4.14

**Table 3.** Performance comparison of wavelet based despeckling method and Wiener filter.

**Figure 6.** Despeckled images using Wiener filter and Wavelet filter

Variance MSE SNR CC PSNR Variance MSE SNR CC PSNR 0.0396 0.0023 16.89 0.991 26.38 0.0482 0.0015 18.94 0.994 28.48 0.0396 0.0025 14.83 0.992 26.02 0.0368 0.0012 19.06 0.995 29.94 0.0370 0.0024 16.04 0.992 26.17 0.0456 0.0016 17.94 0.993 28.56

**Bayes' Shrinkage rule(1)**

216 Advancements and Breakthroughs in Ultrasound Imaging

**Hard (i) Soft (ii) Semi soft (iii)**

**Visu shrinkage rule (2)**

**Hard (i) Soft (ii) Semi soft (iii)**

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

$$\left\|\mathbf{x}\right\|^2 = \sum\_{j=1}^{l} \left\|\mathbf{d}\_j\right\|^2 + \left\|\mathbf{c}\_j\right\|^2 \tag{30}$$

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 ultra‐ sound 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 con‐ trast. For every difference signal of N-level of Laplacian pyramidal decomposition a threshold value is calculated using Bayes' shrinkage rule. Further, thresholding is per‐ formed to reduce speckle. The exponential operation is performed on the filtered output to obtain the despeckled image.

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

**Table 4.** Performance comparison of the LP transform method and the wavelet transform based despeckling method

**Despeckling methods PSNR SNR MSE Variance CC Computational**

28.44 18.23 0.00141 0.0416 0.9972 2.24

29.94 19.06 0.00121 0.0368 0.9954 2.80

**time (in Secs.)**

219

http://dx.doi.org/10.5772/56519

Speckle Noise Reduction in Medical Ultrasound Images

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

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 preserva‐ tion 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 decomposi‐ tion 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

elements that are a variety of elongated shapes with different aspect ratios [37].

follicle detection in an ovarian image [40].

LP transform based despeckling method (LP-HT-

Wavelet transform based despeckling method (WT-L3-

L1)

ST)

**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 des‐ peckling 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.
