Comparison of Reducing the Speckle Noise in Ultrasound Medical Images using Discrete Wavelet Transform

—Speckle noise in ultrasound (US) medical images is the prime factor that undermines its full utilization. This noise is added by the constructive / destructive interference of sound waves travelling through hard-and soft-tissues of a patient. It is therefore generally accepted that the noise is unavoidable. As an alternate researchers have proposed several algorithms to somewhat undermine the effect of speckle noise. The discrete wavelet transform (DWT) has been used by several researchers. However, the performance of only a few transforms has been demonstrated. This paper provides a comparison of several DWT. The algorithm comprises of a pre-processing stage using Wiener filter, and a post-processing stage using Median filter. The processed image is compared with the original image on four metrics: two are based on full-reference (FR) image quality assessment (IQA), and the remaining two are based on no-reference (NR) IQA metrics. The FR-IQA are peak signal-to-noise ratio (PSNR) and mean structurally similarity index measure (MSSIM). The two NR-IQA techniques are blind pseudo-reference image (BPRI), and blind multiple pseudo-reference images (BMPRI). It has been demonstrated that some of these wavelet transforms outperform others by a significant margin.


I. INTRODUCTION
An ultrasound (US) medical image helps in an early diagnosis of kidney stones.These stones cause severe pain in situations where they become large or block the flow of urine.In rare situations, a small stone is stuck in the ureter.The ureter is a small tube connecting kidney and bladder.As per statistics, 1 in 11 persons in USA suffer from kidney stones [1].An early treatment can save someone from severe pain, cost and medical complexities.The US imaging is quick, noninvasive, cost effective, and has no known side effects to the best of our knowledge.The medical complexities involving kidney stones are considered high-risk illnesses.These are life threatening if left untreated for a long time.
The size and location of kidney stone is also important.A patient may experience no symptom at all to severe, incapacitating pain in the loin requiring urgent treatment.At times the patient is complaining burning sensation during passage of urine, blood, or small stone debris in their urine.Stones can block the main outflow of urine from the kidney leading to irreversible kidney damage, disturbances in the biochemical balance of the body and eventually death.Thus a safe, economical and rapid imaging technique requiring no prior preparation is invaluable in saving many lives.
The practical issue of speckle noise is not new.This has been addressed by several researchers, some as early as 1980"s.Jain had initially approximated the speckle noise as multiplicative [2].He suggested to apply homomorphic filter.The pre-and post-processing were performed using Wiener, and the Median filter, respectively.Chien-Min used Bayesian approach for removing the speckle noise [3].Perona and Malik introduced an edge preserving approach of anisotropic diffusion (AD) in [4].The AD filter was subsequently used by several researchers for speckle noise reduction [5]- [8].The estimation of signal and noise using Kuan"s filter has been used for speckle reduction anisotropic diffusion (SRAD) in [5].SRAD using filtering across image contours and principle curvature directions are given in [6].The probabilistic modelbased SRAD is discussed in [7].A comparison of SRAD and several other schemes are presented in [8].
The introduction of wavelets during the early 90"s resulted in several papers on speckle noise.Mallat introduced multichannel decomposition of images using wavelet transform in [9].The wavelet theory for image coding was developed in [10].The application of wavelet packets was presented in [11].The Bayesian maximum a posteriori (MAP) estimator based design were illustrated in [12]- [13].The rational-dilation wavelet transform (RADWT) and non-linear bilateral filter based approaches were proposed in [14].The genetic algorithm based solution was introduced in [15].A quantum-inspired de-speckling method was discussed in [16].The wavelet and fuzzy theory based approach is given in [17].
During the last two decades, several new wavelets have been derived.These wavelets offer numerous benefits.Some of them are preferred over the others in terms of number of stages, linearity, ease of use, etc.The performance of Symlet wavelets has been discussed in [18].A comparison of five wavelets Haar, Daubechies, Symlet, Coiflet and biorthogonal www.ijacsa.thesai.orgwavelets for removing the speckle noise has been given in [19].The applications of various wavelets for identification of bone fracture has been discussed in [20].A comparative study of Birge-Massart strategy for setting a threshold for image compression is given in [21].The identification and classification of colonic polyps using wavelets transforms has been demonstrated in [22].This paper compares the performance of seven wavelets for reducing the speckle noise in US medical images.The selected discrete wavelets are Haar, Daubechies, Symlet, Coiflet, biorthogonal, reverse biorthogonal, and discrete Mayer wavelets.The pre-and post-processing is performed by using Weiner and Median filters, respectively.The introduction is followed by evaluation criteria in Section 2. The description of selected discrete wavelet transforms is given in Section 3. The performances of these wavelets are given in simulations in Section 4. Section 5 concludes this paper.

II. EVALUATION CRITERIA
There are generally two ways of image quality assessment (IQA).The first approach compares the results with the original image.This is termed as full-reference (FR) IQA [23].The second approach is more recent in which the assumption is that no reference image is available.This is referred to as blind or no-reference (NR) IQA [24]- [25].In case of an US image, a truly clean image without speckle noise is not really available.In this regard, the NR-IQA seems to be a better metrics for US images.

A. Full-Reference IQA (FR-IQA)
The performance of various wavelets is tested using two FR-IQA.The first is based on peak signal to noise ratio (PSNR) and the second is mean structurally similarity index measure (MSSIM).The mean square error (MSE) is used as a criterion for comparing the original image with the processed image as given by where, is the pixel value of the original image and ̂ is the estimated value of .The image row and column numbers are given by M, and N, respectively.The peak signal-to-noise ratio in decibels PSNR dB is found by using, √ The superscript B is the number of bits in a pixel.The value of B is taken as 8, resulting in 256 grey shades.Each pixel value therefore varies in the range of 0-255.
The structurally similarity index measure (SSIM) was proposed in 2004 [23].It compares the mean and variance of two images.The two images are considered by x and y.The SSIM is given as, where are the average of images x and y, respectively.
are the variances of images x, and y, respectively.1 = ( 1 ) 2 , 2 = ( 2 ) 2 are two variables to stabilize the division with weak denominator.L is the dynamic range of the pixel value 255.The k 1 and k 2 are equal to 0.01 and 0.03 (taken by default), respectively.Another single parameter used for images is MSSIM given as [23].
The SSIM is the similarity index calculated over a small region of an image.The MSSIM is the mean value of SSIM across all windows.The MATLAB code for MSSIM is available in [26].

B. No-Reference IQA (NR-IQA)
The no-reference (NR) IQA assumes that a clean image is not available, and the image quality is assessed based on the available noisy image.In this paper two NR-IQA techniques are used to test the quality of US images.Both the techniques are essentially based on the use of a pseudo-reference image (PRI).A PRI is generated from the distorted image.In the absence of a clean image, the PRI is used as a reference image for the targeted US image quality assessment.
The first NR assessment technique is blind PRI (BPRI) [24].This technique measures the image quality in terms of block effects, sharpness, and noise.The block effect is found by using pseudo structure similarity (PSS) index.The local binary pattern (LBP) is used for sharpness, and noise measurement.The PRI-based sharpness and noise is used to derive the local structure similarity (LSS) index.The second NR assessment technique is multiple PRI (MPRI) [25].In this technique, the image is distorted with an aggregation of four types of distortions, based on JPEG compression, JPEG2000 compression, Gaussian blur (GB), and white Gaussian noise (WGN).

III. WAVELET SELECTION
The Fourier and Laplace transforms have been used extensively for extracting the significant frequency components of a noisy image.These transforms perform very well by translating the time domain signal into frequency domain.The main limitation of these transforms is that the local information is lost.In an analogue or digital signal transmission system, the location of noise is generally not very critical.This is different in image processing, where the perceived quality of an image depends on the location of noise.As an example, the loss of signal at the edges of an image can be acceptable, but the loss of fine details at the centre of an image, or around the critical regions is quite unacceptable.A significant advantage of wavelet transform over the previously available transforms is that the wavelet not only translates time-domain into frequency-domain, it also preserves the physical location of noise present in an image.This advantage has no parallel in the other transforms.
During the last two decades, several wavelets have been derived with different set of properties.Before proceeding on wavelets, it is important to review few basics of wavelet www.ijacsa.thesai.orgtransforms.An image decomposed by wavelet has two functions.These are scaling function and wavelet function, sometimes refer to as the mother wavelet, and the father wavelet.The scaling component gives lower frequency components corresponding to variations in the grey shades.The wavelet function gives high frequency components like edges.The scaling function, and the wavelet function, are given by, The and are low-pass and high-pass filters, respectively.The n is the periodic shift that implements the filter coefficient index.Both filters are related by, (6) where N is the number of vanishing moments.A wavelet with N vanishing moment has at least a polynomial of order N-1.The vanishing moments represent level of differentiability of a function.A wavelet with vanishing moment N is defined as multi-scale differential of order N. This, in essence, defines the local irregularity of a signal.A smaller value of N is therefore preferred over the larger values.An N vanishing moment corresponds to 2N taps in the filter bank.The filter is implemented as a finite-impulse response (FIR) filter.A smaller value of N, therefore, corresponds to a shorter filter with less number of taps.Haar wavelet is the oldest and the simplest of all wavelets.This is the only orthogonal wavelet having linear phase.Haar wavelet decomposes the discrete signal into two sub-signals of half its length.One sub-signal provide the trend, while the second sub-signal gives difference or fluctuations.The main advantage of Haar wavelet is that it is fast, memory efficient, and conceptually simple to implement.
The biorthogonal wavelets have linear phase.These filters have a pair of scaling functions and an associated scaling filters used for analysis and synthesis.The analysis and synthesis filters can be designed to have different order of vanishing moments.It is possible to have greater number of vanishing moments for sparse representation analysis and a smoother wavelet for reconstruction.In MATLAB notation, the biorthogonal wavelets are designated as "biorN r .N d ".Similarly, the reverse biorthogonal wavelets are represented as "rbioN r .N d ".The "N r " represents the effective number of reconstruction filter, and the "N d " represents effective number of decomposition filters.Table I gives several choices of "N r" and "N d ".
The Daubechies wavelet has several versions represented by vanishing moments, N. In Symlet wavelet, the value of "N" varies from 1 to 10.The Coiflet wavelet has 5 variations represented by the value of "N" that equals 1 to 5. In all the above wavelets, the number of taps in synthesis and analysis filters are same.The discrete Meyer wavelet has a single transform.A comprehensive mathematical analysis of wavelet theory is given in [27].
The Haar and Daubechies wavelets are orthogonal wavelets, while biorthogonal and reverse biorthogonal wavelets are biorthogonal in nature.The discrete Meyer wavelet is simply a discrete version of the continuous Meyer wavelet.A compactly supported wavelet function restricts itself to within certain limits and as a consequence the signal is also restricted to within some limits.The Haar, Daubechies, and reverse orthogonal wavelets have compactly supported functions.Both the Daubechies and reverse biorthogonal wavelets show an arbitrary number of vanishing moments.All the selected wavelets have finite impulse response (FIR).The FIR has an advantage of having only a few non-zero coefficients.

IV. SIMULATIONS
The algorithm is tested on six US medical images of abnormal kidneys with stones.These images are downloaded from the US imaging database [28].The database comprises of large collection of US images that are categorized as fetal, kidney, renal calculi, appendix, urinary bladder, liver, spleen, chest and vascular system.The images in each category consists of high resolution samples.The low resolution images are good for fast processing; however, they are not appropriate as the fine details are lost.Also, the outcome may have little practical value.The high-resolution images provide enough details but they require more processing power as well as time.If these images are being used during an operation, then the speed of processing needs to be sufficiently higher to give real-life response to the ongoing activities, like surgery and other diagnostic treatments.www.ijacsa.thesai.orgThe selected kidney images are US images of patients complaining about minor to severe pain in the left region of their abdomen.All samples are for male patients.The initial diagnosis recommended US imaging for further treatment.The US images clearly showed stones, but the number of stones, and their sizes are not clearly identified as the images contain significant amount of speckle noise.The stones are more identifiable if they are either in larger size or present close enough to form a larger area in concerned region.Unfortunately, in most cases the kidney stones have relatively smaller sizes, and they are scattered across the whole active region of a kidney.In US images, kidney stones usually appear in "white" or lighter grey shades.A distinct feature of the stones is that they always have a long "shadow" that originates from the stone, and spreads out towards the outer edge of the kidney.These shadows are quite visible in most of the US images; however, in few cases the shadow is much lighter and may be overlooked.Image enhancement helps in improving the visibility of these shadows in such situations.In general, there are two primary objectives.The first is to reduce the speckle noise.The second is to improve the contrast level of regions containing significant amount of information like stones.
Fig. 1 gives various steps involving US medical image processing.These are pre-processing, filtering, and postprocessing.In pre-processing, Wiener filter helps in smoothing the image.This has been demonstrated that although this causes sharp edges to be slightly smoothed out, but as a result the overall visibility improves.The preprocessing is followed by filtering, using several discrete wavelet transforms.The selected transforms are Haar, Daubechies, Symlet, Coiflet, biorthogonal, reverse biorthogonal, and discrete Meyer wavelets.The Haar and discrete Meyer wavelets are distinct in nature.The remaining have several combinations.A list of selected wavelets is given in Table I.The post-processing is based on Median filter that helps in restoring the sharp edges.The window size used in the Wiener filter is (3 x 3) pixels.It has been observed that this is a preferred window size than the larger window size of (5 x 5), or (7 x 7) pixels.The two-dimensional Median filter is applied using (5 x 5) pixel window.The subjective tests have supported the use of (5 x 5) pixel windows against the possible (3 x 3) pixel or (7 x 7) pixel windows.
The processing of wavelet is performed only at the first level of decomposition.The coefficients in the horizontal, vertical, and the diagonal directions are filtered using a hardlimiter as, where is the threshold.The value of in the horizontal, vertical, and the diagonal directions are selected as 200.This high value of threshold essentially removes most of the effects present in the above three directions.The low resolution coefficient cut-off is dynamically calculated by taking the mean pixel value of the selected image.It has been observed that a higher threshold value of low resolution coefficients result in severe quality degradation.
The qualitative analysis of the above algorithms is tested on six US images of various granularity levels.The unprocessed, original images are given in the left column of The quantitative analysis of the images are given in Table II through Table V.As mentioned earlier, the performance of various wavelet transforms are analyzed on four criterions broadly categorized as FR-IQA, and NR-IQA.The FR-IQA techniques are PSNR, and MSSIM, while the NR-IQA techniques are BPRI, and BMPRI.The PSNR and MSSIM metrics are given in Table II and Table III, respectively.The PSNR is given in decibels (dB), while the similarity index using MSSIM is given in a value between 0 and 1.A value closer to 1 corresponds to a higher similarity across images.It is clear that the PSNR value depends on image contents.This is the reason that the PSNR of six selected images varies.However, the PSNR of any one image processed through a wavelet has smaller variation in terms of decibels (dB).This is also observed that a higher value of MSSIM corresponds to a higher value of PSNR.The metrics of BPRI, and BMPRI techniques are given in Table IV and Table V, respectively.The BPRI index is in the range of 0 to 1, while the BMPRI is in the range of 20 to 60. www.ijacsa.thesai.orgThe highest PSNR in each of the six images is marked with bold letters.These correspond to Daubechies 24, reverse biorthogonal 3.1, Daubechies 24, biorthogonal 6.8, Daubechies 20, and biorthogonal 3.3 for successive images of Fig. 4(a)-(f).The specific wavelet is mentioned against each image in Fig. 4. The corresponding value of MSSIM, BPRI, and BMPRI are also highlighted with bold letters.The histograms of processed images are given on the right side in Fig. 4. It has been observed that the response to Haar wavelet, Symlet 1 wavelet, and the biorthogonal1.1 wavelet is same.Similarly, the performance of Daubechies 2, and Daubechies 3 are exactly same as Symlet 2 and Symlet 3, respectively.The www.ijacsa.thesai.orgmean PSNR of six images is plotted in Fig. 2 and Fig. 3.The standard deviation of these images is also plotted along with the mean values in Fig. 2 and Fig. 3.The graphs in Fig. 2 corresponds to the Haar wavelet, and Daubechies wavelets.The mean value of the remaining wavelets is given Fig. 3. From Fig. 2, it is clear that variations among the values of Daubechies wavelets are within 1 dB value, but there is a definite increasing trend of PSNR as the value of N is increased from Daubechies 2 to Daubechies 20, and then it decreases slightly afterwards.This is understandable as with more number of taps, the amount of information increases until N becomes equal to 20, and then it slightly reduce as the undesired signal is possibly added in the extracted information.
The mean value of MSSIM is given in Table III.It is clear that the similarity index increases at the higher values of N. The mean values of BPRI, and BMPRI are given in Table IV, and Table V, respectively.The peak index values of BPRI and BMPRI for each image is encircled.It is clear that Daubechies 6 out performs in four out of six images.In the remaining two images, the value of Daubechies 6 closely follows two other wavelets Symlet 7, and Coiflet 2. The BPRI index of image (c) is same for Daubechies 6 and Symlet 7. The maximum value of BMPRI is not consistent, as the peak values of six images correspond to six different wavelets.The possible reason is that BMPRI compares the original image with a noisy image that has been generated by the aggregation of several noisy version of the original image.The image content has a larger role in generating the reference image, and therefore the response to various wavelets does not generate consistent results.

V. CONCLUSIONS
This paper reviews the performance of seven discrete wavelet transforms in reducing the effect of speckle noise of the US medical images.The complete analysis is based on pre-processing, filtering, and post-processing.The preprocessing involves the application of Wiener filter, followed by filtering using seven discrete wavelet transforms.The selected wavelet transforms are Haar, Daubechies, Symlet, Coiflet, biorthogonal, reverse biorthogonal, and discrete Meyer wavelets.The post-processing involves the application of Median filter.The simulation results are based on US images of kidney that have previously been diagnosed for stones.The performance of various wavelets is compared using four metrics.Two are based on full-reference (FR) image quality assessment (IQA), and the remaining two are based on no-reference (NR) IQA.The two FR-IQA techniques are PSNR, and MSSIM, while the two NR-IQA techniques are blind pseudo-reference image (BPRI), and blind multiple PRI (BMPRI) metrics.The images with the highest PSNR are identified, and the other metrics of these images are marked.The output of the selected images with the highest PSNR, along with their histograms, are reproduced for comparison.The qualitative and quantitative analysis clearly show significant improvement in the processed images.In this paper only the first level of wavelet decomposition is considered.An extension of this work by considering multilevel wavelet decomposition has strong potential for better results.It would be interesting to see the effect of using one wavelet function for first level and another wavelet function for the second level.

Fig 1 .
Fig 1. Flow diagram of US image processing.

Fig. 4 .
The location of stones and their shadows are marked with arrows on original images.The processed images are given on the right side.The histogram of original and the processed images are given on right side of the successive image.A quick comparison of the results reveal that a significant amount of speckle noise has been removed.At the same time, the contrast of the images has also improved.The histogram of all images show that the lower pixel values (darker grey) have been filtered out.The grey shades at higher values (lighter grey) have smoothed out, but the general shape has remained same.This has resulted in making the shadow region more prominent.The shadows in original images as shown in Fig.4(a)-(e) are quite visible.These have become clearer in the processed images given on the right side.The shadows in Fig. 4(b), and f are not very clear in the original images, but they have become quite visible in the processed images.

(a) Daubechies 24 ,
Fig 4. (a-f) Original and Processed image along with their histograms processed through selected wavelet transforms.