wKSR-NLM: An Ultrasound Despeckling Filter Based on Patch Ratio and Statistical Similarity

Ultrasound images are affected by the well known speckle phenomenon, that degrades their perceived quality. In recent years, several denoising approaches have been proposed. Among all, those belonging to the non-local (NL) family have shown interesting performance. The main difference among the proposed NL filters is the metric adopted for measuring the similarity between patches. Within this manuscript, a statistical metric based on the ratio between two patches is presented. Compared to other statistical measurements, the proposed one is able to take into account the texture of the patch, to consider a weighting kernel and to limit the computational burden. A comparative analysis with other despeckling filters is presented. The method provided good balance between noise reduction and details preserving both in case of simulated (by means of Field II software) and real (breast tumor) datasets.


I. INTRODUCTION
Ultrasound (US) imaging (or sonography) is a well-known diagnostic medical modality that exploits high-frequency mechanical waves to produce dynamic visual images of organs, tissues or blood flow inside the body. These waves are transmitted to the area to be examined and the returning echoes are captured to provide an image of the investigated area due to the differences in the acoustic properties of the biological tissues.
The great interest in this medical imaging technique is motivated by the fact that it is non-invasive, harmless (due to the use of mechanical waves for the image generation), portable and low cost. Unfortunately, it is affected by multiplicative noise, known as speckle, which corrupts the images, causing a granular pattern that worsens their quality and makes their interpretation difficult. This phenomenon is mainly related to the multi-path interactions of the reflected signals and to the imaging equipment which degrades fine details in the image, making the extraction, analysis, recognition and quantitative measurements complicated and unreliable. Therefore, image despeckling is a challenging issue, The associate editor coordinating the review of this manuscript and approving it for publication was Giulia Matrone . but it represents an essential task to improve the image visual quality and allow enhanced analyses.
In the last few decades, there has been an increasing number of articles describing the fundamentals and the statistical characteristics of the multiplicative speckle noise, especially in the field of synthetic aperture radar [1]- [4] as well as in the biomedical field [5], [6]. Conversely from standard denoising approaches adopted for the additive white Gaussian noise, US imaging requires specific filters due to the peculiar nature of speckle [7].
Among the despeckling methods, a first class is represented by adaptive filters. They are often employed due to their ease of implementation and control. Some examples are represented by Lee's Filter, Frost's filter and Kuan's filter [8]- [11]. In the early 1990s, these filters have been improved by exploiting local image statistics to identify specific areas to be processed further and to adapt spatial filter coefficients better [12].
Partial differential equations (PDE) based approaches have also earned great interest to design speckle removal strategies. The design of such mathematical tools looks at the observed image as initial data for a boundary and initial value problem [13], [14]. In this framework, non-linear diffusion methods have proved to produce good spatially-regularized images, VOLUME 8, 2020 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ but with some undesirable smoothing effect on sharp edges. Speckle reducing anisotropic diffusion filter [15], non-linear coherent diffusion filter [16] and total variation minimization scheme [17] represent some examples. Another class of denoising techniques involves thresholding methods combined with some projection operators, e.g. Wavelet [18], [19] or spatial projections. This multiscale approaches can also exploit a Bayesian framework [20], [21] to enhance the filtering or to improve the recovery performance of standard adaptive filters. Moreover, hybrid approaches which try to combine the advantages of the different methodologies in order to overcome the drawbacks introduced in the recovery procedure have also been proposed [22].
All the aforementioned approaches for speckle reduction are based on a locally adaptive recovery framework, in which the filtering operation is carried out via a local comparison of pixels in the selected patches [23]. Conversely from this paradigm, the so-called non-local means (NLM) approach overcomes the drawback of the local filtering procedure and looks for non-local similar patches in a search area, that could correspond to the whole image [24]. Once the similar patches are identified, the central pixels are fused. Such a procedure is repeated for each pixel in order to obtain the regularized image.
The choice of similarity measurement, i.e the definition of a proper metric to compare the similarity among noise-corrupted patches, plays a fundamental role in such filters [25].
Initially, the NL approaches proposed in literature implemented a similarity metric based on the Euclidean distance, i.e. the energy of the deviation between patches [26]. In recent years, some alternative measures come from synthetic aperture radar (SAR) framework, such as the Pearson distance [27], [28] and stochastic distances [29], [30]. Divergence measures play a major role in statistical inference and discrimination since they are measures of the statistical distance between probability distributions. More in detail, Nascimento et al. in [29] explore the use of stochastic distances in the hypothesis testing of SAR speckled imagery, deriving the expression for the derivation of the stochastic distances of eight well-known divergences. More recent articles [31], [32] exploit the geodesic distance as metric for patch similarity measurement, in which the distance between statistical distributions is used for region discrimination. In the KS-NLM filter [33], specifically developed for US images, a statistic metric based on the Kolmogorov-Smirnov distance between the cumulative distribution functions (CDFs) of the patches has been assumed. As a drawback, in order to effectively measure the KS distance, the KS-NLM method required the acquisition of several images, i.e. several frames of a video. Within this manuscript, we propose an improved similarity metric able to discern among areas characterized by close statistical distributions but different spatial textures, which is the main limit of the KS-NLM. Moreover, it required less pixels for the analysis, allowing the denoising of a single image. More in detail, the proposed approach analyzes the ratio among the patches, and is able to include a weighting kernel in the similarity evaluation. The adoption of a weighting kernel in measuring the distance between patches has shown effective improvements in case of Euclidean distance metric [26]. In brief, such a kernel adjusts the contribution of each pixel depending on its position within the patch [25]. Unfortunately, the adoption of weights based on the position is not trivial in case of KS distance, as such a metric is not able to take into account pixel position. Within this manuscript, we proposed an idea for adopting a weighting kernel in case of statistical distances, such as the KS one. In brief, we propose to apply weights to the computation of the statistical moments (in particular mean and standard deviation), instead of in the estimation of CDF.
The reminder of the paper is organized as follows. In Section II, we provide a general overview of the mathematical formulation of proposed approach, proving the feasibility of the method. Section III presents results of a comparison with other state of art despeckling filters on both numerical as well as real datasets. Finally, some conclusions are drawn in Section IV.

II. MATHEMATICAL FORMULATION A. SPECKLE STATISTICAL CHARACTERIZATION AND NON-LOCAL MEANS
The US signal intensity z (x, r) follows a multiplicative noise model [33]: where y (·) is the noise-free signal, n (·) is the speckle noise, (x, r) are the space indices related to the acquisition geometry defined in the domains (X , R) respectively. Considering a resolution cell with randomly-distributed high-density scatterers and focusing on a single spatial position (x 0 , r 0 ), the the signal intensity z = y · n (i.e. the square modulus of the acquired signal) for this kind of regime is well described by the exponential probability density function (PDF) [6], [34], [35]: in which λ is the so-called rate parameter of the exponential distribution and u(·) is the unitary step function. Eq. (2) states that the statistical behavior of the acquired signal intensity only depends on the noise free value y and on the scale parameter λ, that can be assumed equal to one (λ = 1) [35]. Of course, the PDF of Eq. (2) does not take into account any signal transformation, such as the log compression. As previously stated, the NLM approaches divide the noisy image into patches and evaluate their similarities [26]. Given a position m 0 = (x 0 , r 0 ) and a neighborhood system = ( x , r ), the surrounding patch is defined as: with δ x ∈ − x 2 , x 2 , δ r ∈ − r 2 , r 2 , x and r being the patch sizes along the azimuth and range direction. In order to quantify the level of similarity among two patches, a proper metric has to be introduced, and this step is crucial for the whole procedure. Generally, an Euclidean distance is adopted with a windowing kernel introducing some weights to the patch pixels, but many other strategies have been proposed [25], [26].
Once similar patches are detected, pixels all properly merged in order to produce a regularized filtered image: in which q is a normalization factor, w m 0 , m j represents the weighting coefficients related to the similarity between the patches centered at m 0 and m j , andŷ(m 0 ) is the filtered pixel located in m 0 . A sketch of the whole procedure is reported in Fig. 1. Although NLM approaches are characterized by a strong regularization effect, they have some drawbacks related to typical artifacts like ghosts, i.e. structured signal-like patches that mainly appear in flat areas [36].

B. PROPOSED APPROACH: WEIGHTED KOLMOGOROV-SMIRNOV RATIO NLM
The performance of speckle filtering strongly depends on the kind of metric chosen for the similarity evaluation and fusion step. Recently, several statistical metrics have been proposed in literature [29], [31], [32].
In the KS-NLM method proposed in [33], given two patches, the empirical distributions (i.e. eCDFs) are estimated and compared by means of the well-known Kolmogorov-Smirnov (KS) distance [37]. More in detail, this distance d KS ij is defined as: in whichF Z (m i ) (z; y, λ) andF Z (m y ) (z; y, λ) are the estimated eCDF of the patches centered in m i and m j , respectively. If the statistical distributions coincide, the KS distance will be below a certain threshold and the two patches will be  Since the speckle rate parameter λ is equal to 1, the KS distance of Eq. (5) is only related to the difference of the noise free values y between two considered locations, thus this metric is expected to be an effective similarity measurement.
However, there are some drawbacks related to this similarity evaluation. Firstly, the spatial information of the patch, i.e. its texture, is not taken into account when computing the eCDF, and therefore areas that have the same type and quantity of noisy pixels but different spatial textures cannot be distinguished. This case is reported in Figs. 2(a,c,e), where two patches with different textures are characterized by VOLUME 8, 2020 the same empirical distribution and thus considered similar. Another critical step of the KS-NLM is that the estimation of two different eCDFs is required (i.e.F Z (m i ) andF Z (m j ) ).
Within this manuscript, we propose an evolved similarity metric able to overcome the previous limitations. The reported wKSR-NLM analyzes the patch ratio and implements a weighting kernel to enhance the filtering capability. Given two patches z m i ( ) and z m j ( ), we propose to evaluate the ratio among them R ij = z m i ( ) z m j ( ) , which is an elementwise ratio and takes into account the textures of the patches. Fig. 2 provides a proof of concepts of the proposed approach and of its advantages compared to the standard KS metric. In Fig. 2(a), a reference patch has been considered and corrupted by multiplicative exponential noise. Then, a second patch with the same texture has also been corrupted by another noise realization, as shown in Fig. 2(b), respectively. A third noisy patch, with a different texture (the transpose of the previous one) has also been generated (Fig. 2(c)). The estimated eCDFs are compared with the reference one by means of the KS distance, as reported in Figs. 2(d,e). The similarity measurement is then repeated via the proposed wKSR distance (Figs. 2(f,g)). It is easy to observe that the standard KS distance is low in both cases, as it does not take into account any information about the textures.
Conversely, the proposed wKSR metric is sensitive to the differences in the spatial textures of the considered patches, and it is able to discriminate the two cases.
In case the textures are the same, i.e. y m i ( ) = y m j ( ), the ratio between the noisy patches reduces to the ratio between the noise realizations in which the symbol '' '' refers to the element-wise (Hadamard) product. In case of similar textures, R ij is the ratio of two Exponentially distributed random variables, whose distribution is Log-logistic with scale and shape parameters equal to 1 (see Appendix A): where u(r) is the unitary step function. Once the patches ratio R ij is computed, its empirical CDF (eCDF)F R ij (r) can be estimated and compared to the theoretical one of Eq. (7). In case of good matching, the textures can be assumed to be equal, i.e. high similarity between patches centered in m i and m j will be measured. The Log-logistic distribution is completely defined by the scale and the shape parameters α and β, thus the eCDF estimation and the similarity measurement can be achieved by estimating such parameters. More in detail, the KS distance can be computed in closed form according to: The estimation of the scale and the shape parameters α and β can be easily computed in a transformed domain. By applying a logarithmic transformation to the patch ratio, its statistical distribution becomes logistic, which is completely characterized by the mean and scale parameters µ log(R) and s log(R) , the latter related to the variance of log(r) according to: The Log-logistic parameters α and β can easily obtained from µ log(R) and s log(R) via the transformations: To sum up, the proposed algorithm consists of the following steps: 1) select two patches centered in m i and m j ; 2) compute the patch ratio R ij ; 3) apply the logarithmic transformation to the patch ratio; 4) estimate the mean µ log(R) and the variance σ 2 log(R) ; 5) compute the scale parameter s log(R) via Eq. (9); 6) obtain the α and β parameters according to Eq. (10); 7) evaluate the KSR distance d KSR by solving Eq. (8). In case of other proposed distance metrics, e. g. the Euclidean one, the adoption of a weighting kernel has shown to improve the metric effectiveness, as it add the possibility of increasing the weight of pixels in the center of the patch. It is not possible to include such a kernel to the eCDF estimation step, this is why KS-NLM approach does not consider any weighting. Nevertheless, the proposed approach allows such a possibility by including a weighting kernel while computing the mean and the variance of the ratio patch in the logarithmic domain (fourth processing step). Once the weighted moments have been computed, they are are substituted within the eCDF. Fig. 3 reports the weighting kernel adopted within this work, which is in accordance with the one reported in [25]. Thus, the proposed methodology has been named weighted Kolomogorov-Smirnov Ratio NLM (wKSR-NLM).
Once the distance d KSR ij has been computed, the value of the filtered pixelŷ(m i ) is obtained via Eq. (4) by considering the following weights: where T is a tuning parameter that controls the filter strength.

III. NUMERICAL RESULTS
First, an analysis on the optimal parameters configuration (in terms of patch size, search area and tuning parameter) has been conducted. A magnetic resonance imaging (MRI) kidney image, reported in Fig. 4(a), has been degraded by means of a multiplicative exponential noise, obtaining the noisy image of Fig. 4(b). More in details, the speckle samples have been generated as independent and subsequently a correlation on the horizontal axis has been added by means of a 5 pixels long smoothing window. The structural similarity index (SSIM) [38] and the mean square error (MSE) metrics have been computed as quality indices in case of different configurations. The SSIM is a measurement based on luminance, contrast and structure metrics and its optimal value is 1. The MSE measures the energy of the difference between two images, and in case of perfect denoising it is expected to be 0. Both metrics require the knowledge of the noise free image to be computed. Results are reported in Fig. 5. We found that the best performance was achieved with a patch of 19 × 19 pixels, a search area of 33 × 33 pixels and a tuning parameter T = 0.1. This configuration has been adopted for the wKSR-NLM case for all the results reported in the following.
In order to test the effectiveness of the proposed methodology, a comparison with other widely adopted approaches has been proposed both in case of simulated and real datasets. More in detail, we chose to include the following denoising methodologies: • Block-matching and 3-D filtering (BM3D) [39], which belongs to the NLM family; • Optimized Bayesian non-Local means filter (OBNLM) [27], an NLM methodology specifically developed for despeckling US images; • Speckle-reducing anisotropic diffusion (SRAD) [15], a non-linear despeckling approach based on partial differential equations; • Probabilistic patch-based filter (PPB) [40], which belongs to the NLM family and has been specifically tuned for non Gaussian noise. The PPB filter can be applied iteratively or not. For all these algorithms, the configuration parameters suggested by their authors have been adopted. In the case of absence of a default or suggested configuration, the parameters have been manually found by looking for the optimal trade-off between over-regularization and noise removal effectiveness. The parameters are reported in Table 1.
In order to measure the performance of the approaches, the following quality indices have been implemented: • The speckle suppression index (SSI) [41], defined as: where µ and σ are the mean and the standard deviation of the image before (subscript z) and after (subscriptŷ) denoising computed over an homogeneous region. VOLUME 8, 2020  Low SSI values indicate that the denoising methodology is effective in reducing the speckle noise.
• The speckle suppression and mean preservation index (SMPI) [42], defined as: A low SMPI value indicates good filter performances in terms of mean preservation and noise reduction.
• The Edge Preservation Index (EPI), defined as [43]: where y F andŷ F are the filtered version of y andŷ with a 3 × 3 Laplacian operator and after subtracting the mean value and function (·, ·) is defined as: The EPI value is closed to unity when the estimated image resemble the reference one.

A. SIMULATED DATASETS
Two different simulation typologies are considered. The first one consists in considering a reference image (we chose the MRI kidney image adopted previously) and corrupting it, according to the acquisition model of Eq. (1), with multiplicative exponentially-distributed noise (Eq. 2). In this case, the Cartesian geometry instead of the US one has been considered. The second one consists in generating a realistic dataset by considering all the multi-path interactions of the acoustic waves within the image object. For this approach, the Field 2 simulation software (available at http://field-ii.dk) has been adopted, and two datasets were selected, namely the cyst and the kidney ones, reported in Figures 6(a) and 7(a). More in details, the output of the simulation software has been under-sampled in the range direction up to a resolution of 580 × 128 pixels. Further details about the Field 2 simulation software can be found in [44]. We chose to adopt these two approaches in order to have a less accurate case study but with the reference (i.e. noise free) image available, and a very accurate case study (provided by Field 2), but without any ground truth.
The results related to the first approach (Cartesian geometry) are reported in Fig. 4. More in detail, the noisy image of Fig. 4(b) has been processed by all the considered approaches, achieving the results reported in Figs. 4(c-h). All the regularized images are characterized by a different trade-off between noise reduction and details preservation. The proposed wKSR-NLM methodology is able to combine such two aspects in an effective way. Quantitative results are reported in Table 2. We recall that the first column (the SSI values) reports the effectiveness in reducing noise, the second one (the SMPI values) takes into account also the preservation of the mean value of the image, while the last one (the EPI values) evaluates the effectiveness in preserving edges. Among all, the OBNLM approach is characterized by a strong speckle suppression. Nevertheless, from the SMPI values it can be noticed that the proposed approach is able to better conjugate noise reduction and image preservation. Regarding the edges preservation, the EPI scores of OBNLM and wKSR-NLM are very close, indicating similar performances. As expected from the qualitative results of Fig. 4, the proposed wKSR-NLM filter is very The results related to the cyst dataset provided by Field 2 software are reported in Fig. 6, and the computed SSI and SMPI values are reported in Table 3 for three regions of the image, namely the hypo-intense cysts, the background and the hyper-intense cysts highlighted by the red, the blue and the green squares in Fig. 6(a). Similarly to the previous case study, the proposed wKSR-NLM filter provides very interesting results mainly by looking at the SMPI index, indicating that it is able to combine speckle suppression and signals level preservation. Fig. 6(f) is characterized by a very clean background together with sharp edges for all the cysts across the image. This dataset has also been considered for the evaluation of the computation complexity related to the proposed approach. The same code has been run twice, the first time with an euclidean distance metric, the second one with the proposed algorithm and the computational times have been measured. For this analysis, an Intel R Core TM i7 workstation equipped with 32 GB of RAM has been considered. On the software side, the operative system was the Linux Debian (kernel 4.19) and the computing environment was the Mathworks Matlab R R2018b. An computation time increase of VOLUME 8, 2020  about 18% was found, as the code took 523 seconds to run for the Euclidean metric, and 641 seconds for the proposed one. In order to have a fair comparison, no code optimization was performed, thus it is expected that these time will be greatly reduced by engineering the code.
The results related to the second dataset provided by Field 2 software, i.e. the kidney dataset in the US geometry, are reported in Fig. 7. This is a more challenging image as it is rich in small details and edges. Also in this case the wKSR-NLM filter is able to provide a clean result with a satisfactory level of details.

B. REAL DATASETS
For the real data analysis, two datasets acquired by the Ultrasonic Imaging Laboratory at the University of Illinois at Urbana-Champaign have been considered. The two acquisitions are related to a breast benign (fibroadenoma, Fig. 8(a)) and malignant (invasive ductal carcinoma, Fig. 9(a)) tumors obtained from biopsy-verified studies.
The radio-frequency frames have been acquired using a linear transducer array from the Antares System. The B-mode ultrasound images adopted in this comparison have been formed via the URI offline processing tools (URI-OPT). Each one is composed of 518×180 pixels (lateral size × axial size). The first dataset (benign tumor) is characterized by some structures in the top part of the image, which are useful for evaluating the details preservation capabilities, and a smooth area in the bottom part, that allows to analyze the performance in terms of noise reduction. From the results reported in Fig. 8, the proposed wKSR-NLM methodology produces a high quality image. By looking at the other filters, the PPB and the OBNLM approaches manage to preserve details well, at the cost of lower noise reduction, while BM3D and SRAD balance the two aspects similarly to wKSR-NLM.
The malignant tumor dataset is much harder to be analyzed, as no sharp structures are clearly visible and a severe speckle noise is present. Nevertheless, similar findings than the previous dataset can be found: PPB and OBNLM are characterized by a slightly higher level of residual noise, while BM3D, SRAD and wKSR-NLM regularize the image in a more pronounced way.

IV. CONCLUSION
Within this manuscript, a despeckling filter for ultrasound images, namely the wKSR-NLM, has been presented. The approach belongs to the Non Local family, i.e. identifies similar patches across the image and fuse them in order to produce the regularized image. The wKSR-NLM is able to statistically measure the similarity between patches taking into account the texture. Compared to the standard KS-NLM, the proposed approach is able to take into account the textures of the patches when measuring their similarity. Moreover, the wKSR-NLM adopts a kernel for weighting differently the inner and the outer parts of the patch. The performance of the filter has been tested in both synthetic and real data sets. Moreover, a comparison with other state of the art despeckling filters showed the ability of the proposed approach of combining effective noise reduction with edges and details preservation.

APPENDIX
Let us consider two noise samples n 1 and n 2 . They can be modeled as two independent Exponentially distributed random variables (N 1 and N 2 ) with parameter λ = 1, i.e.: f N 1 (n 1 ) = e −n 1 u(n 1 ), f N 2 (n 2 ) = e −n 2 u(n 2 ) . (A.1) Note that, due to the independence between the variables, the joint pdf is the product of the marginal ones, i.e. f N 1 N 2 (n 1 , n 2 ) = f N 1 (n 1 )f N 2 (n 2 ).
Let us define the two random variables R and S as the ratio and the sum of N 1 and N 2 , respectively: We are interested in deriving the joint statistical distribution of R and S, namely f RS (r, s), starting from f N 1 N 2 (n 1 , n 2 ). Following the rules for transforming random variables, it can be computed as follows: by considering that the integrand is the pdf of the Erlang distribution with shape and rate parameters equal to 1, the integral reduces to 1, so: from which the CDF of R can be derived, obtaining: which is a Burr III, also called Dagum, distribution with unitary scale and shape parameters.