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 ultrasound noise model and the bilateral filter
Thermal noise and speckle model
Artifacts in ultrasound images contain additive electronic or thermal noise and multiplicative speckle. When an image is corrupted by thermal noise, each original pixel value is added with a noise value \( n_{i,j} \) produced from a zero-mean Gaussian distribution. Then the noisy image \( u_{i,j} \) is related to the original image \( f_{i,j} \) by:
$$ u_{i,j} = f_{i,j} + n_{i,j} $$
(1)
Speckle is an interfering phenomenon. It occurs when two or more waves traveling to the probe from the scatters interfere with each other, whose severity depends on the relative phase between two overlapping returning ultrasonic echoes. A reasonable trade-off between accuracy and simplicity is to model the speckle as a multiplicative artifact with Rayleigh distribution [26]. The relationship between the uncorrupted signal \( f_{i,j} \) and the observed signal \( u_{i,j} \) contaminated by speckle is:
$$ u_{i,j} = f_{i,j} *\eta_{i,j} $$
(2)
where \( \eta \) is a Rayleigh probability density function with the scale parameter \( \sigma \):
$$ p(\eta ) = \frac{\eta }{{\sigma^{2} }}e^{{\left( { - \frac{{\eta^{2} }}{{2\sigma^{2} }}} \right)}} $$
(3)
The noisy ultrasound image, \( u_{i,j} \) contaminated by noise and speckle can be modeled as:
$$ u_{i,j} = n_{i,j} + f{}_{i,j} * \eta_{i,j} $$
(4)
The bilateral filter
The Bilateral filter was proposed to remove Gaussian noise while preserving edges [24]. Each pixel is replaced by a weighted average of the intensities in the filtering window, where the pixels near or similar to the target one are assigned a high weight.
In a \( (2N + 1)*(2N + 1) \) window, let \( u_{i,j} \) and \( u_{i + s,j + t} \) represent the value of the central pixel and the pixels surrounding the central one, respectively, where \( (i,j) \) and \( (i + s,j + t) \) indicate the location of \( u_{i,j} \) and \( u_{i + s,j + t} \). Then the filtered result of the bilateral filter \( y_{i,j} \) is defined as:
$$ y_{i,j} = \frac{{\sum\nolimits_{s = - N}^{N} {\sum\nolimits_{t = - N}^{N} {W_{G} (s,t)W_{R} (s,t)u_{i + s,j + t} } } }}{{\sum\nolimits_{s = - N}^{N} {\sum\nolimits_{t = - N}^{N} {W_{G} (s,t)W_{R} (s,t)} } }} $$
(5)
where
$$ W_{G} (s,t) = \exp - \frac{{(i - s)^{2} + (j - t)^{2} }}{{2\sigma_{S}^{2} }} $$
(6)
$$ W_{R} (s,t) = \exp - \frac{{(u_{i,j} - u_{i + s,j + t} )^{2} }}{{2\sigma_{R}^{2} }} $$
(7)
\( \sigma_{S} \) and \( \sigma_{R} \) are the two parameters that control the bilateral filter, whose values are determined based on experiments.
The bilateral filter performs well in suppressing Gaussian noise while keeping the edge, but it is hard to remove ultrasound speckle because speckle is a type of multiplicative noise and it follows Rayleigh distribution. In order to effectively filter noise and speckle in ultrasound image, we propose a Rayleigh-maximum-likelihood bilateral filter (RSBF) with noise/speckle detection scheme, discussed in the following.
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.
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.
-
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.
-
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.
-
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)
-
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.
Take the search window, size of \( (2N + 1) \times (2N + 1) \), and the corresponding subwindows size of \( (N + 1) \times (N + 1) \).
-
2.
Compute \( SQMV \) by Eq. (14) and detect the edge and texture by using Eq. (17).
-
3.
Identify edge feature according to the cluster of \( SQMV \).
-
4.
Calculate \( ref \) by Eq. (18) and \( SQMR \) by Eq. (19).
-
5.
Detect \( u_{i,j} \) as noise, speckle or noise-free pixel.
-
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.
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.
If \( u_{i,j} \) is noise-free then the original intensity is kept unaltered.