- Email: [email protected]
SRBF: Speckle Reducing Bilateral Filtering
Ultrasound in Med. & Biol., Vol. 36, No. 8, pp. 1353–1363, 2010 Copyright 2010 World Federation for Ultrasound in Medicine & Biology Printed in the USA. All rights reserved 0301-5629/$ - see front matter
doi:10.1016/j.ultrasmedbio.2010.05.007
d
Original Contribution SRBF: SPECKLE REDUCING BILATERAL FILTERING SIMONE BALOCCO,* CARLO GATTA,* ORIOL PUJOL,*y JOSEPA MAURI,z and PETIA RADEVA*y * Computer Vision Center, Bellaterra, Spain; y Department Matema`tica Aplicada i Ana`lisi, Universitat de Barcelona, Barcelona, Spain; and z Unitat d’hemodina`mica cardı´aca hospital universitari Germans Trias i Pujol Badalona, Spain (Received 8 October 2009; revised 5 May 2010; in final form 8 May 2010)
Abstract—Speckle noise negatively affects medical ultrasound image shape interpretation and boundary detection. Speckle removal filters are widely used to selectively remove speckle noise without destroying important image features to enhance object boundaries. In this article, a fully automatic bilateral filter tailored to ultrasound images is proposed. The edge preservation property is obtained by embedding noise statistics in the filter framework. Consequently, the filter is able to tackle the multiplicative behavior modulating the smoothing strength with respect to local statistics. The in silico experiments clearly showed that the speckle reducing bilateral filter (SRBF) has superior performances to most of the state of the art filtering methods. The filter is tested on 50 in vivo US images and its influence on a segmentation task is quantified. The results using SRBF filtered data sets show a superior performance to using oriented anisotropic diffusion filtered images. This improvement is due to the adaptive support of SRBF and the embedded noise statistics, yielding a more homogeneous smoothing. SRBF results in a fully automatic, fast and flexible algorithm potentially suitable in wide ranges of speckle noise sizes, for different medical applications (IVUS, B-mode, 3-D matrix array US). (E-mail: [email protected]) 2010 World Federation for Ultrasound in Medicine & Biology. Key Words: Bilateral filter, Ultrasound speckle reduction.
INTRODUCTION
efficient speckle removal filters become an unavoidable need. In literature, de-speckling has been tackled by means of different techniques. Speckle reducing filters based on anisotropic diffusion algorithms were introduced by Perona and Malik (Perona and Malik 1990) and later improved by several authors (Yu and Acton 2002; Black et al. 1998). This filter removes the noise by computing a local weighted average of the central pixel intensity with the ones of its neighbors. The iterative process achieves a balance between averaging (in homogeneous regions) and the identity filter (where edges exist) according to a coefficient proportional to the directional gradient. Yu (Yu and Acton 2002) successively improved Perona’s method by proposing a partial differential equation approach called speckle reducing anisotropic diffusion (SRAD) with the main purpose to enhance edges. Successively, Krissian (Krissian et al. 2007) proposed an extension called oriented speckle reduction anisotropic diffusion (OSRAD) that allows different levels of filtering along the image contours and their principal curvature directions. In 1998, Tomasi (Tomasi and Manduchi 1998) introduced the bilateral filter (BF) framework in which each output pixel value in the image is a Gaussian weighted average of its
Speckle, a form of multiplicative noise, affects imaging applications such as medical ultrasound (US). Speckle is the primary factor that limits the contrast in diagnostic ultrasound imaging, thereby reducing the effective application of image processing and analysis algorithms (i.e., edge detection, segmentation) and two-dimensional (2-D) and three-dimensional (3-D) volume rendering (Aysal and Barner 2007). The effectiveness of a segmentation can be improved when the noise is removed without affecting important image features. Speckle removal filters are often used as a preprocessing step for region-based detection, segmentation and classification algorithms (Munteanu et al. 2008; Li and Acton 2007; Tsai and Watanabe 1998). Their goal is to selectively remove the noise without destroying important image features and enhance object boundaries. Additionally, due to the recent diffusion of 3-D devices and real-time US applications in medical environments, computationally
Address correspondence to: Simone Balocco, Computer Vision Center, Edifici O, Campus UAB, 08193 Bellaterra (Cerdanyola), Barcelona, Spain. E-mail: [email protected] 1353
1354
Ultrasound in Medicine and Biology
neighbors in both space and intensity range. Comaniciu (Comaniciu and Meer 2002) proposed the mean shift approach based on a statistical local modes analysis of the image distribution in the joint spatial-range domain. The relationship between anisotropic diffusion, adaptive smoothing, bilateral filtering and mean shift procedure was then established by Barash and Comaniciu (Barash and Comaniciu 2004) in 2004. More recently, adaptive filters based on local noise statistics were proposed by Guo and Thakur (Guo et al. 2009; Thakur and Anand 2007). Dantas (Dantas and Costa 2007) introduced a filter based on a set of modified Gabor filters. However, most filters are developed independently of the image nature and its noise statistical model. In contrast, Rayleigh noise statistics were embedded in a basic filter framework by Aysal (Aysal and Barner 2007). In this article, a fully automatic speckle reducing bilateral filter (SRBF) tailored to US Images is proposed. The edge preserving feature for US images is obtained by embedding noise statistics in the filter framework. Consequently, the filter is able to tackle the noise multiplicative behavior modulating selectively the smoothing strength with respect to local statistics. Additionally, the filter support is automatically adapted to the speckle size. To obtain a strong reduction of speckle noise, the SRBF is iterated as suggested by (Tomasi and Manduchi 1998). In the Materials and Methods section, the bilateral filter framework tailored to US speckle statistics is introduced. The accuracy and robustness of our approach is compared with other speckle reducing filters. In the Results section, the filter is then tested on in vivo US images and its influence on a segmentation task is quantified. The final section is devoted to the discussion and conclusion. MATERIALS AND METHODS In this section, a brief description of the classical BF is given, followed by the characterization of the speckle noise. Then the BF framework is adapted to the a priori knowledge on speckle noise statistics and estimated speckle size. Bilateral filter framework The BF has been introduced as a technique for edgepreserving image smoothing (Tomasi and Manduchi 1998). It has been successfully used in different image processing and computer graphics applications (Durand et al. 2002; Choudhury and Tumblin 2005). The general BF functional can be expressed as: hðpÞ5G21 ðpÞeUðpÞ f ðxÞcðx; pÞsðf ðxÞ; f ðpÞÞdx
(1)
with the normalization factor: GðpÞ5eUðpÞ cðx; pÞsðf ðxÞ; f ðpÞÞdx
(2)
Volume 36, Number 8, 2010
where f is the input image, h is the output image, U(p) is the spatial neighborhood of the coordinate of a generic pixel p in the image and x is the integration variable representing pixels coordinates in U. The classical BF framework (Tomasi and Manduchi 1998) defines both c and s functions as unbiased isotropic Gaussian functions. ! k p2x k2 (3) cðx; pÞ5exp 2 2s2c
ðf ðpÞ2f ðxÞÞ2 sðf ðxÞ; f ðpÞÞ5exp 2 2s2s
! (4)
where sc represents the standard deviation of the Gaussian on the spatial support and ss is the standard deviation in the range domain. Function c [eqn (3)] spatially weights the Euclidean distance between c and x while s [eqn (4)] operates on the intensity domain. Speckle noise statistics Speckle appears in medical images when the characteristic size of the scatterers is small compared with the wavelength. Under the hypothesis of fully developed speckle, (Jensen 1996), biologic tissue can be modelled by a network of identical discrete scatterers, randomly distributed in a homogeneous media (Chivers 1977), thus the radio-frequency echo signal can be described as a complex Gaussian probability function with zero mean (Wagner et al. 1983). Envelope detection removes the phase component, generating a speckle distribution with statistics described by a Rayleigh probability function (PRL): ! 2 f ðpÞ 2f ðpÞ PRL ðf ðpÞ; aÞ5 2 exp (5) 2a2 a where f(p) is the image pixel intensity at position p and a is the shape parameterpofffiffi PRL related to the mean of the distribution by m5a p2 . It is worth noting that some commercial ultrasound scanners only provide log-compressed data whose noise statistic is not Rayleigh distributed. Ultrasound raw signal can be recovered from DICOM images using nonlinear compression parametric functions (Seabra and Sanches 2008). Given a spatial neighborhood of pixels with uniform intensity of f(p), a can be computed using the maximum likelihood estimator (MLE) (Steven 1993) defined for Rayleigh functions as: sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 X b5 a f ðqÞ2 (6) 2N q˛U
Speckle reducing bilateral filtering d S. BALOCCO et al.
Design of the s function for ultrasound noise statistics. In the case of images affected by additive zero-mean Gaussian noise, the value of ss is a constant and can be estimated directly from the image. However, the zero-mean assumption on the noise statistics, typical of the classical bilateral filter and anisotropic filter frameworks, is no longer valid for US images, where the noise is multiplicative and modelled by Rayleigh distributions. To extend the BF framework to US envelope images, the noise statistic is embedded in the filter by designing s as the probability that f(p) belongs to the Rayleigh distribution that generated f(x). Given a local neighborhood around p characterized by a Rayleigh distribution RLp(a), the probability of the central pixel intensity f(p) which belongs to RLx(a) can be expressed, combining eqns (5) and (6), as follows: ! f ðpÞ f ðpÞ2 P½f ðpÞ˛RLx ðaÞ5 2 exp 2 (7) b 2x bx 2a a b x is approximated by the pixel value f(x) In this article, a and s(f(x), f(p)) is set equal to P[f(p) ˛ RLx(a)]. Since the BF computes the s function on all pixels in the spatial neighborhood, the errors of the Rayleigh estimation are averaged converging to the local MLE. Therefore, when the central pixel in a neighborhood x and its neighbors belong to a similar Rayleigh distribution, function s yields high values and promotes smoothing. On the contrary, s has a small value when a contour is present, thus preserving the edges. Design of speckle size adaptive c function. The value of sc depends on the image resolution and must be set as a trade-off between computation efficiency and accuracy. An oversized support may promote undesired edge smoothing and increase computational time while a too small sc may hinder the smoothing capability of the BF and require more iterations to converge to the optimal solution. In this article, sc is set equal to the size of the speckle noise measured directly on the image. As proposed by Smith and Wagner (Smith and Wagner 1984), the speckle size can be obtained as the duration of the pulse waveform, which can be estimated by computing the full width at half maximum of the noisy autocorrelation coefficient. Given the autocorrelation matrix M is defined as follows: Mði; jÞ5
M21 X N21 X
f ðm; nÞ,f ðm1i; n1jÞ
(8)
m51 n51
IMis defined as the binary image containing the values of M superior to 12max{M}. The size of the adaptive support c corresponds to the maximum spatial width of the positive values of IM.
1355
RESULTS In silico experiments A set of in silico experiments was designed to evaluate the performance of the proposed approach. The edge-preserving performance of our method has been tested by evaluating the sum of square differences (SSD) criterion between the ground truth and the denoised images, which provides local comparison of the filtering results with the ground truth data. Experimental image set A set of US synthetic images, consisting of two regions, of low (R2) and high (R1) scatterer intensity, respectively, was generated using the convolutional approach (Yu and Acton 2002; Bamber and Dickinson 1980) detailed in the appendix. The echogenicity model amplitude in R1 (I0R1) is unitary while in R2 (I0R2) can be tuned using a parameter g 5 I0R1/I0R2 representing the brightness contrast between the two areas. Since the noise is multiplicative, g indirectly indicates the signal to noise ratio (SNR): the higher the contrast between the two areas, the higher the SNR. In the whole set of images, a 5 MHz central frequency probe is simulated and an anisotropic noise distribution is obtained using a point spread function (PSF) with parameters sx 5 1$k and sy 5 1.5$k, where k will be used as noise scaling factor. Figure 1 illustrates a 128 3 128 pixel image of the set with echogenicity model amplitudes I0R1 5 1 and I0R2 5 0.2. The histogram computed in two uniform regions of the image exhibits two Rayleigh distributions (Fig. 1b), with b s R2 50:272 and mR2 5 0.154 (in the central region) and b s R1 51:55 and mR1 5 1.55 (in the upper bright region), respectively. The ratio between the Rayleigh parameter estimators in the bright (R1) and the dark (R2) areas of the envelope image (Fig. 1a) are equal to the ratio of I0R1 over I0R2 since speckle noise is purely multiplicative. SRBF perfomance In this experiment, SRBF is applied to synthetic images generated with increasing intensities of I0R2 (by varying g while keeping k 5 1) to progressively reduce the contrast between the low and high brightness areas, thus increasing the overlapping of the two Rayleigh distributions. Therefore, it becomes more difficult to discriminate between real signal edges and variations due to the noise. Figure 2 shows the denoising result obtained when the noise intensity increases. As can be observed in Figure 2a in presence of high contrast (g 5 5.7), the algorithm converges in each uniform region to the MLE of the Rayleigh speckle distribution. At the same time, the
1356
Ultrasound in Medicine and Biology
Volume 36, Number 8, 2010
Fig. 1. (a) Synthetic ultrasound envelope-detected images B(x, y) (128 3 128 pixels size), simulated using a 5 MHz central frequency and sx 5 2 and sy 5 1.5 PSF parameters. (b) Histogram and Rayleigh curve fitting in BR1and BR2.
reconstructed image strongly correlates with the echogenicity image I0 showing that the edges have been accurately preserved. Figure 2c illustrates that the performance of the filter reduces proportionally with the decrease of the contrast as shown in Figure 2b (g 5 3.6).
SRBF accuracy and efficiency In this section, SRBF is compared with the state-ofthe-art speckle reducing techniques. Classic bilateral filter (CLASSIC BF) (Tomasi and Manduchi 1998), anisotropic diffusion (AD) (Black et al. 1998), Rayleigh-maximumlikelihood filter (RMLF) (Aysal and Barner 2007), adaptive filter based on second order statistics (AF) (Thakur and Anand 2007), modified Gabor filters (MGF) (Dantas and Costa 2007), Guo filter (GUO) (Guo et al. 2009), mean-shift (MEAN-SHIFT) (EDISON: Code for the Edge Detection and Image SegmentatiON system). (Barash and Comaniciu 2004), speckle reduction anisotropic diffusion (SRAD) (Yu and Acton 2002), oriented speckle reduction anisotropic diffusion (OSRAD) (Krissian et al. 2007) and SRBF are applied to a set of 50 synthetic images generated with the same contrast and speckle size (g 5 5 and k 5 1) to qualitatively (Fig. 3) and quantitatively (Table 1) evaluate the respective performances. Table 1 lists the parameters for each method, compares the features and reports the accuracy scores. As it can be observed in Figure 3a and b, the proposed method presents a superior edge preserving behavior, thus outperforming filtering method based on zero-mean noise assumptions (CLASSIC BF and AF). RMLF, AF, MGF, GUO and MEAN-SHIFT methods are unable to fully remove the speckle noise in the uniform area although they exhibit border preservation properties (see Fig. 3c, d, e, f and g). On the other hand, SRAD and OSRAD better smooth the noise close to the lower
edge of the image because of their edge-enhancing feature. Even though SRAD and OSRAD use a fixed spatial support, the latter achieves better performances since it emphasizes the smoothing behavior along the tangential direction to the object boundary (Fig. 3h and i). These results are in agreement with the average SSD similarity score, computed over the 50 synthetic images, reported in Table 1 (last column). In this experiment, the reconstruction accuracy of SRBF and OSRAD are similar, however, it is worth noting that SRBF is faster (about 20 times) and fully automatic, as summarized in the feature description listed in Table 1. SRBF robustness Even though in the previous experiment the denoising efficacy of OSRAD and SRBF has been found to be similar, the robustness of such filters is investigated when the size of the speckle changes. Both filters are applied to synthetic images generated with different speckle sizes (by varying k while keeping g 5 5.7) and the results are illustrated in Figure 4. As shown in Figure 4a, the performance of OSRAD is comparable to SRBF when the noise pattern is small and progressively decreases in presence of big speckle patterns. In this last case, as it can be seen in Figure 4c, OSRAD enlarges the edges. The same behavior can be observed at different contrasts (g). The performance decrease of OSRAD filter is probably due to its constant size neighbor support while SRBF automatically adapts its spatial support depending on the speckle pattern. In conclusion, SRBF exhibits higher robustness because of its fully automatic settings of relevant parameters, the speckle noise statistics and the speckle size. In vivo experiments Despeckle algorithms are generally used as preprocessing step for automatic segmentation of US images.
Speckle reducing bilateral filtering d S. BALOCCO et al.
1357
tested on intravascular ultrasound (IVUS) images because of their challenging properties: the anisotropy of the speckle patterns, the presence of artifacts and the variable echogeniticy of the tissues.
Fig. 2. Intensity profile and denoising result obtained with g 5 5.7 (a) and g 5 3.6 (b), respectively. Average sum of square differences (SSD) between echogenicity and denoised images (c) while g varies from 20 to 2.86. The plots (a) and (b) correspond to the dashed lines in (c), respectively.
IVUS in vivo images Fifty slices, chosen to represent different vascular structures (plaque and vessel shape, presence of stent, lumen area and diameter), were extracted from a database of 3000 in vivo coronary images. The acquisition was performed using an IVUS Galaxy II System with a catheter Atlantis SR Pro 40 MHz (Boston Scientific, Natick, MA, USA). For each slice, the RF signal is available. From all the 50 slices in polar coordinates, different sets of images are computed: the envelope of the RF signal (data set A), the OSRAD filtered envelope (data set B) and the SRBF filtered envelope (data set C). To avoid the discontinuity at 0/360 degrees of the IVUS image, we mirrored the leftmost (0 ) and rightmost (360 ) texture before filtering, copying 30 of both on the opposite side, thus obtaining a continuous signal at the boundaries. The catheter ring-down artifact located at the top of the polar image was suppressed by filling the region with a uniform colour computed as the mean intensity of the underlying stripe, 0.3 mm thick, adjacent to the catheter area. The same stripe, containing only blood, has been defined as uniform region required by the OSRAD filter. All the images were log compressed and converted to Cartesian coordinates prior to segmentation and visualization. A classic level-set active contour (snake) algorithm (Chan and Vese 2001), commonly used to segment US medical images (Michailovich and Tannenbaum 2007), is applied to the three data sets. The initialization shape and the level-set smoothness parameter are fixed and tuned to obtain the best performance on the envelope data set. The same settings are used when the snake is applied to the data sets B and C. The snake propagates from an initial circular contour until stopping at an energy-stable configuration dependent on the image gradients. The snake inner area is obtained from the zero levelset shape at convergence. Ground truth shapes, indicating the expected segmentation result, were manually delineated in each of the IVUS frames by two independent experts. The operator performing the analysis was blinded to the automatic and to the other manual segmentation results.
The performance of the sole filter cannot be validated on in vivo images since reference denoised data is obviously not available. Instead, the accuracy improvement of a segmentation algorithm is quantified when the images are denoised. For this purpose, the algorithm is compared against the best state-of-the-art filtering method (OSRAD) and
Segmentation comparison Figure 5 presents the manual delineated border and the results of the automatic segmentation for seven representative frames of the three data sets A, B and C . Frame 1 illustrates an accurate automatic segmentation obtained in all the three data sets. In frame 2, the snake crosses the vessel border between the stent wires (highly echogenic
1358
Ultrasound in Medicine and Biology
Volume 36, Number 8, 2010
Fig. 3. Normalized intensity profile of denoised images, superimposed on the ground truth shape I0. Pairwise filtering comparison of (a) CLASSIC BF , (b) AD , (c) RMLF, (d) AF, (e) MGF, (f) GUO, (g) MEAN-SHIFT, (h) SRAD and (i) OSRAD vs. SRBF algorithm.
spots) and spreads in the surrounding tissue. On the contrary, the snake propagation stops correctly at the lumen frontier in data sets B and C, since both filters enhance the outline of the vessel membrane by smoothing the noise and by preserving the sharp vessel edges. Similarly, in frame 3, the denoising methods enhance the borders making the segmentation easier. However, in this case (frame 3) the segmentation of data set B is slightly less accurate since the filtering method fuses a small catheter-ringing artifact with the vessel border, creating an artificial boundary. The same behavior can be observed in frame 4 where, in data set B, the noise lying between the catheter artifact and the vessel is not removed but homogenized.
Frame 5 illustrates that the high amount of noise (data set A) may hamper the snake propagation on the envelope data, while an efficient smoothing promote an accurate and automatic segmentation (data sets B and C). However, it is worth noting that, in data set B (frame 5), the image filtered with OSRAD is locally smoothed but radial ripples are still presents since OSRAD has a fixed support. On the other hand, SRBF better homogenizes uniform areas since the spatial support controlled by the function c is wider and has been designed to automatically contain consecutive scatter peaks. For this reason probably, in frame 6, the snake applied to B is trapped in a local minimum and stops propagating after 560 iterations, while in C the snake keeps spreading inside the whole lumen.
Table 1. Features and performance comparison of speckle reducing filters Support
User interaction
Noised image I0 CLASSIC BF AD
– User defined N 4 neighbors
– Fully automatic Noise variance size (s)
RMLF
8 neighbors
MGF
Freq. filter
2 smoothing coefficients (k1,2) 1 coefficient of variation (t) Probe central frequency k0, Lateral and axial Gabor filter orientations sax, slat
AF
Adaptive W
a0,D, CT, Liny, Hiny, Lasm Hasm, Et as described in [13]
GUO MEAN-SHIFT SRAD OSRAD SRBF
24 neighbors User defined JV 4 neighbors 4 neighbors Adaptive
Fully automatic Noise variance (a) Choice of uniform region Choice of uniform region Fully automatic
Experimental parameters – 8 s 51.3 k1 5 0.5 k2 5 0.05 t 5 0.05 k0 54 sax 55 sax 5 10 W 5 2, a0 5 0.04 D 5 2, CT 5 3.2, Et 5 1 Liny 5 04, Lasm 5 0.05 Hiny 5 0.8, Hasm 5 0.2 None N 5 6, s 5 1 top stripe of IR2 top stripe of IR2 None
Computat. time (1283128 pixels)
SSD (mean 6 SD)
6 300
– 0.364 [s] 1.411 [s]
7104 6 180 2736 6 130 3692 6 139
1
0.028 [s]
3138 6 114
1
0.055 [s]
2384 6 190
Iterations n. –
1 1 10 1000 1000 6
47 [s] 3.117 [s] 3.312 [s] 1.57 [s] 6.305 [s] 0.364 [s]
2810 6 148 2622 6 129 2246 6 129 1913 6153 1635 6 168 1288 6 127
CLASSIC BF 5 classic bilateral filter; AD 5 anisotropic diffusion; RMLF 5 Rayleigh-maximum-likelihood filter; MGF 5 modified Gabor filter; AF 5 adaptive filter; GUO 5 Guo filter; MEAN-SHIFT 5 mean-shift; SRAD 5 speckle reduction anisotropic diffusion; OSRAD 5 oriented speckle reduction anisotropic diffusion; SRBF 5 speckle reducing bilateral filter.
Speckle reducing bilateral filtering d S. BALOCCO et al.
Fig. 4. Average sum of square differences (SSD) (a) between the echogenicity image I0 and the image filtered using OSRAD and SRBF when the speckle noise scaling factor k ranges from 0.5 to 3.5. Intensity profile of I0 (b-c), superimposed on the OSRAD and SRBF filtering result when the noise scaling is equal to 1 and to 2.5, respectively.
Finally in frame 7, the snake spans outside the vessel for data sets A and B, but not in C. In this frame, the OSRAD filter enhances the shadow area on the left of the catheter, and creates an undesired gradient attracting the snake in B, while SRBF homogenizes the shadow region with the rest of the lumen in C. In fact, compared with OSRAD, which is a zero mean filter, the SRBF noise model better fits the speckle distribution at low image intensities and discriminate more accurately between noise and tissues. The last column of Figure 5 represents the average distance of the automatic segmentation to the two manual segmentations. In frame 1, the high slope of the curve indicates that the snake converges faster in filtered data sets (B and C) since the noise is accurately removed. In frame 4, a similar behavior can be observed, but the segmentation of B results as the slowest and the less accurate because of the wrinkled lumen region. In frames 2, 3
1359
and 5, the segmentation applied to the nonfiltered image diverges while the two filters show similar convergence speed. These examples illustrate the interest of using an accurate filtering method before the segmentation stage. In frame 6, the snake applied to data set B stops expanding earlier than data set C leading to incorrect vessel segmentation and to lumen area under-estimation. Finally, in frame 7, the SRBF filtered the snake applied to data set A and B diverges and only the snake applied to the data set C converges to the vessel wall. In all the cases, the segmentation on data set C is either the fastest or the only one converging to the vessel boundary, showing the superior performances to SRBF. Figures 6 and 7 show a pair wise error comparison for all the frames of the three data sets (A, B and C). In particular, Figure 6 illustrates the error between the automatic and ground-truth border segmentations while Figure 7 shows the error in the lumen area estimations. In Figure 6a most of the markers are located below the solid line of unitary slope, indicating that the average distance between the automatic and the manual segmentations is smaller in the filtered data sets (B, C) then in the nonfiltered (A). Figure 6b, shows the superior accuracy obtained by filtering the images with SRBF. The average distance over all the images of the data sets A, B and C with respect to the manual segmentations are 0.27 6 0.47 mm, 0.15 6 0.24 mm and 0.12 6 0.14 mm, respectively, confirming that the contours obtained from the data set C exhibit the lowest segmentation error. The average distance over all the images of the interobserver variability is 0.021 mm. To quantify the error in estimating the area for the three data sets the Jaccard coefficient (Jaccard 1912) jFXJj defined as the ratio J5jFWJj is used, where F is the set of pixels of the manual segmentation and J is the set of pixels resulting from automatic segmentation. The score J is one when the areas matches exactly and rapidly decreases when one of the two areas differs from the other. The performance of the automatic segmentation on the filtered data sets B and C, is superior to the non filtered images A (see Fig. 7a), since J is higher than the line with unitary slope for most of the images. Analyzing the scatter plot of Figure 7b and comparing directly the performance of the snake on the two denoised data sets, it can be noticed that the J-score is similar but some outliers, representing unsuccessful segmentations of data set B, are present. In more than 50% of the images , the J-score corresponding to OSRAD filter ranges from 0.12 to 0.33 while for SRBF J varies from 0.06 to 0.24 (Fig. 6b). Analyzing the scatter plot of Figure 7b in detail, it can be noticed that in about 50% of the images, OSRAD J-score ranges from 0.2 to 0.8 while for SRBF J varies from 0.6 to 0.9. The average J-score over all the images of the data sets A, B and C are 0.41 6 0.29, 0.76 6 0.18 and 0.83 6 0.09
1360
Ultrasound in Medicine and Biology
Volume 36, Number 8, 2010
Fig. 5. Snake manual (first column) and automatic (second, third and forth columns) segmentations of seven representative frames from the three data sets A, B and C . The last column illustrates the average distance of the automatic segmentations to the manual segmentations.
Speckle reducing bilateral filtering d S. BALOCCO et al.
1361
Fig. 6. Average distance of the automatic segmentation to the two manual segmentations. (a) comparison of data sets B and C vs. A ; (b) pair wise comparison of data sets B vs. C.
respectively, showing that the areas recovered on the data set C are the most similar to the manually segmented regions. ANOVA and t-test statistical analysis were performed for both segmentation distances and J-scores in data sets A, B and C. In all cases the scores of data sets B and C are significantly different (p values ,1$10–9 for both ANOVA and t-test) from data set A. Comparing data sets B and C, the ANOVA provides p values of 0.0032 and 0.0153, respectively, for the distances and J-scores; and the t-test provides p values of 5.5$10–4 and 0.0039, respectively, for the distances and J-scores. These tests confirm the statistical significance of the results. The J-score over all the images of the interobserver variability is 0.95. DISCUSSION AND CONCLUSIONS In this article, a bilateral filter tailored to US images is proposed. The edge preservation property is obtained by embedding noise statistics in the filter framework. The filter implicitly estimates the Rayleigh distribution parameters on the pixel neighborhood. Moreover, its spatial
support is automatically defined by estimating the size of the speckle noise directly on the image. A series of in silico experiments has been designed with the aim to compare the performances of the proposed approach against other denoising methods on synthetic images corrupted by speckle noise. As expected, SRBF outperforms the filtering methods based on zero mean noise assumptions. Compared with OSRAD, our approach is fully automatic, locally adaptive and the number of iterations needed to reach the optimal solution is smaller. Additionally, SRBF results in a more robust algorithm suitable for a wide ranges of speckle noise sizes. The proposed method has been tested on 50 in vivo IVUS images, by quantifying the influence of the denoising filter on the segmentation accuracy of a generic level-set algorithm. For this purpose, three data sets were created: the first composed of simple envelope images and the other two filtered using OSRAD and SRBF algorithms. The use of a more sophisticated segmentation algorithm, such as active contours with shape priors (Eduard et al. 2008), statistic snakes (Ayache et al. 1993) or including a leakage detection feature would improve the
Fig. 7. Area ratio J representing the likeness between the automatically segmented areas. (a) comparison of data sets B and C vs. A; (b) pair wise comparison of data sets B and C.
1362
Ultrasound in Medicine and Biology
segmentation result and succeed on all the envelope data sets, while a weakly constrained snake performs according to the quality of the original data. As shown in Figures 6 and 7, an accurate denoising method able to enhance the boundary and remove the noise combined with an ordinary snake reduces one third the average distance between the manual and the automatically segmented lumen. Considering that the average diameter of the coronary artery of the used image database is about 3 mm, the segmentation error reduces from about 0.3 mm to 0.1 mm representing 20% and 6% of the arterial radius, respectively. The experiment shows that the accuracy of the snake applied to filtered images is superior to that obtained in envelope data sets since both despeckle algorithms enhance the outline of the vessel membrane and preserve the sharp vessel edges. However, as expected, in some cases (frames 6 and 7) the snake applied to OSRAD filtered images diverges or stops too early, leading to incorrect segmentation results. Indeed, SRBF not only locally smoothes the images, but fills the gaps in uniform areas accurately given its adaptive spatial support. For this reason, the snake propagates faster as demonstrated by the low number of iterations required to converge to the correct vessel contour. The adaptive scaling of the filter, obtained through the estimation of the speckle size, is a feature that can be implemented only in denoising methods allowing a support size adjustment (see Table 1), thus, OSRAD filter cannot accommodate this feature. Moreover, in almost all the 50 tested images, the contours and the lumen areas obtained on the SRBF filtered data set were the closest to the manual delineated contours. Finally, it is worth noting that SRBF is fully automatic and does not require the definition of a user-provided uniform region to tune the intensity of the filtering, thus making it suitable for unsupervised filtering of large amounts of clinical data. In summary, SRBF is a fast and flexible algorithm potentially suitable in different medical applications (IVUS, B-mode, 3-D matrix array US). Future work will include studying a possible extension of the proposed filter to more complex noise models, by replacing the Rayleigh function with a Rician distribution, to extend the filter applicability to other medical modalities (e.g., MRI images). Acknowledgments—The works of S. Balocco and C. Gatta are supported by the Catalan Agency for Administration of University and Research (AGAUR) under a Beatriu de Pinos Fellowship. The authors would like to thank Dr. T. Aysal and Dr. K. Barner for providing the source code of the RMLF filter and Dr. K. Krissian for providing a Windows compiled version of the OSRAD filter and for the fruitful discussions about anisotropic diffusion filters.
REFERENCES Ayache N, Cohen I, Herlin I. Medical image tracking. Cambridge, MA: MIT Press; 1993. p. 285–302.
Volume 36, Number 8, 2010 Aysal TC, Barner KE. Rayleigh-maximum-likelihood filtering for speckle reduction of ultrasound images. IEEE Trans Med Imaging 2007;26:712–727. Bamber JC, Dickinson RJ. Ultrasonic b-scanning: A computer simulation. Phys Med Biol 1980;25:463–479. Barash D, Comaniciu D. A common framework for nonlinear diffusion, adaptive smoothing, bilateral filtering and mean shift. Image Vis Comput 2004;22:73–81. Black M, Sapiro G, Marimont D, Heeger D. Robust anisotropic diffusion. IEEE Trans Image Process 1998;7:421–432. Chan TF, Vese LA. Active contours without edges. IEEE Trans Image Process 2001;10:266–277. Chivers R. The scattering of ultrasound by human tissues, some theorical models. Ultrasound Med Biol 1977;3:1–13. Choudhury P, Tumblin J. The trilateral filter for high contrast images and meshes. International Conference on Computer Graphics and Interactive Techniques ACM SIGGRAPH 2005 Courses Los Angeles, CA: SESSION: Computational photography. Article No.: 5. 2005. Comaniciu D, Meer P. Mean shift: A robust approach toward feature space analysis. IEEE Trans PAMI 2002;24:603–619. Dantas RG, Costa ET. Ultrasound speckle reduction using modified Gabor filters. IEEE Trans Ultrason Ferroelectr Freq Control 2007;54:530–538. Durand F, Dorsey J. Fast bilateral filtering for the display of high-dynamicrange images. New York, NY: ACM Press; 2002. p. 257–266. Guo Y, Cheng HD, Tian J, Zhang Y. A novel approach to speckle reduction in ultrasound imaging. Ultrasound Med Biol 2009;35:628–640. Jaccard P. The distribution of the flora in the alpine zone. New Phytologist 1912;11:37–50. Jensen J. Estimation of blood velocities using ultrasound: A signal processing approach. Cambridge: Cambridge University Press; 1996. Kay SM. Fundamentals of statistical signal processing: Estimation theory. Upper Saddle River, NJ: Prentice-Hall, Inc; 1993. p. 595. Krissian K, Westin CF, Kikinis R, Vosburgh KG. Oriented speckle reducing anisotropic diffusion. IEEE Trans Image Process 2007;16: 1412–1424. Li B, Acton ST. Active contour external force using vector field convolution for image segmentation. IEEE Trans Image Process 2007;16: 2096–2106. Michailovich O, Tannenbaum A. Segmentation of medical ultrasound images using active contours. Proceedings of the IEEE International Conference on Image Processing (ICIP). Sep 2007;513–516. Munteanu C, Morales FC, Fernandez JG, Rosa A, Deniz LG. Enhancing obstetric and gynecology ultrasound images by adaptation of the speckle reducing anisotropic diffusion filter. Artif Intell Med 2008; 43:223–242. Perona P, Malik J. Scale space and edge detection using anisotropic diffusion. IEEE Trans Pattern Anal Machine Intell 1990;12:629–639. Seabra J, Sanches J. Modeling log-compressed ultrasound images for radio frequency signal recovery. Conf Proc IEEE Eng Med Biol Soc 2008;2008:426–429. Smith SF, Wagner RF. Ultrasound speckle size and lesion signal to noise ratio: Verification of theory. Ultrason Imaging 1984;6:174–180. Sojka E, Gaura J, Krumnikl M. Active contours without edges and with simple shape priors source. Lecture Notes In Computer Science; Vol. 5259 archive. Proceedings of the 10th International Conference on Advanced Concepts for Intelligent Vision Systems. Juan-les-Pins, France 2008;730–741. Thakur A, Anand RS. Speckle reduction in ultrasound medical images using adaptive filter based on second order statistics. J Med Eng Technol 2007;31:263–279. Tomasi C, Manduchi R. Bilateral filtering for gray and color images. In Proceedings of the Sixth International Conference on Computer Vision, Bombay, India, January 1998. Tsai DY, Watanabe S. A method for optimization of fuzzy reasoning by genetic algorithms and its application to discrimination of myocardial heart disease. 1998;. 3:1756–1761. Wagner RF, Smith SW, Sandrik JM, Lopez H. Statistics of speckle in ultrasound b-scans. IEEE Trans Sonics Ultrason 1983;30:156–163. Yu Y, Acton S. Speckle reducing anisotropic diffusion. IEEE Trans Image Process 2002;11:1260–1271.
Speckle reducing bilateral filtering d S. BALOCCO et al.
APPENDIX In this study, synthetic images corrupted by speckle noise have been generated using a convolutional approach described by (Yu and Acton 2002; Bamber and Dickinson 1980). The simulation assumes that the imaging system has a linear, space-invariant point spread function (PSF) and the transducer is linear. Let I0(x, y) be an echogenicity model (an image with different intensities corresponding to different tissues of the object being imaged) in which the variables x and y are the lateral and axial coordinates, respectively. First, subresolution variations in object impedance should be introduced by multiplying the echogenicity model by a Gaussian white noise with zero mean and unitary variance: Tðx; yÞ5I0 ðx; yÞ$Gðx; yÞ
(9)
T(x, y), accounts for acoustic impedance inhomogeneities in the object due to density and acoustic speed perturbation that generates the scattering. The V(x, y) ultrasonic radio-frequency (RF) echo data can then be obtained by a convolution: Vðx; yÞ5hðx; yÞ Tðx; yÞ
(10)
1363
where h(x, y) is the point spread function or impulse response of a hypothetical US imaging system. It is assumed (Yu and Acton 2002; Bamber and Dickinson 1980) that h(x, y) is separable, i.e., h(x, y) 5 h1(x) h2(y) where h1 is a Gaussian-weighted sinusoidal function determined by: h1 ðx; sx Þ5sinðk0 xÞ exp x2 =2s2x
(11)
where k0 5 2pf0/v, v is the speed of sound in tissue, f0 is the center frequency, and sx represents the pulse-width of transmitting ultrasonic wave. The spatial response for the transmitting and receiving aperture h2(x) is determined by: h2 ðyÞ5exp y2 =2s2y
(12)
where sy represents the beam-width of transmitting ultrasonic wave. The image of the envelope-detected amplitude, B(x, y) shown in Figure 1a, is given by: Bðx; yÞ5jVðx; yÞ1i Vbðx; yÞj
(13)
where Vbðx; yÞ is the Hilbert transform of V(x, y) and i is the imaginary unit.
doi:10.1016/j.ultrasmedbio.2010.05.007
d
Original Contribution SRBF: SPECKLE REDUCING BILATERAL FILTERING SIMONE BALOCCO,* CARLO GATTA,* ORIOL PUJOL,*y JOSEPA MAURI,z and PETIA RADEVA*y * Computer Vision Center, Bellaterra, Spain; y Department Matema`tica Aplicada i Ana`lisi, Universitat de Barcelona, Barcelona, Spain; and z Unitat d’hemodina`mica cardı´aca hospital universitari Germans Trias i Pujol Badalona, Spain (Received 8 October 2009; revised 5 May 2010; in final form 8 May 2010)
Abstract—Speckle noise negatively affects medical ultrasound image shape interpretation and boundary detection. Speckle removal filters are widely used to selectively remove speckle noise without destroying important image features to enhance object boundaries. In this article, a fully automatic bilateral filter tailored to ultrasound images is proposed. The edge preservation property is obtained by embedding noise statistics in the filter framework. Consequently, the filter is able to tackle the multiplicative behavior modulating the smoothing strength with respect to local statistics. The in silico experiments clearly showed that the speckle reducing bilateral filter (SRBF) has superior performances to most of the state of the art filtering methods. The filter is tested on 50 in vivo US images and its influence on a segmentation task is quantified. The results using SRBF filtered data sets show a superior performance to using oriented anisotropic diffusion filtered images. This improvement is due to the adaptive support of SRBF and the embedded noise statistics, yielding a more homogeneous smoothing. SRBF results in a fully automatic, fast and flexible algorithm potentially suitable in wide ranges of speckle noise sizes, for different medical applications (IVUS, B-mode, 3-D matrix array US). (E-mail: [email protected]) 2010 World Federation for Ultrasound in Medicine & Biology. Key Words: Bilateral filter, Ultrasound speckle reduction.
INTRODUCTION
efficient speckle removal filters become an unavoidable need. In literature, de-speckling has been tackled by means of different techniques. Speckle reducing filters based on anisotropic diffusion algorithms were introduced by Perona and Malik (Perona and Malik 1990) and later improved by several authors (Yu and Acton 2002; Black et al. 1998). This filter removes the noise by computing a local weighted average of the central pixel intensity with the ones of its neighbors. The iterative process achieves a balance between averaging (in homogeneous regions) and the identity filter (where edges exist) according to a coefficient proportional to the directional gradient. Yu (Yu and Acton 2002) successively improved Perona’s method by proposing a partial differential equation approach called speckle reducing anisotropic diffusion (SRAD) with the main purpose to enhance edges. Successively, Krissian (Krissian et al. 2007) proposed an extension called oriented speckle reduction anisotropic diffusion (OSRAD) that allows different levels of filtering along the image contours and their principal curvature directions. In 1998, Tomasi (Tomasi and Manduchi 1998) introduced the bilateral filter (BF) framework in which each output pixel value in the image is a Gaussian weighted average of its
Speckle, a form of multiplicative noise, affects imaging applications such as medical ultrasound (US). Speckle is the primary factor that limits the contrast in diagnostic ultrasound imaging, thereby reducing the effective application of image processing and analysis algorithms (i.e., edge detection, segmentation) and two-dimensional (2-D) and three-dimensional (3-D) volume rendering (Aysal and Barner 2007). The effectiveness of a segmentation can be improved when the noise is removed without affecting important image features. Speckle removal filters are often used as a preprocessing step for region-based detection, segmentation and classification algorithms (Munteanu et al. 2008; Li and Acton 2007; Tsai and Watanabe 1998). Their goal is to selectively remove the noise without destroying important image features and enhance object boundaries. Additionally, due to the recent diffusion of 3-D devices and real-time US applications in medical environments, computationally
Address correspondence to: Simone Balocco, Computer Vision Center, Edifici O, Campus UAB, 08193 Bellaterra (Cerdanyola), Barcelona, Spain. E-mail: [email protected] 1353
1354
Ultrasound in Medicine and Biology
neighbors in both space and intensity range. Comaniciu (Comaniciu and Meer 2002) proposed the mean shift approach based on a statistical local modes analysis of the image distribution in the joint spatial-range domain. The relationship between anisotropic diffusion, adaptive smoothing, bilateral filtering and mean shift procedure was then established by Barash and Comaniciu (Barash and Comaniciu 2004) in 2004. More recently, adaptive filters based on local noise statistics were proposed by Guo and Thakur (Guo et al. 2009; Thakur and Anand 2007). Dantas (Dantas and Costa 2007) introduced a filter based on a set of modified Gabor filters. However, most filters are developed independently of the image nature and its noise statistical model. In contrast, Rayleigh noise statistics were embedded in a basic filter framework by Aysal (Aysal and Barner 2007). In this article, a fully automatic speckle reducing bilateral filter (SRBF) tailored to US Images is proposed. The edge preserving feature for US images is obtained by embedding noise statistics in the filter framework. Consequently, the filter is able to tackle the noise multiplicative behavior modulating selectively the smoothing strength with respect to local statistics. Additionally, the filter support is automatically adapted to the speckle size. To obtain a strong reduction of speckle noise, the SRBF is iterated as suggested by (Tomasi and Manduchi 1998). In the Materials and Methods section, the bilateral filter framework tailored to US speckle statistics is introduced. The accuracy and robustness of our approach is compared with other speckle reducing filters. In the Results section, the filter is then tested on in vivo US images and its influence on a segmentation task is quantified. The final section is devoted to the discussion and conclusion. MATERIALS AND METHODS In this section, a brief description of the classical BF is given, followed by the characterization of the speckle noise. Then the BF framework is adapted to the a priori knowledge on speckle noise statistics and estimated speckle size. Bilateral filter framework The BF has been introduced as a technique for edgepreserving image smoothing (Tomasi and Manduchi 1998). It has been successfully used in different image processing and computer graphics applications (Durand et al. 2002; Choudhury and Tumblin 2005). The general BF functional can be expressed as: hðpÞ5G21 ðpÞeUðpÞ f ðxÞcðx; pÞsðf ðxÞ; f ðpÞÞdx
(1)
with the normalization factor: GðpÞ5eUðpÞ cðx; pÞsðf ðxÞ; f ðpÞÞdx
(2)
Volume 36, Number 8, 2010
where f is the input image, h is the output image, U(p) is the spatial neighborhood of the coordinate of a generic pixel p in the image and x is the integration variable representing pixels coordinates in U. The classical BF framework (Tomasi and Manduchi 1998) defines both c and s functions as unbiased isotropic Gaussian functions. ! k p2x k2 (3) cðx; pÞ5exp 2 2s2c
ðf ðpÞ2f ðxÞÞ2 sðf ðxÞ; f ðpÞÞ5exp 2 2s2s
! (4)
where sc represents the standard deviation of the Gaussian on the spatial support and ss is the standard deviation in the range domain. Function c [eqn (3)] spatially weights the Euclidean distance between c and x while s [eqn (4)] operates on the intensity domain. Speckle noise statistics Speckle appears in medical images when the characteristic size of the scatterers is small compared with the wavelength. Under the hypothesis of fully developed speckle, (Jensen 1996), biologic tissue can be modelled by a network of identical discrete scatterers, randomly distributed in a homogeneous media (Chivers 1977), thus the radio-frequency echo signal can be described as a complex Gaussian probability function with zero mean (Wagner et al. 1983). Envelope detection removes the phase component, generating a speckle distribution with statistics described by a Rayleigh probability function (PRL): ! 2 f ðpÞ 2f ðpÞ PRL ðf ðpÞ; aÞ5 2 exp (5) 2a2 a where f(p) is the image pixel intensity at position p and a is the shape parameterpofffiffi PRL related to the mean of the distribution by m5a p2 . It is worth noting that some commercial ultrasound scanners only provide log-compressed data whose noise statistic is not Rayleigh distributed. Ultrasound raw signal can be recovered from DICOM images using nonlinear compression parametric functions (Seabra and Sanches 2008). Given a spatial neighborhood of pixels with uniform intensity of f(p), a can be computed using the maximum likelihood estimator (MLE) (Steven 1993) defined for Rayleigh functions as: sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 X b5 a f ðqÞ2 (6) 2N q˛U
Speckle reducing bilateral filtering d S. BALOCCO et al.
Design of the s function for ultrasound noise statistics. In the case of images affected by additive zero-mean Gaussian noise, the value of ss is a constant and can be estimated directly from the image. However, the zero-mean assumption on the noise statistics, typical of the classical bilateral filter and anisotropic filter frameworks, is no longer valid for US images, where the noise is multiplicative and modelled by Rayleigh distributions. To extend the BF framework to US envelope images, the noise statistic is embedded in the filter by designing s as the probability that f(p) belongs to the Rayleigh distribution that generated f(x). Given a local neighborhood around p characterized by a Rayleigh distribution RLp(a), the probability of the central pixel intensity f(p) which belongs to RLx(a) can be expressed, combining eqns (5) and (6), as follows: ! f ðpÞ f ðpÞ2 P½f ðpÞ˛RLx ðaÞ5 2 exp 2 (7) b 2x bx 2a a b x is approximated by the pixel value f(x) In this article, a and s(f(x), f(p)) is set equal to P[f(p) ˛ RLx(a)]. Since the BF computes the s function on all pixels in the spatial neighborhood, the errors of the Rayleigh estimation are averaged converging to the local MLE. Therefore, when the central pixel in a neighborhood x and its neighbors belong to a similar Rayleigh distribution, function s yields high values and promotes smoothing. On the contrary, s has a small value when a contour is present, thus preserving the edges. Design of speckle size adaptive c function. The value of sc depends on the image resolution and must be set as a trade-off between computation efficiency and accuracy. An oversized support may promote undesired edge smoothing and increase computational time while a too small sc may hinder the smoothing capability of the BF and require more iterations to converge to the optimal solution. In this article, sc is set equal to the size of the speckle noise measured directly on the image. As proposed by Smith and Wagner (Smith and Wagner 1984), the speckle size can be obtained as the duration of the pulse waveform, which can be estimated by computing the full width at half maximum of the noisy autocorrelation coefficient. Given the autocorrelation matrix M is defined as follows: Mði; jÞ5
M21 X N21 X
f ðm; nÞ,f ðm1i; n1jÞ
(8)
m51 n51
IMis defined as the binary image containing the values of M superior to 12max{M}. The size of the adaptive support c corresponds to the maximum spatial width of the positive values of IM.
1355
RESULTS In silico experiments A set of in silico experiments was designed to evaluate the performance of the proposed approach. The edge-preserving performance of our method has been tested by evaluating the sum of square differences (SSD) criterion between the ground truth and the denoised images, which provides local comparison of the filtering results with the ground truth data. Experimental image set A set of US synthetic images, consisting of two regions, of low (R2) and high (R1) scatterer intensity, respectively, was generated using the convolutional approach (Yu and Acton 2002; Bamber and Dickinson 1980) detailed in the appendix. The echogenicity model amplitude in R1 (I0R1) is unitary while in R2 (I0R2) can be tuned using a parameter g 5 I0R1/I0R2 representing the brightness contrast between the two areas. Since the noise is multiplicative, g indirectly indicates the signal to noise ratio (SNR): the higher the contrast between the two areas, the higher the SNR. In the whole set of images, a 5 MHz central frequency probe is simulated and an anisotropic noise distribution is obtained using a point spread function (PSF) with parameters sx 5 1$k and sy 5 1.5$k, where k will be used as noise scaling factor. Figure 1 illustrates a 128 3 128 pixel image of the set with echogenicity model amplitudes I0R1 5 1 and I0R2 5 0.2. The histogram computed in two uniform regions of the image exhibits two Rayleigh distributions (Fig. 1b), with b s R2 50:272 and mR2 5 0.154 (in the central region) and b s R1 51:55 and mR1 5 1.55 (in the upper bright region), respectively. The ratio between the Rayleigh parameter estimators in the bright (R1) and the dark (R2) areas of the envelope image (Fig. 1a) are equal to the ratio of I0R1 over I0R2 since speckle noise is purely multiplicative. SRBF perfomance In this experiment, SRBF is applied to synthetic images generated with increasing intensities of I0R2 (by varying g while keeping k 5 1) to progressively reduce the contrast between the low and high brightness areas, thus increasing the overlapping of the two Rayleigh distributions. Therefore, it becomes more difficult to discriminate between real signal edges and variations due to the noise. Figure 2 shows the denoising result obtained when the noise intensity increases. As can be observed in Figure 2a in presence of high contrast (g 5 5.7), the algorithm converges in each uniform region to the MLE of the Rayleigh speckle distribution. At the same time, the
1356
Ultrasound in Medicine and Biology
Volume 36, Number 8, 2010
Fig. 1. (a) Synthetic ultrasound envelope-detected images B(x, y) (128 3 128 pixels size), simulated using a 5 MHz central frequency and sx 5 2 and sy 5 1.5 PSF parameters. (b) Histogram and Rayleigh curve fitting in BR1and BR2.
reconstructed image strongly correlates with the echogenicity image I0 showing that the edges have been accurately preserved. Figure 2c illustrates that the performance of the filter reduces proportionally with the decrease of the contrast as shown in Figure 2b (g 5 3.6).
SRBF accuracy and efficiency In this section, SRBF is compared with the state-ofthe-art speckle reducing techniques. Classic bilateral filter (CLASSIC BF) (Tomasi and Manduchi 1998), anisotropic diffusion (AD) (Black et al. 1998), Rayleigh-maximumlikelihood filter (RMLF) (Aysal and Barner 2007), adaptive filter based on second order statistics (AF) (Thakur and Anand 2007), modified Gabor filters (MGF) (Dantas and Costa 2007), Guo filter (GUO) (Guo et al. 2009), mean-shift (MEAN-SHIFT) (EDISON: Code for the Edge Detection and Image SegmentatiON system). (Barash and Comaniciu 2004), speckle reduction anisotropic diffusion (SRAD) (Yu and Acton 2002), oriented speckle reduction anisotropic diffusion (OSRAD) (Krissian et al. 2007) and SRBF are applied to a set of 50 synthetic images generated with the same contrast and speckle size (g 5 5 and k 5 1) to qualitatively (Fig. 3) and quantitatively (Table 1) evaluate the respective performances. Table 1 lists the parameters for each method, compares the features and reports the accuracy scores. As it can be observed in Figure 3a and b, the proposed method presents a superior edge preserving behavior, thus outperforming filtering method based on zero-mean noise assumptions (CLASSIC BF and AF). RMLF, AF, MGF, GUO and MEAN-SHIFT methods are unable to fully remove the speckle noise in the uniform area although they exhibit border preservation properties (see Fig. 3c, d, e, f and g). On the other hand, SRAD and OSRAD better smooth the noise close to the lower
edge of the image because of their edge-enhancing feature. Even though SRAD and OSRAD use a fixed spatial support, the latter achieves better performances since it emphasizes the smoothing behavior along the tangential direction to the object boundary (Fig. 3h and i). These results are in agreement with the average SSD similarity score, computed over the 50 synthetic images, reported in Table 1 (last column). In this experiment, the reconstruction accuracy of SRBF and OSRAD are similar, however, it is worth noting that SRBF is faster (about 20 times) and fully automatic, as summarized in the feature description listed in Table 1. SRBF robustness Even though in the previous experiment the denoising efficacy of OSRAD and SRBF has been found to be similar, the robustness of such filters is investigated when the size of the speckle changes. Both filters are applied to synthetic images generated with different speckle sizes (by varying k while keeping g 5 5.7) and the results are illustrated in Figure 4. As shown in Figure 4a, the performance of OSRAD is comparable to SRBF when the noise pattern is small and progressively decreases in presence of big speckle patterns. In this last case, as it can be seen in Figure 4c, OSRAD enlarges the edges. The same behavior can be observed at different contrasts (g). The performance decrease of OSRAD filter is probably due to its constant size neighbor support while SRBF automatically adapts its spatial support depending on the speckle pattern. In conclusion, SRBF exhibits higher robustness because of its fully automatic settings of relevant parameters, the speckle noise statistics and the speckle size. In vivo experiments Despeckle algorithms are generally used as preprocessing step for automatic segmentation of US images.
Speckle reducing bilateral filtering d S. BALOCCO et al.
1357
tested on intravascular ultrasound (IVUS) images because of their challenging properties: the anisotropy of the speckle patterns, the presence of artifacts and the variable echogeniticy of the tissues.
Fig. 2. Intensity profile and denoising result obtained with g 5 5.7 (a) and g 5 3.6 (b), respectively. Average sum of square differences (SSD) between echogenicity and denoised images (c) while g varies from 20 to 2.86. The plots (a) and (b) correspond to the dashed lines in (c), respectively.
IVUS in vivo images Fifty slices, chosen to represent different vascular structures (plaque and vessel shape, presence of stent, lumen area and diameter), were extracted from a database of 3000 in vivo coronary images. The acquisition was performed using an IVUS Galaxy II System with a catheter Atlantis SR Pro 40 MHz (Boston Scientific, Natick, MA, USA). For each slice, the RF signal is available. From all the 50 slices in polar coordinates, different sets of images are computed: the envelope of the RF signal (data set A), the OSRAD filtered envelope (data set B) and the SRBF filtered envelope (data set C). To avoid the discontinuity at 0/360 degrees of the IVUS image, we mirrored the leftmost (0 ) and rightmost (360 ) texture before filtering, copying 30 of both on the opposite side, thus obtaining a continuous signal at the boundaries. The catheter ring-down artifact located at the top of the polar image was suppressed by filling the region with a uniform colour computed as the mean intensity of the underlying stripe, 0.3 mm thick, adjacent to the catheter area. The same stripe, containing only blood, has been defined as uniform region required by the OSRAD filter. All the images were log compressed and converted to Cartesian coordinates prior to segmentation and visualization. A classic level-set active contour (snake) algorithm (Chan and Vese 2001), commonly used to segment US medical images (Michailovich and Tannenbaum 2007), is applied to the three data sets. The initialization shape and the level-set smoothness parameter are fixed and tuned to obtain the best performance on the envelope data set. The same settings are used when the snake is applied to the data sets B and C. The snake propagates from an initial circular contour until stopping at an energy-stable configuration dependent on the image gradients. The snake inner area is obtained from the zero levelset shape at convergence. Ground truth shapes, indicating the expected segmentation result, were manually delineated in each of the IVUS frames by two independent experts. The operator performing the analysis was blinded to the automatic and to the other manual segmentation results.
The performance of the sole filter cannot be validated on in vivo images since reference denoised data is obviously not available. Instead, the accuracy improvement of a segmentation algorithm is quantified when the images are denoised. For this purpose, the algorithm is compared against the best state-of-the-art filtering method (OSRAD) and
Segmentation comparison Figure 5 presents the manual delineated border and the results of the automatic segmentation for seven representative frames of the three data sets A, B and C . Frame 1 illustrates an accurate automatic segmentation obtained in all the three data sets. In frame 2, the snake crosses the vessel border between the stent wires (highly echogenic
1358
Ultrasound in Medicine and Biology
Volume 36, Number 8, 2010
Fig. 3. Normalized intensity profile of denoised images, superimposed on the ground truth shape I0. Pairwise filtering comparison of (a) CLASSIC BF , (b) AD , (c) RMLF, (d) AF, (e) MGF, (f) GUO, (g) MEAN-SHIFT, (h) SRAD and (i) OSRAD vs. SRBF algorithm.
spots) and spreads in the surrounding tissue. On the contrary, the snake propagation stops correctly at the lumen frontier in data sets B and C, since both filters enhance the outline of the vessel membrane by smoothing the noise and by preserving the sharp vessel edges. Similarly, in frame 3, the denoising methods enhance the borders making the segmentation easier. However, in this case (frame 3) the segmentation of data set B is slightly less accurate since the filtering method fuses a small catheter-ringing artifact with the vessel border, creating an artificial boundary. The same behavior can be observed in frame 4 where, in data set B, the noise lying between the catheter artifact and the vessel is not removed but homogenized.
Frame 5 illustrates that the high amount of noise (data set A) may hamper the snake propagation on the envelope data, while an efficient smoothing promote an accurate and automatic segmentation (data sets B and C). However, it is worth noting that, in data set B (frame 5), the image filtered with OSRAD is locally smoothed but radial ripples are still presents since OSRAD has a fixed support. On the other hand, SRBF better homogenizes uniform areas since the spatial support controlled by the function c is wider and has been designed to automatically contain consecutive scatter peaks. For this reason probably, in frame 6, the snake applied to B is trapped in a local minimum and stops propagating after 560 iterations, while in C the snake keeps spreading inside the whole lumen.
Table 1. Features and performance comparison of speckle reducing filters Support
User interaction
Noised image I0 CLASSIC BF AD
– User defined N 4 neighbors
– Fully automatic Noise variance size (s)
RMLF
8 neighbors
MGF
Freq. filter
2 smoothing coefficients (k1,2) 1 coefficient of variation (t) Probe central frequency k0, Lateral and axial Gabor filter orientations sax, slat
AF
Adaptive W
a0,D, CT, Liny, Hiny, Lasm Hasm, Et as described in [13]
GUO MEAN-SHIFT SRAD OSRAD SRBF
24 neighbors User defined JV 4 neighbors 4 neighbors Adaptive
Fully automatic Noise variance (a) Choice of uniform region Choice of uniform region Fully automatic
Experimental parameters – 8 s 51.3 k1 5 0.5 k2 5 0.05 t 5 0.05 k0 54 sax 55 sax 5 10 W 5 2, a0 5 0.04 D 5 2, CT 5 3.2, Et 5 1 Liny 5 04, Lasm 5 0.05 Hiny 5 0.8, Hasm 5 0.2 None N 5 6, s 5 1 top stripe of IR2 top stripe of IR2 None
Computat. time (1283128 pixels)
SSD (mean 6 SD)
6 300
– 0.364 [s] 1.411 [s]
7104 6 180 2736 6 130 3692 6 139
1
0.028 [s]
3138 6 114
1
0.055 [s]
2384 6 190
Iterations n. –
1 1 10 1000 1000 6
47 [s] 3.117 [s] 3.312 [s] 1.57 [s] 6.305 [s] 0.364 [s]
2810 6 148 2622 6 129 2246 6 129 1913 6153 1635 6 168 1288 6 127
CLASSIC BF 5 classic bilateral filter; AD 5 anisotropic diffusion; RMLF 5 Rayleigh-maximum-likelihood filter; MGF 5 modified Gabor filter; AF 5 adaptive filter; GUO 5 Guo filter; MEAN-SHIFT 5 mean-shift; SRAD 5 speckle reduction anisotropic diffusion; OSRAD 5 oriented speckle reduction anisotropic diffusion; SRBF 5 speckle reducing bilateral filter.
Speckle reducing bilateral filtering d S. BALOCCO et al.
Fig. 4. Average sum of square differences (SSD) (a) between the echogenicity image I0 and the image filtered using OSRAD and SRBF when the speckle noise scaling factor k ranges from 0.5 to 3.5. Intensity profile of I0 (b-c), superimposed on the OSRAD and SRBF filtering result when the noise scaling is equal to 1 and to 2.5, respectively.
Finally in frame 7, the snake spans outside the vessel for data sets A and B, but not in C. In this frame, the OSRAD filter enhances the shadow area on the left of the catheter, and creates an undesired gradient attracting the snake in B, while SRBF homogenizes the shadow region with the rest of the lumen in C. In fact, compared with OSRAD, which is a zero mean filter, the SRBF noise model better fits the speckle distribution at low image intensities and discriminate more accurately between noise and tissues. The last column of Figure 5 represents the average distance of the automatic segmentation to the two manual segmentations. In frame 1, the high slope of the curve indicates that the snake converges faster in filtered data sets (B and C) since the noise is accurately removed. In frame 4, a similar behavior can be observed, but the segmentation of B results as the slowest and the less accurate because of the wrinkled lumen region. In frames 2, 3
1359
and 5, the segmentation applied to the nonfiltered image diverges while the two filters show similar convergence speed. These examples illustrate the interest of using an accurate filtering method before the segmentation stage. In frame 6, the snake applied to data set B stops expanding earlier than data set C leading to incorrect vessel segmentation and to lumen area under-estimation. Finally, in frame 7, the SRBF filtered the snake applied to data set A and B diverges and only the snake applied to the data set C converges to the vessel wall. In all the cases, the segmentation on data set C is either the fastest or the only one converging to the vessel boundary, showing the superior performances to SRBF. Figures 6 and 7 show a pair wise error comparison for all the frames of the three data sets (A, B and C). In particular, Figure 6 illustrates the error between the automatic and ground-truth border segmentations while Figure 7 shows the error in the lumen area estimations. In Figure 6a most of the markers are located below the solid line of unitary slope, indicating that the average distance between the automatic and the manual segmentations is smaller in the filtered data sets (B, C) then in the nonfiltered (A). Figure 6b, shows the superior accuracy obtained by filtering the images with SRBF. The average distance over all the images of the data sets A, B and C with respect to the manual segmentations are 0.27 6 0.47 mm, 0.15 6 0.24 mm and 0.12 6 0.14 mm, respectively, confirming that the contours obtained from the data set C exhibit the lowest segmentation error. The average distance over all the images of the interobserver variability is 0.021 mm. To quantify the error in estimating the area for the three data sets the Jaccard coefficient (Jaccard 1912) jFXJj defined as the ratio J5jFWJj is used, where F is the set of pixels of the manual segmentation and J is the set of pixels resulting from automatic segmentation. The score J is one when the areas matches exactly and rapidly decreases when one of the two areas differs from the other. The performance of the automatic segmentation on the filtered data sets B and C, is superior to the non filtered images A (see Fig. 7a), since J is higher than the line with unitary slope for most of the images. Analyzing the scatter plot of Figure 7b and comparing directly the performance of the snake on the two denoised data sets, it can be noticed that the J-score is similar but some outliers, representing unsuccessful segmentations of data set B, are present. In more than 50% of the images , the J-score corresponding to OSRAD filter ranges from 0.12 to 0.33 while for SRBF J varies from 0.06 to 0.24 (Fig. 6b). Analyzing the scatter plot of Figure 7b in detail, it can be noticed that in about 50% of the images, OSRAD J-score ranges from 0.2 to 0.8 while for SRBF J varies from 0.6 to 0.9. The average J-score over all the images of the data sets A, B and C are 0.41 6 0.29, 0.76 6 0.18 and 0.83 6 0.09
1360
Ultrasound in Medicine and Biology
Volume 36, Number 8, 2010
Fig. 5. Snake manual (first column) and automatic (second, third and forth columns) segmentations of seven representative frames from the three data sets A, B and C . The last column illustrates the average distance of the automatic segmentations to the manual segmentations.
Speckle reducing bilateral filtering d S. BALOCCO et al.
1361
Fig. 6. Average distance of the automatic segmentation to the two manual segmentations. (a) comparison of data sets B and C vs. A ; (b) pair wise comparison of data sets B vs. C.
respectively, showing that the areas recovered on the data set C are the most similar to the manually segmented regions. ANOVA and t-test statistical analysis were performed for both segmentation distances and J-scores in data sets A, B and C. In all cases the scores of data sets B and C are significantly different (p values ,1$10–9 for both ANOVA and t-test) from data set A. Comparing data sets B and C, the ANOVA provides p values of 0.0032 and 0.0153, respectively, for the distances and J-scores; and the t-test provides p values of 5.5$10–4 and 0.0039, respectively, for the distances and J-scores. These tests confirm the statistical significance of the results. The J-score over all the images of the interobserver variability is 0.95. DISCUSSION AND CONCLUSIONS In this article, a bilateral filter tailored to US images is proposed. The edge preservation property is obtained by embedding noise statistics in the filter framework. The filter implicitly estimates the Rayleigh distribution parameters on the pixel neighborhood. Moreover, its spatial
support is automatically defined by estimating the size of the speckle noise directly on the image. A series of in silico experiments has been designed with the aim to compare the performances of the proposed approach against other denoising methods on synthetic images corrupted by speckle noise. As expected, SRBF outperforms the filtering methods based on zero mean noise assumptions. Compared with OSRAD, our approach is fully automatic, locally adaptive and the number of iterations needed to reach the optimal solution is smaller. Additionally, SRBF results in a more robust algorithm suitable for a wide ranges of speckle noise sizes. The proposed method has been tested on 50 in vivo IVUS images, by quantifying the influence of the denoising filter on the segmentation accuracy of a generic level-set algorithm. For this purpose, three data sets were created: the first composed of simple envelope images and the other two filtered using OSRAD and SRBF algorithms. The use of a more sophisticated segmentation algorithm, such as active contours with shape priors (Eduard et al. 2008), statistic snakes (Ayache et al. 1993) or including a leakage detection feature would improve the
Fig. 7. Area ratio J representing the likeness between the automatically segmented areas. (a) comparison of data sets B and C vs. A; (b) pair wise comparison of data sets B and C.
1362
Ultrasound in Medicine and Biology
segmentation result and succeed on all the envelope data sets, while a weakly constrained snake performs according to the quality of the original data. As shown in Figures 6 and 7, an accurate denoising method able to enhance the boundary and remove the noise combined with an ordinary snake reduces one third the average distance between the manual and the automatically segmented lumen. Considering that the average diameter of the coronary artery of the used image database is about 3 mm, the segmentation error reduces from about 0.3 mm to 0.1 mm representing 20% and 6% of the arterial radius, respectively. The experiment shows that the accuracy of the snake applied to filtered images is superior to that obtained in envelope data sets since both despeckle algorithms enhance the outline of the vessel membrane and preserve the sharp vessel edges. However, as expected, in some cases (frames 6 and 7) the snake applied to OSRAD filtered images diverges or stops too early, leading to incorrect segmentation results. Indeed, SRBF not only locally smoothes the images, but fills the gaps in uniform areas accurately given its adaptive spatial support. For this reason, the snake propagates faster as demonstrated by the low number of iterations required to converge to the correct vessel contour. The adaptive scaling of the filter, obtained through the estimation of the speckle size, is a feature that can be implemented only in denoising methods allowing a support size adjustment (see Table 1), thus, OSRAD filter cannot accommodate this feature. Moreover, in almost all the 50 tested images, the contours and the lumen areas obtained on the SRBF filtered data set were the closest to the manual delineated contours. Finally, it is worth noting that SRBF is fully automatic and does not require the definition of a user-provided uniform region to tune the intensity of the filtering, thus making it suitable for unsupervised filtering of large amounts of clinical data. In summary, SRBF is a fast and flexible algorithm potentially suitable in different medical applications (IVUS, B-mode, 3-D matrix array US). Future work will include studying a possible extension of the proposed filter to more complex noise models, by replacing the Rayleigh function with a Rician distribution, to extend the filter applicability to other medical modalities (e.g., MRI images). Acknowledgments—The works of S. Balocco and C. Gatta are supported by the Catalan Agency for Administration of University and Research (AGAUR) under a Beatriu de Pinos Fellowship. The authors would like to thank Dr. T. Aysal and Dr. K. Barner for providing the source code of the RMLF filter and Dr. K. Krissian for providing a Windows compiled version of the OSRAD filter and for the fruitful discussions about anisotropic diffusion filters.
REFERENCES Ayache N, Cohen I, Herlin I. Medical image tracking. Cambridge, MA: MIT Press; 1993. p. 285–302.
Volume 36, Number 8, 2010 Aysal TC, Barner KE. Rayleigh-maximum-likelihood filtering for speckle reduction of ultrasound images. IEEE Trans Med Imaging 2007;26:712–727. Bamber JC, Dickinson RJ. Ultrasonic b-scanning: A computer simulation. Phys Med Biol 1980;25:463–479. Barash D, Comaniciu D. A common framework for nonlinear diffusion, adaptive smoothing, bilateral filtering and mean shift. Image Vis Comput 2004;22:73–81. Black M, Sapiro G, Marimont D, Heeger D. Robust anisotropic diffusion. IEEE Trans Image Process 1998;7:421–432. Chan TF, Vese LA. Active contours without edges. IEEE Trans Image Process 2001;10:266–277. Chivers R. The scattering of ultrasound by human tissues, some theorical models. Ultrasound Med Biol 1977;3:1–13. Choudhury P, Tumblin J. The trilateral filter for high contrast images and meshes. International Conference on Computer Graphics and Interactive Techniques ACM SIGGRAPH 2005 Courses Los Angeles, CA: SESSION: Computational photography. Article No.: 5. 2005. Comaniciu D, Meer P. Mean shift: A robust approach toward feature space analysis. IEEE Trans PAMI 2002;24:603–619. Dantas RG, Costa ET. Ultrasound speckle reduction using modified Gabor filters. IEEE Trans Ultrason Ferroelectr Freq Control 2007;54:530–538. Durand F, Dorsey J. Fast bilateral filtering for the display of high-dynamicrange images. New York, NY: ACM Press; 2002. p. 257–266. Guo Y, Cheng HD, Tian J, Zhang Y. A novel approach to speckle reduction in ultrasound imaging. Ultrasound Med Biol 2009;35:628–640. Jaccard P. The distribution of the flora in the alpine zone. New Phytologist 1912;11:37–50. Jensen J. Estimation of blood velocities using ultrasound: A signal processing approach. Cambridge: Cambridge University Press; 1996. Kay SM. Fundamentals of statistical signal processing: Estimation theory. Upper Saddle River, NJ: Prentice-Hall, Inc; 1993. p. 595. Krissian K, Westin CF, Kikinis R, Vosburgh KG. Oriented speckle reducing anisotropic diffusion. IEEE Trans Image Process 2007;16: 1412–1424. Li B, Acton ST. Active contour external force using vector field convolution for image segmentation. IEEE Trans Image Process 2007;16: 2096–2106. Michailovich O, Tannenbaum A. Segmentation of medical ultrasound images using active contours. Proceedings of the IEEE International Conference on Image Processing (ICIP). Sep 2007;513–516. Munteanu C, Morales FC, Fernandez JG, Rosa A, Deniz LG. Enhancing obstetric and gynecology ultrasound images by adaptation of the speckle reducing anisotropic diffusion filter. Artif Intell Med 2008; 43:223–242. Perona P, Malik J. Scale space and edge detection using anisotropic diffusion. IEEE Trans Pattern Anal Machine Intell 1990;12:629–639. Seabra J, Sanches J. Modeling log-compressed ultrasound images for radio frequency signal recovery. Conf Proc IEEE Eng Med Biol Soc 2008;2008:426–429. Smith SF, Wagner RF. Ultrasound speckle size and lesion signal to noise ratio: Verification of theory. Ultrason Imaging 1984;6:174–180. Sojka E, Gaura J, Krumnikl M. Active contours without edges and with simple shape priors source. Lecture Notes In Computer Science; Vol. 5259 archive. Proceedings of the 10th International Conference on Advanced Concepts for Intelligent Vision Systems. Juan-les-Pins, France 2008;730–741. Thakur A, Anand RS. Speckle reduction in ultrasound medical images using adaptive filter based on second order statistics. J Med Eng Technol 2007;31:263–279. Tomasi C, Manduchi R. Bilateral filtering for gray and color images. In Proceedings of the Sixth International Conference on Computer Vision, Bombay, India, January 1998. Tsai DY, Watanabe S. A method for optimization of fuzzy reasoning by genetic algorithms and its application to discrimination of myocardial heart disease. 1998;. 3:1756–1761. Wagner RF, Smith SW, Sandrik JM, Lopez H. Statistics of speckle in ultrasound b-scans. IEEE Trans Sonics Ultrason 1983;30:156–163. Yu Y, Acton S. Speckle reducing anisotropic diffusion. IEEE Trans Image Process 2002;11:1260–1271.
Speckle reducing bilateral filtering d S. BALOCCO et al.
APPENDIX In this study, synthetic images corrupted by speckle noise have been generated using a convolutional approach described by (Yu and Acton 2002; Bamber and Dickinson 1980). The simulation assumes that the imaging system has a linear, space-invariant point spread function (PSF) and the transducer is linear. Let I0(x, y) be an echogenicity model (an image with different intensities corresponding to different tissues of the object being imaged) in which the variables x and y are the lateral and axial coordinates, respectively. First, subresolution variations in object impedance should be introduced by multiplying the echogenicity model by a Gaussian white noise with zero mean and unitary variance: Tðx; yÞ5I0 ðx; yÞ$Gðx; yÞ
(9)
T(x, y), accounts for acoustic impedance inhomogeneities in the object due to density and acoustic speed perturbation that generates the scattering. The V(x, y) ultrasonic radio-frequency (RF) echo data can then be obtained by a convolution: Vðx; yÞ5hðx; yÞ Tðx; yÞ
(10)
1363
where h(x, y) is the point spread function or impulse response of a hypothetical US imaging system. It is assumed (Yu and Acton 2002; Bamber and Dickinson 1980) that h(x, y) is separable, i.e., h(x, y) 5 h1(x) h2(y) where h1 is a Gaussian-weighted sinusoidal function determined by: h1 ðx; sx Þ5sinðk0 xÞ exp x2 =2s2x
(11)
where k0 5 2pf0/v, v is the speed of sound in tissue, f0 is the center frequency, and sx represents the pulse-width of transmitting ultrasonic wave. The spatial response for the transmitting and receiving aperture h2(x) is determined by: h2 ðyÞ5exp y2 =2s2y
(12)
where sy represents the beam-width of transmitting ultrasonic wave. The image of the envelope-detected amplitude, B(x, y) shown in Figure 1a, is given by: Bðx; yÞ5jVðx; yÞ1i Vbðx; yÞj
(13)
where Vbðx; yÞ is the Hilbert transform of V(x, y) and i is the imaginary unit.