Digital Path Approach Despeckle Filter for Ultrasound Imaging and Video

We propose a novel filtering technique capable of reducing the multiplicative noise in ultrasound images that is an extension of the denoising algorithms based on the concept of digital paths. In this approach, the filter weights are calculated taking into account the similarity between pixel intensities that belongs to the local neighborhood of the processed pixel, which is called a path. The output of the filter is estimated as the weighted average of pixels connected by the paths. The way of creating paths is pivotal and determines the effectiveness and computational complexity of the proposed filtering design. Such procedure can be effective for different types of noise but fail in the presence of multiplicative noise. To increase the filtering efficiency for this type of disturbances, we introduce some improvements of the basic concept and new classes of similarity functions and finally extend our techniques to a spatiotemporal domain. The experimental results prove that the proposed algorithm provides the comparable results with the state-of-the-art techniques for multiplicative noise removal in ultrasound images and it can be applied for real-time image enhancement of video streams.


Introduction
Medical ultrasound is an imaging technique widely used in the diagnosis and assessment of internal body structures, and it plays a key role in treating various diseases. Compared to the other techniques of medical imaging, it is safe, noninvasive, and well tolerated by the patient, and ultrasound images are captured in real time at reasonable price. An accurate analysis of ultrasound images and thus an appropriate diagnosis are difficult due to the fact that the images are contaminated with characteristic granular structures called speckle noise, which deteriorates contrast and hinders the identification of important image details [1]. Although the ultrasound images are subjected to an initial improvement during the acquisition process, their quality is still far from optimal. Therefore, the main aim of image denoising is to remove the noise, while preserving the important details.
Recently, a plethora of methods capable of diminishing speckle noise have been proposed [2]. According to the works presented in [3][4][5], one of the most promising results for ultrasound images was obtained with algorithms based on anisotropic diffusion techniques [6] and the idea of nonlocal means [7][8][9][10][11]. However, majority of them were designed for static images, and not much attempt has been made to video by considering temporal coherence.
Video sequence processing algorithms can take an advantage of high correlation between adjacent frames, exploring spatial and temporal neighborhood. These properties are incorporated using different variants of averaging the pixel intensities in successive video frames. Simple temporal filters, such as temporal Gaussian, efficiently remove noise for videos with slowly moving objects, but averaging operation of subsequent frames can cause "ghosting" artifacts in the output results. Some works indicate that motion compensation allows us to reduce the blurring effect [12,13]. Additionally, "ghosting" effects can be reduced by application of the statistical models that reflects the joint distributions of wavelet coefficients over space and time [14,15].
Ghosting artifacts may be also omitted using temporal version of a bilateral filter that was successfully applied in the adaptive spatiotemporal accumulation (ASTA) filter [16].
In our previous works on video filtering, a spatiotemporal denoising scheme utilized local switching between spatial digital path approaches and temporal Gaussian filtering [17] and a fuzzy spatiotemporal filter that was later described [18]. Another interesting study connected with video denoising introduces a 3D filtering framework that is based on fuzzy set logic to combine the gradient values in different directions between previous and current temporal frames [19].
The main aim of this research is to develop a filter that will efficiently cope with multiplicative noise in ultrasound images and videos. The proposed algorithm is based on the idea of spatial digital paths presented in [20,21]. The original 2D algorithm presented in [20] perfectly removes Gaussian noise and after some modifications impulsive noise but fails in the presence of multiplicative interferences. The new denoising scheme explores spatial or spatiotemporal pixel neighborhood in order to calculate the filter weights for pixels belonging to the processing window. Digital paths in a spatiotemporal domain could be understood as trajectories or object displacements in subsequent frames of a video stream.
The proposed technique utilizes a specific kind of digital paths, the so called escaping paths, and extends this concept from the spatial domain (2D) to the spatiotemporal domain (3D) that allows us to efficiently reduce the speckle noise. Furthermore, to increase the filter denoising ability, the new method of path creation utilizes the extended neighborhood introduced by von Neumann for cellular automata [22]. A detailed description of the proposed algorithm and its extensions are presented in Section 2. In Section 3, we present the results of experiments and the comparison with competitive filters. Finally, the conclusions are drawn in Section 4.

Method
2.1. Speckle Noise Model. Speckle noise is a signal-dependent and non-Gaussian multiplicative image distortion. Such noise is generally more difficult to remove than additive noise [23]. This type of distortion appears in sonar, laser, and synthetic aperture radar (SAR), and it depends on the structure of the material being imaged and various acquisition parameters [24]. Speckle noise is a random process and can be modeled using gamma, Rayleigh, and Fisher-Tippett distribution [25][26][27]. In our work, the simplified speckle noise model has been employed, since it has been applied successfully in many studies [11,24,28,29].
where u x , v x , and n x denote the observed signal, original unknown signal, and zero-mean Gaussian noise, respectively.

General Filter Framework.
Since smoothing is commonly used to decrease the level of random noise, an averaging operation is required in order to replace the noisy pixel, v x , with a suitable pixel representative for the local spatiotemporal neighborhood of point x = x, y, t . So the general form of the fuzzy adaptive filters in this work is defined as the weighted average of inputs inside the spatiotemporal window W that are in neighborhood relation N with a center pixel x [30,31].
where u x i andv x denote filter inputs and outputs, respectively, and μ x, x i is a similarity function computed along digital path starting at window central point x, associated with its neighbor x i , and bounded by the spatiotemporal processing window W. To describe the model of digital paths, a few notions should be introduced: digital lattice ℋ = X, N defined by X, which is the set of all points of the image sequence, and a neighborhood relation N between the lattice points [32].
A digital path P = p i n i=0 defined on the image lattice is a sequence of adjacent points p i , p i+1 ∈ N , and n is a number of path segments. Let L P denote the length of the digital path P p i n i=0 that is calculated as ∑ n−1 i=1 ρ p i , p i+1 , where ρ denotes the spatiotemporal Euclidean distance between two adjacent points of the path.
The connection cost over the single digital path P = p i n i=0 can be determined as a measure of dissimilarity between image pixels p 0 , p 1 , …, p n that forms a specific path linking p 0 and p n [33,34]. In the new approach, the connection cost will be calculated as a combination of the topological length of the path and the cumulative differences of gray levels. Thus, the connection cost for the entire path Λ P can be defined as follows: 3. Similarity Functions. Weights for the filter described by (2) can be defined in many ways; in our case, we use two different approaches based on connection costs calculated along digital paths. We created two kinds of membership functions which leads to the filters with slightly different properties.
2.3.1. DPA 1st . In this approach, similarity functions are defined for all neighbors of the central point x that remain in the neighborhood relation. To define the similarity function μ x, x i between the filtered point x and its ith neighbor, we create digital paths starting at center point x = p 0 , intersecting x i = p 1 and finally terminating at point p n , which may be reached in n steps from p 0 . The similarity function defined for points p 0 and p 1 uses paths exploring the further neighborhood of the central point passing points p 2 , …, p n , so that the filtering result will be better suited to local image structures. An illustration of this idea is presented in Figure 1. This approach will be further denoted as DPA 1st . In this case, the similarity function takes the form as follows: where ω denotes the number of all possible paths P m = p 0 , p 1 , p m 2 , …, p m n with n steps totally included in the processing window W, originating at x = p 0 and crossing x i = p 1 ; m is the index of a specific path; Λ ⋅ is a dissimilarity value along a specific path; and g ⋅ is a smooth function of Λ. In this work, the exponential function is used as g ⋅ [20], so the similarity function takes the form as follows: where β is the filter design parameter.

DPA last .
Another approach is to determine the similarity function between pixels x and x i by all possible paths connecting them ( Figure 2). This approach will be denoted as DPA last , and the similarity function between image points x and x i can be defined as follows: where ω denotes the number of all possible paths P m = p 0 , p m 1 , p m 2 , …, p n with n steps connecting x and x i and totally included in the processing window W, Λ ⋅ is a dissimilarity value along a specific path, and g ⋅ is a smooth function of Λ. Finally, the DPA last similarity function takes the form as follows: . Filter Output and β Normalization. The proposed method can use different path lengths and bit depths, and therefore, to ensure the comparability of the results, it was necessary to rescale the parameter β. In this case, the β parameter will be divided by the maximum possible cost of the single path; thus, in (5), the normalized valueβ will be used.
Maximum path cost can be determined using the formula as follows: where bpp denotes the number of bits per pixel and n the number of path steps. The next stage is a normalization of the similarity function, which can be defined as follows: Figure 1: Illustration of paths created on the 2D image lattice with the DPA 1st approach, used to determine the similarity function between two adjacent points. Let x = p 0 denote the pixel under consideration, with u x i representing the noisy pixel x i ; the filter outputv x is given as follows:  (Figure 3).
The new spatial filters with N 8 neighborhood and the spatiotemporal ones with N 18 and N 26 are very efficient, but they fail for heavily degraded images. Therefore, in the proposed denoising design, we introduced the extended von Neumann neighborhood [22] originally defined for cellular automata. Various neighborhood systems for static 2D images are drawn in Figure 4, while 3D neighborhoods are shown in Figure 5.
Additionally, the efficiency of the proposed denoising framework is strongly connected with the type of digital paths. Different models of paths allow us to suppress a certain type of noise [20].
Previous research has demonstrated that the best results are obtained in the presence of impulsive Gaussian as well as of mixed-type noise for the self-avoiding path (SAP) model, but our experiments suggest that in the case of ultrasound images, the greatest results are achieved for the so called escaping path model (EPM); thus, the new filter will be denoted as the escaping path filter (EPF). In the proposed denoising scheme, the topological distance from the initial point in the following steps must be increased. Exemplary spatial escaping paths created with various neighborhood systems are illustrated in Figure 6, while Figure 7 depicts the spatiotemporal case. Later in this paper, the proposed filters will be marked as EPF2D for the case of two-dimensional filtering (2D) and EPF3D for the spatiotemporal case (3D).
In situations when the images are highly contaminated, we can increase the efficiency of the filter in one of two ways: (1) extend the length of the used paths or (2) apply it in an  iterative manner. The second option is much faster and more accurate and allows us to control the filter strength adjusting the β parameter in subsequent iterations. In this way, β can be updated as follows:

Simulation Results
3.1. Static Images. The commonly used benchmark images goldhill, boats, and artificially generated phantom were chosen to compare efficiency of different filters. Besides, we used a synthetic test image that was a phantom for a 3-month-old fetus denoted later as fetus. This test image was obtained using the Field II applications, which simulate the ultrasound field that is based on linear acoustics using the Tupholme-Stepanishen method for calculating the spatial impulse response [35][36][37]. The proposed ultrasound data seems to be more reliable, due to the fact that it contains artifacts typical for an ultrasound acquisition process and thus it can be used for image quality evaluation. A reference noise-free image was achieved by averaging of 500 simulated images and was depicted in Figure 8. The described filtering design has been compared with the following state-of-the-art methods capable of suppressing a speckle noise: (i) Wiener filter [2] (ii) Speckle reducing anisotropic diffusion (SRAD) [6] (iii) Nonlocal means (NLM) [7,8] (iv) Optimized Bayesian nonlocal means (OBNLM) [11] (v) Probabilistic nonlocal means (PNLM) [9] (vi) Probabilistic patch-based weights (PPBW) [10] (vii) Digital path approach (DPA last ) [20] The source codes of the algorithms used in our comparison were provided by the authors of the respective papers. In order to determine the optimal values of parameters for the reference filters, we tested a wide range of parameters according to the authors' recommendations to obtain the highest possible PSNR value.
The recommended values of the parameters α and β for the proposed technique were adjusted, so that the best PSNR value was achieved. The test images were deteriorated by the multiplicative noise described by (1) with mean μ = 0 and σ 2 = 0 2, 0 4, and 0 6 except for the fetus image, which was contaminated using a more realistic artifacts resulting from the physical model of an ultrasound image acquisition process. The restoration efficiency has been assessed using the peak signal-to-noise ratio (PSNR) and more sophisticated mean structural similarity index (MSSIM) calculated with Gaussian kernel and default parameters [38]. Numerical results obtained for the static images are summarized in Table 1. The conducted simulations revealed that the proposed EPF approach outperforms other techniques for highly deteriorated real images in terms of PSNR metric, while algorithms that are based on nonlocal means provide slightly better results for lower noise contamination level. Additionally, the proposed filter gives better results for synthetic images. A visual comparison of the results achieved for the phantom image is drawn in Figure 9. From this figure, it can be seen that most filters removed multiplicative noise and produces more visually pleasing results, but the outcomes are slightly blurred.
The most significant and valuable results seems to be obtained for simulated fetus image, because it is much similar to the realistic ultrasound images. The visual comparison of the performance of the analyzed filters for the fetus image is presented in Figure 10, while exemplary results for the real ultrasound images of fingers are depicted in Figure 11. The assessment of the achieved results can be also evaluated by segmentation accuracy applying image denoising as the preprocessing step, but currently, it is out of scope in this work. The proposed denoising scheme requires also a much  smaller neighborhood than that used in family of nonlocal means methods; therefore, our approach is much faster and less aggressively blurs the image. It should also be emphasized that for all methods based on the concept of NLM, even a small modification of optimal values of parameters gives a significant decrease in performance, while the EPF framework gives acceptable results for a wide spectrum of parameters.  approach. Additionally more complex techniques, based on nonlocal means [8] and BM3D [39], were added to comparison; however, the computational complexity of those filters limits their use for offline processing. The performance of the following filters was evaluated:

Video
(i) Wiener 2D-a spatially adaptive Wiener filter The video denoising algorithms were tested using publicly available video sequences: foreman, salesman, and tennis, contaminated by the multiplicative noise described in (1) with mean μ = 0 and σ 2 = 0 2, 0 4, and 0 6. To obtain a more realistic comparison, we have prepared a test sequence based on the fetus image. Base fetus sequence consists of 150 frames subjected to different transformations that simulate the possible displacements during the ultrasound acquisition process. Then, all the frames have been subjected to a simulation using the Field II application. Virtually noise-free reference video was obtained by averaging of 200 simulation results (reference videos could be downloaded from http:// dip.aei.polsl.pl/usg/fetus_video.7z). Table 2 summarizes results obtained for the set of test video sequences. Based on synthetic tests only, it is difficult to choose the best filter. In most cases, especially for small levels of noise, the best PSNR results are obtained for the BM3D filter; however, our spatiotemporal solution gives only slightly worse results at much higher processing speed. For heavily corrupted sequences, we can clearly see the advantage of the proposed filtrating technique. It should be noted that for the more realistic ultrasound noise model, obtained using the Field II application, the advantages of our solution is clear (the best results were obtained for the NLM3D filter, but the computational complexity disqualifies it entirely, even for offline processing). In the case of the tennis sequence with minimal disruption, most filtering techniques deteriorate quality ratios. This is due to the specific background that resembles an impulsive noise pattern. Figure 12 presents exemplary filtering results and SSIM maps for the mentioned sequence. It is thus clear that the worst results were obtained for filtering based on the median, which effectively removes elements of background texture. In addition, median 3D filtering introduces some temporal artifacts, such as blurred or erased moving objects. It should be noted that the ranking obtained with the SSIM index is slightly different in favor of the escaping path filters.

Conclusions
In this paper, a new class of fast spatial and spatiotemporal filters was presented. The proposed filtering techniques were designed for multiplicative noise suppression, specifically for ultrasound image and video filtering. The novel approach is based on the special type of digital paths, the so called escaping paths, created on an image lattice-spatial in a 2D case or spatiotemporal for video processing. Additionally, the new extended neighborhood model was introduced, based on von Neumann concept derived from cellular automata theory. The presented methods give comparable or better results to the other methods, both for static image and video sequences. Another beneficial feature of the proposed denoising scheme is its lower computational complexity than that of other state-of-the-art techniques, which allows us to apply it in real-time image processing tasks. The proposed filter is also more stable for a wide range of input parameters and gives satisfactory results in terms of different quality metrics and visual inspection.