09-11-2016, 09:51 AM
1467220608-paper.docx (Size: 360.02 KB / Downloads: 5)
Assuming that insonification and the resulting echo signal absorption have sufficiently obtained appropriate dynamic compensation from the ultrasound imaging system, the final ultrasonic envelope signal obtained consists of two parts: the reflected signal of the human body, which is a useful signal, and the noise itself, which is made up of two components, multiplicative noise and additive noise. Multiplicative noise is associated with the principle of ultrasonic signal imaging, which results from a random scattering phenomenon in imaging cell resolution. Additive noise can be considered system noise. The ultrasonic envelope signal fpre(i,j) is modeled as follows:
f(i,j)=g(i,j)*n(i,j)+ w(i,j)(6.1)
Where (i, j)∈Z2 is the two-dimensional spatial coordinates, and the superscript pre is the preliminary signal obtained by system. gpre(i,j) and fpre(i,j) denote the original signal and the observed signal, respectively. npre(i,j) and wpre(i,j) represent the multiplicative and additive components of the noise, respectively, where the npre(i,j) is the main component of noise.
The effect of additive noise wpre(i,j) on the qualities of the medical ultrasound images is less significant than the multiplicative noise npre(i,j), and in order to simplify the model (6.1), the additive noise wpre(i,j) is generally omitted and the following model is obtained
f(i,j)=g(i,j)*n(i,j)(6.2)
Logarithmic amplification converts the model (6.2) into the classical additive noise model as follows:
log(f(i,j) )=log (g(i,j)+lo g(n(i,j) )(6.3)
Since wavelet transformation is a linear transformation, the following model is obtained after two-dimensional discrete wavelet transformation for model (6.3):
F_(l,k)^j= G_(l,k)^j+N_(l,k)^j(6.4)
Where j = 1, 2. . . J, l, k∈Z2. F_(l,k)^j,G_(l,k)^j and N_(l,k)^jrepresent the wavelet coefficients of noisy images, noise-free images and speckle noise, respectively. The superscript j represents the decomposition layers of wavelet transformation, and the subscripts (l, k) are wavelet domain coordinates. J denotes the largest decomposition layers.
Since the Bayesian maximum a posteriori estimation will be used to develop a new wavelet shrinkage algorithm, and since the prior probability of noise-free signal and speckle noise is the premise of using the Bayesian maximum a posteriori estimation, we consider that the wavelet coefficients of noise-free signal will obey the generalized Laplacian distribution, and the wavelet coefficients of speckle noise will obey zero mean Gaussian distribution.
The wavelet coefficients of noise-free signal G_(l,k)^j obey the generalized Laplacian distribution, and the probability distribution is:
〖 P〗_G (g)= 1 /(2*s) exp(-|□(g/s)|^v ) (6.5)
Where s is scale parameter, and v is shape parameter.
The wavelet coefficients of speckle noise obey zero mean Gaussian distribution
P_N (n)=□(1 /(√2π σ_N ) exp(-□(n^2/(2σ_N^2 )))) (6.6)
Where σN denotes the standard deviation of noise in wavelet domain.
6.2.2 An improved de-speckling method based on wavelet shrinkage algorithm and fast bilateral filter
On the basis of the traditional wavelet de-noising method, three contributions are made and an improved de-speckling method based on the wavelet shrinkage algorithm with bilateral filter is proposed as shown in figure 6.1.
1) According to the statistical properties of medical ultrasound images in the wavelet domain, an improved threshold function based on the universal threshold function is proposed, which is more consistent with the de-noising of medical ultrasound images.
2) According to the statistical models of noise and signal in the wavelet domain, Bayesian maximum a posteriori estimation is applied to design a new shrinkage algorithm.
3) The speckle noise in the low-pass approximation component is filtered by the bilateral filter to ensure that this de-noising algorithm has a stronger ability to suppress speckle noise in the low-pass approximation component.
An improved wavelet threshold function
Donoho and Johnstone designed a universal wavelet threshold function, i.eT=σ_N √(logM ), where M is the number of the wavelet coefficients in the corresponding wavelet domain. However, when M becomes very large, the bigger threshold could smooth out some useful information, and thus this threshold function is ineffective on the noise removal in medical ultrasonic images. The following improvement is designed in order to obtain the desirable effect:
T_j=〖a_j σ〗_N √(logM ) (6.7)
Where j = 1, 2. . . J are the decomposition layers of wavelet transformation, J denotes the largest decomposition layers. aj represents the adaptive parameter of j layer, and is determined experimentally. After wavelet decomposition, the wavelet coefficients in different decomposition layers have different distribution, and thus the selection of aj should be based on the j, and aj is selected as 1/ln(j + 1). But this choice of aj is not the optimal one, and if aj is selected appropriately, the proposed method will reflect more advantageously.
6.2.2.2 An improved wavelet shrinkage algorithm
The wavelet coefficients of noise-free signal G_(l,k)^j and speckle noise N_(l,k)^j obey the generalized Laplasse distribution and Gaussian distribution, respectively. In order to simplify the calculation, v is selected as 1, so Eq. (6.5) can be simplified as follows:
〖 P〗_G (g)= 1 /(2*s) exp(-|□(g/s)| )(6.8)
Bayesian maximum a posteriori estimation is used to obtain estimation of signal in the wavelet domain. In the process of calculating posterior probability, the following Bayesian equation is used:
P_(G|F) (g|f)=1/(P_F (f) ) P_(F|G) (f|g).P_G(g) =1/(P_F (f) ) P_N (f-g).P_G(g) (6.9)
With Equations (6.6) and (6.8), Equation (6.9) can be restated as:
P_(G|F) (g|f)= 1/(P_F (f) ) □(1 /(2√2π σ_N ) exp(-□((〖2σ〗_N^2 |g|-s(f-g)^2)/(2σ_N^2 ))) )(6.10)
In order to get the maximum a posteriori probability, the first-order derivative equal to zero of ln(p_(G|F) (g|f)) with respect to g leads to Eq. (6.11):
g ̂=sign(f)*(|f|-(σ_N^2)/s)(6.11)
Where g ̂ is the estimation of g, f is assumed in phase with noise-free signal g.
Finally, according to the above equation, a new wavelet shrinkage algorithm based on the improved threshold function and shrinkage algorithm is proposed:
0 ; f≤T_j
g ̂=
sign(f)*(|f|-(σ_N^2)/s) ; f>T_j(6.12)
6.2.2.3 The bilateral filter
The main idea of wavelet de-noising method is to retain the wavelet coefficients of the low-pass component (LL), while wavelet coefficients of high-pass component (LH, HL, HH) are shrunk with the wavelet threshold function. However, only using the wavelet transformation de-noising algorithm to suppress speckle noise in medical ultrasonic images does not produce optimal results. In the experiment, it can be found that the wavelet coefficients in the low-pass component still contain some speckle noise. In order to effectively eliminate speckle noise in the low-pass component, the bilateral filter is selected to filter the speckle noise.
The bilateral filter output at each pixel is a weighted average of its neighbors. The weight is assigned to each neighbor, which decreases with both the distance of the image plane (the spatial domain S) and the distance on the intensity axis (the range domain R). When a Gaussian Gσ is used as a decreasing function, the structure of traditional bilateral filter is given as the following:
I_p^b= 1/(w_p^b ) ∑_(q∈s)▒〖G_(σ_s ) (‖p-q‖ ) G_(σ_r ) (‖I_p-I_q ‖ ) I_q 〗(6.13)
w_p^b=∑_(q∈s)▒〖G_(σ_s ) (‖p-q‖)G_(σ_r ) (‖I_p-I_q ‖)〗normalizes the sum of the weights.
Where I is the input image and Ib is the result of the bilateral filter. The parameter σ_s defines the size of the spatial neighborhoodto filter a pixel, and σ_r controls the weight with its adjacentpixels according to the intensity difference.
The bilateral filter is controlled by two parameters: σsand σr.
• As the range parameter σrincreases, the bilateral filter gradually approximates Gaussian convolution more closely because the range Gaussian G_(σ_r )widens and flattens, i.e., is nearly constant over the intensity interval of the image.
• Increasing the spatial parameter σssmooths larger features.
In practice, in the context of denoising, Liu et al. show that adaptingthe range parameter σrto estimates of the local noise level yieldsmore satisfying results. The authors recommend a linear dependenceσr= 1.95.σn, where σnis the local noise level estimate. An important characteristic of bilateral filtering is that the weights are multiplied: if either of the weights is close to zero, no smoothing occurs. As an example, a large spatial Gaussian coupled with narrow range Gaussian achieves limited smoothing despite the large spatial extent. The range weight enforces a strict preservation of the contours.