Gabor-based anisotropic diffusion with lattice Boltzmann method for medical ultrasound despeckling

1.   Introduction

Speckle noise is a granular pattern that is commonly seen in coherent imaging systems, such as ultrasound (US), optical coherence tomography, and synthetic aperture radar imaging. It is a form of locally correlated multiplicative noise caused by constructive and destructive coherent interferences of sound or light echoes from randomly distributed scatters, whose resolution is beyond the capabilities of the imaging systems ^[1]. In medical ultrasonography, effective and efficient despeckling methods are required to guarantee the accuracy of diagnosis. The main purpose of despeckling is to reduce the speckle noise while preserving the details in images for further processing such as image segmentation, feature extraction and quantitative measurement. That is to say, ultrasound image despeckling is an essential preprocessing step for subsequent analysis.

Many different approaches to reducing ultrasound speckle noise have been proposed. Classic ultrasound despeckling methods are based on local statistics ^[2,3,4], such as the adaptive weighted median filter, Frost filter, and Lee filter, which can effectively reduce noise but may blur some details of tissues. Other methods are based on the anisotropic diffusion (AD) ^[5], and they can suppress noise and smooth an image ^[6]. However, whether they can keep or even enhance image details depends on whether the edge detector it uses can efficiently distinguish edges from noise ^[7].

The AD method, also known as the P-M model, uses partial differential equation (PDE) to filter an image step by step by iterated diffusion ^[8]. Unlike the traditional spatial filtering, that does not protect image edges and details, the AD can preserve or even enhance edges while eliminating noise. However, when a medical ultrasound image is seriously corrupted by speckle noise, the AD method cannot effectively distinguish image edges from noise, resulting in the unsatisfactory performance of noise reduction. The edge detectors in PDEs are critical to filtering. An appropriate edge detector for medical ultrasound images should be able to differentiate edges and noise effectively ^[9]. However, speckle noise is a multiplicative noise, which means that the gradient changes caused by noise may be equal to or even larger than that caused by edges. In the AD model, the gradient operator acts as the edge detector, but it is very sensitive to speckle noise; it is also the main reason why the AD cannot obtain acceptable results for ultrasonic image denoising ^[9,10,11].

We have recently proposed a robust edge detector called the Gabor edge detector based on the Gabor transform and embedded it into the AD's PDE to guide the diffusion. This improved method is called the Gabor-based anisotropic diffusion (GAD) ^[12]. Because the Gabor edge detector can distinguish edges from noise effectively, the GAD can not only suppress speckle noise, but also preserve the edges, structures and details in an ultrasound image ^[13,14]. However, the GAD requires the finite difference method to solve the PDE iteratively, and thus it is a time-consuming method, which limits its application in clinical practice.

Indeed, various existing AD filters including the GAD have a common drawback, the poor efficiency, because they use the finite difference discretization depending on an explicit scheme, which requires a very short time step to make the algorithm stable. Hence, many iterations are required and cause the entire procedure be rather time-consuming. Alternatively, the lattice Boltzmann method (LBM) is a numerical scheme for solving the PDE by constructing a simplified discrete microscopic dynamics to simulate the macroscopic model described by the PDE ^[1,15,16,17]. One large superiority of the LBM is that it is suitable for parallel computation ^[18,19]. Recently, the LBM has been used for solving image diffusion problems, including image denoising with the AD ^[20,21,22]. Thus, we speculate that the LBM may also have the potential for the GAD denoising.

In this paper, we propose using the LBM to fast solve the PDE of the GAD and thus derive a despeckling approach named the Gabor-based anisotropic diffusion coupled with the lattice Boltzmann method (GAD-LBM). In the following sections, we will show that the GAD-LBM not only has the advantages of the GAD including effective speckle noise reduction and detail preservation, but also has the advantages of fast computation speed, high parallelization, and good stability attributed to the LBM.

2.   Methods

2.1.   Gabor-based anisotropic diffusion

The two-dimensional Gabor kernels (i.e., functions) can be expressed as ^[23,24]:

where v and θ are scales and directions of Gabor kernels, z = (x, y) is the spatial coordinates; σ is a constant related to the frequency bandwidth of the wavelet and |∙| denotes the norm operator. The wave vector is defined as:

where k_max is the maximum frequency, and f represents the kernel interval factor in the frequency domain. The spatial frequencydomain information at different scales and directions can be extracted by changing v and θ.

In the paper, the Gabor functions are defined at one scale and multiple directions, namely v = 0, and θ = θ_d = d ∙ 2π/D (d = 0, 1, ..., D – 1; D is the number of directions). Two-dimensional Gabor functions can be decomposed into real and imaginary parts. The imaginary parts are often used to compose the edge detector of images ^[12]. The imaginary part of the Gabor function g_b can be expressed as follows:

where k_d = k_max(cosθ_d, sinθ_d). In this paper, we set σ = 2π, k_max = π/2 and D = 16 ^[24].

The Gabor transform of an image I(x, y) can be expressed by convolution between the image and the Gabor kernel:

G_d(x, y) is the dth convoluted image resulted from the Gabor kernel at the direction θ_d. The details in the input image I(x, y) at the specified direction θ_d can be captured through the Gabor transform. Then the edge-enhanced image G_sd(x, y) by using the Gabor edge detector can be determined according to details from multiple directions:

The GAD model can be built as the following PDE, which iteratively eliminates small variations in an image due to speckle noise and to preserve large variations due to edges:

Here g(G_sd) is the diffusion coefficient, which controls the diffusion process that inhibits smoothing across edges. g(G_sd) is large where G_sd is small (namely where the edge is weak), and it can be expressed as:

where k is an edge magnitude threshold; generally, a large k-value produces a smoother result in a homogeneous region than a smaller one ^[25].

2.2.   Lattice Boltzmann method

The LBM is built with a lattice together with the lattice vectors c_α (α = 0, 1, …, q–1) along q directions. On each lattice node, there are a set of particles, consisting of a rest particle (α = 0) and q–1 moving particles ^[26,27]. Each particle has a particle density function I_α (α = 0, 1, …, q–1) ata vector direction c_α ^[28]. Our research is based on a two-dimensional LBM model D2Q9 (2 dimensions and 9 directions) ^[29]. The lattice vectors of the D2Q9 are defined as:

where α = 1, 2, 3, 4 shows the vertical and horizontal directions, α = 5, 6, 7, 8 indicates the directions of 45,135,225, and 315 degrees, respectively.

I_α(x, t) is the particle densityfunction (i.e. the image intensity) at the position x = [x, y] and time t with a velocity of c_α. Its evolution equation is:

where ω represents the relaxation factor, and $I_\alpha ^{eq}$(x, t) is the equilibrium distribution which is the predicted value of I_α(x, t).

The diffusion is conducted under the law of conservation:

The equilibrium distribution function of D2Q9 is:

After setting the relaxation factor and equilibrium distribution function in Eq (9), we can solve the PDE with given initial values and certain boundary conditions.

2.3.   The proposed GAD-LBM algorithm

The GAD is a competent despeckling model for medical ultrasonography which can effectively distinguish edges from speckle noise ^[12]. Moreover, the LBM is an algorithm with advantages of high stability, easy parallel implementation and fast calculation speed. The GAD-LBM model proposed in this paper combines the GAD and the LBM.

First, the relaxation factor ω in Eq (9) is not a constant across the image but a function of the image intensities denoted by ω(I), thus Eq (9) is rewritten as follows:

We denote x = (x₁, x₂) and c_α = (c_α1, c_α2). We then use a Taylor expansion on the left side of Eq (12) ^[26] and obtain

where ε is a small parameter, and i, j∈ {1, 2}. To derive the macroscopic equation, we use the Chapman-Enskog expansion ^[27] of I_α around $I_\alpha ^{eq}$ in powers of ε:

Then we have ^[26]

We substitute (15) into (13) and equate the ε term of both sides of Eq (13) and get

Combining Eq (15) and Eq (16), we obtain

By taking the sum over α on both sides of Eq (13), we have

Here we have used Eqs (10) and (11). By substituting (17) into (18) and equating the ε²-term of both sides in Eq (18), we get

where the diffusion tensor

Equation (19) is equivalent to

with the diffusion coefficient

With regard to the GAD-LBM method, we use the Gabor transform to derive the edge-enhanced image G_sd via Eq (5) and the diffusion coefficient g(G_sd) via Eq (7). Then we obtain the diffusion equation of the D2Q9 for GAD-LBM as follows:

where the relaxation factor is given by:

2.4.   Diffusion control

The diffusion time t is an important parameter in the process of diffusion-based image filtering. If a fixed diffusion time is used to despeckle ultrasound images with different noise levels, the filtering method will lose flexibility, and the filtering efficacy and efficiency cannot be guaranteed. In order to improve the adaptive ability of the GAD-LBM, we automatically set the diffusion time by searching for a fully formed speckle region (FFSR) in the image, which is a homogeneous region containing no edges but only speckle noise ^[12].

Firstly, we divide the edge-detected image with the Gabor edge detector, G_sd(x, y), into several blocks in m_s× n_s size, and then a global search is carried out, from which a block with the smallest mean is selected as the FFSR. During the diffusion iterations in FFSR, the mean absolute error (MAE) of two successive filtered images (e.g., at the current time t and the previous time t–Δt) is computed:

where $I_{{\rm{FFSR}}}^t$ and I_t–ΔtFFSR(i, j) are the gray levels of pixels (i, j) in the FFSR at two successive iterations. When the MAE^t value is less than the threshold of 0.002, the diffusion will be automatically terminated.

2.5.   The procedure of GAD-LBM

The complete procedure of the GAD-LBM model can be summarized as eight steps.

Step 1: We take an initial noisy image I (x, y, t = 0) with a size of M × N to initialize the distribution functions at all directions:

where α = 0, 1, …, 8 denoting 9 directions. We then search for the FFSR region in I (x, y, t = 0).

Step 2: We initialize the relaxation factor according to Eqs (5), (7) and (24).

Step 3: Collision step:

where n is the iteration time, and the x- and y-axis coordinates are i = 1, …, M and j = 1, …, N.

Step 4: Translation step:

where $I^{n+1^{{\prime}}}$ is the difference between Iⁿ⁺¹ and Iⁿ.

Step 5: We deal with the boundary condition:

Step 6: We update the gray values of pixels:

Step 7: We update the edge-enhanced image G_sd in the (n+1)th iteration and the relaxation factor of pixels:

Step 8: If the MAE^t is still greater than the threshold, we change the diffusion time to t + Δt, and then return to Step 3 for the next iteration. Otherwise, we stop the diffusion, and the final output I is the despeckled image.

3.   Experiments and results

In this section, we present the experimental results on synthetic and clinical ultrasound images by using the proposed algorithm GAD-LBM and nine other algorithms for comparison including the Lee filter ^[30], AD ^[5], AD-LBM, speckle reducing anisotropic diffusion (SRAD) ^[31], Laplacian pyramid-based nonlinear diffusion (LPND) ^[32], complex diffusion (CD) ^[33], non-local means (NLM) ^[34], James-Stein typenon-local means (JSNLM) ^[35] and GAD (solved with finite difference) ^[12].

The parameters of all the existing methods were tuned in accordance with the values suggested in the publications. The computer used in our experiments was a server with the Intel Xeon CPU (E5-2620v2 2.1GHz) and 64GB RAM.

Clinical ultrasound images were acquired from human participants. All procedures performed in studies involving human participants were in accordance with the ethical standards of the institutional and/or national research committee and with the 1964 Helsinki declaration and its later amendments or comparable ethical standards. Informed consent was obtained from all individual participants included in the study.

3.1.   Results of synthetic images

We used a multiplicative model to add multiplicative speckle noise n_s(x, y) to a synthetic noise-free image I(x, y) with a size of 256 × 256. The noise was in uniform distribution with a mean of zero and a variance of v_n. Hence we generated an ultrasound image Y(x, y) contaminated by speckle:

Firstly, we changed the v_n-values from 0.03 to 0.21 in an interval of 0.03 in order to obtain synthetic images with various amounts of speckle noise. Simulated images with speckle variances of 0.03 and 0.12 were shown in Figure 1 and Figure 2 (entitled "Noisy"), respectively.

In Figure 1 and Figure 2, the GAD-LBM algorithm had the best filtering result through visual observation of the image. When the noise variance was 0.03, four methods including the AD-LBM, SRAD, GAD, and GAD-LBM yielded satisfactory results (Figure 1). However, when the variance of noise increased to 0.12, only the GAD and GAD-LBM seemed to obtain robust despeckling results, while the performance of the AD-LBM and SRAD degenerated (Figure 2). Compared with the GAD, the GAD-LBM appeared better in speckle noise suppression and edge preservation (Figure 1 and Figure 2).

In addition, we used three quantitative indices, i.e., the peak signal to noise ratio (PSNR), mean structural similarity (MSSIM), and Pratt's figure of merit (FOM), to evaluate the algorithms' capability of noise reduction, structure preservation and edge preservation, respectively ^[36,37,38].

The quantitative indices for synthetic images are listed in Table 1, where the values were calculated when the speckle variance v_n was 0.12. From the table, we can know that the GAD-LBM largely enhanced the image quality measured by the PSNR, MSSIM and FOM, exceeding other nine methods in terms of all three indices.

We changed the variance of noise from 0.03 to 0.21 in an interval of 0.03 when synthesizing the images with the multiplicative model and compared the performance of all ten filtering methods. The despeckling results are illustrated in Figure 3, where we can see that the GAD-LBM and GAD outperformed other methods on almost every noise level in terms of every quantitative index.

Our experimental results show that when the variance of the speckle noise was relatively low (0.03 and 0.06; Figure 3), the GAD, GAD-LBM, AD and AD-LBM all achieved good PSNR values, indicating that the GAD-LBM and the traditional algorithms suppressed the noise well when the noise level was low. When the noise variances increased from 0.03 to 0.21, the PSNR values of the GAD-LBM and GAD were always above 28 dB, demonstrating they maintained excellent noise reduction performance at high noise levels. The FOM values of most methods decreased sharply with the increase of noise levels. Nevertheless, the FOM of the GAD-LBM still reached 0.67 when the noise increased to the highest variance of 0.21. The results of MSSIM show that the GAD-LBM was stable and always higher than 0.92 at a noise variance from 0.03 to 0.21, demonstrating the superiority of our method in structural preservation.

From Figure 4, we see that the computation time of the GAD-LBM was stable from 2 to 4s for denoising simulated noisy images with a variety of noise variances, which was about 1/50 to 1/20 of the GAD, indicating the advantage of the GAD-LBM in efficiency.

3.2.   Results of clinical images

We also compared the GAD-LBM with nine methods on despeckling two clinical B-scan ultrasound images of breast tumor acquired from Huashan Hospital, Fudan University, Shanghai, China. The imaging device was the HI VISION Preirus system (Hitachi Medical, Tokyo, Japan) equipped with a linear array transducer. In Figure 5 and Figure 6, the yellow line in the original image is a sampling line, and we show both the two-dimensional image and the one-dimensional signal sampled from the line.

It can be seen from Figure 5 and Figure 6 that the GAD-LBM performed best and removed noise and preserved edges appropriately. In the denoised images by ten filtering methods, the speckle noise was suppressed to different extents. The Lee and the CD methods retained much noise while the other methods all smoothed the tumor. Among the ten methods, the GAD-LBM seemed to seek the best compromise between the noise reduction and the detail preservation. A sampling line is superimposed in the initial noisy image. From the plots of the intensities on the line, we can see that the GAD-LBM appropriately eliminated noise and smoothed intensities on the line while the other methods either over-smoothed or under-smoothed the intensities.

Table 2 shows the computation time of ten methods. It can be seen that it took 2 to 4s for the GAD-LBM to denoise a clinical ultrasound image, which was comparable to the Lee and LPND methods and much quicker than the GAD method. The computation time of the GAD-LBM was about 1/90 to 1/30 to that of the GAD, indicating the advantage of the GAD-LBM in efficiency.

4.   Discussion

Speckle noise corrupts ultrasound images and limits the signal-to-noise ratio. Moreover, the speckle noise is multiplicative, increasing the difficulty for noise suppression. In order to improve the feasibility and reliability of subsequent image analysis and interpretation, it is necessary to preserve or enhance the edge and texture of tissues while reducing speckle noise. However, the commonly used denoising methods cannot effectively denoise the speckle noise. The recently proposed GAD algorithm has a good performance for despeckling, but it is a time-consuming method and hence has limited applications in clinical practice.

In order to overcome the drawback of traditional algorithms, we combine the GAD with the LBM and propose the GAD-LBM model, leveraging the advantages of the GAD on edge preservation and the advantages of the LBM on fast parallel implementation. The GAD-LBM model achieves excellent results on both simulated and clinical images. Visual observation and quantitative comparison has revealed the improvement of the noise reduction, detail preservation, and calculation speed of the GAD-LBM compared with other nine methods. However, whether the despeckling method can be valuable for disease diagnosis has to be further validated in a future study. In addition, compared with the recently published methods, such as the nonlocal total variation ^[39] and other methods surveyed in ^[40], the GAD-LBM method has the advantage of computation speed on the premise of ensuring the quality of denoising. In our future work, we would combine these recent methods with AD and LBM and expect better results. We would also introduce higher order lattices ^[41,42] into the LBM models in future studies to further improve the despeckling performance.

The advantage of the GAD-LBM over the GAD is the computation speed (Figure 4 and Table 1). For synthetic and clinical images, the computation time of the GAD-LBM is about 1/90 to 1/20 of the GAD, demonstrating the superiority of the GAD-LBM in computational efficiency. From Figure 3, we have also found that although the PSNR and MSSIM are almost the same between the GAD-LBM and GAD, the FOM-values are different. It may be attributed to the stability of different numerical solutions to the diffusion. The LBM may derive more stable results for diffusion thus yield better results in FOMs than the finite difference.

For clinical US images, the GAD-LBM has showed a better balance between noise removal and detail preservation than other methods, and it has also demonstrated advantages in computational speed. Thus compared with other despeckling methods, the GAD-LBM can produce denoised images that could be more suitable for subsequent image analysis and interpretation, such as lesion detection, segmentation, feature extraction and intelligent diagnosis ^[43,44].

5.   Conclusions

In this paper, we propose a new model called the GAD-LBM for despeckling US images. The experiments of filtering simulated and clinical US images have demonstrated that the GAD-LBM not only has excellent ability of noise reduction and edge preservation but also has advantages in computation efficiency.

Acknowledgments

The work was funded by the National Science Foundation of China (Nos. 61671281, 61911530249, 81627804).

Conflict of interest:

The authors declare that they have no conflict of interest.