The raw images, as the output of ultrasonic instruments, always hold characters and irons in a black background around their inside borders. Since these are useless for the image processing here, a morphological mask is designed to acquire the real ultrasonic images of breast tissues, as shown in Fig. (2a, b). The stepped gray variation of the customized mask makes it possible to detect the real margins, and this pre-processing is helpful for the complete automation. Subsequently, given the condition of speckle noises and weak edges, the real ultrasonic images need to be denoised. It has been proven that the traditional denoising algorithms (such as: Gaussian low-pass filtering, Wiener filtering and so forth) may blur the boundaries at the same time in implementing noise reduction. While, the anisotropic diffusion algorithm overcomes this weakness by introducing an edge detection operator to achieve different blur effects in areas nearby and far away from the edge area respectively. Furthermore, according to the characteristics of ultrasonic images, the speckling reducing anisotropic diffusion (SRAD) algorithm [22Yu YJ, Acton ST. Speckle reducing anisotropic diffusion IEEE T Image Process 2002; 11(11 ): 1260-70.] meliorates the edge detection by integrating the performance of gradient operator and Laplace operator, in which way the grey variation caused by the speckle noise could be exclude from the detected edges. Thus, the real ultrasonic images are treated by the SRAD, to reduce speckle noises, meanwhile to preserve and even enhance the edge of tumors, as shown in Fig. (2c).

In this section, the image obtained by the above pre-processing is divided into non-overlapping grids, as shown in Fig. (2c). Here the size of each grid is l_grid×l_grid, and the selection of grid size l_grid will be discussed in the following experiments in detail. Then the local texture, local gray level co-occurrence matrix and position features are extracted to describe each grid.

Local texture features [23Gonzalez RC, Woods RE, Eddins SL. Digital Image Processing Using MATLAB, Pearson Prentice Hall, Upper Saddle River, NJ 2004.] include six measures, which have been formulated as quantitative shape descriptions of the first-order histogram of each grid.

1. Average gray level m, which describes the difference in gray between the ROI and non-ROI, is defined by:
   
   
   (1)  m=∑i−0L−1zip zi
   
   where z, with a range of {z₀, z₁, …, z_L-1}, is the gray level variable of the L-grayscale image, and {p(z_i), i=0,1,∙∙∙,L-1} is the corresponding proportion.
2. Standard deviation σ, which characterizes the texture roughness, is expressed as:
   
   
   (2)  σ=∑i=0L−1zi−m2p zi1/2
3. Normalized smoothness R, which characterizes the texture smoothness, can be calculated by:
   
   
   (3)  R=1−11+σ2/L−12
4. Third moment μ₃, which describes the histogram skewness, is defined as:
   
   
   (4)  μ3=∑i=0L−1zi−m3p zi
5. Uniformity U, which reflects the uniformity degree of the gray distribution, can be written as:
   
   
   (5)  U=∑i=0L−1p2zi
   
   The quantity U depends on the local distribution of the image intensities. Its value reaches the maximum when the local image intensity distribution is uniform.
6. Entropy E_g, a randomness metric which can measure the image’s information content, is defined by:
   
   
   (6)  Eg=−∑i=0L−1pzilog2pzi
   
   When the image gray distribution has the biggest randomness, E_g reaches its maximum.

The Gray Level Co-occurrence Matrix (GLCM) features reflect the comprehensive information of image gray on different directions, adjacent intervals and variations. Thus these can be used to analyze the image primitive and its structure. To measure the local GLCM features, the image after the pre-processing is written as f(x, y), with a size of M×N and a grayscale of N_g. The local GLCM of each sub-block is indicated by P, which is an N_g×N_g matrix. Then P(g₁, g₁ ), the element at row g₁ and column g₂ of P, which meets a certain space relation, can be expressed by:

where g₁, g₂ are the gray level values of the current pixel pair, #{∙} is the number of elements in this set. The normalized local GLCM elements Pˆg1,g2 equals to Pg1,g2/Npair , where Npair is summation of all elements in P. If the distance between points (x₁, y₁) and (x₂, y₂) is d, and their angle with the abscissa axis is θ, Pˆg1,g2,d,θ is defined as the joint probability density of gray levels g₁ and g₂ separated by the distance d and along the direction θ. So it can not only reflect the brightness distribution, but also reflect a tabulation of how often different combinations of pixel brightness values (gray levels) occur in an image. Here, since the texture in ultrasonography has no specific direction, the normalized local GLCM element Pˆa..g1,g2 can be calculated by ∑θPˆg1,g2,d,θ/4, where θ=0,π/4,π/2,3π/4 , and d is set to 2 due to the texture roughness. Thus, three typical statistics can be calculated by the normalized local GLCM Pˆä of each grid [24Baraldi A, Parmiggiani F. An investigation of the textural characteristic associated with gray level cooccurrence matrix statistical parameters IEEE T Geosci Remote 1995; 33(2 ): 293-304.]:

1. Contrast C_n, which presents to be higher when there are more pixel pairs with high gray level difference, is defined as:
   
   
   (8)  Cn=∑g1=0Ng−1∑g2=0Ng−1||g1−g22Pˆä2g1,g2
2. Correlation C_r, a measure of the similarity degree of GLCM elements along the vertical or horizontal direction, is written as:
   
   
   (9)  Cr=∑g1=0Ng−1∑g2=0Ng−1g1g2pˆä2g1,g2−μhorμverσhorσver
   
   where
   μhor=∑g1=0Ng−1g1∑g2=0Ng−1P∧äg1,g2,μver=∑g2=0Ng−1g2∑g1=0Ng−1P∧äg1,g2,σhor2=∑g1=0Ng−1g1−μhor2∑g2=0Ng−1P∧äg1,g2,andσver2=∑g2=0Ng−1g2−μver2∑g1=0Ng−1P∧äg1,g2.
3. Energy E_t, namely the quadratic sum of the GLCM elements, which reflects the uniformity of the image’s gray and texture distribution, is expressed as:
   
   
   (10)  Et=∑g1=0Ng−1∑g2=0Ng−1P∧a..2g1,g2
   
   The image with coarser texture holds a higher E_t.

By dividing the ultrasonic image into grids, each sub-block corresponds to a coordinate value (x_b, y_b). In human’s visual system, the neighborhood property may influence the judgment, which comes from the similarity measured by other features. In this way, two sub-blocks, which are far away from each other, still hold a low possibility belonging to the same class, even if their other characteristics are quite similar. Meanwhile, it is indicated that sub-blocks locating in the center of the image is more likely to be part of the ROI. Here, the introduced position features (x_b, y_b) are multiplied by the weight w_pos to limit its impact. Considering that the position information is only used to support the distinction of objects with a large distance between them, the value of w_pos generally ranges between 0.01 and 0.1. Here it is empirically selected as 0.05.

Due to the correlation between the above 11 features, an effective feature subset needs to be selected to avoid dimension disasters. In this study, the breast tumor ROI and non-ROI background marked by the physician are taken as the training set. The feature selection tools such as the classic Principle Component Analysis (PCA) and the Sequential Floating Forward Selection (SFFS) are used to obtain the optimal feature vector and minimize the system time. Finally the achieved optimal features are: five basic local texture features (average gray, standard deviation, smoothness, third moment and uniformity), two local GLCM features (contrast and energy), and two position characteristics.

The feature vector of a grid, which consists of nine features, is normalized to [0, 1], and the position features are scaled by the weighting coefficient w_pos. Then, taking the nine dimensional vector of each grid as the input of the classifier, the ROI and non-ROI background in ultrasonic breast tumor images can be distinguished.

In this section, the self-organizing feature map neural network (SOMNN) [25Kohonen T. Self-organizing maps. New York, NY: Springer-Verlag 1997.], proposed by Kohonen according to characteristics of the nervous system, is used for the classification of these grids to detect the ROI automatically. The SOM is trained by the unsupervised learning to produce a low-dimensional discretized representation of the input space of the training samples, and it uses a neighborhood function to preserve the topological properties of the input space, which makes it different from other artificial neural networks.

In this way, the SOM is able to classify the sub-blocks into n clusters only by the unsupervised learning. In our study, the output layer number n is set to 2 (ROI and non-ROI output) and the feature vector extracted from each grid is taken as the input. To map the sorted sub-blocks back to the size of the original ultrasonic images, the initial result of ROI detection can be represented by a binary image as shown in Fig. (2d), where white areas indicate the ROI candidates while black areas the non-ROIs.

Due the influence of the speckle noise, there will exist some isolated points in the binary image after pixel classification. Here, a morphological opening is used to remove isolated dots, burrs and tiny bridges, and a morphological closing is used to fill holes and small cracks, while keeping the general position and shape unchanged.

Since the hypoechoic areas also comprise the subcutaneous fat and muscle region, which scatter around the upper and lower bound of the image, the following rules are designed to remove their undesirable impact and detect the real ROI, according to the general size and position of a real ROI in an ultrasonic image. Firstly, the largest connected region is found out as B₁ and the second largest connected region B₂. Then,

1. if B₁ does not reach the image bound, take B₁ as ROI, else continue;
2. if B₂ reaches the image bound, take B₁ as ROI, else continue;
3. if the area of B₂<1/25 of the original image, take B₁ as ROI, else take B₂.

Furthermore, the pre-processing using the SRAD and the adjustable sub-block size l_grid instead of a constant value in [6Liu B, Cheng HD, Huang JH, Tian JW, Liu JF, Tang XL. Automated segmentation of ultrasonic breast lesions using statistical texture classification and active contour based on probability distance Ultrasound Med Bio 2009; 35(8 ): 1309-24.], are proven to be beneficial to remove the adjacent shadowing artifacts which share similar features as the real ROI does to some extent.

After the grids selection, an irregularly shaped white area R_cand is identified as the candidate ROI. The enclosing rectangle of R_cand is found, and its border is then expanded outwards by l_grid×4 pixels (if it exceeds the bound of the original image, take the bound of image as the final result) to gain the final rectangular ROI used for the subsequent segmentation. The magnification coefficient 4 is obtained by the experiments as described later.

As mentioned above, the implementation flow of the automatic ROI generation is shown in Fig. (2), and Fig. (2e, f) is the final detected ROI.

The basic principle of clustering is to minimize the similarity among different clusters whereas to maximize the similarity inside each cluster itself. In this study, the graph theory-based clustering algorithm is applied for image segmentation. Firstly, the ROI generated above is processed by the SRAD [22Yu YJ, Acton ST. Speckle reducing anisotropic diffusion IEEE T Image Process 2002; 11(11 ): 1260-70.], whilst the resultant I_S is taken as the object image to be segmented. Then a weighted undigraph G = (V, E) is obtained, where V is the node set (composed of all pixels in I_S) and E is the similarity between different pixels.

Given an image with a size of N × N, its adjacency matrix can be noted as W ∈ N₂ × N₂, which is regarded as the similarity matrix as well. By dividing the image into two disjoint parts A and B (A ∪ B = V ), we can define the sum of weight between A and B as WcutWcut A,B=∑i∈A,j∈BW i,j . Then the optimum segmentation of the image can be achieved by minimizing the value of W_cut. However, there are always some clusters which only consist of individual nodes in the result of Minimum Cut [26Wu Z, Leahy R. An optimal graph theoretic approach to data clustering: Theory and its application to image segmentation IEEE T Pattern Anal 1993; 15(11 ): 101-3.]. This is an unfavorable phenomenon for image segmentation. Fortunately, it can be overcome by Normalized Cut (Ncut) [27Shi JB, Malik J. Normalized cuts and image segmentation IEEE T Pattern Anal 2000; 22(8 ): 888-905.]. Ncut is targeted to find a normalized minimum cut W_cut (A,B) by firstly respectively dividing W_cut (A,B) by assoc (A,V) and assoc (B,V), then adding up the ratios, where assoc (A,V) represents the total connections between A and V, and so it is the same with assoc (B, V).

where assocA,V=∑i∈A,j∈VWi,j .

The optimum solution is to minimize the value of W_Ncut between A and B, which is a complex NP-complete. In [27Shi JB, Malik J. Normalized cuts and image segmentation IEEE T Pattern Anal 2000; 22(8 ): 888-905.], Shi and Malik defined X as a matrix of N²×N_C, where N_C is the number of clusters. X_l is the l^th column vector of X. If a node i belongs to C_l , x_{i ,l} which is the i^th row and l^th column element of X is set to 1, otherwise x_i.l = 0. They also defined the connectivity of the i^th node as di=∑jWi,j, and a diagonal matrix D (D(i,i) = d_i). Based on a series of linear algebra derivation, the search of the minimum WNcutNC can be transformed to the optimization of X.

The constraint condition is: &KHgr;∈|0,1N2×NC,&KHgr;1NC=1N2 where , 1NC and 1N2 are vectors of all ones with sizes of N_C×1 and N²×1 respectively. Let &Zgr;=&KHgr;&KHgr;TDX−1/2 , Equation (12) can be expressed by:

Here, the constraint condition can be rewritten as ZTDZ=INC . In (13), tr represents the trace of a matrix, and INC is an identity matrix of N_C×N_C. By extending the range of Z to the set of real numbers, the discrete problem will be transformed to a continuous optimization problem, which is much easier to be handled.

In this paper, by taking pixels in ultrasound images as nodes in the graph theory, the similarity matrix W can be calculated based on the different position and gray level of each pixel.

where L_C (i) represents the position of node i, or equivalently its coordinate inside I_S. F(i) describes the gray level information of node i. r is the set effective range, σF and σLC are the control parameters. In our study, the brightness of i’s neighborhood nodes is also considered when F (i) is defined, due to the ultrasonic image characteristics.

In [28Qiu PH. Image processing and jump regression analysis. Hoboken NJ: John Wiley & Sons 2005.], Qiu PH demonstrates that a proper local smoothing of the observed image intensities would benefit edge detection, because neighboring image intensities often carry helpful information about the edge structure at a given pixel. Thus, a mask of 5×5 is designed here as shown in Fig. (3). The gray level of the image center is g (i,i), with its weight set to 1. The weights of its neighborhood pixels are also marked. The further away a node is to the image center, the smaller its weight will be, and the darker its brightness is equivalently. Then in the center of an image, F (i) can be expressed by ∑kFik/2,k=0,1,2,3,4,5 , where Fik can be defined as follow:

where a should satisfy the condition

4×a+a2+a2+a22+8×a5=1

that is a = 0.0724. By substituting the value of a into the definition of F (i), the pixel i’s gray level information considering weighted neighborhood gray level is obtained. Then by substituting the value of F (i) into Equation (14), the similarity matrix W used in NCut fully automatic segmentation can be calculated.

In ultrasonic images, due to the severe presence of speckle noise, low contrast, as well as the irregular shape and area of objects, it is hard for 2-way Ncut to segment the tumor and the background directly. So we partition the ROI image into N_C (N_C≧2) clusters. According to the gray information and space distribution of each cluster, the ones belong to the tumor area can be easily screen out from the background tissue. For instance, the clustering result of the modified Ncut with N_C set to 5 is shown in Fig. (4a, b). Firstly, by sorting the average gray of each cluster in ascending or descending order, the kink of the gray distribution curve can be found. Then, the clusters with a higher average gray than that of the kink point (e.g. Cluster a2, a4, a5 and b1, b3, b4 in Fig. (4a, b)) are removed. Subsequently, due to the fact that tumors always locate centrally, an elongated mask across the center of the ROI image is taken to eliminate clusters which have no intersects with this mask ((e.g. Cluster a3 in Fig. (4c))). Lastly, a morphological post-processing is used to eliminate the minute crevices. In this way, the boundary C_I of the checked regions I_C is the result of N_C-way Ncut automatic segmentation, just as shown in Fig. (4c, d).

It is proven that a quite different result of partition is obtained with a different N_C parameter. So it is of great significance to choose an appropriate N_C. Traditionally the Ncut algorithm is invoked recursively to over-segment the image, and then the region merging is used to get the lesions area [29Liu X, Huo ZM, Zhang JW. Automated segmentaion of breat lesions in ultraound images In: Proceddings of the 27 th Annual Conference of IEEE Engineering in Medicine and Biology Society; 2006 Jan 17-18; Shanghai, China. 2005.]. However, our study proves that the over-segmentation is unnecessary with an automatic searching algorithm of N_C. Here, by setting N_C to different values ranges from 2 to 10, it is found that the modified Ncut results in an ideal segmentation with N_C belongs to {5, 6, 7}. So the optimal N_C value is searched from 5 to 7, with the stopping criterion of area(I_C(t))-area(I_C(t+1))< (1/16)×area(I_C(t)) (t is the trial times). Then I_C(t) could be found as the final checked regions I_C.

For natural images, the Ncut can result in an ideal segmentation. However, for some poor quality and edge-missing ultrasonic breast tumor images, the boundary C_I gained by the method in section 2.2.1, may have a certain Hausdorff distance and mean distance error (MDE) from the real contour (i.e. gold standard, the boundary marked by physicians). To overcome this problem, we employ the active contour model containing a local gray level item in the energy function. So this model may refine the imprecise cases to much more accurate breast tumor segmentation.

The concept of Region-scalable Fitting Active Contours (RFAC) is proposed in [30Li CM, Kao CY, Gore JC, Zhao HD. Minimization of region-scalable fitting energy for image segmentation IEEE T Image Process 2008; 17(10 ): 1940-9.]. For each node i∈Ω (the image domain), the regional energy item ϵiFit is defined as:

where k=1 or 2 denotes the node i outside or inside the closed curve respectively, f₁, f₂ represent the regional gray energy of the node outside and inside the close curve, λk is the corresponding weighting coefficient, and the kernel function K (i-j) indicates the weight assigned to each intensity I (j) at j. In level set methods, a contour ξ is represented by the zero level set of a Lipschitz function ϕ, then the following energy function is defined for ϕ:

where H is the Heaviside Function, and M1φ=Hφ,M2φ=1−Hφ . Taking C_I as its initial contour, RFAC will continuously evolve the closed curve ξ, to reduce the value of E, until it finally reaches the real edge of the target. At this time, the optimal f1, f2 have been selected to stop the whole process. Then the comprehensive consideration of global and regional energy makes it possible to reach the precise segmentation.

Pixels of an ultrasound image hold different gray level intensities. Different lesions always have different textures. In this study, T_σtumor the standard deviation of the gray level intensity in tumors, T_σedge the standard deviation of the gray level intensity in the annular boundary area, and T_RltRrt the ratio of the gray level intensity in tumors to the value in the annular boundary area are exploited as texture characteristics. Their definitions are given in Table 1.

Table 1

Texture and Morphologic Features Extracted

Table 2

The Results of ROI Extraction

Table 3

Efficiency Comparison of Three Segmentation Methods

Table 4

Performances of Different Classifiers

In a morphological way, a physician clinically differentiate between benign and malignant breast tumors based on the principle that benign ones usually hold smooth shapes whilst malignant ones tend to have irregular borders. To describe this hypothesis using the mathematical expression, five morphologic characteristics including M_compt, M_spicu, M_lobus, M_rectlike and M_eccent can be selected to describe the pathological region [31Joo SY, Yang YS, Moon WK, Kim HC. Computer-Aaided Diagnosis of Solid Brest Nodules: Use of n Aartificil Neural Network Based on Multiple Sonographic Feature IEEE T Med Imag 2004; 23(10 ): 1292-300.-32Chang RF, Wu WJ, Moon WK, Chen DR. Automatic ultrasound segmentation and morphology based diagnosis of solid breast tumors Breast Cancer Res Treatment 2005; 89(2 ): 179-85.]. They are also defined in Table 1.

In Table 1, Area and Perimeter are the area and perimeter of the tumor respectively. The polar coordinates (r, θ) of boundary pixels are calculated with the origin at their center of mass. R, (ω)) is the Fourier transform of the curve-fitted radial distance (r, θ), and the #{·}means the element number of this set. M_lobus is defined as the number of local extrema in the curve-fitted (r, θ). Stavros AT, et al. proved that the spiculation shows higher odds ratio for malignant ultrasonographic features versus benign findings in [33Stavros AT, Thickman D, Rapp CL, Dennis MA, Parker SH, Sisney GA. Solid breast nodules: use of sonography to distinguishbetween benign and malignant lesions Radiology 1995; 196: 123-34.], and this could be reflected by M_lobus. The area of the tumor versus the area of its minimum enclosing rectangle is represented as M_rectlike. L_long means the length of the longest chord connected by a pair of boundary points, and L_short the maximum of a chord set perpendicular to the longest chord.

In this way, three texture features and five morphologic features form a 1×8 vector to figure a single tumor. Lastly, it is necessary to normalize each element in the feature matrix of our ultrasonic breast tumor database into [0,1].

The Affinity Propagation (AP) clustering [20Brendan JF, Delbert D. Clustering by passing messanges between data points Science 2007; 315(5814 ): 972-7.] is a novel clustering algorithm with good performance and high speed. This algorithm avoids the influence of initialization by treating all data points as potential exemplars at its initial stage. Here, instead of traditional methods, the AP clustering is presented as an effective classifier to differentiate benignancies from malignancies in an existing database.

The main iteration procedure of the AP clustering is the updating and transmission of two kinds of message: “responsibility” r(i, k) and “availability” a(i, k). Here, r(i, k) denotes the accumulated evidence for how well-suited the point k is to be the exemplar for the point i, while a(i, k) denotes the accumulated evidence of the appropriateness if the point i selects the point k as its exemplar. In order to calculate these two kinds of message, the matrix denotes similarity between each data point is introduced. The “similarity” s(i, k) shows how well the data point k is suited to be the exemplar for the data point i. The “preference” values s(k, k) all k data points are set to a same negative real number. Consequently, the input of AP clustering is a collection of real-valued similarities between each pair of data points. Taking the whole dataset as a network of data points, these two kinds of message are exchanged along edges of this network. For a point i, if the value of k maximizes a(i, k)+ r(i, k), k is an exemplar for the point i. The values of responsibilities and availabilities are combined to monitor the exemplar decisions, and thus the iteration procedure will be stopped when any termination condition is satisfied.

In this section, the AP clustering is used to classify the ultrasonic breast tumor database without any training process. Here each tumor’s feature vector is taken as a node in the network. Due to the fact that different initialized value of “preference” results in different number of identified exemplars, an empirical “preference” is chosen to access the most rational result. With the “preference” set to the minimum of input similarities, the AP classifies the breast tumor database into five clusters firstly. Each cluster has specific and similar features and tends to hold a certain malignant probability. Then, referring to the feature distribution of breast tumors, a further identification is implemented on these five clusters to differentiate benign tumors from malignant ones.

Actually, the AP clustering can be modified by introducing a constraint in the process of message passing, so as to realize an effective and specified K-cluster AP, which can directly generate K clusters as user defined just as Zhang XL, et al. introduced in [34Zhang XL, Wang W, Nørvag K, Sebag M. K-AP: Generating specified K clusters by effvient affinity propagation: In: Proceedings of the 10th IEEE International Conference on Data Mining; 2010 Dec 14-17; Sydney, Australia. 2010.]. However, K-cluster (K=2) AP fails to reach an ideal result, because of the complex clustering structure in our database. On the contrary, the original AP clustering with the empirical “preference” can finish the superior assortment based on the similarity between each tumor case.