Rayleigh-maximum-likelihood bilateral filter for ultrasound image enhancement

Ultrasound imaging has been used as one of the most prevalent diagnostic techniques due to its advantage of being non-invasive, portable and cost-effective. However, ultrasound images are affected by many types of artifacts, therefore it is hard for an observer to interpret the images and obtain quantitative information from them. Noise in ultrasound can be modeled as the combined effect of two components: one is additive, such as electronic and thermal noise, and the other is multiplicative, called “speckle”. Speckle is the result of the constructive and destructive coherent summation of ultrasound echoes when ultrasound pulses randomly interfere with objects of comparable size to the sound wavelength and then the superposition of acoustical echoes produces an intricate interference pattern [1]. Noise and speckle, considered as undesirable consequence of the image formation process in coherent imaging, directly impact the visualization of the ultrasound image by the physician, deteriorate the quality and the perceivable resolution of diagnostically important features and thus lead to inaccuracy in clinical diagnosis. Therefore, it is essential to remove noise and speckle in the ultrasound images without compromising important image details.

Prior to noise and speckle elimination, the major differences between noise and speckle in ultrasound image should be investigated. Besides the additive nature of noise and multiplicative nature of speckle, their statistical characteristics are completely different in that noise and speckle are subject to Gaussian distribution and Rayleigh distribution, respectively [1]. Figure 1 illustrates the Gaussian distribution whose mean is set as 5 and the variance is set as 1 as well as the Rayleigh distribution whose variance is set as 1. The difference should be considered for effective enhancement.

To alleviate the negative effects of noise and speckle, many efforts have been done to enhance ultrasound images and made the images more valuable after they had been generated and digitized. The state-of-the-art denoising methods for ultrasound images include the median filter [2], adaptive filters such as the Lee filter [3], the Frost filter [4], the Kuan filter [5] as well as the non-local means filter and the anisotropic diffusion method [6].

The median filtering [2] and numerous improved versions [7–9] are often effective for noise reduction. The median intensity of a properly sized and shaped filtering window surrounding the central pixel is used as the output of the target pixel. It thus can eliminate impulsive artifacts whose size is less than a half of the filtering window. Since the amount of smoothing performed by the median filter is determined only by the size of the filtering window, the median filter removes some of the high frequency signal while it results in obscuring the edges. Moreover, the median filter is ineffective when the speckle size is larger than a half of the window size since the filtering scheme does not take into account the statistical characteristics of speckle. Therefore, it undermines the despeckling effectiveness.

Statistical adaptive filters replace the regions within the image whose statistical characteristics similar to those of speckles with the local mean value while keeping the regions with property least similar to speckle unaltered. Therefore they remain effective for speckle suppression. Two of these works are the filters proposed by Lee [3] and Frost [4]. They were first applied on synthetic aperture radar images and then were used on ultrasound images. Lee used the minimum mean square error (MMSE) method to design the filter for additive noise, multiplicative noise and a mixed of the two, its output was estimated by a weighted average based on the mean and the variance of the sub-regions. The frost filter was designed by using an exponentially damped convolution kernel adapted to image fine details and the output was calculated based on local statistics. The Kuan filter [6] is closely similar to the Lee filter but with a different weighting function. Although these filters perform well for removing speckle, they have to compromise between averaging in homogeneous regions and preserving sharp features in the original image.

Non-local means (NLM), a weighted Gaussian filter was proposed for denoising by making use of the high degree of redundancy in the original ultrasound image [10]. NLM performs well in Gaussian noise suppression and sharp edge reservation since it uses the region comparison instead of pixel comparison, which the pattern redundancy is not restricted to be local. However, NL-means cannot be directly applied to ultrasound images since speckle differs from Gaussian noise significantly and is subject to Rayleigh distribution. To extend the application of the NLM method to speckle reduction, relevant NLM-based methods have been studied and proposed. Because these methods determine pixel similarity based on the noisy image patches, speckle in US images will lead to inaccurate computation of the similarity and thus lead to fine detail distortion [11–15]. Furthermore, additional time is required for computing the scale and shape parameters of the distribution of speckle [16].

The anisotropic diffusion was first proposed to reduce noise in images by smoothing in homogeneous regions without blurring the edges [17]. Thereafter, Yu and Acton [6] analyzed the statistical methods for speckle suppression and developed speckle reducing anisotropic diffusion (SRAD), a non-linear and space-variance filter. The SRAD approach reduced speckle by applying isotropic diffusion in homogeneous regions and enhanced edges by inhibiting diffusion across edges, which achieved a balance between despeckling and edge preservation. Flores et al. [18] extended the SRAD to a Log-Gabor guided anisotropic diffusion (ADLG), handling the trade-off between smoothing level and preservation of lesion contour details. Thereafter, a lot of work has been done with anisotropic diffusion equations in such a way that the important structural information can be retained in the denoised images [19–21]. But these SRAD based methods often produce a visually disappointing outputs when they are applied to filter the primary noise contained in ultrasound images, which is subject to Gaussian distribution.

Based on the assumption of Rayleigh distribution of speckle, Aysal [22] first proposed a Rayleigh-trimmed filter for speckle reduction in medical images. For ensuring high efficiency of removing the primary noise and speckle in ultrasound images, two filters are applied to the original image. An alpha-trimmed mean filter was used for suppressing the primary noise and the anisotropic diffusion was subsequently used to further reduce speckle. In addition, Deng et al. [23] proposed a Rayleigh-trimmed anisotropic diffusion filter for speckle reduction in ultrasound images. A Rayleigh-trimmed filter was first applied to estimate the relative standard deviations of local signals and then the anisotropic diffusion was utilized to reduce speckle. However, fine details were simultaneously removed by each filter because speckle detection was not performed before removing.

Elad proposed a new denoising method, named bilateral filter (BF) [24]. It is a non-linear weighted Gaussian filter, taking advantage of adaptive weights based upon spatial and radiometric similarity. Compared with the above-mentioned denoising methods, the BF replaced each pixel by a weighted average of the intensities in the window, which the weighting function gave high weight to those pixels near or similar to the central pixel. Hence it performed well in Gaussian noise reduction and sharp edges preservation. Lin [25] proposed a switching BF (SBF), where the BF could classify a pixel as Gaussian noise, impulse noise or noise-free one. And then the improved SBF switched between the Gaussian and impulse mode depending on the classification result to effectively remove Gaussian noise and impulse noise. However, these BF methods suffered from the drawback that they became ineffective when denoising speckle since the speckle model in ultrasonic images is subject to Rayleigh distribution.

Though these denoising and despeckling approaches operate well in some situations, they have several limitations. Firstly, some approaches do not take into account the statistical characteristics of noise or speckle, undermining the denoising effectiveness. Secondly, most algorithms, for example, the above mentioned approaches, except SBF, do not identify pixel property, such as noise, speckle or edge, before denoising, so they cannot balance effectively in enhancing edges and small structure while reducing noise and speckle, especially when the quality of the original image is poor. Finally, these methods can effectively suppress Gaussian noise or speckle but their performances are not satisfactory in the case of enhancing ultrasound images since ultrasound image is contaminated by addictive background noise and multiplicative speckle.

In order to effectively remove the primary noise and speckle contained in ultrasound images while preserving fine edges and details, a novel and robust method, named as Rayleigh-maximum-likelihood switching bilateral filter (RSBF), which performs the “detect and replace” mechanism before filtering, is proposed. To detect noise, a reference median [25] in the filtering window is first calculated based on the property of the edge in an image, and then the target pixel is identified as noise, speckle or noise-free texture according to the absolute difference between the target pixel and the reference median. Subsequently, noise is removed by the bilateral filter and speckle is suppressed by the Rayleigh-maximum-likelihood filter while noise-free pixels are kept unaltered. The performance and effectiveness of the proposed approach are demonstrated by experiments by using both simulated and clinical ultrasound images.

The remainder of this paper is organized as follows. "Methods" introduces the noise speckle model and the bilateral filter, followed by the detailed description of the proposed RSBF. "Experiments" gives a brief introduction to the experiment environment. "Results and discussion" presents the simulation of the ultrasound images and experiment results on synthetic ultrasound images and real clinical cases in which the proposed method is compared with six state-of-the-art methods for ultrasound speckle reduction. Finally, some conclusions are drawn in "Conclusions".

In this section, the characteristics of noise and speckle is first analyzed. Based on the analysis, a sorted quadrant median vector (SQMV) scheme is performed to classify the target pixel as noise, speckle or noise-free. Thereafter, the bilateral filter is applied to remove noise based on the Gaussian distribution of noise. And the Rayleigh-maximum-likelihood filter is proposed to suppress speckle based on the Rayleigh distribution of speckle. Noise free pixels are kept unchanged in order to enhance images and meanwhile preserve the edge details.

The Rayleigh-maximum-likelihood switching bilateral filter(RSBF)

The Rayleigh-maximum-likelihood bilateral filter (RSBF) consists of two steps. Firstly, the target pixel is identified as noise, speckle or noise-free by using a sorted quadrant median vector [25]. Subsequently, noise and speckle are suppressed by the bilateral filter and the Rayleigh-maximum-likelihood filter, respectively, while noise-free pixels are left unaltered. Figure 2 shows the coarse structure of the proposed method.

Noise/texture detection with the sorted quadrant median vector [25]

1. 1.

   The Sorted Quadrant Median Vector (SQMV)

    

SQMV was proposed for noise detection to avoid false noise detection or blurring image details with inappropriate window size [25]. A large window with the size of \( (2N + 1) \times (2N + 1) \) is divided into four subwindows of size \( (N + 1) \times (N + 1) \), where the central pixel is the corner pixel in the four subwindows, shown in Fig. 3 (\( N = 2 \)).

Let \( C_{i,j} \) be the intensity of the central pixel in the large window, then the set of pixels within the window, denoted as \( W \), is defined as:

$$ W = \{ C_{i + x,j + y} : - N \le x \le N, - N \le y \le N\} $$

(8)

The set of pixels in the four subwindows, denoted as \( W_{1} ,W_{2} ,W_{3} \) and \( W_{4} \), can be expressed as

$$ W_{1} = \{ C_{i + x,j + y} :0 \le x \le N,0 \le y \le N\} $$

(9)

$$ W_{2} = \{ C_{i + x,j + y} : - N \le x \le 0,0 \le y \le N\} $$

(10)

$$ W_{3} = \{ C_{i + x,j + y} : - N \le x \le 0, - N \le y \le 0\} $$

(11)

$$ W_{4} = \{ C_{i + x,j + y} :0 \le x \le N, - N \le y \le 0\} $$

(12)

In case that \( N = 2 \) then the window size is \( 5 \times 5 \) and there are four sub-windows with size of \( 3 \times 3 \). The sequence of the sub-windows is shown in Fig. 3. The median value of each sub-window is denoted as:

$$ m_{j} = median\{ W_{j} \} ,\quad j = 1\;to\;4 $$

(13)

The \( SQMV \) is defined as the four median values in the subwindows in a sorted increasing order:

$$ SQMV = [SQM_{1} ,SQM_{2} ,SQM_{3} ,SQM_{4} ] $$

(14)

where \( SQM_{i} (i = 1,2,3,4) \) is the median of the four subwindows sorted in an increasing order.

1. 2.

   Edge/texture detection based on the \( SQMV \) clusters

    

The edge and texture pixels in a window can be detected based on the difference between the maximal and the minimal \( SQMV \) values:

$$ SQM_{\text{max} } = (SQM_{4} - SQM_{1} ) $$

(15)

$$ SQM_{\text{min} } = (SQM_{3} - SQM_{2} ) $$

(16)

where \( SQM_{\text{max} } \) and \( SQM_{\text{min} } \) represent the maximal and the minimal value of the \( SQMV \), respectively.

If \( SQM_{\text{max} } \) is small, then the pixel intensities in the window are close to each other, so there is no edge. In case that \( SQM_{\text{max} } \) is large but \( SQM_{\hbox{min} } \) is small, then the window contains a weak edge. If \( SQM_{\hbox{max} } \) and \( SQM_{\hbox{min} } \) are both large, then the pixel intensities are different, so the target pixel is identified as strong edge or texture. Based on the analysis, the edge and texture can be detected by:

$$ Edge /texture = \left\{ {\begin{array}{llll} {without \,edge} & & {SQM_{\text{max} } \le T} \\ {weak\, edge} & \hfill & {SQM_{\text{max} } \ge T} \hfill & {\& \, SQM_{\text{min} } \le T} \hfill \\ {strong\, edge/texture} & \hfill & {SQM_{\text{min} } \ge T} \hfill \\ \end{array} } \right. $$

(17)

where T is a threshold.

1. 3.

   Features of edge/texture

    

The cluster of \( SQM_{i} \) and the order of four medians \( m_{j} \) mapping to the clusters implicate the edge features in the window. (1) If \( SQM_{i} \) is clustered as two unequal classes, such as \( \{ (SQM_{1} ,SQM_{2} ,SQM_{3} ,),(SQM_{4} )\} \) or \( \{ (SQM_{1} )(SQM_{2} ,SQM_{3} ,SQM_{4} )\} \) then a diagonal edge is contained in the window. (2) If \( SQM_{i} \) is clustered as two equal classes, \( \{ (SQM_{1} ,SQM_{2} )(SQM_{3} SQM_{4} )\} \), furthermore, \( \{ (m_{2} ,m_{3} )\} \) and \( \{ (m_{1} ,m_{4} )\} \) are associated to \( \{ SQM_{1} ,SQM_{2} \} \) and \( \{ SQM_{3} ,SQM_{4} \} \), respectively, then there is a vertical edge within the window. (3) If \( SQM_{i} \) is clustered as two equal classes,\( \{ (SQM_{1} ,SQM_{2} )(SQM_{3} SQM_{4} )\} \), furthermore, \( \{ (m_{1} ,m_{2} )\} \) and \( \{ (m_{3} ,m_{4} )\} \) are associated to \( \{ (SQM_{1} ,SQM_{2} )\} \) and \( \{ (SQM_{3} ,SQM_{4} )\} \), respectively, then the window contains a horizontal edge.

1. 4.

   Reference median

    

If “without edge” or “weak edge” are detected by Eq. (17), the majority feature of the window is denoted by the medians in the major cluster. Thus the reference median is defined as the average of \( SQM_{2} \) and \( SQM_{3} \). When there is edge or texture in the window, a directional average, used to determine which subwindow the target pixel is more similar to, is calculated by the average of the four pixels in the major pattern depending on the feature of edge or texture, expressed by \( ref \).

$$ ref = \left\{ {\begin{array}{*{20}l} {(SQM_{2} + SQM_{3} )/2} \hfill & {without\, edge\quad \& \;weak\, edge} \hfill \\ {\left( {\sum\limits_{i = 1}^{4} {x_{i} } } \right)/4} \hfill & {strong\, edge /texture} \hfill \\ \end{array} } \right. $$

(18)

where \( x_{i} (i = 1,2,3,4) \) are the four pixels in the major pattern, defined according to the feature of edge or texture, shown in Fig. 4.

If \( ref \) is close to \( (SQM_{1} ,SQM_{2} ) \) then the reference median is defined as \( SQM_{2} \). If \( ref \) is close to \( (SQM_{3} ,SQM_{4} ) \) then the reference median is defined as \( SQM_{3} \). The reference median (\( SQMR \)) can be expressed as:

$$ SQMR\begin{array}{*{20}c} = & {\left\{ {\begin{array}{*{20}l} {(SQM_{2} + SQM_{3} ) /2} \hfill & {SQM_{\text{min} } \le T} \hfill &{} \hfill \\ {SQM_{2} } \hfill & {SQM_{\text{min} } \ge T} \hfill &{\& \, ref \in (SQM_{1} ,SQM_{2} )} \hfill \\ {SQM_{3} } \hfill & {SQM_{\text{min} } \ge T} \hfill & {\& \, ref \in (SQM_{3} ,SQM_{4} )} \hfill \\ \end{array} } \right.} \\ \end{array} $$

(19)

1. 5.

   Noise and edge/texture identification

    

The noise, speckle and edge/texture detection mechanism is implemented by the reference median and the threshold \( T \). When the target pixel is greatly different from the reference median then it is detected as speckle generated by ultrasound echoes [1]. When the target pixel is close to the reference median then it may be a Gaussian background noise or edge/texture. Then the noise and edge/texture is further identified based on the clusters of \( SQMV \) by using Eq. (17). The noise/speckle and edge/texture detection algorithm is shown as follows:

$$ noise /speckle\;detector = \left\{ {\begin{array}{lll} {speckle} \hfill & {\left| {x_{i,j} - SQMR} \right| \ge T_{1} } \hfill \\ {noise} \hfill & {\left| {x_{i,j} - SQMR} \right| \ge T_{2} } \hfill \\ {noise{ - }free} \hfill &\quad{otherwise} \hfill \\ \end{array} } \right. $$

(20)

where \( T_{1} \) and \( T_{2} \) are thresholds determined by experiments.

The Rayleigh-maximum-likelihood filter

Speckle in ultrasound images is modeled as a multiplicative form defined by Eq. (2). The estimation of uncorrupted image \( f \) can be adopted by the robust maximum likelihood (ML) approach [22]. \( f(i,j) \) is assumed to be a constant in the observation window \( \varOmega \) from which \( f(i,j) \) is to be estimated [22], i.e., \( f(i,j) \approx \beta \) for \( \forall (i,j) \in \varOmega \). Hence, for \( \{ u(i,j):(i,j) \in \varOmega \} \) defined in Eq. (2) can be calculated by using:

$$ u(i,j) = \beta *\eta (i,j) $$

(21)

where \( (i,j) \) denotes the coordinates of a pixel in the ultrasound image. \( u \), denoting the image corrupted by speckle, is a Rayleigh distribution with parameter \( \sigma_{I} \). The ML estimation of \( \sigma_{I} \), represented as \( \hat{\sigma } \), is given by

$$ \hat{\sigma }_{I} = \mathop {\arg \hbox{max} }\limits_{\beta } \psi (\vec{u}|\sigma_{I} ) $$

(22)

where \( \vec{u} \) denotes the pixel values of \( u(i,j) \) in the observation window \( \varOmega \) arranged in a vector format. The solution to Eq. (22) is

$$ \hat{\sigma }_{I} = \left( {\frac{1}{2\left\| \varOmega \right\|}\sum\limits_{(i,j) \in \varOmega } {u^{2} (i,j)} } \right)^{1/2} $$

(23)

where \( \left\| \varOmega \right\| \) denotes the cardinality of \( \varOmega \). For example, \( \left\| \varOmega \right\| \) = 9 in a \( 3 \times 3 \) window. Then the ML estimation of the uncorrupted image \( f(i,j) \), denoted as \( \hat{f}(i,j) \), can be calculated as

$$ \begin{aligned} \hat{f}(i,j) & = \hat{\beta } = \hat{\sigma }_{I} (\sigma )^{ - 1} \\ & = (\sigma )^{ - 1} \left[ {\frac{1}{2\left\| \varOmega \right\|}\sum\limits_{(i,j) \in \varOmega } {u^{2} (i,j)} } \right]^{1/2} \\ \end{aligned} $$

(24)

where \( (i,j) \) denotes the pixel location of the estimated original image.

The Rayleigh-maximum-likelihood bilateral filter

Inspired by the “noise/speckle detection and reduction” scheme and statistics of noise and speckle, we combine the bilateral filter and Rayleigh-maximum-likelihood filter together after performing noise/texture identification. Firstly, the target pixel is classified as noise, speckle or noise-free one. Subsequently, noise is removed by using the bilateral filter and speckle is suppressed by using the Rayleigh-maximum-likelihood filter while the noise free pixels are kept unaltered. The detailed algorithm is described below, which can be directly implemented by MATLAB or C language.

For each pixel \( u_{i,j} \) in the noisy image:

1. 1.

   Take the search window, size of \( (2N + 1) \times (2N + 1) \), and the corresponding subwindows size of \( (N + 1) \times (N + 1) \).

    
2. 2.

   Compute \( SQMV \) by Eq. (14) and detect the edge and texture by using Eq. (17).

    
3. 3.

   Identify edge feature according to the cluster of \( SQMV \).

    
4. 4.

   Calculate \( ref \) by Eq. (18) and \( SQMR \) by Eq. (19).

    
5. 5.

   Detect \( u_{i,j} \) as noise, speckle or noise-free pixel.

    
6. 6.

   If \( u_{i,j} \) is classified as noise, replace the target pixel with a weighted average intensities in the window by using Eq. (5).

    
7. 7.

   If \( u_{i,j} \) is detected as speckle, filter the target pixel with the Rayleigh-maximum-likelihood filter in the window by using Eq. (24).

    
8. 8.

   If \( u_{i,j} \) is noise-free then the original intensity is kept unaltered.

    

In the experimental study, synthetic and clinical ultrasound images are used as test sources to evaluate the performance of the proposed RSBF by comparing it with those of six previously proposed filters, including the switching bilateral filter (SBF), the median filter, the SRAD, the NL-means filter, the Lee filter and the Kuan filter.

The proposed RSBF is implemented by using Matlab 7.1 and the experiments are performed on a PC with 2.66 GHz Intel Core 2 processor. Here the parameter \( \sigma_{S} \) and \( \sigma_{R} \) of the bilateral filter are set as 35 and 40 based on the experiments, respectively. The parameter \( \sigma \) of the Rayleigh-maximum-likelihood filter is set as 1. The size of the observation window is \( 5 \times 5 \).

The higher value of SNR and PSNR indicate better performance of denoising.

The less the RMSE, the better the image quality.

A novel and robust method, RSBF, has been proposed to remove speckle and background noise in ultrasound images by implementing a “detect and replace” two-step mechanism. Firstly, each central pixel in the observation window is classified as noise, speckle or noise-free texture according to the absolute difference between the target pixel and the reference median, which is calculated based on the property of the edge in an image. Subsequently, a Rayleigh-maximum-likelihood filter and a bilateral filter are switched to eliminate speckle, assumed to be Rayleigh distributed, and noise, subjecting to Gaussian distribution, respectively, while keeping the noise-free pixel unaltered. Experiments are performed on synthetic and clinical ultrasound images by comparing seven different despeckling methods. Visual evaluation and three numerical indices are applied to evaluate the performance of the proposed method. Results show that the proposed method performs effectively in speckle and noise suppression as well as edge preservation, and is superior to some well-accepted state-of-the-art filters in despeckling, especially when the image quality is poor.

The first conclusion is that the proposed method yields excellent noise/speckle attenuation and edge enhancement because it detects the target pixel as speckle, noise and noise-free before performing filter. Besides, our method achieves robust performance at various image quality levels, especially when the image is greatly deteriorated because the switched filters take into account the statistics of speckle and noise.

In the application of the RSBF, the parameter \( \sigma_{S} \) and \( \sigma_{R} \) of the bilateral filter and the parameter \( \sigma \) of the Rayleigh-maximum-likelihood filter are set as fixed. It is not adaptive in processing the real clinical ultrasound images. Thus, how to optimize the parameters needs to be further researched.

In addition, even though the proposed method is theoretically suitable for ultrasound images of various organs, such as abdominal ultrasound images, the experiment only performs on one dataset of breast ultrasound images due to the limitation of the data. Therefore, more clinic ultrasound images of various organs should be utilized to test the performance of the proposed method in future research.