Combination of Spatial Domain Filters for Speckle Noise Reduction in Ultrasound Medical Images

The occurrence of speckle noise in medical Ultrasound (US) images poses a big challenge to medical practitioners over last several years. Speckle noise reduces the fine details present in the images and hence make it more difficult to diagnose. In this paper, a novel method based on the combination of three spatial domain filters is presented. The output of these filters is combined on the basis of an Intensity Classifier Map (ICF) formed using Coefficient of Dispersion (CoD) parameter. Experiments were conducted on synthetic and real US images. Quantitative and qualitative analysis of the proposed method is carried out in comparison to other six existing methods. It has been found from the obtained results that proposed method delivers observable improvement in all quantitative parameters undertaken for synthetic US image and in MVR value for real US images. Also, the effectiveness of the proposed method is found to be consistent with qualitative assessment of the denoised images.


Introduction
Ultrasound (US) imaging is very popular among the various imaging technologies used for the clinical diagnosis of a class of diseases.The popularity of this imaging modality lays over the advantages, such as low resolution, non-invasiveness, low cost, portability and real time operation [1].US images are generated by taking projections of received echo signal by transmitting US waves inside the human body [2].In spite of having many advantages, US image quality is deterio-rated due to a multiplicative type of noise called speckle noise.Speckle is having a granular type of pattern and is very obvious in medical US images.This noise is generated due to interference between transmitted and reflected echo signal.Speckle presents a challenge to the medical practitioners in making subjective assessment of ultrasound findings difficult.Therefore, image processing methods are used to mitigate the effect of this noise in US images.These methods are generally divided in two classes: (i) spatial domain techniques and (ii) transform domain techniques.A large number of spatial domain techniques are available in literature for speckle noise reduction in US images.Lee [3] and Frost [4] filters preserve the image edges while denoising by choosing different masks in homogeneous and non-homogeneous regions.Kuan filter [5] is similar to Lee filter but was presented with different diffusion equations.Although above mentioned filters have good denoising capabilities, performance of these filters is limited to size and shape of the used mask.
For preserving edges while denoising, an isotropic diffusion technique was first proposed by Perona and Malik, referred to as Perona Malik Anisotropic Diffusion filter (PMAD) [6].However this method works well with multiplicative speckle noise, it delivers significantly good results for images with additive noise.Speckle Reducing Anisotropic Diffusion (SRAD) [7] is a diffusion based filter which uses Partial Differential Equation (PDE) and Minimum Mean Square Error (MMSE) related to Lee filter.In SRAD filter, edges are preserved by allowing diffusion in homogeneous areas while restricting it across edges.Oriented SRAD filter [8] uses the local directional variance of the image gray levels for improving the performance of SRAD filter.Although diffusion filters are good in preserving edges while denoising, they employ iterative process for their operation and hence optimal number of iterations is required to preserve fine details present in the image.Non-Local Means (NLM) [9] filters are another class of important filters which perform denoising by weighted averaging of all the pixels in the noisy image.Coupé et al. [10] introduced a new NLM based filter in Bayesian framework for speckle reduction called Optimized Bayesian Nonlocal Mean filter (OBNLM).Some other popular spatial domain despeckling techniques are bilateral filtering [11] and detail preserving anisotropic diffusion (DPAD) filtering [12].Over last few years, due to availability of high speed processors, hybrid techniques are gaining prominence.Although these techniques are not computationally efficient, they are becoming popular because of their ability to provide better quality of denoised US images while preserving edges.A method based on the combination of SRAD and OBNLM filter is presented in [13].
In this paper, a hybrid method based on the combination of bilateral, DPAD and OBNLM filters is presented.The output obtained from the filters is combined on the basis of information obtained from an Intensity Classifier Map (ICM).ICM is formed by using statistical parameter, Coefficient of Dispersion (CoD).CoD is the ratio of variance to mean.The major contribution of the paper lies in the development of framework to obtain ICM which classifies three different regions viz.homogeneous, detail and edge region.ICM is generated by comparing local CoD (T(i,j)) with the estimated values of threshold T1 and T2.T1 and T2 are boundary threshold values for homogeneous-detail region and detail-edge region respectively.ICM so obtained is used to selectively combine the output of three different spatial domain filters.Section 2. describes the noise model and spatial domain filters used for the proposed method.Section 3. describes the proposed method which is a hybrid method based on CoD.In Sec. 4. experimental results are illustrated for proposed method and are compared with six other state of art methods.The comparison is done on the basis of quantitative and qualitative performance evaluation.Finally, conclusion of the proposed work is presented in Sec. 5.

Noise Model
The speckle noise corrupting the US image can be treated as multiplicative in nature and can be expressed as in Eq. (1) given below [14].
where g(x,y) is the noisy version of original image f (x,y) corrupted by noise n(x,y).The multiplicative speckle noise can be converted in additive noise by taking log transformation of the noisy image and can be represented as Eq.(2).
where g' (x,y), f ' (x,y) and n' (x,y) are the log transformed version of g(x,y), f (x,y) and n(x,y) respectively.

Filters Used
Filters used in the proposed work are presented in this section for the sake of illustration:

1) Bilateral Filter
Bilateral filter [11] performs in spatial domain and is used for non-linear filtering of image pixels.This filter is mainly used because of its ability to smooth images while preserving edges.In Bilateral filter, a combination of two Gaussian filters (i.e.domain and range filters) is used to obtain the output.The response of the filter at location 'a' can be obtained as given in below equations: and where N s (a) is the "spatial neighborhood" and 'A' is the "normalization constant".Also, σ d and σ r are geometric spread and photometric spread parameters respectively.The choice of parameter values σ d and σ r is very crucial and chosen as a function of noise variance.

2) Detail Preserving Anisotropic Diffusion (DPAD)
In DPAD filter [12], the diffusion equations are formed according to Kuan filter.The diffusion equations in DPAD filter are given as: and where K i,j,t 3 is the gain and the parameter | ηs | is the no. of pixels in the neighborhood.C i,j is local coefficient of variation of image, c 2 (.) is the diffusion coefficient, and C 2 u is the coefficient of variation of noise.

3) Optimized Bayesian Non Local Mean Filter (OBNLM)
OBNLM [10] filter is mainly used to reduce the computational complexity involved with classical NL-mean method.This approach can be divided in the following steps: • Firstly, a "search area" of size (2M + 1) 2 is formed which is centered at the current pixel x i .
• The image is divided into blocks with overlapping supports of size V = (2D + 1) 2 .
• NL-means filter is applied on each block.
• The block B ik is restored according to the equations described below: and where Z ik is the normalization factor.The restored block of current pixel x i is represented as NL(u)(B ik ) and h is the smoothing parameter which controls the decay of exponential function.
The above process is repeated for all the pixels of speckled image.Finally, the restored pixel x i is obtained by taking the average of the values stored in the vector G i presented as:

Proposed Work
In the proposed work, a method based on the combination of bilateral, DPAD and OBNLM filters is presented.This combination works according to a map called Intensity Classifier Map (ICM).The proposed method is derived in three steps: (i) formation of ICM (ii) obtainment of three reconstructed images using bilateral, DPAD and OBNLM filters (iii) combination of three reconstructed output images using ICM.

Formation of Intensity Classifier Map (ICM)
Intensity classifier Map (ICM) is formed using CoD.
The detailed process of formation of ICM is as follows: • Firstly, local CoD (T(i,j)) is obtained by using a moving window of size W×W in image neighborhood.
• For pixel classification in different regions viz.homogeneous, detail and edge region, two thresholds T1 and T2 are defined.T1 is used to separate the homogenous and detail region pixels in the noisy image while T2 is used for separating detail and edge region pixels.
• T1 is the CoD value obtained from a chosen homogeneous region in the noisy image.
• T2 is the CoD value obtained from the computation of an edge image.The edge image is obtained by the reconstruction of single level detail band coefficients obtained from single level 2D wavelet decomposition.
• A pixel is classified as homogeneous region pixel if the value T(i,j) is less than T1.If the value T(i,j) is greater than T2, then pixel lies in edge region.Otherwise, the pixel lies in intermediate region or detail region.ICM consists of black, gray and white pixels showing homogeneous, detail and edge regions respectively.Figure 1 shows the noisy image and corresponding ICM obtained after applying the abovementioned procedure.The flow chart of the proposed method is shown in Fig. 2. Three filters used in proposed method are Bilateral, DPAD and OBNLM which are used for speckle denoising of homogeneous, detail and edge regions respectively.

Quality Assessment Parameter
In this paper, for performance evaluation of the proposed method and six other methods, parameters such as Peak Signal to Noise Ratio (PSNR) [15], Signal to Noise Ratio (SNR) [15], Edge Keeping Index (EKI) [16], Structure Similarity Index Measures (SSIM) [17] and Pratt's Figure of Merit (FoM) [18] are used for synthetic US image.However, for real US images, Mean to Variance Ratio (MVR) [19] parameter is used for quality assessment.

Experimental Results
The performance of the proposed method is analysed over both synthetic and real US images quantitatively and qualitatively.The quantitative evaluation is done on the basis of image quality assessment parameters listed in Subsec.4.1.The comparative analysis of the proposed method is carried out against six existing methods viz.PMAD, SRAD, OBNLM, bilateral filter, DPAD and a combination of SRAD & OBNLM filters [13].

1) Experiment on Synthetic US Image
Experiments are performed on synthetic US image of kidney generated using field program II simulator [20].Speckle noise at different noise variances (σ = 0.1, 0.2 and 0.3) is introduced artificially in synthetic US image using speckle model given in MATLAB.
For PMAD filter, results are obtained for 10 iterations (σ = 0.1) and 20 iterations (σ =0.2 & 0.3).The chosen parameters for SRAD filter for experiment are 200 iterations each for σ = 0.1, 0.2 & 0.3 for time step size ( t) = 0.7, 1 and 1.2 respectively.For OBNLM filter smoothing parameter (h) is selected as 1.7, 1.9 and 2.2 for σ = 0.1, 0.2 and 0.3 with search area and patch size of each 5×5 with 100 iterations.In case of bilateral filter the values of σ d = 2, 2.2 and 2.6 and σ r = 2 √ σ are chosen for σ = 0.1, 0.2 and 0.3 respectively.Also window size of 3×3, 5×5 and 5×5 are chosen for σ = 0.1, 0.2 and 0.3.For, DPAD filter, time step size ( t) = 0.2, 0.75 and 0.8 is chosen for 100 iterations at σ = 0.1, 0.2 and 0.3 respectively.For method presented in [13] the parameter sf is selected as 0.3 and the window size of 11×11 is chosen.Filter parameters chosen for all the methods described above are found to be optimal values when these filters are used standalone or in combination.These optimal values are obtained after performing several experiments on synthetic images.A window size of 11×11 is chosen in proposed method for calculating the local CoD.
Table 1 highlights the values of image quality metrics viz.PSNR, SNR, EKI, SSIM and FOM for synthetic kidney US image corrupted with noise variance (σ = 0.1, 0.2 and 0.3).It may be noted from Tab. 1 that the performance of the proposed method is superior to all other methods in terms of all quality assessment parameters over a significant range of noise variance.

Methods
Noise     3(h) presents synthetic kidney images for qualitative analysis.It is evident from visual assessment of denoised images that quality of the image obtained using proposed method is more pleasing in terms of detail and edge preservation.Also, the obtained reconstructed image is very similar to original noise free image.The image obtained using PMAD and SRAD filter is the worst among all images and the images obtained from bilateral filter and OBNLM filter are of moderate quality.Images generated using DPAD and SRAD-OBNLM combination also give good quality images as compared to PMAD, SRAD, bilateral and OBNLM filter but not comparable to the one obtained from proposed method.

2) Experiments on Real US Images
Experiments are also conducted on a set of 50 real US images obtained from an online database [21].Qualitative evaluation is generally done for real US images.ance ratio of that region [19].A higher value of MVR is desirable for a superior denoising method.
Figure 5(a) shows a US image of real pancreas and Fig. 5(b) shows that of real gallbladder.Three different regions are highlighted in red, green and blue color in both images which may be abbreviated as Region 1, Region 2 and Region 3 respectively.MVR values are calculated for all three highlighted regions of abovementioned two images and are abbreviated as MVR1, MVR2 and MVR3 respectively.These values are calculated for all methods and are compared (Fig. 6).From Fig. 6(a) and Fig. 6(b), it is worth highlighting that a higher value of MVR is obtained for the proposed method as compared to other methods in all three regions for both images.Moreover, a total of 100 MVR measurements are taken from different regions of 50 different US images for evaluating the performance of proposed method.
In Tab. 2, MVR mean values along with standard deviation for 100 measurements for all the methods are presented.The results of Tab. 2 show the superiority of the proposed method over other methods in terms of MVR values obtained for each method.
Tab. 2: Comparison of MVR values for different methods.

Conclusion
In this paper, a new despeckling method based on the combination of bilateral, DPAD and OBNLM filters is presented.CoD parameter is used for image pixel classification in different areas of noisy images.ICM is a map used to combine the output of three images.
From the quantitative results obtained, it can be seen that proposed method gives the high value of quality assessment parameters for synthetic US image.Moreover, it also shows the superiority in terms of MVR value obtained for real US images.The potential of the proposed method can also be seen from visual inspection of denoised images.Therefore, the presented work gives better performance than existing methods in terms of speckle noise suppression and edge preservation.
a and b are spatial distances.Coefficients corresponding to domain filter are proportional to the spatial distance (b -a) and for range filter coefficients are given by (s(b) -s(a)).D(a,b) and R(a,b) denotes the domain and range filter components.D(a,b) and R(a,b) are represented as in given equations: (a) Noisy Image (σ = 0.1).(b) Intensity classifier map.

1 :
DIGITAL IMAGE PROCESSING AND COMPUTER GRAPHICS VOLUME: 15 | NUMBER: 5 | 2017 | DECEMBER Tab.Comparison of various image quality metric for synthetic US image of kidney for noise variance (σ = 0.1, 0.2, 0.3).

Figure 3 (
Figure 3(a), Fig. 3(b), Fig. 3(c), Fig. 3(d), Fig. 3(e), Fig. 3(f), Fig. 3(g) and Fig.3(h) presents synthetic kidney images for qualitative analysis.It is evident from visual assessment of denoised images that quality of the image obtained using proposed method is more pleasing in terms of detail and edge preservation.Also, the obtained reconstructed image is very similar to original noise free image.The image obtained using PMAD and SRAD filter is the worst among all images and the images obtained from bilateral filter and OBNLM filter are of moderate quality.Images generated using DPAD and SRAD-OBNLM combination also give good quality images as compared to PMAD, SRAD, bilateral and OBNLM filter but not comparable to the one obtained from proposed method.

Figure 4 (
a), Fig. 4(b), Fig. 4(c), Fig. 4(d), Fig. 4(e), Fig. 4(f), Fig. 4(g) and Fig. 4(h) presents the real breast US images for visual quality assessment.It is worth noting that reconstructed image obtained using proposed method is the best in terms of simultaneous denoising and fine edge preservation.Quantitative evaluation is difficult in real US images as no noise free (reference) image is available.However, MVR is considered as a reasonable parameter which gives good estimation of speckle level in denoised images.MVR can be calculated in different regions of an image by selecting the region and calculating the mean to vari-(a) Noisy Image.(b) PMAD.(c) SRAD.(d) OBNLM.(e) Bilateral.(f) DPAD.(g)SRAD+OBNLM.(h)Proposed method.

Fig. 4 :
Fig. 4: Real breast US reconstructed images obtained from various denoising methods.
(a) Pancreas real US image.(b) Gallbladder real US image.

Fig. 5 :
Fig. 5: MVR calculation in three regions of real US images.
MVR plot for gallbladder real US image.

Fig. 6 :
Fig. 6: MVR calculation in three regions of real US images.