Semi-automatic Breast Ultrasound Image Segmentation Based on Mean Shift and Graph Cuts

Abstract

Computerized tumor segmentation on breast ultrasound (BUS) images remains a challenging task. In this paper, we proposed a new method for semi-automatic tumor segmentation on BUS images using Gaussian filtering, histogram equalization, mean shift, and graph cuts. The only interaction required was to select two diagonal points to determine a region of interest (ROI) on an input image. The ROI image was shrunken by a factor of 2 using bicubic interpolation to reduce computation time. The shrunken image was smoothed by a Gaussian filter and then contrast-enhanced by histogram equalization. Next, the enhanced image was filtered by pyramid mean shift to improve homogeneity. The object and background seeds for graph cuts were automatically generated on the filtered image. Using these seeds, the filtered image was then segmented by graph cuts into a binary image containing the object and background. Finally, the binary image was expanded by a factor of 2 using bicubic interpolation, and the expanded image was processed by morphological opening and closing to refine the tumor contour. The method was implemented with OpenCV 2.4.3 and Visual Studio 2010 and tested for 38 BUS images with benign tumors and 31 BUS images with malignant tumors from different ultrasound scanners. Experimental results showed that our method had a true positive rate (TP) of 91.7%, a false positive (FP) rate of 11.9%, and a similarity (SI) rate of 85.6%. The mean run time on Intel Core 2.66 GHz CPU and 4 GB RAM was 0.49 ± 0.36 s. The experimental results indicate that the proposed method may be useful in BUS image segmentation.

Keywords

breast ultrasound semi-automatic segmentation mean shift graph cuts OpenCV

Introduction

Breast cancer is the most commonly diagnosed cancer in women worldwide and is increasing particularly in developing countries where the majority of cases are diagnosed in late stages.¹ Currently, mammography is the most effective method for detecting and diagnosing breast tumors.² However, many unnecessary biopsies are performed because of the low specificity of mammography.³ Moreover, the ionizing radiation of mammography might cause health risks for patients and radiologists.² Because it is cost-effective, noninvasive, and performed in real time, ultrasound (US) has become a useful adjunct to mammography.⁴ However, certain inherent characteristics of breast ultrasound (BUS) images, including low contrast, speckle noise, and tissue-related textures,^9,15 may cause difficulties for radiologists. Moreover, inter- or intraobserver differences in the understanding and diagnosis of BUS images might occur. Hence, developing a computer-aided detection/diagnosis (CAD) system for analyzing BUS images is crucial.^5,24-27 The output of a CAD system can provide a second opinion for radiologists in the detection and diagnosis of breast tumors.⁵

Image segmentation plays an essential role in CAD systems. Many techniques have been developed for US image segmentation. An extensive review of US image segmentation methods was performed by Noble and Boukerroui,⁶ and a systematic review of BUS image segmentation and classification methods was conducted by Cheng et al.¹⁷ BUS image segmentation techniques include histogram thresholding,⁷ region growing,^7,8,14 model-based (active contour models or snakes,^8,9 level set,^10,28,29,32 Markov random fields^11,30), graph-based,^12,13 neural network,^7,31 and watershed.^15,16 Table 1 summarizes some selected BUS image segmentation methods in the following aspects: number of cases (images), accuracy, level of automation, and time efficiency.

Table 1.

Overview of Selected Breast Ultrasound Image Segmentation Methods.

Study	Year	Method	n	Accuracy	Auto	Time (s)	CPU (GHz)
Shan et al.⁷	2012	ADF + AT + RG + NN	120	TP = 92.8%, FP = 12.0%, SI = 83.1%	Auto	9.5	3.0
Madabhushi and Metaxas⁸	2003	Low-level + High-level + Empirical	42	NTPAO = 75.1%, AMBE = 6.6 p	Auto	22.1	3.0
Chen et al.⁹	2003	ACM	8	SI = 89.2%	Semi	302.8	0.866
Huang et al.¹⁰	2007	AT + LS	118	SI = 88.6%, OL = 89.3%, EF = 12.5%	Auto	—	—
Xiao et al.¹¹	2002	MAP + MRF	—	—	Auto	—	—
Huang et al.¹²	2012	ADF + Graph-based + ACM	20	ARE = 9.0	Semi	9.3	—
Gao et al.¹³	2012	HP-NCut	100	SI = 85.6%, HD = 8.9 p, MED = 1.8 p	Semi	191.8	2.0
Drukker et al.¹⁴	2002	RGI + ARG RG	757	TP = 90%, FP = 61%	Auto	60.0	0.4
Huang et al.¹⁵	2004	NN + WS	60	PR = 81.7%, MR = 94.7%	Auto	—	—
Gómez et al.¹⁶	2010	CLAHE + ADF + WS + ARD	50	OL = 86.0%, NRV = 0.16, PD = 6.6%	Semi	25.0	2.2
Shan et al.¹⁸	2012	ADF + PMO + NLM	122	TP = 92.4%, FP = 7.2%, SI = 86.3%, HD = 22.5 p, MAD = 4.8 p	Auto	9.8	3.0
Liu et al.²⁸	2010	LS	79	TP = 93.9%, FP = 7.0%, SI = 88.1%	Semi	60.4	3.0
Gao et al.²⁹	2012	LS	20	TP = 90.8%, FP = 4.2%, SI = 86.3%	Semi	—	—
Takemura et al.³⁰	2010	AB + MRF + ACM	400	JI = 93.4%	Auto	280.0	—
Jiang et al.³¹	2012	AB + K-means + SVM + RW	112	OL = 83.5%	Auto	—	—
Chang et al.³²	2005	ADF + Stick + AT + LS	210	—	Auto	32.6	—

n = number of images; p = pixels; Auto = automatic; Semi = semi-automatic; ADF = anisotropic diffusion filtering; AT = adaptive thresholding; RG = region growing; NN = neural network; TP = true positive; FP = false positive; SI = similarity; NTPAO = normalized true positive area overlap; AMBE = average mean boundary error; ACM = active contour model; LS = level set; OL = overlap; EF = extra fraction; MAP = maximum a posteriori; MRF = Markov random field; ARE = averaged radial error; HP-NCut = normalized cut based on homogeneous patches; HD = Hausdorff distance; MED = minimum Euclidean distance; RGI = radial gradient index; ARG = average radial gradient; WS = watershed; PR = precision ratio; MR = match rate; CLAHE = contrast-limited adaptive histogram equalization; ARD = average radial derivative; NRV = normalized residual value; PD = proportional distance; PMO = phase in max-energy orientation; NLM = neutrosophic l-means; MAD = mean absolute distance; JI = Jaccard index; AB = AdaBoost; SVM = support vector machine; RW = random walk.

Improving the accuracy, time efficiency, and level of automation is the research goal in the area of BUS image segmentation. However, these three aspects are always contradictory; for instance, refining the segmentation to obtain a higher accuracy generally requires a longer run time. Histogram thresholding and region growing are easy to implement and fast, but they cannot get precise segmentation due to the inherent characteristics of BUS images. Model-based methods can be used for high-level and refined segmentation, but they are sensitive to the initial contours and require an iterative process that is time consuming. Graph-based methods have been used in a variety of applications of computer vision. However, if applied directly to BUS images without some preprocessing steps, these methods may not produce satisfying results, because when the graph vertices are constructed at the image pixels and the segmentation energy is based on the pixel intensities, the speckle noises and low contrast of BUS images may bias the segmentation. Neural network-based approaches can yield reasonable tumor contours, but their training process is time consuming. Watershed-based methods may be used for fast segmentation, but they face the problem of over-segmentation and are sensitive to user-specified markers. Many automatic BUS image segmentation methods have been proposed.^{7,8,10,11,14,15,18,30-32} Automatic segmentation methods can reduce intra- and interoperator variance. In general, these techniques have two steps: (a) region of interest (ROI) detection, for instance, using machine learning techniques and (b) ROI segmentation, that is, extracting the tumor contour from the ROI determined in Step 1. This makes the segmentation accuracy reliant on the sensitivity and specificity of ROI detection. Moreover, the time efficiency of existing automatic segmentation methods still needs to be improved (Table 1). In addition, few of current BUS image segmentation techniques have been validated for the images collected from different US scanners or manufacturers.

To obtain desired accuracy and time efficiency as well as level of automation in BUS image segmentation is still challenging due to the low-contrast, speckle noise, and tissue-related textures of BUS images. This inspired us to develop a method that is practical and easy to implement. A large number of computer vision techniques have been employed in medical image processing and analysis.^37,38 Mean shift and graph cuts are two techniques that may be used for fast segmentation. The former was proposed by Comaniciu and Meer,²⁰ and the latter by Boykov and Jolly.²¹ We considered that it is possible to employ a hybrid method combining mean shift and graph cuts to obtain desired accuracy and computational complexity. Tao et al. proposed a hybrid method for color image segmentation, which incorporates the advantages of the mean shift segmentation and the normalized cut partitioning methods and requires low computational complexity and is, therefore, very feasible for real-time image segmentation processing.³³ A hybrid method incorporating the mean shift filtering and the normalized cuts segmentation methods was employed by Zhang et al. for automated biometric measurement of early gestational sac in US images.³⁴ Ye et al. proposed a hybrid method incorporating the mean shift superpixels and the graph cuts segmentation methods for automatic segmentation of lesions in computer tomography (CT) images.³⁵ Recently, Yun et al. proposed a hybrid method combining the mean shift segmentation and graph cuts segmentation methods for vascular US image segmentation.³⁶ These hybrid methods have demonstrated their benefits to some extent. However, because of the low contrast and speckle noises of BUS images, a hybrid method combining mean shift and graph cuts may not produce desired BUS image segmentation accuracy, so we considered incorporating filtering and contrast enhancement prior to the hybrid method. Besides, the graph cuts segmentation always requires a manual selection of object and background seeds, and we considered automatic seed generation. Notice that the implementation language also affects the run time of a BUS image segmentation method. OpenCV is an open source computer vision library and has been widely used worldwide. Open source software can be very helpful. Ince et al. argued that anything less than the release of source programs is intolerable for results that depend on computation.³⁹ We considered using an OpenCV implementation with the C++ language.

In this paper, we propose a hybrid method for semi-automatic BUS image segmentation combining Gaussian filtering, histogram equalization, pyramid mean shift,²⁰ and graph cuts.²¹ The flowchart of the proposed method is shown in Figure 1. The method was implemented with OpenCV 2.4.3 and Visual Studio 2010 using C++. To evaluate the effectiveness of this method, 69 BUS images collected from different US scanners were tested. There were 38 benign cases and 31 malignant cases. In the next sections, we describe the materials and methods and present experimental results that allow a discussion of the proposed method.

Figure 1.

Flow chart of the proposed method.

Materials and Method

Data Acquisition

This study was approved by the Institutional Review Board of West China Hospital of Sichuan University and the patients signed informed consent forms. BUS images were collected with three commercial US scanners: a HDI 3500 (Philips Medical Systems, Bothell, Washington) and a HDI 5000 (Philips Medical Systems, Bothell, Washington), both with the 5-12 MHz linear probe L12-5, and a LOGIQ 500 PRO (GE Healthcare, Madison, Wisconsin) with the 6-13 MHz linear probe LA39. A sonographer performed the US scanning, and breast tumors were identified as benign or malignant according to biopsy reports. There were 38 benign and 31 malignant tumors from 69 female patients. The image size was 768 × 576 pixels. The size of the tumors ranged from 7.8 to 47.9 mm and the median size was 20.6 mm. For each BUS image, the manual contouring of the tumor by an experienced radiologist is taken as the reference standard.

Image Preprocessing

The proposed method requires selecting two diagonal points to determine a ROI on an input BUS image I (Figure 1). The preprocessing of the ROI image I_ROI includes the following three steps:

Shrink I_ROI by a factor of 2 using bicubic interpolation to reduce computation time. The shrunken image $I_{ROI}^{1 \ 2}$ was obtained. The bicubic interpolation was conducted by fitting a cubic spline between the 4 × 4 surrounding pixels in I_ROI and then reading off the corresponding destination value from the fitted spline.¹⁹

Smooth $I_{ROI}^{1 \ 2}$ by a Gaussian filter to reduce high-frequency noises. The smoothed image I_G was obtained. Gaussian filtering was done by convolving each pixel in $I_{ROI}^{1 \ 2}$ with a Gaussian kernel and then summing to produce the output image I_G. The scale of the Gaussian kernel was 3 × 3 pixels.

Enhance I_G by histogram equalization to improve the contrast of BUS images. The enhanced image I_HE was obtained. Histogram equalization was performed by mapping one distribution (the given histogram of intensity values in I_G) to another distribution (a wider and, ideally, uniform distribution of intensity values in I_HE). The classical histogram equalization method was used.⁴⁴

Pyramid Mean Shift Filtering

After the contrast-enhanced image I_HE was obtained, we improved its homogeneity by pyramid mean shift filtering. Pyramid mean shift filtering is based on mean shift clustering over color²⁰ and consists of the following four steps:

Decompose I_HE into L layers (L ≥ 2) by the Gaussian pyramid: I₁, . . ., I_L, using

I_{k + 1} = DOWN (I_{k} \otimes G) .

where the operator DOWN(.) downsizes the original image by a factor of 2 by removing every even-numbered row and column, the ⊗ symbol denotes convolution, G is a 5 × 5 Gaussian kernel, and 2 ≤ k ≤ L.

Conduct mean shift filtering for I_L to obtain $I_{L}^{m s}$ .

Reconstruct an upsampled image $I_{L - 1}^{u}$ from $I_{L}^{m s}$ using the Gaussian pyramid:

I_{L - 1}^{u} UP (I_{L}^{ms}) \otimes G,

where the operator UP(.) upsizes the original image by a factor of 2 with the new (even) rows and columns filled with 0’s and maps each pixel in location (x, y) in the original image to pixel (2x + 1, 2y + 1) in the destination image. Then, conduct mean shift filtering for $I_{L - 1}^{u}$ to get $I_{L - 1}^{ms}$ .

Repeat Step 3 until getting $I_{L}^{ms}$ , which is the pyramid mean shift filtered image I_f.

Mean shift filtering consists of two steps: discontinuity preserving smoothing and mean shift clustering. Let ${x_{j}} j = 1, . . ., n$ be the set of pixels of the original image and ${z_{j}} j = 1, . . ., n$ be the resulting convergence point of each pixel. And let s and r denote the spatial and range parts of the vectors, respectively. The discontinuity preserving smoothing procedure is as follows²⁰:

For each j = 1, . . ., n,

initialize k = 1 and y_k = x_j;

compute $y_{k + 1} = \frac{1}{n_{k}} \sum x_{i} \in S (y_{k})$ , $k \leftarrow k + 1$ till convergence; and

assign $z_{j} = (x_{j}^{s}, y_{conv}^{r})$

The last assignment specifies that the filtered data at the spatial location of x_j will be assigned to the range components of the point of convergence y_conv. The window S(y_k) is centered on y_k and has n_k points. The original data are normalized with σ _s for the spatial part and σ _r for the range (color) part.

Let ${L_{j}} j = 1, . . ., n$ be a set of labels (scalars). The mean shift clustering procedure is as follows²⁰:

For each j = 1, . . ., n, run the discontinuity preserving smoothing procedure for x_j and store the convergence point in z_j.

Indentify clusters ${C_{p}} p = 1, . . ., m$ of convergence points by linking together all z_j, which are closer than 0.5 from each other in the joint domain. The number of clusters m is controlled by (σ _s , σ _r ).

For each j = 1, . . ., n, assign $L_{j} = {p | z_{j} \in C_{p}}$ .

It should be noted that when L = 1, pyramid mean shift turns into conventional mean shift. We used pyramid mean shift rather than mean shift because the former can better improve the homogeneity of BUS images and suppress the speckle noise and tissue-related textures, as shown in Figure 2. Figure 2(a) is the preprocessed image same as shown in Figure 1; Figure 2(b) is the pyramid mean shift filtered image, using the parameters listed in Table 3; and Figure 2(c) is the mean shift filtered image, using the EDISON code^40,41 (Edge Detection and Image SegmentatiON, http://coewww.rutgers.edu/riul/research/code/EDISON/) and the default parameters provided by EDISON: spatial bandwidth = 7, color bandwidth = 6.5, and minimum region = 20.

Figure 2.

(a) Preprocessed image. (b) Pyramid mean shift filtered image. (c) Mean shift filtered image.

Graph Cuts Segmentation

The pyramid mean shift filtered image I_f was segmented by interactive graph cuts to obtain a binary image I_b. The interactive graph cuts approach²¹ generally requires manual selection of seed pixels in the object region O and the background region B using a few strokes of an image brush. To reduce operator-dependence and speed up BUS image segmentation, the object and background seeds were automatically generated on I_f (Figure 1). Let w_I_f and h_I_f denote the width and height of I_f . The crosshair as the object seeds is centered at $(w_{i_{f}} / 2, h_{i_{f}} / 2)$ and has a width of lh wIf and a height of lh hIf, where λ_w and λ_h are two coefficients controlling the weights of the width and height of the crosshair and 0 < λ_w < 1, 0 < λ_h < 1, The background seeds include eight lines: four lines near to the corners and four lines near to the edges. The end points of the four lines near to the corners are

top left corner: $[(λ_{l} w_{I_{f}}, 0), (0, λ_{l} h_{I_{f}})]$ ;

top right corner: $[(w_{I_{f}} - 1 - λ_{l} w_{I_{f}}, 0), (0, λ_{l} h_{I_{f}})]$ ;

bottom right corner: $[(w_{I_{f}} - 1 - λ_{l} w_{I_{f}}, 0), (0, h_{I_{f}} - 1 - λ_{l} h_{I_{f}})]$ ; and

bottom left corner: $[(λ_{l} w_{I_{f}}, 0), (0, h_{I_{f}} - 1 - λ_{l} h_{I_{f}})]$ ,

where λ_l is a coefficient controlling the positions of the four lines and 0 < λ_l < 0.5. The four lines near to the edges include two lines parallel to the top and bottom edges, Line _t and Line _b , and two lines parallel to the left and right edges, Line _t and Line _r . The width of Line _t and Line _b is equal to that of the crosshair, and the height of Line _t and Line _r is equal to that of the crosshair. The distances between Line _t , Line _b and the top, bottom edges, and the distances between Line _t , Line _r and the left, right edges are all λ_offset.

Let P denote the set of pixels of an input image and N denote the set of adjacent pixels. Let $A = (A_{1}, . . ., A_{p}, . . ., A_{| P |})$ , be a segmentation of the image. A_p belongs to either object (A_p = OBJ) or background (A_p = BKG), where OBJ and BKG are labels for an object pixel and a background pixel, respectively. And let E(A), R(A), and B(A) denote the energy function, region energy, and boundary energy, respectively.

E (A) = λ \cdot R (A) + B (A),

where

R (A) = \sum_{p \in P} R_{p} (A_{p});

B (A) = \sum_{{p, q} \in N} B_{{p, q}} \cdot δ (A_{p}, A_{q});

and

δ (A_{p}, A_{q}) = {\begin{cases} 1 if A_{p} \neq A_{q} \\ 0 if A_{p} = A_{q} \end{cases} .

The parameter λ reflects the weight of the region energy R(A). R(A) represents the individual penalties for assigning the pixel p to OBJ and BKG, that is, R_p(OBJ) and R_p(BKG). R_p(OBJ) is large when p is similar to the user-specified object, and R_p(BKG) is large when p is similar to the user-specified background. The coefficient B_{p,q} serves as a penalty for a discontinuity between p and q. B_{p,q} is large when p and q are similar.

The goal of segmentation is to solve the following optimization problem:

\min E (A) = λ \cdot R (A) + B (A);

s . t . A_{p} = {\begin{cases} OBJ if p \in O \\ BKG if p \in B \end{cases} .

Graph cuts can serve as a good solution to this problem. To this end, we define a graph $G = 〈 V, E 〉$ , where the vertex set $V = P \cup {S, T}$ and E is the set of edges connecting these vertices. S is the OBJ terminal and T is the BKG terminal. The edge set is $E = N \cup {{p, S}, {p, T}}$ . The weights of edges in E is computed by

{p, q} = B_{{p, q}} if {p, q} \in N;

{p, S} = {\begin{cases} » \cdot R_{p} (BKG) if p \in P - (O \cup B) \\ K if p \in O \\ 0 if p \in B \end{cases};

{p, T} = {\begin{cases} » \cdot R_{p} (OBJ) if p \in P - (O \cup B) \\ 0 if p \in O \\ K if p \in B \end{cases},

where $K = 1 + \max_{p \in P} \sum_{q : {p, q} \in N} B_{{p, q}} .$

It can be proven that the minimum cost cut $\hat{C}$ on graph G corresponds to the optimum segmentation of the image.²¹ If ${p, S} \in \hat{C}$ , the pixel p will be segmented as the background; otherwise, there must be ${p, T} \in \hat{C}$ , and p will be segmented as the object.

Image Post-processing

After the binary image I_b containing the object (tumor) and the background (tissue) was obtained, it was expanded to the original size. Morphological processing was conducted to refine the boundary, with opening to reduce the spicules and closing to fill the holes. The post-processing includes the following three steps:

Expand I_b by a factor of 2 using bicubic interpolation. The expanded image $I_{b}^{2}$ was obtained.

Conduct opening and closing to $I_{b}^{2}$ . The processed image I_m was obtained. Both the opening and closing used a 9 × 9 elliptical kernel.

Find the contour of the object on I_m.

Implementation of the Proposed Method

The proposed method was implemented with the C++ language using OpenCV 2.4.3 and Visual Studio 2010. OpenCV is an open source computer vision library (http://opencv.org/), which has aroused increasing interest from users and developers worldwide.¹⁹ We used OpenCV to shorten the development cycle of BUS image segmentation software. To solve the minimum cut problem in graph cuts, we used the max-flow/min-cut code, which was kindly provided by Y.Y. Boykov (see http://vision.csd.uwo.ca/code/).

Table 2 shows the OpenCV functions called by the proposed method. Table 3 shows the values of some important parameters. These values were obtained by trial and error on 58 BUS images collected from the Model T2000 US scanner (Terason Ultrasound, Burlington, Massachusetts) for the sake of both segmentation precision and time efficiency. These images were provided by Symbow Medical Technology Co., Ltd. (Beijing, China). The feasibility of these values for other manufacturers’ BUS images will be demonstrated.

Table 2.

OpenCV Functions Called by the Proposed Method.

Purposes	OpenCV Functions
Image resizing	cvResize()
Gaussian filtering	cvSmooth()
Histogram equalization	cvEqualizeHist()
Pyramid mean shift filtering	cvPyrMeanShiftFiltering()
Morphological processing	cvMorphologyEx()
Contour finding	cvFindContours()

Table 3.

Values of Some Important Parameters.

Parameters	Values
Layers of the Gaussian pyramid: L	3
The spatial part: σ_s	20
The range (color) part: σ_r	40
The coefficient controlling the weight of the width of the crosshair: λw	0.4
The coefficient controlling the weight of the height of the crosshair: λh	0.5
The coefficient controlling the positions of the four lines near to the corners: λl	0.2
The coefficient controlling the distances between the lines and the edges: λoffset	4

Quantitative Evaluation of the Proposed Method

To quantitatively evaluate the segmentation precision of the proposed method, we used area error metrics that can evaluate the percentage of how much the tumor region is correctly covered by the segmented tumor region and how much is wrongly covered.⁷ The TP area ratio, FP area ratio, false negative (FN) area ratio, and SI rate were computed.

TP Area Ratio = \frac{| A_{i} \cap A_{m} |}{| A_{m} |};

FP Area Ratio = \frac{| A_{i} \cup A_{m} - A_{m} |}{| A_{m} |};

FN Area Ratio = 1 - TP Area Ratio = \frac{| A_{i} \cup A_{m} - A_{i} |}{| A_{m} |};

SI = \frac{| A_{i} \cap A_{m} |}{| A_{i} \cup A_{m} |},

where A_i is the pixel set of the segmented tumor region by the proposed method and A_m is the pixel set of the manually delineated tumor region by the radiologist. Since FN can be computed from TP, we evaluated three metrics: TP, FP, and SI.

Results

Three examples of benign tumor segmentation are shown in Figures 3 to 5. Another three examples of malignant tumors are shown in Figures 6 to 8. For each figure (a) is the original image; (b) is the ROI image; (c) is the preprocessed image; (d) is the mean shift filtered image; (e) is the graph cuts segmented and post-processed image; (f) is the segmentation result image; (g) is the manual delineation by the radiologist; (h) is the result image using only graph cuts; (i) is the level set segmented image, using the DRLSE (Distance Regularized Level Set Evolution) code⁴² (http://www.engr.uconn.edu/~cmli/) and the default parameters therein, with the initial zero-level contour set as a 50 × 50 square centered at the center of the ROI image (number of iterations = 510); (j) is the K-means segmented image, K = 2, and (k) is the graph-based segmentation result image, using the Efficient Graph-Based Image Segmentation code⁴³ (http://cs.brown.edu/~pff/segment/) and the following parameters: sigma (to smooth the image) = 0.8, k (constant for threshold function) = 100, and min (minimum component size, enforced by post-processing stage) = 40. The original images in Figures 3 to 5 and Figures 6 to 8 are from HDI 5000, HDI 3500, and LOGIQ 500 PRO, respectively. Figure 3 shows a case of a large benign tumor with small noisy regions inside it and near to the contour. The proposed method can precisely segment the tumor without including any non-tumor region. Figure 4 shows a case of a small benign tumor located in a noisy, low-contrast image. The method can successfully find the tumor contour, while the defect is that the segmented contour is not as smooth as the manual delineation and there is a slight under-segmentation region at the left of the tumor. Figure 5 shows a case of a medium benign tumor, which is noisy, low contrast, and has posterior shadow effect at the bottom left and bottom right of the tumor. The method is able to handle with this case well, but the contour is a little less smooth, and there is a small over-segmentation at the bottom right of the tumor. The reason to such over-segmentation may be due to the shadow effect. Note that the bottom left part of the tumor is segmented successfully. Figure 6 shows a case of a medium malignant tumor with angular margins and noises inside it. The segmentation result is capable to keep these angular shapes. Figure 7 shows a case of a large malignant tumor that is noisy, blurry, and low contrast. Our method manages to detect the real contour of the tumor, with a trivial drawback that the bottom right corner of the tumor is a little under-segmented. Figure 8 shows a case of a small malignant tumor in a noisy, low-contrast image with a large noisy region inside the tumor. Our method is insensitive to these factors and succeeds in contouring the tumor precisely. The experiments demonstrate that the proposed method can segment the tumors of different sizes and shapes accurately; even when the background is quite complicated, noisy, or low contrast, when the tumor boundary has many angular margins, and when the image quality is poor. Figures 3(h) to 8(h) show that using solely graphs cuts is insufficient to extract the tumor contour. Figures 3(i) to 8(i) show that using the level set evolution technique cannot find the true tumor contour. Figures 3(j), 3(k) to 8(j), 8(k) show that using K-means and graph-based segmentation methods fail to detect the tumor region precisely.

Figure 3.

(a) Original image. (b) Region of interest image. (c) Preprocessed image. (d) Mean shift filtered image. (e) Graph cuts segmented and post-processed image. (f) Segmentation result image. (g) Manually segmented image. (h) Result image using only graph cuts. (i) Level set segmented image. (j) K-means segmented image, K = 2. (k) Graph-based segmentation result image.

Figure 4.

Figure 5.

Figure 6.

Figure 7.

Figure 8.

Figure 9 shows the segmentation results of the methods by different combinations. Table 4 shows the quantitative evaluation of these combinations (n = 69). It can be seen that using Graph cuts alone did not yield a satisfying result (53.7% TP, 26.5% FP, 41.8% SI). When using Gaussian filtering + graph cuts not only TP and SI were increased but FP was also increased (79.4% TP, 32.3% FP, 54.5% SI). Combining histogram equalization not only improved SI and reduced FP but also decreased TP (72.8% TP, 19.6% FP, 63.2% SI). Incorporating mean shift can improve TP and SI and reduce FP to some extent (91.3% TP, 11.4% FP, 85.1% SI). Combining morphological processing not only increased slightly TP and SI but also increased FP a little. The proposed method had a TP rate of 91.7%, an FP rate of 11.9%, and a SI rate of 85.6%. These values indicate that our method may be used for BUS image segmentation. Table 5 shows the quantitative evaluation of the proposed method for different US scanners. It can be found that HDI 5000 (n = 22, 93.1% TP, 9.7% FP, 89.6% SI) has the best performance and LOGIQ 500 PRO (n = 20, 91.5% TP, 11.3% FP, 85.8% SI) is a little better than HDI 3500 (n = 27, 90.7% TP, 14.1% FP, 82.2% SI). This performance difference among different scanners may be attributed to the difference in their imaging quality, which is related to such factors as number of channels, A/D converter, beamforming, and signal/image processing.

Figure 9.

Segmentation results by using (a) graph cuts, (b) Gaussian filtering + graph cuts, (c) Gaussian filtering + histogram equalization + graph cuts, (d) Gaussian filtering + histogram equalization + mean shift + graph cuts, (e) Gaussian filtering + histogram equalization + mean shift + graph cuts+ morphological processing (the proposed method), and (f) manually segmented image.

Table 4.

Quantitative Evaluation of the Methods by Different Combinations (n = 69).

	TP (%)	FP (%)	SI (%)
GC	53.7	26.5	41.8
GF + GC	79.4	32.3	54.5
GF + HE + GC	72.8	19.6	63.2
GF + HE + MS + GC	91.2	11.4	85.1
GF + HE + MS + GC + MP	91.7	11.9	85.6

TP = true positive; FP = false positive; SI = similarity; GC = graph cuts; GF = Gaussian filtering; HE = histogram equalization; MS = mean shift; MP = morphological processing.

Table 5.

Quantitative Evaluation of the Proposed Method for Different Ultrasound Scanners.

	n	True positive (%)	False positive (%)	Similarity (%)
LOGIQ 500 PRO	20	91.5	11.3	85.8
HDI 3500	27	90.7	14.1	82.2
HDI 5000	22	93.1	9.7	89.6

The run time of each case was recorded for statistical analysis. The computer configuration was Intel Core 2.66 GHz CPU and 4 GB RAM. The mean run time was computed as M ± SD, where SD is the standard deviation. The mean run time was 0.49 ± 0.36 s, which is the time required for the segmentation process, but not including the time for interaction, that is, selecting the two end points for a ROI. The run time of our method demonstrates that the method is fast. This fast speed is owing to the OpenCV implementation, fast max-flow/min-cut solution, and ROI image shrinking by a factor of 2. It can be seen that most of the methods listed in Table 1 have a run time greater than 9 s. Note that the run time of a BUS image segmentation method is also related to the CPU and implementation language in addition to its computational complexity.

Discussion

Graph cuts have been widely used in computer vision applications. However, when directly applied to BUS images, graph cuts cannot produce satisfying segmentation results, as shown in Figures 3(h) to 8(h). This is because BUS images have inherent characteristics including low-contrast, speckle noise, and tissue-related textures. To overcome these problems, we combine Gaussian filtering, histogram equalization, and pyramid mean shift filtering with graph cuts. Gaussian filtering is used to filter high-frequency noises. Histogram equalization solves the low-contrast problem. Pyramid mean shift filtering enhances the homogeneity of BUS images at the same time of suppressing speckle noise and tissue-related textures (Figure 2). By applying these processing steps before graph cuts, the tumor is successfully segmented, as shown in Figures 3(f) to 8(f). The segmentation result by our method is close to the manual delineation by an experienced radiologist as shown in Figures 3(g) to 8(g), with the only defect that the detected boundary of the tumor is a little less smooth when the BUS image is quite noisy and blurry. This combination has endowed our method with a promising accuracy, especially when segmenting low-quality BUS images, whereas using solely the K-means, graph-based, or level set methods is not easy to produce accurate segmentation results, as shown in Figures 3(i), 3(k) to 8(i), 8(k). In addition, the feasibility of the proposed method for the segmentation of BUS images from different US scanners has been demonstrated (Table 5).

Besides its potential applicability for a BUS CAD system, the proposed method might be used for semi-automatic BUS biometrics or measurements. To date, the majority of radiologists still prefers to use manual contouring for biometrics or measurements on a US scanner, either during real-time scanning or after the desired image is frozen. Although manual segmentation can produce satisfying contours, it is time consuming and operator-dependent. Therefore, the proposed semi-automatic segmentation technique may be helpful in performing biometrics or measurements. In addition, the proposed method may be used for segmentation of other kinds of tumors on US images, such as the liver tumor.

One limitation of our method is that the smoothness of the detected tumor contour is a little lower when the original boundary is noisy or blurry. Another limitation is that the interaction procedure of selecting the two diagonal points for a ROI affects the final segmentation accuracy to some extent because the seeds for object and background in graph cuts are determined by the width and height of the ROI. We suggest that a ROI be selected, which is 1.1 to 1.3 times the rectangular closely encompassing the tumor. As the ROI can be automatically detected by using a nonparametric model based on order statistics²² or a method combing hybrid filtering, multifractal processing, and thresholding segmentation,²³ and so forth, the proposed method may be extended to fully automated BUS image segmentation in the future.

Conclusion

This paper presented a new semi-automatic tumor segmentation method for BUS images combining Gaussian filtering, histogram equalization, pyramid mean shift, and graph cuts. First, the ROI image was shrunken by a factor of 2 to reduce computation time. After smoothed by Gaussian filter and contrast-enhanced by histogram equalization, the ROI image was filtered by pyramid mean shift and segmented by graph cuts. The object and background seeds for graph cuts were automatically generated according to the ROI size. The method was implemented with OpenCV 2.4.3 and Visual Studio 2010 using the C++ language. Sixty-nine BUS images collected from different US scanners were used to validate the effectiveness of the proposed method. The method had a TP rate of 91.7%, a FP rate of 11.9%, and a SI rate of 85.6%. The mean run time on Intel Core 2.66 GHz CPU and 4 GB RAM was 0.49 ± 0.36 s. These results indicate that the proposed method may be helpful for BUS image segmentation. Future work can be conducted on fully automatic BUS image segmentation based on the proposed method. One potential way is to automatically detect the ROI.

Footnotes

Acknowledgements

The authors would like to thank the anonymous reviewers for their valuable comments and suggestions. They thank Y. Luo from West China Hospital of Sichuan University. They would like to thank Y.Y. Boykov for kindly providing the max-flow/min-cut code.

Declaration of Conflicting Interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work was supported by the Doctoral Fund of Innovation of Beijing University of Technology and in part by the National Natural Science Foundation of China (Grant 81127006).

References

World Health Organization. Breast cancer: prevention and control. Available from http://www.who.int/cancer/detection/breastcancer/en/index1.html

Cheng

Shi

Min

Cai

. Approaches for automated detection and classification of masses in mammograms. Pattern Recogn. 2006;39:646-68.

Jesneck

Baker

. Breast mass lesions: computer-aided diagnosis models with mammographic and sonographic descriptors. Radiol. 2007;244:390-8.

McCavert

O’Donnell

Aroori

Badger

Sharif

Crothers

. Ultrasound is a useful adjunct to mammography in the assessment of breast tumours in all patients. Int J Clin Pract. 2009;63:1589-94.

Giger

Karssemeijer

Schnabel

. Breast image analysis for risk assessment, detection, diagnosis, and treatment of cancer. Annu Rev Biomed Eng. 2013;15:327-57.

Noble

Boukerroui

. Ultrasound image segmentation: a survey. IEEE Trans Med Imaging. 2006;25:987-1010.

Shan

Cheng

Wang

. Completely automated segmentation approach for breast ultrasound images using multiple-domain features. Ultrasound Med Biol. 2012;38:262-75.

Madabhushi

Metaxas

. Combining low-, high-level and empirical domain knowledge for automated segmentation of ultrasonic breast lesions. IEEE Trans Med Imaging. 2003;22:155-69.

Chen

Chang

Moon

. 3-D breast ultrasound segmentation using active contour model. Ultrasound Med Biol. 2003;29:1017-26.

10.

Huang

Jiang

Chen

Moon

. Level set contouring for breast tumor in sonography. J Digit Imaging. 2007;20:238-47.

11.

Xiao

Brady

Noble

Zhang

. Segmentation of ultrasound B-mode images with intensity in homogeneity correction. IEEE Trans Med Imaging. 2002;21:48-57.

12.

Huang

Lee

Liu

Jin

. A robust graph-based segmentation method for breast tumors in ultrasound images. Ultrasonics. 2012;52:266-75.

13.

Gao

Yang

Liao

Liu

Feng

Chen

. Segmentation of ultrasonic breast tumors based on homogeneous patch. Med Phys. 2012;39:3299-318.

14.

Drukker

Giger

Horsch

Kupinski

Vyborny

Mendelson

. Computerized lesion detection on breast ultrasound. Med Phys. 2002;29:1438-46.

15.

Huang

Chen

. Watershed segmentation for breast tumor in 2-D sonography. Ultrasound Med Biol. 2004;30:625-32.

16.

Gómez

Leija

Alvarenga

Infantosi

Pereira

. Computerized lesion segmentation of breast ultrasound based on marker-controlled watershed transformation. Med Phys. 2010;37:82-95.

17.

Cheng

Shan

Guo

Zhang

. Automated breast cancer detection and classification using ultrasound images: a survey. Pattern Recogn. 2010;43:299-317.

18.

Shan

Cheng

Wang

. A novel segmentation method for breast ultrasound images based on neutrosophic l-means clustering. Med Phys. 2012;39:5669-82.

19.

Bradski

Kaehler

. Learning OpenCV: Computer Vision with the OpenCV Library. Sebastopol, CA: O’reilly Media; 2008.

20.

Comaniciu

Meer

. Mean shift analysis and applications. In: 7th IEEE International Conference on Computer Vision, Kerkyra, Greece, 20-27 September 1999, pp. 1197-203.

21.

Boykov

Jolly

. Interactive graph cuts for optimal boundary & region segmentation of objects in ND images. In: 8th IEEE International Conference on Computer Vision, Vancouver, British Columbia, Canada, 7-14 July 2001, pp. 105-12.

22.

Mogatadakala

Donohue

Piccoli

Forsberg

. Detection of breast lesion regions in ultrasound images using wavelets and order statistics. Med Phys. 2006;33:840-9.

23.

Yap

Edirisinghe

Bez

. A novel algorithm for initial lesion detection in ultrasound breast images. J Appl Clin Med Phys. 2008;9:181-99.

24.

Shen

Huang

Chang

. Computer-aided multiview tumor detection for automated whole breast ultrasound. Ultrason Imaging. 2014;36:3-17.

25.

Varghese

Jiang

Zagzebski

. In vivo classification of breast masses using features derived from axial-strain and axial-shear images. Ultrason Imaging. 2012;34:222-36.

26.

Liao

Yeh

. Texture feature analysis for breast ultrasound image enhancement. Ultrason Imaging. 2011;33:264-78.

27.

Alam

Feleppa

Rondeau

Kalisz

Garra

. Ultrasonic multi-feature analysis procedure for computer-aided diagnosis of solid breast lesions. Ultrason Imaging. 2011;33:17-38.

28.

Liu

Cheng

Huang

Tian

Tang

Liu

. Probability density difference-based active contour for ultrasound image segmentation. Pattern Recogn. 2010;43:2028-42.

29.

Gao

Liu

Chen

. Phase-and GVF-based level set segmentation of ultrasonic breast tumors. J App Math. 2012;2012:Article 810805.

30.

Takemura

Shimizu

Hamamoto

. A cost-sensitive extension of AdaBoost with markov random field priors for automated segmentation of breast tumors in ultrasonic images. Int J Comput Assist Radiol Surg. 2010;5:537-47.

31.

Jiang

Peng

Zhang

Cheng

Megalooikonomou

Ling

. Learning-based automatic breast tumor detection and segmentation in ultrasound images. In: 9th IEEE International Symposium on Biomedical Imaging (ISBI), Barcelona, Spain, 2-5 May 2012, pp. 1587-90.

32.

Chang

Moon

Chen

. Automatic ultrasound segmentation and morphology based diagnosis of solid breast tumors. Breast Cancer Res Treat. 2005;89:179-85.

33.

Tao

Jin

Zhang

. Color image segmentation based on mean shift and normalized cuts. IEEE Trans Syst Man Cybern Part B-Cybern. 2007;37:1382-9.

34.

Zhang

Chen

Chin

Wang

. Intelligent scanning: automated standard plane selection and biometric measurement of early gestational sac in routine ultrasound examination. Med Phys. 2012;39:5015-27.

35.

Beddoe

Slabaugh

. Automatic graph cut segmentation of lesions in CT using mean shift superpixels. Int J Biomed Imaging. 2010;2010:Article 983963.

36.

Yun

Gao

Wang

. Ultrasound image segmentation based on the mean-shift and graph cuts theory. Telkomnika. 2013;11:5600-8.

37.

Szeliski

. Computer Vision: Algorithms and Applications. Berlin, Germany: Springer; 2011.

38.

Sonka

Hlavac

Boyle

. Image Processing, Analysis, and Machine Vision. Boston, MA: PWS Publishing; 1999.

39.

Ince

Hatton

Graham-Cumming

. The case for open computer programs. Nature. 2012;482:485-8.

40.

Meer

Georgescu

. Edge detection with embedded confidence. IEEE Trans Pattern Anal Mach Intell. 2001;23:1351-65.

41.

Comaniciu

Meer

. Mean shift: a robust approach toward feature space analysis. IEEE Trans Pattern Anal Mach Intell. 2002;24:603-19.

42.

Gui

Fox

. Distance regularized level set evolution and its application to image segmentation. IEEE Trans Image Process. 2010;19:3243-54.

43.

Felzenszwalb

Huttenlocher

. Efficient graph-based image segmentation. Int J Comput Vision. 2004;59:167-81.

44.

Gonzalez

Woods

. Digital Image Processing. 2nd ed. Upper Saddle River, NJ: Prentice Hall; 2001.