- Email: [email protected]
A New Speckle Reducing Anisotropic Diffusion for Ultrasonic Speckle
A New Speckle Reducing Anisotropic Diffusion for Ultrasonic Speckle LI Can-Fei1
WANG Yao-Nan1
XIAO Chang-Yan1
LU Xiao1
Abstract In the area of speckle noise suppression for ultrasound image, a new speckle reducing anisotropic diffuse (NSRAD) is proposed by introducing a sigmoid function to the diffusion coefficient in the traditional speckle reducing anisotropic diffuse (SRAD). The basic route of NSRAD is to combine the image characteristics of different regions using segmented diffusion coefficient, in which the diffusion coefficient approaches to a constant number one in the homogeneous regions, and declines rapidly in the transition regions, then approaches to zero in the strong boundary regions. Benifited by the sigmoid function in the new diffusion function, the NSRAD has better noise-suppressing ability, as well as the abilities of retaining details and weak edges, even sharping the strong boundary, compared to the traditional SRAD. Moreover, the control ratio of the homogeneous regions and the speed control coefficient are adjustable, which well meet different application requirements by solving the problem of different homogeneous regions scales in the ultrasound image. The performances of the proposed model are verified in experiments with synthetic and clinical images. Especially, a quantitative comparison between NSRAD and SRAD is made upon simulation data, where we demonstrate that NSRAD suppresses speckle more efficiently in homogenous regions hereby has less distortion, while in inhomogeneous regions it enhances details and sharpens the edges more effectively than SRAD. Key words Anisotropic diffusion (AD), speckle reduction, ultrasonic image, local statistics Citation Li Can-Fei, Wang Yao-Nan, Xiao Chang-Yan, Lu Xiao. A new speckle reducing anisotropic diffusion for ultrasonic speckle. Acta Automatica Sinica, 2012, 38(3): 412−418 DOI 10.1016/S1874-1029(11)60301-7
Medical ultrasound (US) has been widely employed for imaging human organs and tissues such as the heart, kidney, breast and abdomen, since it is nonradioactive, noninvasive and less expensive, and also allows real-time imaging. Despite of these merits, the diagnostic usage of ultrasound imaging is still limited by its low quality of image[1] . The low image quality mainly results from the presence of signal-dependent noise known as speckle. Speckle typically shows up as noise, and it reduces image contrast and obscures image details. Besides, speckle also affects human interpretation of the acquired ultrasound images and degrades the speed and accuracy of ultrasound image processing tasks. To improve the quality, many efforts have been made to despeckle the image. For instance, various despeckle filtering methods based on local statistics have been proposed by Lee[2−4] , Frost et al.[5] , and Kuan et al.[6−7] . Essentially, both Lee and Kuan filters form an output image by computing a linear combination of the pixel intensity within a moving window. They are sensitive to the size and shape of the window, and different window sizes may greatly affect the quality of the processed images. The response of the Frost filter varies locally with the coefficient of variation. In case of low coefficient of variation, the filter is more averaging-like, and in cases of high coefficient of variation, the filter attempts to sharpen features by preserving the original intensity. Recently, several despeckling methods, including anisotropic diffusion[8−13] , speckle reducing anisotropic diffusion (SRAD)[14] and coherence anisotropic diffusion[15] , were proposed with nonlinear filtering techniques. Those Manuscript received January 28, 2011; accepted September 27, 2011 Supported by National High Technology Research and Development Program of China (863 Program) (2007AA04Z244, 2008AA04Z214), National Natural Science Foundation of China (61172160, 60835004, 60872130), the Fundamental Research Funds for the Central Universities of China (53110704032), and Hunan Province Natural Science Foundation (09J13118) Recommended by Associate Editor LIU Yi-Jun 1. College of Electrical and Information Engineering, Hunan University, Changsha 410082, China
diffusion methods allow for simultaneously performing contrast enhancement and noise reduction by using the coefficient of variation[14] . They can preserve or even enhance prominent edges when removing speckle. The SRAD exploits the instantaneous coefficient of variation to reduce speckle, which is based on the relation between Lee and Frost filters and the anisotropic diffusion (AD) equation. The SRAD can eliminate speckle without distorting useful image information and without destroying the important image edges. Moreover, the SRAD provides superior performance in comparison to the conventional anisotropic diffusion, the enhanced Lee filter and the enhanced Frost filter in terms of smoothing uniform regions and preserving edges and features. Despite the fact that SRAD achieved great progress, it causes blocks due to an imperfect design of the diffusion coefficient. To remedy the shortcoming of conventional SRAD, we propose a new speckle reducing anisotropic diffuse (NSRAD). NSRAD is developed from SRAD whereas a different diffusion coefficient is suggested. We show that NSRAD suppresses speckle more efficiently while preserving sharp edges and detailed features. Of particular note is that the detail preserving anisotropic diffusion (DPAD) in [16] and the oriented SRAD in [17] are the extentions of SRAD, which improve a lot in preserving edges and detailed features. The remainder of this paper is organized as follows. In Section 1 we briefly review the speckle models and SRAD. Then, we introduce NSRAD in Section 2, followed by the experimental results in Section 3, in which we quantify the performance of our model on simulation and real datasets. In Section 4 we conclude the paper.
1
Speckle reducing anisotropic diffusion
A generalized model of the speckle imaging, as proposed in [18], can be written as g = fu + v
(1)
LI Can-Fei, et al.: A New Speckle Reducing Anisotropic Diffusion for · · ·
No. 7
where g is the observed image, u indicates the original image, f and v are the multiplicative and additional noise, respectively. Yu et al.[14] compared the Lee filter with the anisotropic diffusion filter proposed by Rerona et al.[13] , leading to a system of nonlinear partial differential equations for smoothing speckle image g = f u: ⎧ ∂u ⎪ ⎪ ⎪ ⎨ ∂t = div[c(q)u(x, y; t)] (2) u(t = 0) = u0 = g ⎪ ⎪ ∂u ⎪ ⎩ = 0, u ∈ ∂Ω n ∂n
413
depressed with the increase of z, and this will be unfavorable to the preserving of the weak edges and to distinguish the details for there are usually a larger z in the weak edges and the details than in the homogeneous regions.
where ∂Ω denotes the border of Ω, n is the outer normal to the ∂Ω, and c(q) can be taken as c1 (q) =
or
1 q 2 (x, y; t) − q02 (t) 1+ 2 q0 (t)(1 + q02 (t))
2 q (x, y; t) − q02 (t) c2 (q) = exp − 2 q0 (t)(1 + q02 (t))
(3)
(4)
where q0 (t) is the speckle scale function, q(x, y; t) is the instantaneous coefficient of variation determined by 2 2 2 1 |u| 1 u − 2 2 u 4 u (5) q(x, y; t) =
2 2
1 u 1+ 4 u In SRAD, the instantaneous coefficient of variation q(x, y; t) serves as the edge detector in speckled image, and this function presents large values at edges or in the area with high-contrast features while has small values in homogeneous regions. To illustrate the diffuse effect, we carry out analysis on both diffusion coefficients. First, we introduce z = q 2 and z0 = q02 , then c(z) in (3) and (4) can be respectively simplified to 1 (6) c1 (z) = z(x, y; t) − z0 (t) 1+ z0 (t)(1 + z0 (t)) and
z(x, y; t) − z0 (t) c2 (z) = exp − z0 (t)(1 + z0 (t))
(7)
When z is small, it refers to the homogenous region, and when z is large, it stands for the edge, while in between, it is the transition area. For simplicity, we further normalize above diffusion coefficients by dividing the corresponding maximum values, and the final normalized diffusion coefficients are plotted in Fig. 1. z0 stands for noise, which is generally small in a real image, thus we set z0 = 0.002 in this plot. In Fig. 1, the dashed curve is for c1 (z) and the one marked by circle is for c2 (z). Both curves show that c1 (z) and c2 (z) are decreasing and concave functions. At least this will lead to two disadvantages: 1) The diffuse coefficient is very sensitive to z in the small value, and this means that the diffuse coefficient c(z) varies greatly when z has a small deviation, that is, it will magnify the deviation. For the reason that the small z is usually corresponding to the homogeneous regions, and the diffuse coefficient is related to diffuse speed, so it will affect block in homogeneous regions. 2) The sensitivity of the diffuse coefficient will be
Fig. 1 Comparison of SRAD and NSRAD as functions of the diffuse coefficient c(z) with z0 = 0.002, where K = 40, and β = 0.8 in NSRAD
Next, we consider different values of z0 (t). From (3) and (4), one can find that when z 2 (x, y, t) = z02 (t), we have c(x, y, t) = 1, and when z 2 (x, y, t) < z02 (t) we have c(x, y, t) > 1; otherwise c(x, y, t) < 1. Thus z02 (t) provides a reference point for the diffusion effect within various regions. Clearly, the value of reference point z0 (t) is differently obtained by different methods, and such differences could be very large, that is, z0 (t) obtained in SRAD is unstable and insecure, so it is very important to estimate noise accurately. In order to estimate noise particularly for the recursive filtering, considering the fact that noise varies from one iteration to another[19] and converges finally, AjaFernandez et al.[16, 19] proposed a method as follows: z02 = 0.5 ∗ (C ∗ median(M AD(∇ log(u))))2
(8)
where C = 1.4826, median is the median function, which is to return the median value, and M AD(∇ log(u)) = |∇ log(u) − median(∇ log(u))|
(9)
Even though the reference point is exact and reliable, i.e., we presume z0 is very accurate, the mode SRAD is also worth discussing. To illustrate, we consider the relation between the diffuse coefficient c(z) and the reference point z0 . In Fig. 2 (a) is the origin image for plotting a profile in line 65, as highlighted with the vertical line. Fig. 2 (b) is the filtered result by SRAD. Fig. 2 (c) is the gray value corresponding to the profile. The instantaneous coefficient of variation z of the section plane is shown in Fig. 2 (d), where z0 computed by (8) is also shown for comparison. Fig. 2 confirms that, in the homogeneous regions, where z ≤ z0 , as analyzed before, the diffuse coefficient is very sensitive to z in such regions, and this means that the diffuse coefficient c(z) varies greatly when z has a small deviation, that is, it will magnify the deviation. Since the diffuse coefficient is related to the diffuse speed, so it will affect block in the homogeneous regions. Moreover, in the region of the weak edges and details, where z > z0 , as we pointed out before, the diffuse coefficient is less sensitive,
414
ACTA AUTOMATICA SINICA
and this shows that it is disadvantageous for preserving the weak edges and to distinguish the details. We show these in Fig. 2 (b), which is the diffused image by SRAD. There exist apparent blocks in homogeneous regions, as labeled by arrows. Moreover, the details and weak edges are blurred even for the strong edges.
Vol. 38
value to get
c(q) = 1 − 2.2
1 1 + exp(−k(q 2 − βq02 ))
1 + exp(kβq02 ) 2 exp(kβq0 ) (12)
The diffuse coefficient
To illustrate the diffuse effect, we carry out analysis on the diffusion coefficients. If we introduce z = q 2 and z0 = q02 , then c(z) in (12) simplifies to
1 + exp(kβz0 ) 1 c(z) = 1 − 1 + exp(−k(z − βz0 )) exp(kβz0 ) (13) And its derivative reads c (z) =
(b)
(a)
−k exp(−k(z − βz0 ))[1 + exp(kβz0 )] [1 + exp(−k(z − βz0 ))]2 exp(kβz0 )
(14)
As seen from (14), the derivative of c(z) depends on k, i.e., k is the variable to control the drop speed. Moreover, c(z) is also associated with variables z and βz0 . In some way, z0 may be regard as an average of z. So we choose parameter z as a self-correcting parameters to noise estimation z0 as follows: K (15) k= z0 (c)
(d)
Fig. 2 An example analysis for SRAD ((a) Original live ultrasound image; (b) Filtered images by SRAD, Δt = 0.05, 500 iterations; (c) The gray value in the profile; (d) The instantaneous coefficient of variation z of the section plane and the reference point z0 .)
2 2.1
The new speckle reducing anisotropic diffusion The new speckle reducing anisotropic diffusion
To make up the shortcoming of SRAD as analyzed before, an ideal diffuse coefficient should act so as to diffuse homogeneous regions with speeds changing stolidly but the transition regions with speeds changing sensitively. This can be realized by a sigmoid function. So a new speckle reducing anisotropic diffusion (NSRAD) can be written as follows: ⎧ ∂u ⎪ ⎪ ⎪ ⎨ ∂t = div[c(q)u(x, y; t)] (10) u(t = 0) = u0 = g ⎪ ⎪ ∂u ⎪ ⎩ = 0, u ∈ ∂Ω n ∂n with a new diffuse coefficient c(q) = 1 −
1 1 + exp(−k(q 2 − βq02 ))
(11)
where k > 0, it is a tunable parameter to control the dropping speed of the diffusion coefficient. In (11), β is the rate of the homogeneous regions to the noise estimation q02 (t), and we name it as the homogeneous region control coefficient. Note that the homogeneous region in an image is not unanimous to different observers, thus the homogeneous region scale is also different to different observers, which can be achieved by the homogeneous region control coefficient β. Similarly, we normalize c(q) by dividing its maximum
where K is a suitable constant. The new diffuse coefficient is plotted in Fig. 1. Here, we set z0 = 0.002 as before in SRAD, and K = 40, β = 10. The new diffuse coefficient is a sigmoid function and there are three parts in it, which correspond to three regions. The first part is z < βz0 , which corresponds to the homogeneous region. The second part is z ≈ βz0 , which corresponds to the transition region. The last part is z > βz0 , which corresponds to the edge region. In the first part (z < βz0 ), i.e., in the homogeneous region, the new diffuse coefficient c(z) ≈ 1. This means that the diffuse is almost uniform in the homogeneous region. Moreover, the diffuse coefficient c(z) changes slowly in this part. It means that c(z) varies small when z has a small deviation, which will not magnify the deviation of z in this region. So it will not affect block in the homogeneous region. In the second part (z ≈ βz0 ), i.e., in the transition region, the diffuse coefficient decreases rapidly. It means that c(z) varies greatly when z has a small deviation in this part, which will magnify the deviation of z in this region. So it will enhance the resolution power in the transition regions and improves the sharp of weak edges and details. In the third one (z > βz0 ), i.e., in the edge region, the diffuse coefficient c(z) tends to 0 fast, which will strength the sharp and the strong of the edges. 2.3
Analysis of the scale control coefficient β
The homogeneous region scale control coefficient β plays a crucial role in NSRAD. If we ignore the impact of k to β, the diffuse is almost uniform under scale β. The block will not occur in the homogeneous region for there is almost a uniform diffuse speed under scale β. It is well know that one homogeneous region under larger scale (larger β) will be divided into one new homogeneous region and one transition region under smaller scale (smaller β). So it is useful to preserve more details with a smaller β. On the contrary, one homogeneous region and one transition region under smaller scale (smaller β) will combine into one new homogeneous region under larger scale (larger β), this will blur details with a larger β. So, a different result will be gotten by a different scale β, that is, the diffusing can be adjusted
LI Can-Fei, et al.: A New Speckle Reducing Anisotropic Diffusion for · · ·
No. 7
by β, which is useful to fit different application. This will be further explainded by the function image of the diffuse coefficient. The diffuse coefficient function with different values of β (β = 1.0 and β = 1.2 , respectively) are plotted in Fig. 3 (b). From Fig. 3 (b), the homogeneous region under scale β = 1.2 can be seen as a new homogeneous region and a transition region under scale β = 1.0. This means that the different diffusion result will be got in the region if use different scales of β. Moreover, from Fig. 3 (b), the curve of the diffusion coefficient c(z) is steeper in the transition region when β is larger, so the diffusion coefficient c(z) is more sensitive to z in this one, which will improve the sharpness of the weak edges and enhance the details. On the contrary, the curve of the diffusion coefficient c(z) is more gently in the transition region when β is smaller, so the diffusion coefficient c(z) is more stolid to z in this one, which will blur the weak edges and details. All of these will affect the diffusion result. 2.4
(a)
n qi,j
Fig. 3 Analysis of the parameters in the diffuse coefficient for NSRAD ((a) The function image of c(z) in NSRAD with different values of K: solid line for K = 10 and dotted line for K = 40, respectively, in which β = 1.0, and z0 = 0.002; (b) The function image of c(z) in NSRAD with different values of β: solid line for β = 1.0 and dotted line for β = 1.2, respectively, where K = 40, and z0 = 0.002.)
Discretization
The partial differential (10) is solved numerically using a Jacobi iterative method[14] . Choosing a time-step Δt and a grid-size h in both x and y directions, discretize the time and space coordinates as t = nt, n = 0, 1, · · · , x = ih, y = jh, i = 0, 1, · · · , M − 1, j = 0, 1, · · · , N − 1, where M h × N h is the area of the image domain. Let un i,j = u(ih, jh; nt); an numerical approximation to (10) is given by 1 tdn i,j 4
2 2 n 2 1 ∇ ui,j |∇un 1 i,j | − 2 2 un 4 un i,j i,j =
2 ∇2 un 1 i,j
1+ 4 un i,j
(16)
(18)
and q02 = 0.5 ∗ (C ∗ median(M AD(∇ log(u))))2 ∇ R un i,j =
(b)
n un+1 i,j = ui,j +
1 n n n n n [ci+1,j (un i+1,j − ui,j ) + ci,j (ui−1,j − ui,j ) + h2 n n n n n cn (17) i,j+1 (ui,j+1 − ui,j ) + ci,j (ui,j−1 − ui,j )]
with
Analysis of the rate of the parameter
In the NSRAD, for k = K/z0 , to explain the influence of k is equal to explain the influence of K. The diffuse function images with different values of K (K = 10 and K = 40, respectively) are plotted in Fig. 3 (a), in which β = 1.0 and z0 = 0.002. As seen from Fig. 3, parameter k not only affects the rate of decline function but also affects the actual role of the homogeneous region scale control coefficient β. When K is smaller, it will reduce the value of the actual role of the scale control coefficient β. Moreover, the diffuse function drops more slowerly, which is disadvantageous to the resolution over the transition region. On the contrary, when K is larger, it will enlarge the value of the actual role of the scale control coefficient β. Moreover, the diffuse function drops faster, which is advantageous to the resolution over the transition region.
2.5
dn i,j =
415
∇ L un i,j =
∇ 2 un i,j =
n n un un i+1,j − ui,j i,j+1 − ui,j − h h
n n n un i,j − ui−1,j ui,j − ui,j−1 , h h
(19)
(20)
n n n n un i+1,j + ui−1,j + ui,j+1 + ui,j−1 − 4ui,j h
(21)
(22)
Moreover, to avoid distortion at the image boundaries, symmetric boundary condition is required, that is, the image intensity function has equal values at both sides of the boundary.
3
Experiments
To evaluate the performance of NSRAD, simulation and in vivo studies were conducted. The simulation studies may give quantitative performance analysis, and the in vivo study may demonstrate the feasibility and usefulness of NSRAD in real applications. In all the experiments, c(q) were computed by (3) in SRAD, and by (11) in NSRAD, respectively. In order to eliminate the patchy of the gray further (when needed), in experiments, a Gaussian smoothing was made: q (23) q = u 3.1
Results from simulation study
We ran the SRAD and NSRAD on Lena images corrupted by a multiplicative noise; see Fig. 4 (b), and Fig. 4 (a) is the reference image. Here, we use (23) to eliminate the patchy of the gray. To assure stability in the diffusion process, we took a time step t = 0.1 and ran 150 iterations, and the results for SRAD and NSRAD with K = 2 and β = 1 are shown in Fig. 4 (c) and Fig. 4 (d), respectively. As seem from Fig. 4, NSRAD gives a better diffusion but SRAD makes very strong block affect. At the same time, NSRAD preserves the sharpness of the edges. We repeated this experiment for greater values of Δt (see Fig. 5), such as Δt = 0.5 and Δt = 1, and Δt = 1 only for NSRAD. The result that seems not to match the statement in [14], and the fact that the diffusion coefficients observed with this filter may at some point be greater than could explain this effect[16] . But NSRAD is stable, even for Δt = 1.
416
ACTA AUTOMATICA SINICA
(a)
(b)
(c)
(d)
Vol. 38
different factors: loss of correlation, luminance distortion, and contrast distortion. The structure similarity (SSIM)[22] evaluates the overall processing quality. We show the MSE (Fig. 6 (a)), the Q (Fig. 6 (b)), the SNR (Fig. 6 (c)), and the SSIM (Fig. 6 (d)) evolutions with time for Δt = 0.1. The index MSE in NSRAD reduced as compared to SRAD, but index SNR increased which shows that the method NSRAD suppresses speckle more effectively. Moreover, note that NSRAD shows better results in other indexes than SRAD, such as SSIM and Q. These advantages even become more prominent with the increasing of iterations, which shows that the filtered image by NSRAD has less distortion and well structure similarity.
(a) MSE
(b) Q
(c) SNR
(d) SSIM
Fig. 4 Simulation experiment ((a) Original image; (b) Noisy image; (c) Filtered by SRAD: ∇t = 0.1, 150 iterations; (d) Filtered by NSRAD: K = 2, β = 1, ∇t = 0.1, 150 iterations.)
Fig. 6
3.2 (a)
(b)
(c)
Fig. 5 Simulation experiment with different ∇t values, 150 iterations ((a) Filtered by SRAD: ∇t = 0.5; (b) Filtered by NSRAD: ∇t = 0.5, K = 2, β = 1; (c) Filtered by NSRAD: ∇t = 1.0, K = 2, β = 1.)
To quantify the achieved performance improvements, four quality index metrics were computed on the filtered image and reference image. Mean squared error (MSE) measures the quality change between the original and processed images. Signal-to-noise ratio (SNR)[20] measures image fidelity. The mathematically defined universal quality index Q[21] measures any distortion as a combination of three
Index for 150 iterations (solid line) NSRAD and (dashed line) SRAD
Results from in vivo study
The SRAD and NSRAD were further tested with two real clinical images. Here, both the time steps were set to Δt = 0.05 and the same 500 iterations were run. Fig. 7 is a live ultrasound image filtering, the original image is Fig. 2 (a). Fig. 7 (a) and Fig. 7 (b) are the filtered image obtained from NSRAD with K = 400, β = 1.0 and K = 500, β = 0.8, respectively. The filtered image obtained from SRAD is shown in Fig. 2 (b). In Fig. 2 (b), SRAD makes very strong block effect, as labeled with arrows. But the proposed NSRAD suppresses speckle noises more effectively and the fine characters (labeled with circles) are better preserved compared with the traditional SRAD. The undesired block effects (labeled with arrows) in the SRAD can not be seen in NSRAD filtered image. Moreover, the fine characters (labeled with circles) in Fig. 7 (b) are increased with a larger K and a smaller β than the ones in Fig. 7 (a) in the proposed NSRAD filtered image. As shown in Fig. 8 (a), the second experiment was conducted by using a prostate ultrsound image, which is different from the previous data in that the objects appear to be irregular and more fine structures are included due to the inherent anatomical composition of soft tissues. Therefore, the speckle suppression of prostate image remains a challenge to many traditional methods. However, from our filtering results in Fig. 8 (c) and Fig. 8 (d), it is obvious that the majority of fine structures have been well preserved as
LI Can-Fei, et al.: A New Speckle Reducing Anisotropic Diffusion for · · ·
No. 7
marked with circles, and the fine characters (labeled with circles) in Fig. 8 (d) are increased with a larger K and a smaller β than the ones in Fig. 8 (c). In addition, the frequently seen block effects of SRAD as labeled by arrows in Fig. 8 (a) and Fig. 8 (b) are totally eliminated in our results.
417
2 Lee J S. Speckle analysis and smoothing of synthetic aperture radar images. Computer Graphics and Image Processing, 1981, 17(1): 24−32 3 Lee J S. Digital image enhancement and noise filtering by use of local statistics. IEEE Transactions on Pattern Analysis and Machine Intelligence, 1980, 2(2): 165−168 4 Lee J S. Refined filtering of image noise using local statistics. Computer Graphics and Image Processing, 1981, 15(4): 380−389
(a) K = 400, β = 1.0
Fig. 7
(b) K = 500, β = 0.8
The live ultrasound experiment by NSRAD with ∇t = 0.05, 500 iterations
5 Frost V S, Stiles J A, Shanmuggam K S, Holtzman J C. A model for radar images and its application to adaptive digital filtering of multiplicative noise. IEEE Transactions on Pattern Analysis and Machine Intelligence, 1982, 4(2): 157−165 6 Kuan D T, Sawchuk A A, Strand T C, Chavel P. Adaptive restoration of images with speckle. IEEE Transactions on Acoustics Speech and Signal Processing, 1987, 35(3): 373− 383 7 Kuan D T, Sawchuk A A, Strand T C, Chavel P. Adaptive noise smoothing filter for images with signal dependent noise. IEEE Transactions on Pattern Analysis and Machine Intelligence, 1985, 7(2): 165−177 8 Wang Yi, Niu Rui-Qing, Yu Xin, Shen Huan-Feng. Time dependent robust anisotropic diffusion processes. Acta Automatica Sinica, 2009, 35(9): 1253−1256 (in Chinese)
(a) Original ultrasound image
(b) Filtered images by SRAD
9 Jin J S, Wang Y, Hiller J. An adaptive nonlinear diffusion algorithm for filtering medical images. IEEE Transactions on Information Technology in Biomedicine, 2000, 4(4): 298−305 10 Meng Xiang-Lin, Wang Zheng-Zhi. Image diffusion based on visual masking effect. Acta Automatica Sinica, 2011, 37(1): 21−27 (in Chinese) 11 Ying Shi-Hui, Peng Ji-Gen, Zheng Kai-Jie, Qiao Hong. Lie group method for data set registration problem with anisotropic scale deformation. Acta Automatica Sinica, 2009, 35(7): 867−874 (in Chinese)
(c) Filtered images by NSRAD: K = 400, β = 0.9
Figure. 8
4
(d) Filtered images by NSRAD: K = 500, β = 0.7
A pancreas ultrasound experiment by SRAD and NSRAD with ∇t = 0.05, 500 iterations
Conclusion
A new anisotropic diffusion model with a sigmoid function based diffusion coefficient definition for ultrasonic speckle reduction has been investigated. The new method is able to smooth homogeneous region while efficiently avoiding block effects. Moreover, the new method strengthens the sensitivity of the transition region to preserver more details and weak edges. Experiments in simulation data and in vivo ultrasound show its better capacity in diffusing homogeneous region and preserving image edges and details than SRAD. References 1 Lamont D, Parker L, White M, Unwin N, Bennett S M A, Cohen M, Dickinson H O, Adamson A, Alberti K G M M, Craft A W. Risk of cardiovascular disease measured by carotid intima-media thickness at age 49-51: life course study. British Medical Journal, 2000, 320(7230): 273−278
12 Black M J, Sapiro G, Marimont D H, Heeger D. Robust anisotropic diffusion. IEEE Transactions on Image Processing, 1998, 7(3): 421−432 13 Rerona P, Malik J. Scale-space and edge detection using anisotropic diffusion. IEEE Transactions on Pattern Analysis and Machine Intelligence, 1990, 12(7): 629−639 14 Yu Y, Acton S T. Speckle reducing anisotropic diffusion. IEEE Transactions on Image Processing, 2002, 11(11): 1260−1270 15 Abd-Elmoniem K Z, Youssef A B, Kadah Y. Real-time speckle reduction and coherence enhancement in ultrasound imaging via nonlinear anisotropic diffusion. IEEE Transactions on Biomedical Engineering, 49(9): 997−1014 16 Aja-Fernandez S, Alberola-Lopez C. On the estimation of the coefficient of variation for anisotropic diffusion speckle filtering. IEEE Transactions on Image Processing, 2006, 15(9): 2694−2701 17 Krissian K, Westin C F, Kikinis R, Vosburgh K G. Oriented speckle reducing anisotropic diffusion. IEEE Transactions on Image Processing, 2007, 16(5): 1412−1424
418
ACTA AUTOMATICA SINICA
18 Jain A K. Fundamental of Digital Image Processing. New Jersey: Prentice Hall, 1989 19 Aja-Fernandez S, Vegas-Sanchez-Ferrero G, MartinFernandez M, Alberola-Lopez C. Automatic noise estimation in images using local statistics: additive and multiplicative cases. Image and Vision Computing, 2009, 27(6): 756−770 20 Sakrison D. On the role of observer and a distortion measure in image transmission. IEEE Transactions on Communications, 1977, 25(11): 1251−1267 21 Wang Z, Bovik A. A universal image quality index. IEEE Signal Processing Letters, 2002, 9(3): 81−84 22 Wang Z, Bovik A C, Sheikh H R, Simoncelli E P. Image quality assessment: from error visibility to structural similarity. IEEE Transactions on Image Processing, 2004, 13(4): 600−612
LI Can-Fei Ph. D. candidate at the College of Electrical and Information Engineering, Hunan University. She received her B. S. and M. E. degrees from the College of Electrical and Information Engineering, Hunan University in 2001 and 2005, respectively. Her research interest covers image recognition, computer vision, and medical image processing. Corresponding author of this paper. E-mail: [email protected]
Vol. 38 WANG Yao-Nan Professor at the College of Electrical and Information Engineering, Hunan University. He received his Ph. D. degree from Hunan University. He was a postdoctor at National University of Defense Technology and an alexander von Humboldt Stiftung. His research interest covers intelligent control, intelligent image processing, and intelligent robotics. E-mail: [email protected] XIAO Chang-Yan Received his B. S. and master degrees from National University of Defense Technology, China in 1994 and 1997, and his Ph. D. degree from Shanghai Jiao Tong University, China in 2004. From 2008 to 2009, he worked as a postdoctor in Leiden University, the Netherlands. Now he is an associate professor at the College of Electrical and Information Engineering, Hunan University, China. E-mail: [email protected] LU Xiao Ph. D. candidate at the College of Electrical and Information Engineering, Hunan University. She received her B. S. degree from the College of Electrical and Information Engineering, Hunan University in 2007, and her master degree from the College of Automation, Southeast University in 2010, respectively. Her research interest covers computer vision, and robot navigation. E-mail: xlu [email protected]
WANG Yao-Nan1
XIAO Chang-Yan1
LU Xiao1
Abstract In the area of speckle noise suppression for ultrasound image, a new speckle reducing anisotropic diffuse (NSRAD) is proposed by introducing a sigmoid function to the diffusion coefficient in the traditional speckle reducing anisotropic diffuse (SRAD). The basic route of NSRAD is to combine the image characteristics of different regions using segmented diffusion coefficient, in which the diffusion coefficient approaches to a constant number one in the homogeneous regions, and declines rapidly in the transition regions, then approaches to zero in the strong boundary regions. Benifited by the sigmoid function in the new diffusion function, the NSRAD has better noise-suppressing ability, as well as the abilities of retaining details and weak edges, even sharping the strong boundary, compared to the traditional SRAD. Moreover, the control ratio of the homogeneous regions and the speed control coefficient are adjustable, which well meet different application requirements by solving the problem of different homogeneous regions scales in the ultrasound image. The performances of the proposed model are verified in experiments with synthetic and clinical images. Especially, a quantitative comparison between NSRAD and SRAD is made upon simulation data, where we demonstrate that NSRAD suppresses speckle more efficiently in homogenous regions hereby has less distortion, while in inhomogeneous regions it enhances details and sharpens the edges more effectively than SRAD. Key words Anisotropic diffusion (AD), speckle reduction, ultrasonic image, local statistics Citation Li Can-Fei, Wang Yao-Nan, Xiao Chang-Yan, Lu Xiao. A new speckle reducing anisotropic diffusion for ultrasonic speckle. Acta Automatica Sinica, 2012, 38(3): 412−418 DOI 10.1016/S1874-1029(11)60301-7
Medical ultrasound (US) has been widely employed for imaging human organs and tissues such as the heart, kidney, breast and abdomen, since it is nonradioactive, noninvasive and less expensive, and also allows real-time imaging. Despite of these merits, the diagnostic usage of ultrasound imaging is still limited by its low quality of image[1] . The low image quality mainly results from the presence of signal-dependent noise known as speckle. Speckle typically shows up as noise, and it reduces image contrast and obscures image details. Besides, speckle also affects human interpretation of the acquired ultrasound images and degrades the speed and accuracy of ultrasound image processing tasks. To improve the quality, many efforts have been made to despeckle the image. For instance, various despeckle filtering methods based on local statistics have been proposed by Lee[2−4] , Frost et al.[5] , and Kuan et al.[6−7] . Essentially, both Lee and Kuan filters form an output image by computing a linear combination of the pixel intensity within a moving window. They are sensitive to the size and shape of the window, and different window sizes may greatly affect the quality of the processed images. The response of the Frost filter varies locally with the coefficient of variation. In case of low coefficient of variation, the filter is more averaging-like, and in cases of high coefficient of variation, the filter attempts to sharpen features by preserving the original intensity. Recently, several despeckling methods, including anisotropic diffusion[8−13] , speckle reducing anisotropic diffusion (SRAD)[14] and coherence anisotropic diffusion[15] , were proposed with nonlinear filtering techniques. Those Manuscript received January 28, 2011; accepted September 27, 2011 Supported by National High Technology Research and Development Program of China (863 Program) (2007AA04Z244, 2008AA04Z214), National Natural Science Foundation of China (61172160, 60835004, 60872130), the Fundamental Research Funds for the Central Universities of China (53110704032), and Hunan Province Natural Science Foundation (09J13118) Recommended by Associate Editor LIU Yi-Jun 1. College of Electrical and Information Engineering, Hunan University, Changsha 410082, China
diffusion methods allow for simultaneously performing contrast enhancement and noise reduction by using the coefficient of variation[14] . They can preserve or even enhance prominent edges when removing speckle. The SRAD exploits the instantaneous coefficient of variation to reduce speckle, which is based on the relation between Lee and Frost filters and the anisotropic diffusion (AD) equation. The SRAD can eliminate speckle without distorting useful image information and without destroying the important image edges. Moreover, the SRAD provides superior performance in comparison to the conventional anisotropic diffusion, the enhanced Lee filter and the enhanced Frost filter in terms of smoothing uniform regions and preserving edges and features. Despite the fact that SRAD achieved great progress, it causes blocks due to an imperfect design of the diffusion coefficient. To remedy the shortcoming of conventional SRAD, we propose a new speckle reducing anisotropic diffuse (NSRAD). NSRAD is developed from SRAD whereas a different diffusion coefficient is suggested. We show that NSRAD suppresses speckle more efficiently while preserving sharp edges and detailed features. Of particular note is that the detail preserving anisotropic diffusion (DPAD) in [16] and the oriented SRAD in [17] are the extentions of SRAD, which improve a lot in preserving edges and detailed features. The remainder of this paper is organized as follows. In Section 1 we briefly review the speckle models and SRAD. Then, we introduce NSRAD in Section 2, followed by the experimental results in Section 3, in which we quantify the performance of our model on simulation and real datasets. In Section 4 we conclude the paper.
1
Speckle reducing anisotropic diffusion
A generalized model of the speckle imaging, as proposed in [18], can be written as g = fu + v
(1)
LI Can-Fei, et al.: A New Speckle Reducing Anisotropic Diffusion for · · ·
No. 7
where g is the observed image, u indicates the original image, f and v are the multiplicative and additional noise, respectively. Yu et al.[14] compared the Lee filter with the anisotropic diffusion filter proposed by Rerona et al.[13] , leading to a system of nonlinear partial differential equations for smoothing speckle image g = f u: ⎧ ∂u ⎪ ⎪ ⎪ ⎨ ∂t = div[c(q)u(x, y; t)] (2) u(t = 0) = u0 = g ⎪ ⎪ ∂u ⎪ ⎩ = 0, u ∈ ∂Ω n ∂n
413
depressed with the increase of z, and this will be unfavorable to the preserving of the weak edges and to distinguish the details for there are usually a larger z in the weak edges and the details than in the homogeneous regions.
where ∂Ω denotes the border of Ω, n is the outer normal to the ∂Ω, and c(q) can be taken as c1 (q) =
or
1 q 2 (x, y; t) − q02 (t) 1+ 2 q0 (t)(1 + q02 (t))
2 q (x, y; t) − q02 (t) c2 (q) = exp − 2 q0 (t)(1 + q02 (t))
(3)
(4)
where q0 (t) is the speckle scale function, q(x, y; t) is the instantaneous coefficient of variation determined by 2 2 2 1 |u| 1 u − 2 2 u 4 u (5) q(x, y; t) =
2 2
1 u 1+ 4 u In SRAD, the instantaneous coefficient of variation q(x, y; t) serves as the edge detector in speckled image, and this function presents large values at edges or in the area with high-contrast features while has small values in homogeneous regions. To illustrate the diffuse effect, we carry out analysis on both diffusion coefficients. First, we introduce z = q 2 and z0 = q02 , then c(z) in (3) and (4) can be respectively simplified to 1 (6) c1 (z) = z(x, y; t) − z0 (t) 1+ z0 (t)(1 + z0 (t)) and
z(x, y; t) − z0 (t) c2 (z) = exp − z0 (t)(1 + z0 (t))
(7)
When z is small, it refers to the homogenous region, and when z is large, it stands for the edge, while in between, it is the transition area. For simplicity, we further normalize above diffusion coefficients by dividing the corresponding maximum values, and the final normalized diffusion coefficients are plotted in Fig. 1. z0 stands for noise, which is generally small in a real image, thus we set z0 = 0.002 in this plot. In Fig. 1, the dashed curve is for c1 (z) and the one marked by circle is for c2 (z). Both curves show that c1 (z) and c2 (z) are decreasing and concave functions. At least this will lead to two disadvantages: 1) The diffuse coefficient is very sensitive to z in the small value, and this means that the diffuse coefficient c(z) varies greatly when z has a small deviation, that is, it will magnify the deviation. For the reason that the small z is usually corresponding to the homogeneous regions, and the diffuse coefficient is related to diffuse speed, so it will affect block in homogeneous regions. 2) The sensitivity of the diffuse coefficient will be
Fig. 1 Comparison of SRAD and NSRAD as functions of the diffuse coefficient c(z) with z0 = 0.002, where K = 40, and β = 0.8 in NSRAD
Next, we consider different values of z0 (t). From (3) and (4), one can find that when z 2 (x, y, t) = z02 (t), we have c(x, y, t) = 1, and when z 2 (x, y, t) < z02 (t) we have c(x, y, t) > 1; otherwise c(x, y, t) < 1. Thus z02 (t) provides a reference point for the diffusion effect within various regions. Clearly, the value of reference point z0 (t) is differently obtained by different methods, and such differences could be very large, that is, z0 (t) obtained in SRAD is unstable and insecure, so it is very important to estimate noise accurately. In order to estimate noise particularly for the recursive filtering, considering the fact that noise varies from one iteration to another[19] and converges finally, AjaFernandez et al.[16, 19] proposed a method as follows: z02 = 0.5 ∗ (C ∗ median(M AD(∇ log(u))))2
(8)
where C = 1.4826, median is the median function, which is to return the median value, and M AD(∇ log(u)) = |∇ log(u) − median(∇ log(u))|
(9)
Even though the reference point is exact and reliable, i.e., we presume z0 is very accurate, the mode SRAD is also worth discussing. To illustrate, we consider the relation between the diffuse coefficient c(z) and the reference point z0 . In Fig. 2 (a) is the origin image for plotting a profile in line 65, as highlighted with the vertical line. Fig. 2 (b) is the filtered result by SRAD. Fig. 2 (c) is the gray value corresponding to the profile. The instantaneous coefficient of variation z of the section plane is shown in Fig. 2 (d), where z0 computed by (8) is also shown for comparison. Fig. 2 confirms that, in the homogeneous regions, where z ≤ z0 , as analyzed before, the diffuse coefficient is very sensitive to z in such regions, and this means that the diffuse coefficient c(z) varies greatly when z has a small deviation, that is, it will magnify the deviation. Since the diffuse coefficient is related to the diffuse speed, so it will affect block in the homogeneous regions. Moreover, in the region of the weak edges and details, where z > z0 , as we pointed out before, the diffuse coefficient is less sensitive,
414
ACTA AUTOMATICA SINICA
and this shows that it is disadvantageous for preserving the weak edges and to distinguish the details. We show these in Fig. 2 (b), which is the diffused image by SRAD. There exist apparent blocks in homogeneous regions, as labeled by arrows. Moreover, the details and weak edges are blurred even for the strong edges.
Vol. 38
value to get
c(q) = 1 − 2.2
1 1 + exp(−k(q 2 − βq02 ))
1 + exp(kβq02 ) 2 exp(kβq0 ) (12)
The diffuse coefficient
To illustrate the diffuse effect, we carry out analysis on the diffusion coefficients. If we introduce z = q 2 and z0 = q02 , then c(z) in (12) simplifies to
1 + exp(kβz0 ) 1 c(z) = 1 − 1 + exp(−k(z − βz0 )) exp(kβz0 ) (13) And its derivative reads c (z) =
(b)
(a)
−k exp(−k(z − βz0 ))[1 + exp(kβz0 )] [1 + exp(−k(z − βz0 ))]2 exp(kβz0 )
(14)
As seen from (14), the derivative of c(z) depends on k, i.e., k is the variable to control the drop speed. Moreover, c(z) is also associated with variables z and βz0 . In some way, z0 may be regard as an average of z. So we choose parameter z as a self-correcting parameters to noise estimation z0 as follows: K (15) k= z0 (c)
(d)
Fig. 2 An example analysis for SRAD ((a) Original live ultrasound image; (b) Filtered images by SRAD, Δt = 0.05, 500 iterations; (c) The gray value in the profile; (d) The instantaneous coefficient of variation z of the section plane and the reference point z0 .)
2 2.1
The new speckle reducing anisotropic diffusion The new speckle reducing anisotropic diffusion
To make up the shortcoming of SRAD as analyzed before, an ideal diffuse coefficient should act so as to diffuse homogeneous regions with speeds changing stolidly but the transition regions with speeds changing sensitively. This can be realized by a sigmoid function. So a new speckle reducing anisotropic diffusion (NSRAD) can be written as follows: ⎧ ∂u ⎪ ⎪ ⎪ ⎨ ∂t = div[c(q)u(x, y; t)] (10) u(t = 0) = u0 = g ⎪ ⎪ ∂u ⎪ ⎩ = 0, u ∈ ∂Ω n ∂n with a new diffuse coefficient c(q) = 1 −
1 1 + exp(−k(q 2 − βq02 ))
(11)
where k > 0, it is a tunable parameter to control the dropping speed of the diffusion coefficient. In (11), β is the rate of the homogeneous regions to the noise estimation q02 (t), and we name it as the homogeneous region control coefficient. Note that the homogeneous region in an image is not unanimous to different observers, thus the homogeneous region scale is also different to different observers, which can be achieved by the homogeneous region control coefficient β. Similarly, we normalize c(q) by dividing its maximum
where K is a suitable constant. The new diffuse coefficient is plotted in Fig. 1. Here, we set z0 = 0.002 as before in SRAD, and K = 40, β = 10. The new diffuse coefficient is a sigmoid function and there are three parts in it, which correspond to three regions. The first part is z < βz0 , which corresponds to the homogeneous region. The second part is z ≈ βz0 , which corresponds to the transition region. The last part is z > βz0 , which corresponds to the edge region. In the first part (z < βz0 ), i.e., in the homogeneous region, the new diffuse coefficient c(z) ≈ 1. This means that the diffuse is almost uniform in the homogeneous region. Moreover, the diffuse coefficient c(z) changes slowly in this part. It means that c(z) varies small when z has a small deviation, which will not magnify the deviation of z in this region. So it will not affect block in the homogeneous region. In the second part (z ≈ βz0 ), i.e., in the transition region, the diffuse coefficient decreases rapidly. It means that c(z) varies greatly when z has a small deviation in this part, which will magnify the deviation of z in this region. So it will enhance the resolution power in the transition regions and improves the sharp of weak edges and details. In the third one (z > βz0 ), i.e., in the edge region, the diffuse coefficient c(z) tends to 0 fast, which will strength the sharp and the strong of the edges. 2.3
Analysis of the scale control coefficient β
The homogeneous region scale control coefficient β plays a crucial role in NSRAD. If we ignore the impact of k to β, the diffuse is almost uniform under scale β. The block will not occur in the homogeneous region for there is almost a uniform diffuse speed under scale β. It is well know that one homogeneous region under larger scale (larger β) will be divided into one new homogeneous region and one transition region under smaller scale (smaller β). So it is useful to preserve more details with a smaller β. On the contrary, one homogeneous region and one transition region under smaller scale (smaller β) will combine into one new homogeneous region under larger scale (larger β), this will blur details with a larger β. So, a different result will be gotten by a different scale β, that is, the diffusing can be adjusted
LI Can-Fei, et al.: A New Speckle Reducing Anisotropic Diffusion for · · ·
No. 7
by β, which is useful to fit different application. This will be further explainded by the function image of the diffuse coefficient. The diffuse coefficient function with different values of β (β = 1.0 and β = 1.2 , respectively) are plotted in Fig. 3 (b). From Fig. 3 (b), the homogeneous region under scale β = 1.2 can be seen as a new homogeneous region and a transition region under scale β = 1.0. This means that the different diffusion result will be got in the region if use different scales of β. Moreover, from Fig. 3 (b), the curve of the diffusion coefficient c(z) is steeper in the transition region when β is larger, so the diffusion coefficient c(z) is more sensitive to z in this one, which will improve the sharpness of the weak edges and enhance the details. On the contrary, the curve of the diffusion coefficient c(z) is more gently in the transition region when β is smaller, so the diffusion coefficient c(z) is more stolid to z in this one, which will blur the weak edges and details. All of these will affect the diffusion result. 2.4
(a)
n qi,j
Fig. 3 Analysis of the parameters in the diffuse coefficient for NSRAD ((a) The function image of c(z) in NSRAD with different values of K: solid line for K = 10 and dotted line for K = 40, respectively, in which β = 1.0, and z0 = 0.002; (b) The function image of c(z) in NSRAD with different values of β: solid line for β = 1.0 and dotted line for β = 1.2, respectively, where K = 40, and z0 = 0.002.)
Discretization
The partial differential (10) is solved numerically using a Jacobi iterative method[14] . Choosing a time-step Δt and a grid-size h in both x and y directions, discretize the time and space coordinates as t = nt, n = 0, 1, · · · , x = ih, y = jh, i = 0, 1, · · · , M − 1, j = 0, 1, · · · , N − 1, where M h × N h is the area of the image domain. Let un i,j = u(ih, jh; nt); an numerical approximation to (10) is given by 1 tdn i,j 4
2 2 n 2 1 ∇ ui,j |∇un 1 i,j | − 2 2 un 4 un i,j i,j =
2 ∇2 un 1 i,j
1+ 4 un i,j
(16)
(18)
and q02 = 0.5 ∗ (C ∗ median(M AD(∇ log(u))))2 ∇ R un i,j =
(b)
n un+1 i,j = ui,j +
1 n n n n n [ci+1,j (un i+1,j − ui,j ) + ci,j (ui−1,j − ui,j ) + h2 n n n n n cn (17) i,j+1 (ui,j+1 − ui,j ) + ci,j (ui,j−1 − ui,j )]
with
Analysis of the rate of the parameter
In the NSRAD, for k = K/z0 , to explain the influence of k is equal to explain the influence of K. The diffuse function images with different values of K (K = 10 and K = 40, respectively) are plotted in Fig. 3 (a), in which β = 1.0 and z0 = 0.002. As seen from Fig. 3, parameter k not only affects the rate of decline function but also affects the actual role of the homogeneous region scale control coefficient β. When K is smaller, it will reduce the value of the actual role of the scale control coefficient β. Moreover, the diffuse function drops more slowerly, which is disadvantageous to the resolution over the transition region. On the contrary, when K is larger, it will enlarge the value of the actual role of the scale control coefficient β. Moreover, the diffuse function drops faster, which is advantageous to the resolution over the transition region.
2.5
dn i,j =
415
∇ L un i,j =
∇ 2 un i,j =
n n un un i+1,j − ui,j i,j+1 − ui,j − h h
n n n un i,j − ui−1,j ui,j − ui,j−1 , h h
(19)
(20)
n n n n un i+1,j + ui−1,j + ui,j+1 + ui,j−1 − 4ui,j h
(21)
(22)
Moreover, to avoid distortion at the image boundaries, symmetric boundary condition is required, that is, the image intensity function has equal values at both sides of the boundary.
3
Experiments
To evaluate the performance of NSRAD, simulation and in vivo studies were conducted. The simulation studies may give quantitative performance analysis, and the in vivo study may demonstrate the feasibility and usefulness of NSRAD in real applications. In all the experiments, c(q) were computed by (3) in SRAD, and by (11) in NSRAD, respectively. In order to eliminate the patchy of the gray further (when needed), in experiments, a Gaussian smoothing was made: q (23) q = u 3.1
Results from simulation study
We ran the SRAD and NSRAD on Lena images corrupted by a multiplicative noise; see Fig. 4 (b), and Fig. 4 (a) is the reference image. Here, we use (23) to eliminate the patchy of the gray. To assure stability in the diffusion process, we took a time step t = 0.1 and ran 150 iterations, and the results for SRAD and NSRAD with K = 2 and β = 1 are shown in Fig. 4 (c) and Fig. 4 (d), respectively. As seem from Fig. 4, NSRAD gives a better diffusion but SRAD makes very strong block affect. At the same time, NSRAD preserves the sharpness of the edges. We repeated this experiment for greater values of Δt (see Fig. 5), such as Δt = 0.5 and Δt = 1, and Δt = 1 only for NSRAD. The result that seems not to match the statement in [14], and the fact that the diffusion coefficients observed with this filter may at some point be greater than could explain this effect[16] . But NSRAD is stable, even for Δt = 1.
416
ACTA AUTOMATICA SINICA
(a)
(b)
(c)
(d)
Vol. 38
different factors: loss of correlation, luminance distortion, and contrast distortion. The structure similarity (SSIM)[22] evaluates the overall processing quality. We show the MSE (Fig. 6 (a)), the Q (Fig. 6 (b)), the SNR (Fig. 6 (c)), and the SSIM (Fig. 6 (d)) evolutions with time for Δt = 0.1. The index MSE in NSRAD reduced as compared to SRAD, but index SNR increased which shows that the method NSRAD suppresses speckle more effectively. Moreover, note that NSRAD shows better results in other indexes than SRAD, such as SSIM and Q. These advantages even become more prominent with the increasing of iterations, which shows that the filtered image by NSRAD has less distortion and well structure similarity.
(a) MSE
(b) Q
(c) SNR
(d) SSIM
Fig. 4 Simulation experiment ((a) Original image; (b) Noisy image; (c) Filtered by SRAD: ∇t = 0.1, 150 iterations; (d) Filtered by NSRAD: K = 2, β = 1, ∇t = 0.1, 150 iterations.)
Fig. 6
3.2 (a)
(b)
(c)
Fig. 5 Simulation experiment with different ∇t values, 150 iterations ((a) Filtered by SRAD: ∇t = 0.5; (b) Filtered by NSRAD: ∇t = 0.5, K = 2, β = 1; (c) Filtered by NSRAD: ∇t = 1.0, K = 2, β = 1.)
To quantify the achieved performance improvements, four quality index metrics were computed on the filtered image and reference image. Mean squared error (MSE) measures the quality change between the original and processed images. Signal-to-noise ratio (SNR)[20] measures image fidelity. The mathematically defined universal quality index Q[21] measures any distortion as a combination of three
Index for 150 iterations (solid line) NSRAD and (dashed line) SRAD
Results from in vivo study
The SRAD and NSRAD were further tested with two real clinical images. Here, both the time steps were set to Δt = 0.05 and the same 500 iterations were run. Fig. 7 is a live ultrasound image filtering, the original image is Fig. 2 (a). Fig. 7 (a) and Fig. 7 (b) are the filtered image obtained from NSRAD with K = 400, β = 1.0 and K = 500, β = 0.8, respectively. The filtered image obtained from SRAD is shown in Fig. 2 (b). In Fig. 2 (b), SRAD makes very strong block effect, as labeled with arrows. But the proposed NSRAD suppresses speckle noises more effectively and the fine characters (labeled with circles) are better preserved compared with the traditional SRAD. The undesired block effects (labeled with arrows) in the SRAD can not be seen in NSRAD filtered image. Moreover, the fine characters (labeled with circles) in Fig. 7 (b) are increased with a larger K and a smaller β than the ones in Fig. 7 (a) in the proposed NSRAD filtered image. As shown in Fig. 8 (a), the second experiment was conducted by using a prostate ultrsound image, which is different from the previous data in that the objects appear to be irregular and more fine structures are included due to the inherent anatomical composition of soft tissues. Therefore, the speckle suppression of prostate image remains a challenge to many traditional methods. However, from our filtering results in Fig. 8 (c) and Fig. 8 (d), it is obvious that the majority of fine structures have been well preserved as
LI Can-Fei, et al.: A New Speckle Reducing Anisotropic Diffusion for · · ·
No. 7
marked with circles, and the fine characters (labeled with circles) in Fig. 8 (d) are increased with a larger K and a smaller β than the ones in Fig. 8 (c). In addition, the frequently seen block effects of SRAD as labeled by arrows in Fig. 8 (a) and Fig. 8 (b) are totally eliminated in our results.
417
2 Lee J S. Speckle analysis and smoothing of synthetic aperture radar images. Computer Graphics and Image Processing, 1981, 17(1): 24−32 3 Lee J S. Digital image enhancement and noise filtering by use of local statistics. IEEE Transactions on Pattern Analysis and Machine Intelligence, 1980, 2(2): 165−168 4 Lee J S. Refined filtering of image noise using local statistics. Computer Graphics and Image Processing, 1981, 15(4): 380−389
(a) K = 400, β = 1.0
Fig. 7
(b) K = 500, β = 0.8
The live ultrasound experiment by NSRAD with ∇t = 0.05, 500 iterations
5 Frost V S, Stiles J A, Shanmuggam K S, Holtzman J C. A model for radar images and its application to adaptive digital filtering of multiplicative noise. IEEE Transactions on Pattern Analysis and Machine Intelligence, 1982, 4(2): 157−165 6 Kuan D T, Sawchuk A A, Strand T C, Chavel P. Adaptive restoration of images with speckle. IEEE Transactions on Acoustics Speech and Signal Processing, 1987, 35(3): 373− 383 7 Kuan D T, Sawchuk A A, Strand T C, Chavel P. Adaptive noise smoothing filter for images with signal dependent noise. IEEE Transactions on Pattern Analysis and Machine Intelligence, 1985, 7(2): 165−177 8 Wang Yi, Niu Rui-Qing, Yu Xin, Shen Huan-Feng. Time dependent robust anisotropic diffusion processes. Acta Automatica Sinica, 2009, 35(9): 1253−1256 (in Chinese)
(a) Original ultrasound image
(b) Filtered images by SRAD
9 Jin J S, Wang Y, Hiller J. An adaptive nonlinear diffusion algorithm for filtering medical images. IEEE Transactions on Information Technology in Biomedicine, 2000, 4(4): 298−305 10 Meng Xiang-Lin, Wang Zheng-Zhi. Image diffusion based on visual masking effect. Acta Automatica Sinica, 2011, 37(1): 21−27 (in Chinese) 11 Ying Shi-Hui, Peng Ji-Gen, Zheng Kai-Jie, Qiao Hong. Lie group method for data set registration problem with anisotropic scale deformation. Acta Automatica Sinica, 2009, 35(7): 867−874 (in Chinese)
(c) Filtered images by NSRAD: K = 400, β = 0.9
Figure. 8
4
(d) Filtered images by NSRAD: K = 500, β = 0.7
A pancreas ultrasound experiment by SRAD and NSRAD with ∇t = 0.05, 500 iterations
Conclusion
A new anisotropic diffusion model with a sigmoid function based diffusion coefficient definition for ultrasonic speckle reduction has been investigated. The new method is able to smooth homogeneous region while efficiently avoiding block effects. Moreover, the new method strengthens the sensitivity of the transition region to preserver more details and weak edges. Experiments in simulation data and in vivo ultrasound show its better capacity in diffusing homogeneous region and preserving image edges and details than SRAD. References 1 Lamont D, Parker L, White M, Unwin N, Bennett S M A, Cohen M, Dickinson H O, Adamson A, Alberti K G M M, Craft A W. Risk of cardiovascular disease measured by carotid intima-media thickness at age 49-51: life course study. British Medical Journal, 2000, 320(7230): 273−278
12 Black M J, Sapiro G, Marimont D H, Heeger D. Robust anisotropic diffusion. IEEE Transactions on Image Processing, 1998, 7(3): 421−432 13 Rerona P, Malik J. Scale-space and edge detection using anisotropic diffusion. IEEE Transactions on Pattern Analysis and Machine Intelligence, 1990, 12(7): 629−639 14 Yu Y, Acton S T. Speckle reducing anisotropic diffusion. IEEE Transactions on Image Processing, 2002, 11(11): 1260−1270 15 Abd-Elmoniem K Z, Youssef A B, Kadah Y. Real-time speckle reduction and coherence enhancement in ultrasound imaging via nonlinear anisotropic diffusion. IEEE Transactions on Biomedical Engineering, 49(9): 997−1014 16 Aja-Fernandez S, Alberola-Lopez C. On the estimation of the coefficient of variation for anisotropic diffusion speckle filtering. IEEE Transactions on Image Processing, 2006, 15(9): 2694−2701 17 Krissian K, Westin C F, Kikinis R, Vosburgh K G. Oriented speckle reducing anisotropic diffusion. IEEE Transactions on Image Processing, 2007, 16(5): 1412−1424
418
ACTA AUTOMATICA SINICA
18 Jain A K. Fundamental of Digital Image Processing. New Jersey: Prentice Hall, 1989 19 Aja-Fernandez S, Vegas-Sanchez-Ferrero G, MartinFernandez M, Alberola-Lopez C. Automatic noise estimation in images using local statistics: additive and multiplicative cases. Image and Vision Computing, 2009, 27(6): 756−770 20 Sakrison D. On the role of observer and a distortion measure in image transmission. IEEE Transactions on Communications, 1977, 25(11): 1251−1267 21 Wang Z, Bovik A. A universal image quality index. IEEE Signal Processing Letters, 2002, 9(3): 81−84 22 Wang Z, Bovik A C, Sheikh H R, Simoncelli E P. Image quality assessment: from error visibility to structural similarity. IEEE Transactions on Image Processing, 2004, 13(4): 600−612
LI Can-Fei Ph. D. candidate at the College of Electrical and Information Engineering, Hunan University. She received her B. S. and M. E. degrees from the College of Electrical and Information Engineering, Hunan University in 2001 and 2005, respectively. Her research interest covers image recognition, computer vision, and medical image processing. Corresponding author of this paper. E-mail: [email protected]
Vol. 38 WANG Yao-Nan Professor at the College of Electrical and Information Engineering, Hunan University. He received his Ph. D. degree from Hunan University. He was a postdoctor at National University of Defense Technology and an alexander von Humboldt Stiftung. His research interest covers intelligent control, intelligent image processing, and intelligent robotics. E-mail: [email protected] XIAO Chang-Yan Received his B. S. and master degrees from National University of Defense Technology, China in 1994 and 1997, and his Ph. D. degree from Shanghai Jiao Tong University, China in 2004. From 2008 to 2009, he worked as a postdoctor in Leiden University, the Netherlands. Now he is an associate professor at the College of Electrical and Information Engineering, Hunan University, China. E-mail: [email protected] LU Xiao Ph. D. candidate at the College of Electrical and Information Engineering, Hunan University. She received her B. S. degree from the College of Electrical and Information Engineering, Hunan University in 2007, and her master degree from the College of Automation, Southeast University in 2010, respectively. Her research interest covers computer vision, and robot navigation. E-mail: xlu [email protected]