Abstract

Automatically extracting breast tumor boundaries in ultrasound images is a difficult task due to the speckle noise, the low image contrast, the variance in shapes, and the local changes of image intensity. In this paper, an improved edge-based active contour model in a variational level set formulation is proposed for semi-automatically capturing ultrasonic breast tumor boundaries. First, we apply the phase asymmetry approach to enhance the edges, and then we define a new edge stopping function, which can increase the robustness to the intensity inhomogeneities. To extend the capture range of the method and provide good convergence to boundary concavities, we use the phase information to obtain an improved edge map, which can be used to calculate the gradient vector flow (GVF). Combining the edge stopping term and the improved GVF in the level set framework, the proposed method can robustly cope with noise, and it can extract the low contrast and/or concave boundaries well. Experiments on breast ultrasound images show that the proposed method outperforms the state-of-art methods.

1. Introduction

Breast cancer is one of the major causes of death among women [1]. The early detection is necessary and helpful to the diagnosis and treatment of breast cancer [2]. Both X-ray mammography and ultrasonic examination can be used for this purpose. The advantages of ultrasonic examination over mammography are its inexpensive, noninvasive, and cost-effect nature. Segmentation of breast ultrasound (US) images is essential for quantitative measurement of tumor characteristics. Manual segmentation is time consuming and it increases inter- and intraobserver variability. Fully automatic segmentation of US image is still a challenging task for two main reasons [3, 4]. First, US images have speckle noise, low contrast, and local changes of intensity [3]. Second, breast tumors have a great variation of their shapes, sizes and locations [4]. Particularly, when the tumors are malignant, they generally have acute concave shapes. Therefore, developing a semiautomatic algorithm for segmentation of breast tumor in US images becomes urgent in clinical application.

Active contour models (also called snakes) have been proved to be very useful for image segmentation [519]. The desirable advantage of this semiautomatic method lies in its ability to facilitate the segmentation process with a priori knowledge and user interaction. Active contour models can be divided into two major classes: parametric [510] and geometric [1119].

Parametric active contour models use parametric equation to explicitly represent the evolving curve. The problems of the traditional parametric active contour model are its limited capture range and its inability to handle concave boundaries. Xu and Prince [6, 7] proposed a gradient vector flow (GVF) field model to increase capture range and encourage convergence to boundary concavities. GVF snake has received a lot of attention in US image segmentation [810]. For example, Rodtook and Makhanov [8] developed an improved GVF approach based on the continuous force field analysis. They used the resulting vector field to generate a modified edge map and then used it to calculate a novel gradient-based source term. A comparison is made between their method and the conventional GVF type methods by using the reference segmentation results outlined by expert physicians. In their experiments, they claimed that their method can successfully solve the active contour convergence problem. They also obtained more accurate tumor segmentation result.

Geometric active contour models implicitly represent evolving curve as the zero-level set of a higher dimensional function, called level set function. Commonly, this function is defined by computing the closest distances between pixels and a given closed curve in an image domain. The points that have positive distances are inside the curve, and ones that have negative distances are outside the curve [11]. The curve is evolved using the partial differential equation [14] derived from the energy function. When the evolving curve maintains a steady status, that is, the curve stops evolving, the segmentation result is obtained as the zero-level set of the function. The partial differential equation [14], which sometimes can be derived from the energy function [13], describes a curve smoothing process. When the curve evolution stops, the zero-level set of the energy function refers to the segmentation result. Compared to the parametric models, level set methods can represent contour of complex topology and it can also deal naturally with topological changes by using the evolution of the level set function rather than the parametric curves. These properties lend themselves well in medical image segmentation where the anatomic structures split and merge through the pixels. One can see that this is especially useful when segmenting the breast tumors, which generally have complex and varying shapes. The level set method was first introduced by Osher and Sethian [11] and successfully applied to image segmentation by Malladi et al. [12] and Caselles et al. [13]. Caselles et al. developed the well-known geodesic active contour (GAC) model in [13]. However, in these methods, reinitialization of the level set function to a signed distance function is often required. Re-initialization is performed by periodically solving a partial differential equation, which is a computationally expensive operation. To address this problem, Li et al. [15] presented a variational level set method to eliminate the need for reinitialization by incorporating a penalty term into the energy functional. Unfortunately, this penalty term may cause an undesirable side effect on the level set function in some circumstances, which may affect the numerical accuracy. Recently, Li et al. [16] developed a more general variational level set formulation by incorporating a distance regularization term into the energy functional. Particularly, they investigated a double-well potential related to the distance regularization term to maintain a desired shape of the level set function. Accordingly, the undesirable side effect produced by the penalty term can be completely avoided. They have demonstrated the capability of segmentation on artificial, natural and real medical images.

Using an additional edge stopping function (ESF) allows the active contour to stop at the edges [1116]. In these models, ESF is related to a convolution of the intensity image with a Gaussian filter. However, the traditional ESF has at least two disadvantages. The first one is that the ESF is never equal to zero, and boundary leakage may occur, that is, the active contour may cross the object’s boundaries. The second one is that the edge detection based on gradient is not appropriate for speckle noise, although it is generally effective for additive noise. In breast US images corrupted by multiplicative noise, using only the gradient information may detect false contours caused by noise, and consequently too many trivial contours may be obtained. An alternative approach to circumvent these problems is to use a stopping term based on the coefficient of variation [17]. This stopping term equals to zero at the edges and accommodates the speckle noise. An important advantage of using this stopping term is the capability of preventing the active contour from crossing the object’s boundaries. Alternatively, Chan and Vese [18] proposed a method called C-V model, where their stopping term is based on a particular image segmentation result obtained by Mumford–Shah technique [19] rather than the image gradients. This segmentation model uses the intensity information of different region to minimize the energy function rather than using the edge information. Liu et al. [4] also proposed a model that considers the differences between the actual probability densities of the intensities in different regions and the corresponding estimated probability densities. Their model shows advantages over the C-V model [18] on breast US images with weak boundaries and low contrast.

In this paper, a novel segmentation method of ultrasonic breast tumor is proposed within the level set framework. First, we present a novel ESF, which is independent of the intensity gradient of image, as most of the active contour models do, but instead is related to local phase information obtained by using the phase asymmetry approach. The use of this ESF is more robust than the traditional gradient magnitude-based ESF because local phase is theoretically invariant to intensity magnitude. Subsequently, in order to increase the capture range of the method and its ability to move into acute concave boundaries, the improved gradient vector flow (GVF) field is adopted as external force and it is added to the proposed model. The results show that the proposed method can obtain better segmentation performance in comparison with some state-of-art methods [11, 16].

The rest of this paper is organized as follows. In Section 2, we first describe the traditional level set methods and the distance regularized level set method. Section 3 details the proposed method. Section 4 shows qualitative and quantitative experimental results on clinical breast US images. Section 5 summarizes our work.

2. Background

2.1. Traditional Level Set Methods

The basic idea of the level set methods is to implicitly represent a closed curve 𝐶(𝑡) at time 𝑡 by the zero-level set of function 𝜙, that is, 𝐶(𝑡)={(𝑥,𝑦)𝜙(𝑥,𝑦,𝑡)=0},with𝜙(𝑥,𝑦,0)=𝐶0,(2.1) where 𝐶0 is the initial curve.

In the level set framework, the general evolution equation of the level set function 𝜙 is 𝜕𝜙||||𝜕𝑡+𝐹𝜙=0,(2.2) where 𝐹 is the speed function. For image segmentation, 𝐹 is the external force which is dependent on some image information, such as image gradient and image intensity. This external force 𝐹 attracts the curves toward the edges. 𝐹 is commonly defined as follows: 𝐹=𝐹image+𝐹geometry.(2.3) The first term 𝐹image is the image force generated from image information. The second term 𝐹geometry is the geometry force to smooth the curve.

2.2. Distance Regularized Level Set Method

In most of the traditional level set methods, reinitialization is used as a numerical remedy for keeping sound curve evolution and guaranteeing reliable results. The main problem with this procedure is its high computational cost. To avoid this problem, Li et al. [15] proposed a level set method without reinitialization by integrating a penalty term into the energy functional. However, this penalty term leads to undesirable side effect. More recently, by adding a distance regularization term into energy functional, Li et al. [16] presented a novel distance regularized level set evolution method, that is, DRLSE model. This model can not only completely eliminate the need for reinitialization but also avoid the undesirable side effect.

Let Ω𝑅2 be a bounded Lipschitz image domain. The energy functional 𝐸() is written as 𝐸(𝜙)=𝜇𝑅𝑝(𝜙)+𝐸ext(𝜙),(2.4) where 𝜙 is a level set function, 𝜇>0 is a constant, and 𝑅𝑝(𝜙) is the distance regularization term defined by 𝑅𝑝(𝜙)=Ω𝑝||||𝜙𝑑𝑥,(2.5) where 𝑝 is the potential function for 𝑅𝑝(𝜙) [16], that is, 1𝑝(𝑠)=(2𝜋)21(1cos(2𝜋𝑠))if𝑠1,2(𝑠1)2if𝑠1.(2.6) In (2.4), 𝐸ext(𝜙) is the external energy, which is borrowed from the GAC model [13], that is, 𝐸ext(𝜙)=𝜆𝐿𝑔(𝜙)+𝛼𝐴𝑔(𝜙)=𝜆Ω||||𝑔𝛿(𝜙)𝜙𝑑𝑥+𝛼Ω𝑔𝐻(𝜙)𝑑𝑥,(2.7) where 𝜆 and 𝛼 are constants, 𝐿𝑔(𝜙) is the length of the interface, 𝐴𝑔(𝜙) is the area of the region that the zero-level curve encloses, 𝛿() is the Dirac delta function, and 𝐻() is the Heaviside function. 𝑔 is an ESF such that lim𝑡𝑔(𝑡)=0, 1𝑔=||1+𝐺𝜎||𝐼2,(2.8) where 𝐺𝜎 is a Gaussian function with the standard deviation 𝜎, 𝐼 is a given 2D image, symbol represents convolution, is the gradient operator, and || is the modulus of the smoothed image gradients.

As can be seen, the level set function 𝜙 will stop moving when 𝑔 approaches zero. In other words, the resulting final contour will halt at where the image gray gradient is large enough to make 𝑔0. A small value of 𝜎 may make the proposed method sensitive to noise, and the evolution will not be stable. On the other hand, a large value of 𝜎 may lead to boundary leakage, and the extracted boundary may not be accurate. The choice of 𝜎 is a dilemma, especially when the images have speckle noise and weak boundaries. This property makes the traditional ESF in (2.8) inappropriate for US images.

The associated Euler-Lagrange equations, obtained by minimizing function equation (2.4) with respect to 𝜙, is 𝜕𝜙𝑑𝜕𝑡=𝜇div𝑝||||𝑔𝜙𝜙+𝜆𝛿(𝜙)div𝜙||||𝜙+𝑎𝑔𝛿𝜀(𝜙),𝜙(𝑥,𝑦,0)=𝜙0(𝑥,𝑦),(2.9) where div() is the divergence operator and 𝑑𝑝 is a function given by 𝑑𝑝𝑝(𝑠)=(𝑠)𝑠,(2.10) that satisfies |𝑑𝑝(𝑠)|<1, for all 𝑠(𝑜,) and lim𝑠0𝑑𝑝(𝑠)=lim𝑠𝑑𝑝(𝑠)=1.

In (2.9), the initial level set function 𝜙0 being a binary step function is given by 𝜙0(𝑥,𝑦)=𝑑if(𝑥,𝑦)Ω0,𝑑,otherwise,(2.11) where 𝑑>0 is a constant and Ω0 is a subset in the image domain Ω.

In practice, the Dirac delta function 𝛿() in (2.9) is replaced by the following function 𝛿𝜀():  𝛿𝜀1(𝑥)=2𝜀1+cos𝜋𝑥𝜀,|𝑥|𝜀,0,|𝑥|>𝜀,(2.12) where 𝜀 is a positive constant. As shown in Figure 1, if 𝜀 is too small, the values of 𝛿𝜀() are prone to approximately zero to make its effective range small, hence the energy functional can be easily trapped into a local minimum. On the contrary, if 𝜀 is large, the final contour location may be inaccurate, even if 𝛿𝜀() tends to achieve a global minimum. For all experiments in this paper, the parameter 𝜀 is set to 1.5.

3. Proposed Method

The goal of our method is to extract a breast tumor from a given US image. This method consists of three stages: ESF definition, improved GVF construction, and model generation. Figure 2 shows the block diagram of the proposed method. In the first stage, we first apply phase asymmetry approach to enhance the image edges and then define a novel phase-based ESF. At the second stage, we generate an edge map and use it to compute the GVF. At the last stage, we construct the model first by using a phase-based ESF and then by incorporating the improved GVF.

3.1. Design of the ESF

The sine and cosine functions with specific amplitudes (energies) can represent any discrete signal. In the time domain, these functions cause a set of scaled waves, synthesizing the original signal. Morrone and Owens [20] gave a general definition of how to describe and detect signal features using the phase of the Fourier components.

Kovesi [21] proposed a low-level feature detector, named phase congruency model [20], in terms of the Fourier analysis. The phase congruency model is based on the assumption that features are perceived at points in an image where the Fourier components are maximally in phase, rather than assume a feature is a point of maximal intensity gradient. In other words, the phase of the sinusoidal components varies and a low phase congruency is achieved at all other nonfeature points. In practical applications, the local frequency, particularly, the phase information, is estimated via a bank of Cauchy filter tuned to various spatial frequencies (scales) rather than Log-Gabor filter. The main reason for this is that Log-Gabor filters are likely not a good choice in the event of feature detection [22, 23].

The measure of phase congruency provides inaccurate localization, noise sensitivity, and high computational cost. Motivated by the properties of the phase congruency, Kovesi considered that symmetries and asymmetries cause special phase patterns in the values of the image intensity [21]. At symmetric points, all Fourier series should have either maximum or minimum value. Kovesi proposed a modified measure to detect symmetry and asymmetry. The symmetric measure is [𝑛]=𝑃𝑆𝑠𝐴𝑠[𝑛]||𝜙cos𝑠[𝑛]||||𝜙sin𝑠[𝑛]||𝑇𝑠𝑠𝐴𝑠[𝑛]=+𝜀𝑠||even𝑠[𝑛]||||odd𝑠[𝑛]||𝑇𝑠𝑠𝐴𝑠[𝑛]+𝜀,(3.1) and the asymmetric one is[𝑛]=𝑃𝐴𝑠||odd𝑠[𝑛]||||even𝑠[𝑛]||𝑇𝑠𝑠𝐴𝑠[𝑛]+𝜀,(3.2) where even𝑠 and odd𝑠 are the responses of the even-symmetric (cosine) and odd-symmetric (sine) filters at scale 𝑠. At a point of symmetry, the absolute value of even𝑠 is large while that of odd𝑠 is small. At a point of asymmetry, the opposite is true. 𝐴𝑠[𝑛] represents the magnitude of the responses of each quadrature pair of filters at a given scale 𝑠 and point 𝑛 in the image, that is, 𝐴𝑠[𝑛]=(odd𝑠[𝑛])2+(even𝑠[𝑛])2. Symbol denotes zeroing of negative values, 𝜀 is a positive constant to avoid division by zero, and 𝑇𝑠 is the scale specific noise threshold.

Extension to two-dimensional (2D) images uses the 1D analysis over several orientations and then combines the results to offer a single measure of the feature significance [24]. For a given breast US image, we first convolve a breast US image with a bank of Cauchy kernel filter in order to calculate the 2D phase asymmetry. In this paper, we yield the asymmetric measures over six orientations and four/five scales at each pixel. The values of 𝑃𝐴 range from 0 to 1, close to one in homogeneous region and close to zero in the vicinity of boundaries. Clearly, 𝑃𝐴 is a dimensionless quantity.

Based on the analysis mentioned above, one can see that 𝑃𝐴[0,1] is the feature asymmetry measure defined by (3.2) and can offer a better control on the edge detection quality. The main difference between this edge detection and the traditional edge detections is that it is invariant to the changes in image intensity or contrast. According to our knowledge of the breast US images, they are often intensity inhomogeneous and have high amounts of speckle noise. Thus, we apply 𝑃𝐴 to define ESF with sufficient robustness and accuracy in this paper. Since ESF is also an inverse edge indicator function, an appropriate definition of ESF is given by 1𝑔=1+(𝑃𝐴)𝛼,(3.3) where 𝛼[0,1] is a hyperparameter. The 𝑔 takes values in [0,1], close to zero in the vicinity of boundaries and close to one in homogeneous regions. Particular advantage of this ESF over the gradient-based ESF (e.g., (2.8)) is its insensitivity to intensity inhomogeneity and speckle noise. Figure 3 illustrates this statement.

As can be seen from Figure 3, two situations are described. The first case, where there exits a single homogeneous region, is shown in Figure 3(a), and the second case, where two distinct homogeneous regions are separated by a sharp edge, is shown in Figure 3(d). The gradient-based ESF and phase-based ESF are calculated on these two speckled images. Figures 3(b) and 3(e) show that the variance of the gradient-based ESF increases with the intensity in homogeneous regions. The multiplicative nature of speckle noise leads to this expected behavior. In contrast, the phase-based ESF has more homogeneous values, as shown in Figures 3(c) and 3(f). Consequently, the detected edges obtained using the phase-based ESF are better than those obtained using the gradient-based one. The very low values of ESF easily make curve stop at spurious edges which are caused by the speckle noise. Aiming at precise segmentation of US images, our proposed model prefers to use the phase-based ESF than the conventional gradient-based one.

3.2. Improved GVF Construction

Xu and Prince [6, 7] introduced the GVF as a bidirectional external force to capture the object’s boundary from either sides and to deal with concave boundary. The construction of the GVF can be divided into two stages: the generation of edge map 𝑓 from the image and the computation of the GVF from the edge map 𝑓. Clearly, the quality of edge map 𝑓 is a crucial factor in the construction of GVF field. A classical choice is defined as ||𝑓=𝐺𝜎||𝐼2.(3.4)

An obvious problem with this edge map is that the traditional GVF could be attracted to strong edges which are caused by noises in an inhomogeneous region with high gray values.

In order to eliminate this shortcoming, we define an edge map as 𝑓(𝑥,𝑦)=𝑔(𝑥,𝑦),(3.5) where 𝑔 is given in (3.3).

The GVF field is defined by the vector field 𝑉GGVF(𝑥,𝑦)=(𝑢(𝑥,𝑦),𝑣(𝑥,𝑦)),(3.6) that minimizes the following objective function [6, 7]: 𝐸GGVF(1𝑢,𝑣)=21||||𝑢𝑓2𝑥+𝑢2𝑦+𝑣2𝑥+𝑣2𝑦+2||||||(||𝑓𝑢,𝑣)𝑓2𝑑𝑥𝑑𝑦,(3.7) where 𝑢𝑥,𝑢𝑦,𝑣𝑥,𝑣𝑦 are the first derivatives of the vector field, 1() and 2() are the weight functions depending on the gradient of the edge map 𝑓, and they are defined as 1||||||||𝑓=exp𝑓𝐾2||||𝑓=11||||,𝑓,(3.8) where 𝐾 is a nonnegative parameter which is used to control the smoothness of the resulting vector field 𝑉GGVF. 1() and 2() should be monotonically decreasing and increasing functions of |𝑓|, respectively.

The force vector field 𝑉GGVF(𝑥,𝑦) can be solved from the following Euler equations: 12𝑢2𝑢𝑓𝑥=0,12𝑣2𝑣𝑓𝑦=0,(3.9) where 𝑓𝑥 and 𝑓𝑦 are the derivatives of 𝑓 with respect to 𝑥 and 𝑦.

For an example of breast US image, Figures 4(b) and 4(c) show the detailed original GVF field and improved GVF field of the window marked in Figure 4(a). It can be seen that more gradient vectors point to the true boundary obviously in the improved GVF field, with compared to the original GVF field.

The major difference between the traditional GVF and the improved GVF is that the former edge map is dependent on gradient magnitude while the latter is dependent on phase congruency. The advantage of the latter over the former is that it makes use of the phase information in the cure convolution and, thus, achieves better edge detection results for objects with intensity inhomogeneities or speckle noise.

3.3. Model Generation

Although most of the traditional level set methods require reinitialization, Li’s method [16] does not require reinitialization and provides stable level set evolution. Therefore, based on the model in (2.9) presented by Li et al. [16], we propose an improved edge-based level set model in this paper. The main difference of the proposed model from the DRLSE model in (2.9) is that it uses the phased-based ESF previously constructed and simultaneously integrates the improved GVF.

The proposed model yields a level set evolution, called phase- and GVF-based level set evolution (PGBLSE). The PGBLSE model can be described by the following differential equations: 𝜕𝜙𝑑𝜕𝑡=𝜇div𝑝||||𝜙𝜙+𝑎𝑔𝛿𝜀(𝜙)+𝜆𝛿𝜀𝑔(𝜙)div𝜙||||[̂𝑣]𝜙̂𝑢,,(3.10) where 𝜇>0, 𝑎𝑅, 𝜆𝑅 are the constants as the weights of the distance regularization term, the area term, and the length term, 𝑑𝑝() is a function defined in (2.10), 𝛿𝜀() is the Dirac delta function in (2.12), 𝑔 is the ESF defined in (3.3), 𝜙 is the level set function, t represents time variable, div() is the divergence operator, and is the gradient operator. ̂[̂𝑢,𝑣] is the normalized GVF field where ̂𝑢=𝑢/𝑢2+𝑣2 and ̂𝑣=𝑣/𝑢2+𝑣2. The field normalization is responsible for obviating the dependency on the absolute value of the image intensity and providing sufficient information to guide the curve evolution.

Now assuming that 𝜙 is negative inside the zero-level set and positive outside, the inward normal vector to the curve can be expressed as 𝑁=𝜙/|𝜙|. The curve motion (shrinkage/expansion) is governed by the two vectors. This curve is expanded if the curve’s normal aligns with the normalized GVF, or it is shrunk if the curve’s normal and the normalized GVF are in the opposite directions and is not deformed if the curve’s normal is orthogonal to the normalized GVF. Obviously, this flow increases the shrinking and expanding capability of curve evolution. This property is useful in concave boundaries.

4. Results and Discussions

4.1. Data Acquisition

This work was approved by Shanghai Sixth People’s Hospital in China. Our method was developed using MATLAB and evaluated using 20 US images of breast tumors. Among these images, 10 were malignant tumors and 10 were benign tumors. The patients’ ages are in the range from 18 and 75 years old. The images were collected by using three kinds of scanners with the linear transducers (7–14 MHz). The ultrasonic subimage of region of interest (ROI) was chosen by an expert radiologist. This area only was then used for image segmentation. An example of an extracted ROI subimage is illustrated in Figure 5.

4.2. Quantitative Measure

To evaluate the accuracy of the segmentation algorithms, we compared them with the gold standard, which was manually produced by an experienced radiologist. In addition to the golden standard segmentation, an error measure is required for evaluating object segmentation algorithms. Measures can be divided into region-based measures and boundary-based ones.

Region-based error measures are made by comparing the regions inside the contours. The overlapping area error [4] is represented by the three parameters including true positive (TP), false positive (FP), and Jaccard similarity (JS). For a given image, let 𝑆𝑎 and 𝑆𝑚 be the obtained segmentation and golden standard, respectively. These three metrics are defined as ||𝑆TP=𝑚𝑆𝑎||||𝑆𝑚||,||𝑆FP=𝑚𝑆𝑎𝑆𝑚||||𝑆𝑚||,||𝑆JS=𝑚𝑆𝑎||||𝑆𝑚𝑆𝑎||,(4.1) where is the set union operator, is the set intersection operator, and || represents the number of pixels in the corresponding pixel set. The higher the TP value is, the more actual tumor regions the segmented tumor regions cover. The lower the FP value is, the fewer normal tissue regions the segmented tumor regions cover. The JS value ranges from 0 to 1. The closer the JS value is to 1, the better the segmentation is.

The Dice similarity coefficient [25] has also been widely used for comparing existing methods: DSC=2JS=2||𝑆1+JS𝑚𝑆𝑎||||𝑆𝑚||+||𝑆𝑎||.(4.2) The closer the value of DSC is to 1, the better the segmentation is.

Boundary-based error measures evaluate segmentations based on the distance between two contours. The boundary error is defined as the average or max distance between two contours. We denote the golden standard boundary by 𝐴={𝑎1,,𝑎𝑚} and the obtained boundary as 𝐵={𝑏1,,𝑏𝑛}, where each element of 𝐴 or 𝐵 is a point on the corresponding contour. For every point 𝑎𝑖 on boundary 𝐴, a distance to boundary 𝐵, denoted as MD(𝑎𝑖,𝐵), can be defined as the minimum distance of point 𝑎𝑖 to all points in 𝐵: 𝑎MD𝑖,𝐵=min𝑗{1,,𝑛}𝑎𝑖𝑏𝑗,(4.3) where |||| denotes the Euclidean distance.

The average minimum Euclidean distance (AMED) [26] is the average distance between two contours: AMED(𝐴,𝐵)=𝑚𝑖=1𝑎MD𝑖,𝐵+2𝑚𝑛𝑗=1𝑏MD𝑗,𝐴2𝑛.(4.4)

The Hausdorff distance (HD) [26], on the other hand, measures the max distance between two contours: HD(𝐴,𝐵)=maxmax𝑖{1,,𝑚}𝑎MD𝑖,𝐵,max𝑗{1,,𝑛}𝑏MD𝑗,𝐴.(4.5) For a perfect segmentation, the distances HD and AMED are equal to zero.

4.3. Comparison with Previous Methods

We compare our method, PGBLSE model (3.10), with two edge-based active contour models: GAC model [11] and DRLSE model [16]. The level set evolution in GAC model is 𝜕𝜙||||𝜕𝑡=𝑔𝜙div𝜙||||||||𝜙+𝑎𝑔𝜙+𝑔𝜙,(4.6) where 𝑎 is a constant, which plays the similar role as in the PGBLSE model (3.10) and 𝑔 is the ESF defined by (2.8).

The same initial contours are used in these methods. We use the following default setting of the parameters in our method: 𝜇=0.2, 𝑎=1.1, 𝜆=0.2. The parameters 𝛼 and 𝑠 are not set to the same values in all experiments. If we have to detect tumors of weak boundaries or concave shapes, then parameters 𝛼 and 𝑠 should be small. Conversely, if we have to detect tumors of strong boundaries and regular shapes, then parameters 𝛼 and 𝑠 should be large. For the GAC and DRLSE models, we tuned the parameters for the best segmentation results for all images. Experiments  1–5 illustrate the comparison results on five US images of breast tumors.

Experiment  1 applies these three models to an US image which contains a malignant breast tumor with obvious intensity inhomogeneity and highly concave boundary as shown in Figure 6(a). The golden standard image was shown in Figure 6(b). The results obtained with the GAC [11], DRLSE [16], and the proposed PGBLSE models are illustrated in Figures 6(c), 6(d) and 6(e), respectively. The GAC model fails in concave region since there is no force that can pull the contour towards the concave portion of the tumor boundary. Although the DRLSE model produces a smoother contour, it also fails in a similar fashion. Obviously, without gradient diffusion, the ability to capture concave boundary is limited. With the help of diffusion process, the proposed PGBLS method (𝛼=0.1,𝑠=4) successfully extracts the concave part of the tumor. In addition, gradient diffusion is capable of removing the boundary effect of weak edges. This is the reason why some parts of the contours from the GAC and DRLSE models leak past weak boundary gradients while the PGBLSE model does not. This also supports the efficiency of the improved GVF.

In experiment  2, these three models are applied to an US image of the malignant tumor in Figure 7(a). This image is extremely challenging for edge-based active contour methods since it contains a complex-shaped tumor. Besides, the boundaries between the tumor region and the surrounding normal tissue region are very weak, and the contrast is also low. The results of the GAC model [11], DRLSE model [16], and the proposed method are illustrated in Figures 7(c), 7(d), and 7(e), respectively. The GAC and DRLSE models capture the entire tumor when parts of the normal tissues are also identified as tumor. Since tumor region is surrounded by normal tissue regions with similar intensities and has local changes of intensity, the gradient-based stopping terms of GAC and DRLSE models are heavily affected and can be easily trapped into unsuitable local minima. In contrast, the proposed PGBLSE model (𝛼=0.1,𝑠=4) uses a novel phase-based edge term; therefore, it can handle intensity inhomogeneity well. In addition, due to the efficiency of the improved GVF mentioned before (the ability to converge the concave boundaries is increased with sufficient gradient diffusion), the PGBLSE model can prevent the leakage through weak edges and extract the concave boundaries of tumor.

In experiment  3, Figure 8 illustrates the results using the GAC [11], DRLSE [16], and PGBLSE models in the segmentation of malignant breast tumor in an US image. The GAC model fails to segment the tumor in that the surrounding normal tissues of the tumor are included as shown in Figure 8(c). It is also difficult to use the DRLSE model to successfully extract the tumor; Figure 8(d) illustrates that some parts of the contour leak past weak boundary gradients while other parts are confined inside the tumor. The PGBLSE model (𝛼=0.5,𝑠=5) separates the real tumor from normal tissue regions, as illustrated in Figure 8(e). That is, a successful segmentation is obtained.

Experiment  4 applies these three methods on an US image of a benign breast tumor, as shown in Figure 9(a). It can be seen that the tumor image has gradually changing intensity and the boundaries between the tumor region and its surrounding normal tissue are very weak. The GAC model cannot deal with weak boundary, and the resulting contour leaks at where there are relatively weak boundary gradients, as depicted in Figure 9(c). The DRLSE model also cannot handle weak boundaries well, and the contour leaks into the surrounding normal tissues, as shown in Figure 9(d). In contrast, the proposed model (𝛼=0.1,𝑠=4) uses the edge map based on phase information to calculate the gradients, which helps to prevent the contour from leaking into normal tissues. Therefore, the proposed model is more immune to leaking effect and the generated tumor contour is very close to the golden standard, as shown in Figure 9(e).

Experiment  5 illustrates the application of the proposed method to a relatively smooth benign breast tumor as shown in Figure 10(a). We analyze the GAC and DRLSE models. The US image suffers from speckle noise. The classical GAC model cannot yield smooth contour even for the round shape, as shown in Figure 10(c). In addition, the GAC contour is confined inside by small clusters of noise with large intensity magnitudes. Even though the DRLSE model generates smoother contour, the resulting contour is still confined inside the tumor region as illustrated in Figure 10(d). On the contrary, Figure 10(e) shows that the proposed PGBLSE model (𝛼=0.5,𝑠=5) is the most accurate among these three methods.

Clearly, the GAC and DRLSE models are sensitive to the noise which seriously affects the gradient-based stopping terms. Therefore, these two models yield many comparatively strong gradients throughout the whole image including the homogeneous regions, which distracts the evolution contours from the real boundaries. This problem can be solved by using the phase-based stopping term because this term is theoretically intensity invariant.

Experiments  1–5 illustrate the comparison results on five US images of breast tumors in comparison with the GAC model [11], the DRLSE model [16], and the golden standard. The results suggest that the tumor contours extracted from our proposed PGBLSE model are better than those from the other two models [11, 16], and they are very close to the golden standard. This conclusion is also supported by the details shown in Figure 11.

Figure 11 shows that the proposed PGBLSE model obtains better accuracy than the GAC and DRLSE models in terms of the JS and DSC measures. Moreover, the PGBLSE model is superior to the GAC and DRLSE models in terms of the AMED and HD measures. Particularly, the PGBLSE model is almost 2 times better than the GAC model in terms of the HD measure.

Figures 610 demonstrate that the proposed model is the most accurate comparing with the other models in the study. Finally, we evaluated our approach on the dataset of 20 US images of breast tumors. For these studied images, Figure 12 shows six box plot graphics using the TP, FP, JS, DSC, AMED, and HD measures for the three segmentation models.

From Figure 12(a), a phenomenon is observed that the GAC model and the DRLSE model performed slightly better than the proposed PGBLSE model in terms of the TP measure. This advantage has been visualized in the box plot showing the large median values of the TP. However, their segmentation results are not accurate, and the higher values of FP also imply this. As shown in Figure 12(b), the median FP values of the GAC and DRLSE models are very high (17.74% and 19.74%, resp.). On the contrary, the proposed PGBLSE model obtains a much lower median value of FP (4.19%), which demonstrates that there are fewer normal tissue regions mistakenly included in tumor regions and thus the segmentation is more reliable. Like the cases in Figures 7 and 9, although the GAC and DRLSE models approximately include the entire real tumor regions, they simultaneously cover a lot of unsatisfied normal tissue regions. Figure 12(c) illustrates that the median values of JS for these models are, respectively, 79.73%, 80.08%, and 86.30%. It is observed that the proposed PGBLSE model presents the best median value, while the GAC and DRLSE models have lower median values, which suggests worst segmentation performance. The GAC model (from 81.63% to 85.1%) and the DRLSE model (from 81.51% to 85.95%) have the upper quartile in a lower position when compared to the position of the same quartile of the proposed PGBLSE model (from 87.74% to 90.04%), which suggests a lower performance of the GAC and DRLSE models. This is also verified by the DSC measure (Figure 12(d)).

As can be seen from Figure 12(f), the median value of HD for our PGBLSE model is around 7 pixels, whereas the median values for the GAC and DRLSE models are all above 16 pixels. High HD value means high disagreement between the segmented contour and the golden standard. Besides, it is observed that our model offers a smaller dispersion range while the GAC and DRLSE models provide large dispersion ranges. Although the median values of AD for the three models are not very different from each other, the GAC and DRLSE models have more outliers and larger dispersion range of values as illustrated in Figure 12(e).

As can be seen from Figure 13, the results are then assessed on 20 clinical US images of breast tumors using the tumor areas that tumor contour enclosed. They are displayed in Bland-Altman [27] plots. The 𝑥-axis is golden standard area and the 𝑦-axis represents the tumor area difference between the three segmentation models and golden standard. For the 20 breast US images, using the GAC and DRLSE models, 95% of the point pairs (19 images) are with ±2SD (95% confidence intervals). Using the proposed PGBLSE model, 100% of the point pairs (20 images) are with ±2SD (100% confidence intervals). The mean differences are 568, 369, and −202.4 pixels for the GAC, DRLSE, and proposed PGBLSE models, respectively.

The goal of this particular experiment is twofold. On the one hand, we demonstrate that tumor areas of the GAC and DRLSE models are larger than those of the proposed PGBLSE model, and, therefore, large values of TP are obtained. On the other hand, some nontumor regions are also incorrectly covered by the GAC and DRLSE models, and, therefore, large values of FP are achieved.

In our experiments, all figures (i.e., Figures 613) show, among the three models, that the best agreement is between the proposed PGBLSE model and the golden standard.

5. Conclusions

We have presented in this paper a new approach for the segmentation of the ultrasonic breast tumors. For the first time, phase asymmetry approach, which can enhance the boundaries, is used to segment the ultrasonic breast tumors. In a level set framework, we integrate the use of a novel ESF and the improved GVF, both constructed using the output of phase asymmetry. This model shows significant improvements, particularly, in robustness against the speckle noise, as well as in handling intensity inhomogeneities and capturing concave boundaries.

The performance of the proposed PGBLSE model is demonstrated on clinical US images of breast tumors. Qualitative and quantitative results show that the PGBLSE model outperforms the classical intensity-based GAC and DRLSE models. The further work will focus on developing a hybrid level set active contour model by investigating the addition of a region-based term, in order to improve the performance of the method.

Acknowledgment

This work was supported by National Basic Research Program of China, the 973 Program (Grant no. 2010CB732501).