Breast Ultrasound Despeckling Using Anisotropic Diffusion Guided by Texture Descriptors

Ultrasound in Med. & Biol., Vol. -, No. -, pp. 1–13, 2014 Copyright
2014 World Federation for Ultrasound in Medicine & Biology Printed in the USA. All rights reserved 0301-5629/$ - see front matter

http://dx.doi.org/10.1016/j.ultrasmedbio.2014.06.005

d

Original Contribution BREAST ULTRASOUND DESPECKLING USING ANISOTROPIC DIFFUSION GUIDED BY TEXTURE DESCRIPTORS y OMEZ FLORES,* WAGNER COELHO DE ALBUQUERQUE PEREIRA, WILFRIDO G y and ANTONIO FERNANDO CATELLI INFANTOSI

* Technology Information Laboratory, Center for Research and Advanced Studies of the National Polytechnic Institute, Ciudad Victoria, Tamaulipas, Mexico; and y Biomedical Engineering Program/COPPE, Federal University of Rio de Janeiro, Rio de Janeiro, Brazil (Received 4 October 2013; revised 22 May 2014; in final form 4 June 2014)

Abstract—Breast ultrasound (BUS) is considered the most important adjunct method to mammography for diagnosing cancer. However, this image modality suffers from an intrinsic artifact called speckle noise, which degrades spatial and contrast resolution and obscures the screened anatomy. Hence, it is necessary to reduce speckle artifacts before performing image analysis by means of computer-aided diagnosis systems, for example. In addition, the trade-off between smoothing level and preservation of lesion contour details should be addressed by speckle reduction schemes. In this scenario, we propose a BUS despeckling method based on anisotropic diffusion guided by Log–Gabor filters (ADLG). Because we assume that different breast tissues have distinct textures, in our approach we perform a multichannel decomposition of the BUS image using Log–Gabor filters. Next, the conduction coefficient of anisotropic diffusion filtering is computed using texture responses instead of intensity values as stated originally. The proposed algorithm is validated using both synthetic and real breast data sets, with 900 and 50 images, respectively. The performance measures are compared with four existing speckle reduction schemes based on anisotropic diffusion: conventional anisotropic diffusion filtering (CADF), specklereducing anisotropic diffusion (SRAD), texture-oriented anisotropic diffusion (TOAD), and interference-based speckle filtering followed by anisotropic diffusion (ISFAD). The validity metrics are the Pratt’s figure of merit, for synthetic images, and the mean radial distance (in pixels), for real sonographies. Figure of merit and mean radial distance indices should tend toward ‘1’ and ‘0’, respectively, to indicate adequate edge preservation. The results suggest that ADLG outperforms the four speckle removal filters compared with respect to simulated and real BUS images. For each method—ADLG, CADF, SRAD, TOAD and ISFAD—the figure of merit median values are 0.83, 0.40, 0.39, 0.51 and 0.59, and the mean radial distance median results are 4.19, 6.29, 6.39, 6.43 and 5.88. (E-mail: [email protected])
2014 World Federation for Ultrasound in Medicine & Biology. Key Words: Breast ultrasound, Speckle filtering, Anisotropic diffusion, Log–Gabor filters, Texture features.

Ultrasound image quality is affected mainly by an inherent imaging artifact called speckle, which results from interference effects between returning echoes produced by discontinuities of tissue below ultrasonic beam resolution (unresolved scatterers) (Thijssen 2003). Speckle can be interpreted as a locally correlated noise (or random granular texture) that degrades the US image by concealing fine structures and reducing the signal-tonoise ratio (SNR). Further, speckle tends to reduce the image contrast and to obscure and blur image details (Michailovich and Tannenbaum 2006). With respect to BUS images, the speckle could make human interpretation difficult and, consequently, influence inter-/intra-observer variations. Moreover, computer-aided diagnosis (CAD) systems commonly extract shape or contour features from segmented breast lesions to classify them as benign or malignant (Cheng

INTRODUCTION Currently, breast ultrasound (BUS) is the most important adjunct to mammography for patients with palpable masses and normal or inconclusive mammograms. Additionally, BUS has the ability to visualize hidden lesions in women with dense breast tissue (Corsetti et al. 2011) and is particularly useful in distinguishing cystic from solid lesions, with an accuracy of nearly 100% (Jackson 1990). BUS is also used to differentiate between benign and malignant tumors, which can be characterized by their shapes, borders, internal echo features and posterior acoustic behavior (Rahbar et al. 1999). Address correspondence to: Wilfrido Gomez-Flores, Technology Information Laboratory, CINVESTAV-IPN, Science and Technology Park TECNOTAM, Km. 5.5 Ciudad Victoria-Soto La Marina Road, ZIP 87130, Ciudad Victoria, Tamaulipas, Mexico. E-mail: [email protected] 1

2

Ultrasound in Medicine and Biology

et al. 2010). However, the performance of the segmentation stage depends not only on its technical strategy, but also on adequate image preprocessing schemes for reducing the speckle. Several filtering approaches aimed at reducing speckle in US images while preserving edge details have been proposed in the literature. Popular despeckling techniques based on local statistics include median filter (Horsch et al. 2001), truncated median filter (Evans and Nixon 1993), adaptive weighted median filter (Loupas et al. 1989) and the well-known Frost, Lee and Kuan filters (Lopes et al. 1990). In addition, the use of directive filtering, such as modified Gabor filters (Dantas and Costa 2007) or the ‘‘stick’’ method (Czerwinski et al. 1999), to reduce ultrasonic speckle has been investigated. Morphologic filters have also been used in BUS images to preserve lesion contour anfractuosities (Alvarenga et al. 2009; Infantosi et al. 2008). Furthermore, non-linear anisotropic diffusion filtering (ADF), proposed by Perona and Malik (1990), has attracted much attention because it is capable of reducing noise in images without blurring the boundaries between homogeneous regions. Several variants of ADF have been developed for US images, including speckle-reducing anisotropic diffusion (SRAD) (Yu and Acton 2002), median-guided anisotropic diffusion (MADF) (Yang and Fox 2004), texture-oriented anisotropic diffusion (TOAD) (Aleman-Flores et al. 2007) and interference-based speckle filtering followed by anisotropic diffusion (ISFAD) (Cardoso et al. 2012). The main advantage of ADF is the computation of the gradient-based conduction coefficient to stop the diffusion process through ‘‘strong’’ or significant edges. Thus, the speckle is reduced efficiently within homogeneous regions while important edges details are preserved (Yu and Acton 2002). This characteristic is desired when filtering BUS images to avoid overblurring lesion contours for further lesion segmentation and extraction of useful shape and contour features by means of CAD systems. However, the reduction of speckle in BUS images is a difficult task because of the large variance in lesions shapes and low contrasts produced by shadows, echo features and blurry or ill-defined boundaries. Therefore, if the local contrast between breast lesion and adjacent tissues is poor, the diffusion process will pass through the edges. To overcome this drawback, texture descriptors have been employed to guide the diffusion process in the ADF approach by using the responses of Gabor filters (Aleman-Flores et al. 2007). In this context, texture can be defined as the spatial variation of pixel intensities at scales smaller than the scales of interest (Petrou and Garcıa-Sevilla 2006). The main idea is to perform a multichannel decomposition of the BUS image to find

Volume -, Number -, 2014

strong edges between tissues with different textures. Therefore, the conduction coefficient is computed by using texture responses instead of pure intensity values. Despite finding strong edges between distinct textured regions could increase the filter performance, in terms of boundary preservation, it is important to consider a suitable value of the edge magnitude parameter in the conduction coefficient, denoted as k, to control the diffusion extension of the ADF process. Such a parameter is commonly fixed either heuristically by the user or by computing a ‘‘noise estimator’’ (Perona and Malik 1990). Hence, depending on its value the filter should perform as an all-pass filter, on strong edges, or as Gaussian smoothing, on homogeneous regions. Also, the k value is considered global; that is, all pixels in the image are treated equally. Thus, the k parameter should be adequately chosen to cope with some characteristics of the image such as noise power and local contrast. However, these characteristics could vary among BUS images, depending on the ultrasound equipment, operator skills, scanned tissues or inherent artifacts (Feldman et al. 2009). The goal of the work described in this article was to reduce the speckle within regions with similar textures while avoiding overblurring of tissue-structure edges. In this sense, we were interested in preserving breast lesion shape for CAD purposes to enhance further tasks such as lesion segmentation. We proposed to adapt the conduction coefficient parameter, k, for each pixel in the BUS image by using 2-D Log–Gabor filters to depict textures in specific directions. Also, we compared the performance of the proposed technique with that of four ADF-based techniques used in US images. METHODS Anisotropic diffusion filtering The 4- or 8-nearest neighbor (N 4 or N 8 ) discretization of the non-linear partial differential equation for the ADF approach is expressed as (Perona and Malik 1990) X ½gðjVd IjÞVd Iti;j (1) Ii;jt11 5 Ii;jt 1t d˛N

d

where t is the iteration step; Ii;jt is the noisy pixel at iteration t; the pair i, j is the pixel location; 0 , t # ¼ for the numerical approximation to be stable; jxj denotes the magnitude; N d indicates the set of d-directions for the nearest-neighbor difference, N 4 5 {N, S, W, E} or N 8 5 {N 4 , NW, SW, NE, SE}, denoted by the symbol V: VN I 5 Ii;j21 2Ii;j ; VNW I 5 Ii21;j21 2Ii;j VS I 5 Ii;j11 2Ii;j ; VSW I 5 Ii21;j11 2Ii;j VW I 5 Ii21;j 2Ii;j ; VNE I 5 Ii11;j21 2Ii;j VE I 5 Ii11;j 2Ii;j ; VSE I 5 Ii11;j11 2Ii;j

(2)

Breast US despeckling d W. G OMEZ FLORES et al.

3

Fig. 1. Effect of anisotropic diffusion filtering with a global k value, where the upper half is the noisy image and the lower half is the filtered image: (a) original free-noise image; (b, c) low- and high-power noise, respectively; (d, e) low- and high-power noise with an illumination gradient, respectively. For all examples, the expression in eqn (4) is used as the conduction coefficient, where its parameters are k 5 5, t 5 1=4, 500 iterations and 4-nearest neighbors (Weickert 2008).

The function g(x) in eqn (1) is the conduction coefficient updated at every iteration defined as   gðVIÞ 5 exp 2ðjVIj=kÞ2 (3) or as gðVIÞ 5

1 11ðjVIj=kÞ2

(4)

where k is an edge magnitude parameter that controls the diffusion extension. Equation (3) enhances high-contrast edges over low-contrast ones, whereas eqn (4) enhances wide regions over smaller ones. The gradient magnitude jVIj detects the image edges or boundaries as a steep discontinuity in intensity. Hence, it guides the diffusion process through the conduction coefficient. If jVIj .. k, then g(jVIj) / 0, resulting in an all-pass filter. On the other hand, if jVIj 5 k, then g(jVIj) / 1, implying an isotropic diffusion process (Gaussian filtering) (Yu and Acton 2002). Therefore, the anisotropic diffusion has the property of blurring small discontinuities and sharpening edges. The selection of a suitable k value then depends strongly on the noise power and regional contrast, and it should be adapted for each pixel in the image instead of using a global k parameter. Figure 1 illustrates the effect of the ADF approach with a global k value for distinct noise power and local contrasts in a speckled image. Figure 1 (a) is the freenoise simulated image (256 3 256 pixels) with three constant intensity regions: dark (zero), medium (127) and bright (230). The images in Figure 1 (b, c) are corrupted with multiplicative speckle noise with variances equal to 0.1 and 0.5, respectively, for simulating ‘‘low’’ and ‘‘high’’ power noise. The images in Figure 1 (d, e) correspond to Figure 1 (b, c), respectively, but include a linear illumination gradient, in the range [0.5, 1.0], for varying the local contrast. The advantages of the ADF approach include intraregion smoothing and edge preservation, as illustrated in Figure 1 (b). However, as observed in Figure 1 (c, e), selection of k the parameter is critical; depending on its

value, the filter should perform as an all-pass filter, on strong edges, or as Gaussian smoothing, on homogeneous regions. Therefore, it should be set in accordance with both the amount of noise and the local contrast. In this work, we use texture features, based on Log–Gabor filters, to compute the local value of the k parameter for preserving lesion edges in BUS images.

Log–Gabor filters Aleman-Flores et al. (2007) proposed the TOAD approach, which is based on the responses of Gabor filters for reducing speckle in BUS images. In this context, a texture feature associated with some textured region depicts locally some texture-relevant characteristic of such a region. Hence, the assumption is that different tissues have distinct textures, which could be represented by a vector of texture features. In TOAD, the gradient of the intensity image in the conduction coefficient (eqn 3) is changed by the gradient of texture responses, depicted by a bank of spatial Gabor filters built as  2 2   x 1y (5) cos kx x1ky y Gskx ;ky ðx; yÞ 5 exp 2 2 2s where s is the scale factor of a spatial (x, y) region; and k is the spatial frequency. The TOAD approach uses a set of S scales, N frequencies within each scale and two orientations (vertical and horizontal). Then, for a fixed scale, s0.0, and frequency, k0 .0, the texture responses can be defined as the convolution 0 FnSx ;ny ðx; yÞ 5 Iðx; yÞGss nx k0 ;ny k0 ðx; yÞ

(6)

where 1 # s # S, –N # nx # N, 0 # ny # N and I (x, y) is the original noisy image. The parameters of the Gabor kernels are determined heuristically as k0 5 0:05, s0 5 3, S 5 3 and N 5 4; the k parameter is undefined. Then, the gradient magnitude in eqn (3) is calculated by using the distinct texture responses from eqn (6). In the frequency domain, the Gabor transfer function is constructed as the sum of two Gaussians around the

4

Ultrasound in Medicine and Biology

center frequency. If the bandwidth of these Gaussians increase more than one-third of the center frequency, their tails will overlap at the origin, resulting in a nonzero direct current (DC) component. Hence, it is difficult to construct Gabor functions of arbitrary wide bandwidth (.1 octave) and still maintain a reasonably small DC component (Kovesi 1997). To cope with this limitation, we used Log–Gabor filters that can be constructed with arbitrary bandwidth, minimal overlap and practically without a DC component, which contributes to improving the contrast between different textures (Field 1987; Wang et al. 2008). The frequency response of the oriented Log–Gabor function is expressed as (Field 1987; Wang et al. 2008)



 Gðu; qÞ 5 exp 2logðu=u0Þ2 2s2u exp 2ðq2q0 Þ2 2s2q (7) where (u, q) denotes the polar coordinates; u0 is the filter center frequency; q0 is the filter orientation angle; su defines the frequency bandwidth; and sq determines the angular bandwidth. The ratio sq/su defines the filter bandwidth; a value of 0.75 will result in a filter bandwidth of approximately one octave, 0.55 will result in two octaves and 0.41 will produce three octaves (Field 1987). The filtered image by a Log–Gabor function in the frequency domain can thus be expressed as 0

F ðu; vÞ 5 Fðu; vÞ$Gðu; vÞ

(8)

where Gðu; vÞ is the expression of the Log–Gabor filter in a Cartesian coordinate system; and Fðu; vÞ and F 0 ðu; vÞ are the original and filtered spectra, respectively. After a bank of Log–Gabor filters complete the filtering, the real part of the inverse Fourier transform of each filtered spectrum is obtained to recover the spatial domain data which contain texture information. Herein, the Log–Gabor filters are constructed in terms of both radial and angular components (Kovesi 1997). The former controls the frequency band and requires the number of scales, Sn , whereas the latter controls the orientation selectivity and requires the number of orientations, qn. In addition, for computing the radial component, the smallest-scale filter, lmin, and the scaling factor between frequencies, S, are also considered. The set values of these parameters are presented elsewhere in this article.

Proposed despeckling approach With the aim of stopping the diffusion process through boundaries that separate regions with different textures and distinct regional contrasts, the local response of the conduction coefficient should be enhanced. Because the k value determines if some noisy pixel is

Volume -, Number -, 2014

blurred or not, we proposed estimating a directionallocal kdi;j parameter for each pixel in the original image, 0 , based on texture information given by Log–Gabor Ii;j functions. The conduction coefficient of eqn (4) then becomes  s g Vd Ii;jt ; kdi;j 5  (9)

.

t 2 d s1 Vd Ii;j ki;j where d is a direction in the set N 8 . Because Ii;jt is normalized to the range ½0; 1, the constant s 5 0.01 for stability. Let Rri;j denote the response of the rth Log–Gabor fil0 (i.e., inverse Fourier ter applied to the original image Ii;j d transform of eqn 8). Then, ki;j is calculated in two stages. First, the N 8 differences of each texture response are determined according to eqn (2), resulting in the vector VN 8 Rri;j 5 ðVN Rri;j ; .; VSE Rri;j Þ. Second, kdi;j is computed by integrating all texture-gradient magnitudes with specific direction d as sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

2ffi X

r d ki;j 5 (10)

Vd Ri;j 1#r#R

Once each kdi;j has been determined, the anisotropic diffusion is performed during tmax iterations. The N 8 dift ferences in eqn (2) are obtained from the current Ii;j filtered image, for 0 # t # tmax. These discrete differences are grouped into eight orthogonal pairs—p1 5 (W, N), p2 5 (N, E), p3 5 (E, S), p4 5 (S, W), p5 5 (NW, NE), p6 5 (NE, SE), p7 5 (SE, SE) and p8 5 (SW, NW)—to form the following gradient vectors (Aleman-Flores et al. 2007): Vp1 Ii;jt 5 ðVW I; VN IÞti;j ; Vp5 Ii;jt 5 ðVNW I; VNE IÞti;j t t Vp2 Ii;jt 5 ðVN I; VE IÞi;j ; Vp6 Ii;jt 5 ðVNE I; VSE IÞi;j t t t t Vp3 Ii;j 5 ðVE I; VS IÞi;j ; Vp7 Ii;j 5 ðVSE I; VSW IÞi;j t Vp4 Ii;jt 5 ðVS I; VW IÞi;j ; Vp8 Ii;jt 5 ðVSW I; VNW IÞti;j (11) For the sake of clarity, the eight orthogonal pairs are grouped in the set D 5 {p1, p2, ., p8}. Similarly, the distinct directions of kdi;j are grouped in agreement to the set of orthogonal pairs D . For instance, considering the orthogonal pair p1 5 ðW; NÞ, we have the N vector kpi;j1 5 ðkW i;j ; ki;j Þ. Thereafter, the conduction coefficient is determined by using eqn (9) to build the vector t t t N ; kpi;j1 Þ 5 ½gðVW Ii;j ; kW gðVp1 Ii;j i;j Þ; gðVN Ii;j ; ki;j Þ. To simplify t the notation of this expression, we use gp1 5 ðgtW ; gtN Þ for each location ði; jÞ. Obviously, this operation is applied to every orthogonal pair in the set D . Next, the product of the conduction coefficient vector and the image difference vector is determined. For example,

Breast US despeckling d W. G OMEZ FLORES et al.

with respect to the orthogonal directions N and W, we t t t obtain the vector gtp1 $Vp1 Ii;j 5 ðgtW $VW Ii;j ; gtN $VN Ii;j Þ. Subsequently, the divergence operator is applied on this resultant vector by summing their elements as t t t Þ 5 gtW ,VW Ii;j 1gtN ,VN Ii;j . Therefore, at divðgtp1 ,Vp1 Ii;j each iteration t, we obtain eight divergences corresponding to the orthogonal pairs in set D . Finally, the solution is updated as Ii;jt11 5 Ii;jt 1t

8 h X q51

it  div gpq $Vpq I

(12)

i;j

where Ii;jt11 ˛½0; 1. In the following, this proposed approach will be denoted as ADLG, that is, anisotropic diffusion filtering guided by Log–Gabor filters. Other ADF-based approaches The performance of the ADLG approach is compared with that of the four ADF-based techniques, which have been used to reduce speckle in US images, and a brief description of each is provided here. Conventional ADF (CADF). This technique is the original anisotropic diffusion approach proposed by Perona and Malik (1990). This technique applies the N 4 -discretization of the non-linear partial differential equation expressed in eqn (1). Also, in this study, the diffusion coefficient in eqn (4) is used. Speckle-reducing anisotropic diffusion (SRAD). This method, proposed by Yu and Acton (2002), uses eqn (1) to reduce speckle iteratively in ultrasound data. In this approach, the conduction coefficient depends on the coefficient of variation that is determined as a function of the ratio of the local variance to the local squared mean of the pixels in a homogeneous region. Thus, the coefficient of variation determines whether a pixel should be smoothed or left intact. However, this algorithm requires a homogeneous area inside the image being processed, and the authors did not suggest any methodology to automatically select this region. We propose dividing the input image with a lattice of nonoverlapping partitions of 20 3 20 pixels. Next, the coefficient of variation of each block is computed as s2i =m2i , where si and mi are the standard deviation and mean value, respectively, of pixels within the ith block. Finally, the block with the lowest value of the coefficient of variation is selected as the homogeneous region. Texture-oriented anisotropic diffusion (TOAD). This approach, proposed by Aleman-Flores et al. (2007), uses texture responses of Gabor filters in eqn (6) to guide the conduction coefficient in eqn (3). As depicted previously, texture-gradient magnitudes are computed, kVRk, to find certain edges despite low-intensity contrast between

5

distinct tissues. Therefore, the diffusion is inhibited at large values of kVRk, that is, at points where there is a rapid transition of the texture characteristics of the image. Interference filtering followed by ADF (ISFAD). Cardoso et al. (2012) proposed the ISFAD technique, which performs the speckle reduction in two stages: (i) interference-based speckle filtering (ISF) is applied to the original image to produce an image with well-defined edges and minimized speckle, and (ii) CADF is applied to the output image of ISF. Specifically, ISF is divided into three steps: median filtering with a circular window, destructive interference suppression (pointwise maximum between the original image and median-filtered image) and constructive interference suppression (median filtering with a small circular window). Image data sets The despeckling methods are evaluated by using synthetic and real BUS images. The simulations are used to compare the filtered edge map with groundtruth information to assess their similarity objectively. On the other hand, real images are employed to determine the effectiveness of the filter in preserving breast lesion boundaries determined manually by radiologists. Both types of images are detailed next. Synthetic images. The simulation process consists of two basic steps: (i) generating BUS image characteristics (geometries and breast tissue echogenicities), and (ii) simulating the physical process of image acquisition based on assumptions regarding the imaging device. The hypothetical ultrasound equipment has the following characteristics: transducer frequency of 7.5 MHz, beam resolution of 0.5 mm, spatial resolution of 0.33 mm/pixel, 8-bit grayscale images. Six typical tumor shapes are simulated: round, oval, lobulated, polygonal, angular and spiculated (Chou et al. 2001). Also, the sizes of the lesion shapes are defined between 1.0 and 3.0 cm (for the largest diameter). Any selected lesion shape is randomly placed on a basis background, which includes common structures found in breast tissue (Fig. 2a). Then, the echogenicity map is created by setting specific values of the ultrasonic attenuation coefficient and speed of sound for each different structure. Table 1 summarizes the simulated acoustic properties obtained from the literature (Madsen et al. 2006a, 2006b). Also, three kinds of lesion histology are mimicked: cyst, fibroadenoma and carcinoma. A backscattered signal is generated as Bðx; yÞ 5 Sðx; yÞh1 ðx; yÞh2 ðx; yÞ, where the operator * denotes spatial convolution; Sðx; yÞ represents the acoustic impedance inhomogeneities in the tissue map where each simulated tissue has a specific attenuation; h1 ðxÞ 5 sinðk0 xÞexp½2x2 =ð2s2x Þ is the transmitting pulse,

6

Ultrasound in Medicine and Biology

Volume -, Number -, 2014

Fig. 2. (a) Breast tissue map used to simulate data set SD-1. Common structures that are found in the female breast, as well as typical lesions, are considered. In this particular example, an angular lesion is shown. The variables w and h denote the lesion’s width and height, respectively. (b) Resultant simulated breast ultrasound image and its region of interest (ROI) (white dashed line).

and h2 ðyÞ 5 sinc2 ðyÞ is the beam aperture (Anderson and Trahey 2006). The constant k0 5 2pf0 =c is the wavenumber; c is the speed of sound in tissue; f0 is the transducer frequency; and sx is the pulse width of the transmitting acoustic wave. For simulating the speckle artifact, the model Sðx; yÞ 5 Tðx; yÞ$Gðx; yÞ is employed, where the centered dot denotes a pointwise product; Tðx; yÞ is the echogenicity model, that is, the breast tissue map with acoustic properties; and Gðx; yÞ is a Gaussian white noise field with zero mean and some variance (Yu and Acton 2002). Next, the echo magnitude is calculated as

b yÞ , where Bðx; b yÞ is the Hilbert Aðx; yÞ 5 Bðx; yÞ1j Bðx; transform of Bðx; yÞ with respect to x. Finally, Aðx; yÞ is log-compressed and scaled to the gray-level range 0– 255 (Yu and Acton 2002), and the resultant image is processed with the proposed despeckling approach. For each combination of shape and histology, 50 different synthetic images are built to increase versatility in terms of lesion arbitrary location and speckle generaTable 1. Simulated breast tissues and their acoustic properties Properties Simulated tissue

Attenuation coefficient, a (dB cm–1 MHz–1)

Speed of sound, c (m s–1)

Skin Subcutaneous fat Muscle Retromammary fat Milk duct Cooper’s ligaments Glandular parenchyma Cyst Fibroadenoma Cancer

0.60 0.47 0.50 0.38 0.10 1.20 0.70 0.11 0.29 0.39

1570 1456 1550 1483 1568 1570 1568 1572 1556 1570

tion. Hence, 900 images constitute the synthetic data set. Also, a region of interest (ROI) is automatically cropped around the lesion (dashed square in Fig. 2b) to reduce the filtering computation time. In the following, this simulated data set will be referred to as SD-1. Additionally, a second simulated data set is built, where round lesion-like objects of different diameters are placed on a homogeneous background. Also, the relative object/background contrast is controlled within the range –15 (high) to –2 (low) dB, where the attenuation coefficient of the background is fixed at 0.5 dB cm–1 MHz–1 and the speed of sound at 1540 m s–1. Fifty simulated images are generated for every contrast value. This synthetic data set will be denoted as SD-2. Real images. Fifty BUS images were acquired during routine breast diagnostic procedures at the National Cancer Institute (INCa) of Rio de Janeiro, Brazil. The INCa Research Ethics Committee approved this study (Protocol 38/2001). Patients were informed of the purpose of the study prior to consenting to participate. These BUS images were obtained with Sonoline Sienna (Siemens, Erlangen, Germany) equipment, using a 7.5-MHz linear array B-mode 40-mm ultrasound transducer, with an image resolution of about 0.33 mm/pixel. The images were captured directly from the 8-bit video signal (i.e., 256 gray levels) and saved in the TIFF format. Note that those ultrasonic scanner characteristics are used for BUS simulations. Two senior radiologists cropped each BUS image containing a single breast lesion to determine the respective ROIs. Also, every lesion was manually delineated to quantitatively evaluate the edge preservation of the aforementioned ADF-based approaches in relation to specialists’ perception.

Breast US despeckling d W. G OMEZ FLORES et al.

Edge preservation indices In the context of CAD systems for BUS images, it is important to determine the capability of the despeckling technique to preserve edge details among breast tissues, for further extraction of useful shape and contour features for tumor segmentation and classification. Thus, the aforementioned ADF-based approaches are evaluated in terms of edge preservation, by quantifying the agreement between both the filtered and the reference edge maps. For the synthetic data sets, we use Pratt’s figure of merit (FOM), which is defined as (Yu and Acton 2002) FOM 5

b N X 1 1  2 max Nb; Nideal i 5 1 11d$Di 

(13)

where Nb and Nideal are the numbers of filtered and ideal edge pixels, respectively; Di is the Euclidean distance between the ith filtered edge pixel and the nearest ideal edge pixel; and d is a constant typically set at 1/9. The FOM index ranges between 0 and 1, where unity indicates perfect edge preservation. The Canny edge detector (s 5 0:1) is used to obtain the binary edge maps from the filtered and the simulated ground-truth images (Canny 1986; Yu and Acton 2002). Note that all the detected edge pixels are involved in the calculation of the FOM value. On the other hand, in the absence of a unique absolute ground-truth edge map, edge preservation could be assessed by using several human interpretations of an image (Unnikrishnan et al. 2007). To perform such comparison in real sonographies, we propose the mean radial distance (MRD) index, which gives the average distance in the radial direction, from the lesion’s centroid, between the filtered edge map (by Canny’s detector) and the radiologist’s lesion delineation. Let M 5 {m1, m2, ., m p} denote the set of p points on the manually delineated lesion and C 5 {c1, c2, ., cq} denote the edge map with q points of the filtered image. Then, the distance from mi to the contour C in the same radial direction from the lesion centroid is defined as b i 2b c k jjÞ; cmi ˛M (14) dðmi ; C Þ 5 minðjjmi 2ck jj$jj m k

b i and cbk are the where k$k is the Euclidean distance; and m unit vectors directed radially outward. Hence, from the contour C , the subset S 5 {s1, s2, ., sp} with p points is extracted, that is, the set of closest points to M that are located in the same radial direction. Then, the MRD is computed as MRD 5

p 1X jjmi 2si jj p i51

(15)

7

If S presents more discrepancies in relation to M , then the value of MRD tends to increase. Therefore, for acceptable edge preservation performance, MRD should tend toward zero pixels. Statistical analysis Two robust statistics for estimating location (sample median, MD) and dispersion (Qn estimator) of FOM and MRD results are used. These statistics are capable of coping with outliers and non-normal distributions. As is well known, the sample median is the middle value in an ordered sequence of data values, that is, the 50th percentile. On the other hand, the Qn estimator computes all possible distances between pairs of observations, sorts them and chooses the smallest distance, corresponding to the first quartile of all interpoint distances (Rousseeuw and Croux 1993). Additionally, statistical significance is analyzed by using the Kruskal–Wallis test (a 5 0.05) to evaluate whether the medians of the approaches compared differ under the assumption that the shapes of the underlying distributions are the same. Also, the correction for multiple testing on the basis of the same data is made by the Bonferroni method. RESULTS Determining the ADLG configuration Because the performance of the ADLG approach depends on kdi;j in the conduction coefficient, it is worth evaluating the effect of different numbers of scales (Sn ) and orientations (qn ) in Log-Gabor filters. Both values are free parameters, which are set before applying the ADLG technique to BUS images. Also, for more accurate localization of texture boundaries, Jain and Farrokhnia (1991) proposed smoothing the texture information with a small window for homogenizing similar textural regions. Herein, the input BUS image is first smoothed with a 9 3 9 hybrid median window (Davies 2004). A search grid is used to test several values of Sn and qn to obtain an appropriate FOM value. A combination of five scales, Sn 5 f2; 3; 4; 5; 6g, five orientations, qn 5 {2, 4, 8, 12, 24}, and two octaves, sq/su 5 0.55, is carried out in our experiments. Hence, 25 distinct ADLG configurations are applied to the synthetic data set SD-1. These configurations provide even and thorough coverage of the whole spectrum (Shan et al. 2012). Also, the number of iterations of the anisotropic diffusion process is set to tmax 5 1000 and the time step t 5 1/4. For the sake of clarity, we use the notation Sn/ qn to indicate a specific ADLG configuration. The robust statistics, MDyQn, are computed by using the maximum FOM value determined over all the iterations of each simulated image in the data set.

8

Ultrasound in Medicine and Biology

Volume -, Number -, 2014

Fig. 3. Responses of ADLG approach through iterations t, where gt is the conduction coefficient (eqn 9) and I t is the despeckled image (eqn 12); the squares contain the corresponding figure of merit results. Also, the input and ideal images are shown. ADLG 5 anisotropic diffusion guided by Log–Gabor filters.

The experimental results point out that the configuration 3/24 reaches the best performance with 0:83y0:03. In addition, the Kruskal2Wallis test, with Bonferroni correction, performed 300 pairwise comparisons, that is, all combinations of 25 configurations are taken two at a time. This statistical analysis revealed that the configurations 4/24, 5/24 and 6/24 did not statistically differ from the best one, 3/24; that is, their edge preservation performances are quite similar. For the remaining experiments, we used the ADLG configuration 3/24, whose set of Log–Gabor functions comprises 72 filters, whereas the configurations with similar performance consider larger banks. Figure 3 illustrates some outcomes (I t ) of the proposed approach (configuration 3/24) applied to a synthetic image through 1000 iterations and the respective FOM values. Also, the corresponding directional conduction coefficient responses (gt ) are shown, where

dark pixels indicate locations that should not be filtered, and bright pixels point indicate regions that should be blurred. It is noticeable that the edges between different textures are preserved in relation to the ideal image.

Stop criterion of ADLG approach Until now, the experimental results considered the maximum FOM value reached after tmax iterations (i.e., the best performance for a single image), because the ground-truth information is available. However, in real data, it is impossible to determine a target FOM value to stop the diffusion process of the ADLG approach. In the SRAD technique, Yu and Acton (2002) proposed measurement of the mean square error (MSE) of images between two iterations as MSE 5

Fig. 4. Joint frequency distribution between mean square error (MSE) and maximum figure of merit (FOM) values.

M 21 X N 21 h i2 1 X Ii;jt 2Ii;jt2Dt MN i 5 0 j 5 0

(16)

where M and N are the image width and height, respectively; t is the current iteration; and Dt denotes the iteration step. If MSE is smaller than 0.01, the diffusion process is stopped automatically. However, the authors did not justify this stop criterion analytically or experimentally. In this sense, we associate the MSE values with maximum FOM results by using the joint frequency of both indices to determine which MSE range appears recurrently. This experiment considered the synthetic data set SD-1 and Dt 5 10. In Figure 4, it is noticeable that the largest peak is located in the MSE range 0.120.5 for FOM values between 0.8 and 0.9, which represents the 70% of all occurrences. Also, if the whole range of FOM values (0.620.9) is considered, 90% of occurrences fall in the MSE range 0.120.5. Thus, on the basis of these experimental results,

Breast US despeckling d W. G OMEZ FLORES et al.

Table 2. Anisotropic diffusion-based filters compared in this study and their respective initial parameters Technique

Parameter k 5 1, t 5 1/8 Coordinates of homogeneous region, t 5 1/8 k 5 5, S 5 3, s0 5 3, N 5 4, k0 5 0:05, t 5 1/8 k 5 1, wradii 5 10 and 3, t 5 1/8 Sn 5 3, qn 5 24, S 5 3, lmin 5 1, 2 octaves, t5¼

CADF SRAD TOAD ISFAD ADLG

CADF 5 conventional anisotropic diffusion filtering; SRAD 5 speckle-reducing anisotropic diffusion; TOAD 5 texture-oriented anisotropic diffusion; ISFAD 5 interference-based speckle filtering followed by anisotropic diffusion; ADLG 5 anisotropic diffusion guided by Log–Gabor filters.

we can automatically stop the diffusion process of the ADLG approach when MSE,0:5. Comparison of ADF-based approaches The performance of the despeckling approaches is compared by using both the synthetic data set SD-1 and the real BUS data set. Table 2 summarizes the set of ADF-based techniques including the parameters used herein (Aleman-Flores et al. 2007; Cardoso et al. 2012; Shan et al. 2012). For the experiments with synthetic data, the number of iterations is set at 1000, and the FOM index assess the edge preservation performance. Table 3 summarizes the statistical results for each ADF-based technique, which are computed by using the maximum FOM values determined over all the iterations, that is, the best result for a particular image. Note that the ADLG approach reaches the best edge preservation with a median FOM value equal to 0.83. In addition, the Kruskal2Wallis test, with Bonferroni correction, confirmed that there exist statistically significant differences regarding all pairwise comparisons that involved the ADLG approach. Table 3. Figure of merit results for each anisotropic diffusion-based approach with respect to the synthetic data set* Technique

MD

Qn

Max

Min

ADLG ISFAD TOAD CADFy SRADy

0.83 0.59 0.51 0.40 0.39

0.04 0.04 0.06 0.04 0.06

0.92 0.72 0.68 0.52 0.72

0.66 0.47 0.33 0.28 0.27

CADF 5 conventional anisotropic diffusion filtering; SRAD 5 speckle-reducing anisotropic diffusion; TOAD 5 texture-oriented anisotropic diffusion; ISFAD 5 interference-based speckle filtering followed by anisotropic diffusion; ADLG 5 anisotropic diffusion guided by Log–Gabor filters. * The techniques are ordered from best to worst performance. y Non-statistically significant differences.

9

Figure 5 illustrates an example of a simulated carcinoma with angular lesion shape filtered by each of the ADF-based approaches evaluated. Also, the FOM results associated with each technique are illustrated; the ADLG method yielded the largest value. In the experiments with the real BUS data set, the parameters of the despeckling approaches are set as indicated in Table 1. However, the number of iterations for each technique is fixed as Cardoso et al. (2012) defined: CADF, 600; SRAD, 400; TOAD, 600; ISFAD, 500. In the case of ADLG, we used the diffusion stop criterion defined at MSE,0:5. Also, the MRD index evaluates lesion boundary preservation. Table 4 outlines the statistical results for each ADFbased technique, with respect to both radiologists’ delineations, where the ADLG technique reached the lowest median MRD value with 4.19 pixels. Additionally, the Kruskal2Wallis test, with Bonferroni correction, determined that all the pairwise comparisons that involved the ADLG approach are statistically significant; that is, our proposed technique performed best. Figure 6 illustrates an example of a real BUS image with distinct contrasts. Note that for the original image, the edge map corresponds to the radiologists’ manual delineation of the lesion. Also, the MRD results associated with each approach are illustrated, where the ADLG method yielded the smallest value. Edge preservation with contrast variations The previous experiments indicated that the ADLG approach is capable of preserving adequately the tissue edges better than the other four ADF-based methods. However, it is important to evaluate the behavior of the ADLG technique in distinct scenarios with controlled contrast variations. Also, the performance of the four aforementioned ADF-based methods is compared with that of our technique in terms of the FOM index. In this experiment, the synthetic data set SD-2 is considered. Figure 7 illustrates the experimental results, where each point in the graph represents the median of the maximum FOM value (after tmax 5 1000) calculated from the entire set of simulated images. These results point out that ADLG approach performed best for distinct contrasts in terms of edge preservation. Additionally, Figure 8 illustrates the ADLG behavior for different contrast values, where the left side corresponds to the filtered image, and the right side, to the original image. One can note that the proposed approach is capable of preserving the lesion-like regions relatively well separated from the background. However, small structures could remain on the object’s contour as the contrast becomes lower, although the general shape is preserved.

10

Ultrasound in Medicine and Biology

Volume -, Number -, 2014

Fig. 5. Simulated breast ultrasound image filtered by anisotropic diffusion-based approaches. Top row: Outputs of evaluated filters. Bottom row: Edge maps obtained with Canny’s detector. The squares contain the corresponding figure of merit results. ADLG 5 anisotropic diffusion guided by Log–Gabor filters; CADF 5 conventional anisotropic diffusion filtering; SRAD 5 speckle-reducing anisotropic diffusion; TOAD 5 texture-oriented anisotropic diffusion; ISFAD 5 interference-based speckle filtering followed by anisotropic diffusion.

DISCUSSION In BUS images, speckle not only reduces image quality, but also results in contrast variations between different tissues. Hence, a filtering approach should consider reducing efficiently the speckle while preserving breast lesion shape characteristics. To accomplish this goal, we proposed the ADLG technique, which adapts the conduction coefficient parameter, k, for each image pixel by using local texture features. In this sense, 2-D Gabor filters have been widely used to depict texture features, as their parameters (such as frequency, orientation, and bandwidth) can be easily tuned. Thus, if we assume that distinct BUS textures represent different tissues, then Gabor functions could be used to distinguish them. Aleman-Flores et al. (2007) built a bank of Gabor filters by setting two free parameters: the scale factor, s, and the spatial frequency, k. Table 4. Mean radial distance (in pixels) results for each anisotropic diffusion-based approach with respect to the real data set* Technique

MD

Qn

Max

Min

ADLG ISFADy CADFy SRADy TOADy

4.19 5.88 6.29 6.39 6.43

1.51 1.51 1.45 1.84 1.81

8.76 10.83 12.36 23.31 14.82

1.48 2.68 3.20 2.02 2.79

CADF 5 conventional anisotropic diffusion filtering; SRAD 5 speckle-reducing anisotropic diffusion; TOAD 5 texture-oriented anisotropic diffusion; ISFAD 5 interference-based speckle filtering followed by anisotropic diffusion; ADLG 5 anisotropic diffusion guided by Log–Gabor filters. * The techniques are ordered from best to worst performance. y Non-statistically significant differences.

However, both parameters are not independent, and tuning them separately could be inadequate for the following reasons. First, the wavelength, l, of the cosine term determines the preferred spatial frequency, that is, k 5 1/l, and it depends on image width (cycles/image width) (Jain and Farrokhnia 1991), which varies between BUS images. Second, s is related to l by the ratio s/l 5 0.56, when the half-response spatial frequency bandwidth is approximately one octave (Kruizinga and Petkov 1999). In addition, the Gabor filters present a major drawback: the filter averaging is not null, and consequently, the DC component influences intermediate bands. Hence, we used the Log–Gabor functions to overcome this bandwidth limitation (Field 1987). The ADLG approach extracts textures features by using Log–Gabor filters to adapt the k parameter. The experimental results, using synthetic data, indicated that the bank of Log–Gabor filters that attained the best edge preservation had 24 orientations and three scales, with median FOM 5 0.83. Contrarily, the filter configuration with two orientations and two scales yielded the worst performance, with median FOM 5 0.55. Therefore, it is obvious that the greater the number of filter orientations, the better is the edge preservation. The simulated data set was also used to compare the performance of the ADLG approach with that of four ADF-based techniques commonly used for speckle reduction. The statistical results pointed out that our approach yielded the largest FOM value. On the other hand, the lowest performance was attained by the CADF and SRAD techniques, both with median FOM z 0.40. These findings suggest that it is convenient to adapt the k parameter locally for each pixel to reach adequate edge preservation while speckle is reduced

Breast US despeckling d W. G OMEZ FLORES et al.

11

Fig. 6. Real breast ultrasound image filtered by anisotropic diffusion-based approaches. Top row: Outcomes of evaluated filters. Middle row: edge maps obtained with Canny’s detector. Bottom row: black lines correspond to radiologist’s contour, M , and gray points are the nearest filtered edge pixels, C , at the same radial direction from the lesion centroid. The squares contain the corresponding mean radial distance results. ADLG 5 anisotropic diffusion guided by Log–Gabor filters; CADF 5 conventional anisotropic diffusion filtering; SRAD 5 speckle-reducing anisotropic diffusion; TOAD 5 texture-oriented anisotropic diffusion; ISFAD 5 interference-based speckle filtering followed by anisotropic diffusion.

efficiently in homogeneous textured regions. Additionally, a simulated experiment with controlled contrast variations revealed that the ADLG technique outperformed the four aforementioned ADF-based techniques. The

Fig. 7. Edge preservation performance (in terms of figure of merit [FOM] index) of distinct anisotropic diffusion-based approaches as a function of relative contrast variations. ADLG 5 anisotropic diffusion guided by Log–Gabor filters; CADF 5 conventional anisotropic diffusion filtering; SRAD 5 speckle-reducing anisotropic diffusion; TOAD 5 texture-oriented anisotropic diffusion; ISFAD 5 interferencebased speckle filtering followed by anisotropic diffusion.

results suggested that our approach is capable of maintaining acceptable edge preservation with distinct contrast values. Real BUS images were also included in this study. The image data set consisted of 50 ultrasonographies; each one contained a single breast lesion delineated manually by two senior radiologists. The goal was to assess the performance of the above-mentioned ADFbased approaches in preserving lesion contour details. We proposed measuring the MRD between the specialist lesion delineation and the filtered edge map generated by Canny’s edge detector. The results indicated that the ADLG approach performed best, with median MRD 5 4.19 pixels. On the contrary, the lowest performance was that of the TOAD approach, with median MRD 5 6.43 pixels. Also, the statistical analysis determined that all pairwise comparisons that involved the ADLG approach revealed significant differences in terms of the MRD index. Hence, our approach performs better than the other ADF-based approaches tested within a 95% confidence interval. We defined a stop criterion for automatically halting the anisotropic diffusion of the ADLG approach by using the synthetic data set. We registered the value of the MSE reached at the same iteration of the maximum FOM value for each filtered image. It was observed that 90% of the filtered images reached the maximum FOM in the MSE range 0.1–0.5. Hence, on the basis of these results, we proposed stopping the ADLG algorithm when

12

Ultrasound in Medicine and Biology

Volume -, Number -, 2014

Fig. 8. Behavior of ADLG approach with different contrast values in simulated images from data set SD-2. The left side represents the filtered speckled image, and the right side, the original input image. ADLG 5 anisotropic diffusion guided by Log–Gabor filters.

MSE,0:5. Note that this stop criterion was used in evaluation of the ADLG approach with real images, and the results reveal that our approach performs best in preserving breast lesion edges. CONCLUSIONS Encouraged by the results of this study, we consider that the ADLG approach (three scales and 24 orientations for Log–Gabor filters) could potentially be used in CAD systems to enhance procedures such as breast lesion segmentation. Both, the synthetic data set and the MATLAB (The MathWorks, Natick, MA, USA) implementation of the ADLG approach are available from the authors on request. REFERENCES  Aleman-Flores M, Alvarez L, Caselles V. Texture-oriented anisotropic filtering and geodesic active contours in breast tumor ultrasound segmentation. J Math Imaging Vis 2007;28:81–97. Alvarenga AV, Infantosi AFC, Pereira WCA, Azevedo CM. Assessing the performance of morphological parameters in distinguishing breast tumors on ultrasound images. Med Eng Phys 2009;32:49–56. Anderson ME, Trahey GE. A seminar on k-space applied to medical ultrasound. Department of Biomedical Engineering, Duke University; 2006. Canny J. A computational approach to edge detection. IEEE Trans Pattern Anal Mach Intell. PAMI-8 1986;679–698. Cardoso FM, Matsumoto MMS, Furuie SS. Edge-preserving speckle texture removal by interference-based speckle filtering followed by anisotropic diffusion. Ultrasound Med Biol 2012;38:1414–1428.

Cheng HD, Shan J, Ju W, Guo Y, Zhang L. Automated breast cancer detection and classification using ultrasound images: A survey. Pattern Recogn 2010;43:299–317. Chou YH, Tiu CM, Hung GS, Wu SC, Chang TY, Chiang HK. Stepwise logistic regression analysis of tumor contour features for breast ultrasound diagnosis. Ultrasound Med Biol 2001;27:1493–1498. Corsetti V, Houssami N, Ghirardi M, Ferrari A, Speziani M, Bellarosa S, Remida G, Gasparotti C, Galligioni E, Ciatto S. Evidence of the effect of adjunct ultrasound screening in women with mammographynegative dense breasts: Interval breast cancers at 1 year follow-up. Eur J Cancer 2011;47:1021–1026. Czerwinski RN, Jones DL, O’Brien WD. Detection of lines and boundaries in speckle images—Application to medical ultrasound. IEEE Trans Med Imaging 1999;18:126–136. Dantas RG, Costa ET. Ultrasound speckle reduction using modified Gabor filters. IEEE Trans Ultrason Ferroelectr Freq Control 2007;54: 530–538. Davies ER. Computer and machine vision: Theory, algorithms, practicalities. Waltham, MA: Academic Press; 2004. Evans AN, Nixon MS. Speckle filtering in ultrasound images for feature extraction. In: Proceedings International Conference on Acoustic Sensing and Imaging, London, UK, 29–30 March 1993:44–49. Feldman MK, Katyal S, Blackwood MS. US Artifacts. Radiographics 2009;29:1179–1189. Field DJ. Relations between the statistics of natural images and the response properties of cortical cells. J Opt Soc Am A 1987;4: 2379–2394. Horsch K, Giger ML, Venta LA, Vyborny CJ. Automatic segmentation of breast lesions on ultrasound. Med Phys 2001;28:1652–1659. Infantosi AFC, Luz LMS, Pereira WCA, Alvarenga AV. Breast ultrasound segmentation using morphologic operators and a Gaussian function constraint. In: Proceedings, 14th Nordic-Baltic Conference on Biomedical Engineering and Medical Physics. Berlin/New York: Springer; 2008. p. 520–523. Jackson VP. The role of US in breast imaging. Radiology 1990;177: 305–311. Jain AK, Farrokhnia F. Unsupervised texture segmentation using Gabor filters. Pattern Recogn 1991;24:1167–1186.

OMEZ FLORES et al. Breast US despeckling d W. G

Kovesi P. Symmetry and asymmetry from local phase. In: Proceedings, 10th Australian Joint Conference on Artificial Intelligence. Cambridge, MA: MIT Press; 1997. p. 185–190. Kruizinga P, Petkov N. Nonlinear operator for oriented texture. IEEE Trans Image Process 1999;8:1395–1407. Lopes A, Touzi R, Nezry E. Adaptive speckle filters and scene heterogeneity. IEEE Trans Geosci Remote Sens 1990;28:992–1000. Loupas T, McDicken WN, Allan PL. An adaptive weighted median filter for speckle suppression in medical ultrasonic images. IEEE Trans Circuits Syst 1989;36:129–135. Madsen EL, Berg WA, Mendelson EB, Frank GR. Anthropomorphic breast phantoms for qualification of investigators for ACRIN Protocol 6666. Radiology 2006a;239:869–874. Madsen EL, Hobson MA, Frank GR, Shi H, Jiang J, Hall TJ, Varghese T, Boley MM, Weaver JB. Antropomorphic breast phantoms for testing elastography systems. Ultrasound Med Biol 2006b;32:857–874. Michailovich OV, Tannenbaum A. Despeckling of medical ultrasound images. IEEE Trans Ultrason Ferroelectr Freq Control 2006;53: 64–78. Perona P, Malik J. Scale-space and edge detection using anisotropic diffusion. IEEE Trans Pattern Anal Mach Intell 1990;35:629–639. Petrou M, Garcıa-Sevilla P. Image processing. Chichester: Wiley; 2006.

13

Rahbar G, Sie AC, Hansen GC, Prince JS, Melany ML, Reynolds HE, Jackson VP, Sayre JW, Bassett LW. Benign versus malignant solid breast masses: US differentiation. Radiology 1999;213:888–894. Rousseeuw PJ, Croux C. Alternatives to the median absolute deviation. J Am Stat Assoc 1993;88:1273–1283. Shan J, Cheng HD, Wang Y. Completely automated segmentation approach for breast ultrasound images using multiple-domain features. Ultrasound Med Biol 2012;38:262–275. Thijssen JM. Ultrasonic speckle formation, analysis and processing applied to tissue characterization. Pattern Recogn Lett 2003;24: 659–675. Unnikrishnan R, Pantofaru C, Hebert M. Toward objective evaluation of image segmentation algorithms. IEEE Trans Pattern Anal Mach Intell 2007;29:929–944. Wang W, Li J, Huang F, Feng H. Design and implementation of LogGabor filter in fingerprint image enhancement. Pattern Recogn Lett 2008;29:301–308. Weickert J. Anisotropic diffusion in image processing. 2nd ed. Stuttgart: Teubner; 2008. Yang Z, Fox MD. Speckle reduction and structure enhancement by multichannel median boosted anisotropic diffusion. EURASIP J Adv Signal Process 2004;2004:2492–2502. Yu Y, Acton ST. Speckle reducing anisotropic diffusion. IEEE Trans Image Process 2002;11:1260–1270.