Semi-automatic Segmentation of Brain Tumors Using Population and Individual Information

Framework of Graph Cuts Algorithm

In graph cuts-based segmentation framework, segmentation problem is formulated in terms of energy minimization. Let X denote the set of image pixels, and N is a set of neighboring pixels. l = {l_p|l_p ∊ L} denotes a labeling, and the L = {0,1} is the segmentation label set, where 0 corresponds to the background and 1 corresponds to the foreground. The cost function is defined as following:

1

On the right-hand side of Eq. (1), the first term is the region term, which measures the sum of penalty for assigning each individual pixel p to l_p. In the standard graph cuts algorithm [16], the region term is defined as:

2

where l^obj is the label of object, l^bkg is the label of background, x_p is the image features, and F_obj, F_bkg are probability density functions of the object and background, respectively. P(l_p = l^obj | x_{p ,}F_obj) denotes the conditional probability of assigning pixel p to l^obj and P(l_p = l^bkg | x_{p ,}F_bkg) denotes the conditional probability of assigning pixel p to l^bkg.

The second term in Eq. (1) is the boundary term, which can be interpreted as a penalty for a discontinuity between p and q. As introduced in paper [16], this term is defined as

3

where G(x_p, x_q) denotes the distance between features of pixel p and q. dist(p,q) is the Euclidean distance between the locations of p and q, Parameter σ is a constant relating to the gray level which measures noise standard deviation. This penalty term is introduced when pixels p and q have different labels.

Overview of the Proposed Algorithm

We introduce high-dimensional features into graph cuts segmentation framework. NCA is proposed to learn two optimal distance metrics, which contain population and patient-specific information, respectively. The probability of each pixel belonging to the foreground (tumor) and the background is estimated by the kNN classifier under the learned optimal distance metrics. A new cost function for segmentation is constructed through these probabilities and is optimized using graph cuts. Construction of the graph is based on Kolmogorov et al. [24], and minimization of the cost function is conducted using the min-cut/max-flow algorithm [25]. Finally, some morphological operations are performed to improve the achieved segmentation results. Figure 2 shows a flowchart of the proposed brain tumor segmentation algorithm.

Image Feature Extraction

The texture features which captures the spatial correlation among adjacent pixels such as gray level co-occurrence matrix (GLCM) [26] and the moments of the gray-level histogram (MGH) [27] which represent the local structure of the image have achieved satisfactory results [28, 29]. In our paper, these two methods are used to extract the image features.

Before extracting the image features, we should normalize the intensity value. Normalization is a vital step in MR image analysis, since the intensity values of MR images are sensitive to acquisition conditions. For each slice, intensity values at the 1 and 99 % quantiles are computed for the brain region. Then, these two values are used to normalize intensity to [0, 1] through the min–max method. After intensity normalization, MGH- and GLCM-based features are extracted.

1. MGH-Based Features

Gray level histogram is a summary of the statistical information. MGH is based on the statistical properties of intensity distribution in a local region, and its effectiveness for the classification of tissues has been evaluated in [27, 28]. The feature extraction method adopted here is similar to that of [28]. Nine statistics, including the mean, variance, smoothness, uniformity, entropy, and so on, are calculated and compose the image features. The following is the list of these features:

4

5

6

7

8

9

10

where z_i indicates the pixel intensity, p(z_i) is the gray level histogram in a region, which is calculated from a W × W window with the pixel of interest in the center. H is the number of possible intensity levels. m is the mean intensity of the pixels in the W × W window, and σ² is the variance. They are defined as

11

• 2.

  GLCM-Based Features

Only using intensity distribution-based features is not sufficient because they do not take spatial information into consideration. Features extracting from the GLCM are based on the joint probability distribution of pairs of pixels, which have been proved to be effective for texture image classification [29].

The dimension of the co-occurrence matrix is equal to the number of gray levels. The gray levels in medical images are usually reduced to M discrete levels to save computational time. GLCM is a W × W dimensional matrix. Distance d and angle θ are used to calculate the joint probability distribution between pixels. The (i, j)th entry p(i, j) in the matrix gives the number of times that the gray level j follows the gray level i for a distance with an angle.

Fourteen well-known coefficients were introduced as the GLCM-based features [26]. But only six of these features are the most relevant [29]. These are energy, entropy, contrast, variance, correlation, and inverse difference moment. Here, we extract these six effective features. The following equations define these features:

12

13

14

15

16

17

where M is the dimension of GLCM. Coefficientsμ_x, μ_y and σ_x, σ_y denote the mean and standard deviations of the row and column sums of the matrix, respectively, and μ is the mean of μ_x and μ_y. They are defined as follows:

18

Coefficient p(i, j) is the (i, j)th entry in GLCM. After calculating six GLCM-based features for each p(i, j), the mean and variance of the four values (for θ = 0°, 45°, 90°, 135°) were taken as the final features, and then, 12 GLCM-based features are present.

kNN-Based Probability Model

In high-dimensional feature space, it is hard to estimate the probability density of the tumor and the background. kNN is a probability estimating method, which can be used for classification. This algorithm is implemented by computing the distances of feature vectors from the test sample to all training samples. In this model, a test sample is classified into the majority class among its k-nearest neighbors in the training dataset. It is simple to implement without any assumption on the shape of the class distributions. In this paper, a kNN-based probability model is introduced to construct the region term, estimating the probabilities of assigning each pixel to tumor or background.

Let S = {x_i, l_i | i = 1,…,n} denote a training set, x_i is the feature vector of training sample i, and l_i is the corresponding label. Under a distance metric G, the K-nearest training samples of a test pixel are determined. C_obj, C_bkg denote the number of these K training samples belonging to the object and background, respectively. The probability of pixel p belonging to the object and background can be denoted as

19

20

Given the training set S and distance metric G, kNN(x_p, S, G) denotes the K-nearest training samples of x_p. P(l_p = l^obj | kNN(x_p, S, G)) is the conditional probability of assigning pixel p to l^obj and P(l_p = l^bkg | kNN(x_p, S, G)) denotes the conditional probability of assigning pixel p to l^bkg. Using the kNN model to predict the probability, the gray level histograms need not be constructed. In our method, the region term in Eq. (1) is defined as

21

Optimal Distance Metric Learning

In the current paper, kNN classifier is introduced to construct the region term, and it can be significantly improved by learning a distance metric from labeled seeds via NCA [23]. NCA is a method for learning a Mahalanobis distance measure for kNN, and it can also reduce high-dimensional linear space that can be used for immediate classification. In this section, we briefly describe the basic principle of NCA algorithm.

Let S = {x_i, l_i | i = 1,…, n} denote a training set with sample x_i ∊ Rⁿ¹and its corresponding label l_i ∊ {0,1}. The goal of NCA is to find a transformation matrix , which projects the samples to space Rⁿ². In this new space, kNN classifier can reach the highest classification accuracy:

22

where O(T) is the sum of classification accuracies for all training samples in S.

In this new space, the distance between transformed samples x_i and x_j is measured by

23

And the similarity between x_i,x_j can be defined as

24

For one sample x_i, the probability of being correctly classified according to kNN classification rule can be written as

25

Maximizing the leave one out error on the training dataset, a transformation matrix T can be achieved. The objective function of NCA is defined as

26

Let

27

Then, the objective function of NCA can be denoted as

28

The sum of the correct classification probabilities of all training samples is calculated via Eq. (26). Using the gradient-based approaches, the maximum of the objective function can be obtained, thus, achieves the distance transformation matrix T and the optimal distance metric G (Eq.(23)) is finally used to construct the region term (Eq. (21)) and boundary term (Eq. (3)).

The Combination Model

Using graph cuts segmentation framework, we need to provide “object” and “background” seeds through manual interactions, and these training samples are very important for the construction of kNN model and the optimal distance metric learning.

There are two drawbacks in traditional interactive modes. First, the interactive users often obtain the seeds away from the boundary of the tumor, in order to avoid that the brush straddles the edge, thus neglecting the information near the edge. The kNN classifier may have difficulty near the edge of tumor because of lacking seeds in this margin region. Second, in traditional interactive modes, the “object” and “background” seeds are only taken from one patient-specific image with manual interaction, and sometimes the tumor is hard to distinguish from other tissues if the seeds are not sufficient.

To solve these problems, two groups of training sets are constructed here: (1) 68 brain MR images of other patients are used to build the population training set, denoted as S_Pop. The manual segmentations of the tumor in each image are also provided. To solve the first problem mentioned above, which is that the interactive users often obtain the seeds away from the boundary of the tumor thus neglecting the information near the edge, this population dataset is built by extracting seeds near the tumor boundaries, where the edge information can be emphasized; (2) The individual training set S_Ind is made up of the manual interactive seeds on the test image.

Based on the population and individual sets, kNN_Pop and kNN_Ind classifiers are constructed, and two transformation matrics G_Pop, G_Ind are learned from via NCA, respectively. Two segmentation models based on population and individual information are constructed, and their cost functions are denoted as E_Pop, E_Ind, respectively. In our method, a weighted combination of these two cost functions is finally used:

29

The weight of E_Pop and E_Ind is dynamically adjusted by factor α according to the position of each test pixel during the segmentation process. The choice of α is based on the following assumption: If the test pixel is closer to the manual interactive seeds, the image information between the test pixel and the individual set is more similar, thus increasing the weight of E_Ind. The weight of E_Pop increases as the test pixel gets far to the individual seeds. To restrict the weight factor α between 0 and 1, it is defined as:

30

where d_i is the shortest distance from the test pixel i to the manual interactive seeds. The value of parameter δ is determined through experiments.

Data Acquisition and Workstation

Our dataset consists of 21 patients (10 for training and 11 for testing) with total 137 T1C MR images, including 68 for training and 69 for testing. All these images were extracted from 42 scans with brain gliomas. The data were acquired at Nanfang Hospital, Guangzhou, China and General Hospital, Tianjin Medical University, China from September 2005 to October 2010 using spin echo weighted images with 512 × 512 matrix. The manual segmentations of the tumors in each image were also provided to evaluate the performance of the segmentation.

A PC with Intel P4 3.0 GHz processor and 4 GB RAM was used as the workstation for the segmentation. The software environment was a Visual Studio C++2008 MFC under Windows XP.

Parameter Choices

Many parameters exist in the current paper. In this section, we discuss how to choose the suitable values and a summary of the settings is listed in Table 1.

The coefficient β in Eq.(1) denotes the relevance of the region term to the boundary term in the cost function, computing the ratio of the region term to the boundary term, and setting the reciprocal of the ratio to β. Normally, β ranges from 5 to 10 in our dataset.

Parameter σ in the boundary term (Eq. (3)) measures the noise standard deviation. To set this parameter properly, σ was varied from 1 to 1,000 with 50 for interval in our experiments. Finally, σ was set as 500.

Parameter H in Eq. (4) is the number of possible intensity levels. However, in brain MR images, the intensity levels vary a lot from different images. The maximal gray value even exceeds 3,000 in many cases. We decrease this highest level and set H as 256 after several trials.

Parameter d, θ, and M are used to calculate features based on GLCM which is a classical feature extracting method. Normally, only the distances d = 1 and 2 with angles θ = 0°, 45°, 90°, and 135° are used. To get a fine textural segmentation, d was set to a small value (d = 1) in our experiments. As suggested in paper [30], we set M as 16 in Eq. (12) to compute GLCM.

MGH and GLCM features were computed in a W × W window with test pixel at the center. A 15 × 15 window was finally chosen after several trials.

Parameter K is the number of nearest training samples of each test pixel in kNN model. Larger value of K can reduce the effect of noise, but make boundaries between tumor and background less distinct. A good K can be selected by cross-validation. It is helpful to choose K to be an odd number in segmentation problem, and it was set as 7 in our experiments.

Parameter δ in Eq.(30) is used to adjust factor α. To set this parameter properly, δ was varied from 1 to 1,000 with 50 for interval in our experiments. Finally, δ was set as 200.

Effectiveness of the Combination Model

Some experiments were carried out to evaluate the effectiveness of the combined segmentation model. Dice similarity coefficient (DSC) was used to measure the similarity between segmentation results and ground truth.

31

where A_G represents the area of the ground truth of target and A_A represents the area of segmented target by the testing method. DSC is expressed as mean ± standard deviation. Higher mean DSC value denotes better segmentation.

Table 2 displays the performance results of 69 test images using different segmentation models. Compared to the model using Euclidean distance metric (using kNN without NCA), the mean DSC of our method is improved from 83.5 to 94.1 % with 10.6 % increase (t test p value less than 0.0001), which shows the necessity of using NCA to learn a new distance metric. Compared to population and individual segmentation models, the mean DSC of our approach increases by 6.7 % (t test p value less than 0.0001) and 1.8 % (t test p value of 0.0003), respectively. These results indicate that the proposed method using the combination model is effective.

All segmentation results of the 69 images are shown in Fig. 3. After investigating all the results, we find that even the minimum DSC of our method is 79.6 %, indicating that the combined segmentation model is robust. These results can be explained as follows. First, as introduced above, kNN classifier is used to predict the probability density of the tumors and the background. When the feature scales are not consistent with their importance, the accuracy of kNN classifier decreases using Euclidean distance metric. Thus, using NCA to learn a Mahalanobis distance metric is effective. Second, if we use individual model, taking seeds from one patient-specific image, it might be hard to distinguish the tumors from other tissues if the seeds are few. Additionally, the interactive users often obtain the seeds away from the boundary of the tumor and thus neglect the information near the edge, which leads to difficulty for kNN classifier near the boundary of the tumor. Therefore, the combination model can use the advantage of both population and individual information, thus providing better segmentation results. Two samples are shown in Fig. 4.

Effectiveness of the Combination Features

The performance of the GLCM-MGH combination features is also evaluated on the 69 test images in our experiments. Table 3 shows the DSC computed between the segmentation results and the ground truth.

Compared to method using gray value as image feature, the mean DSC of the GLCM-MGH features is improved from 86.9 to 94.1 %. This is because just using intensity changes of foreground and background is sometimes insufficient to identify the tumors from medical images. From Fig. 5, we can see the high-dimensional features based on GLCM and MGH produce more accurate segmentation results.

Comparison with Other Semi-automatic Segmentation Methods

In the introduction, many semi-automatic methods are mentioned. To demonstrate the performance of our method, other five semi-automatic methods were tested in our experiments. Figure 6 shows the segmentation results of six methods: (a) our method, (b) Random walker [13], (c) grow cut [14], (d) graph cuts [19], (e) GrabCut [17] and (f) SLIC superpixels [31].

Superpixel is another novel segmentation method. In superpixel-based methods, an image is divided into many adjacent patches (see Fig. 6f). A superpixel is an image patch where the pixel intensity is well aligned with each other. The paper [31] supplies a novel technique which can output a desired number of superpixels with regular shapes and sizes, called simple linear iterative clustering (SLIC) superpixels. Unlike other semi-automatic methods, this approach need not provide seeds, but the user should choose the patches which look like the desired object with mouse. Because manual interactions exist in this approach, we also regarded it as a semi-automatic method and compared it with other algorithms.

From Fig. 6, it can be seen that our method and GrabCut appear to outperform the other methods. In GrabCut algorithm, if we only use a single initial user interaction (rectangle in Fig. 6e), the segmentation result is not satisfying. Additionally, some brain tumors such as glioma often have high gray value near the tumor edges. This specific characteristic increases the difficulty for GrabCut to succeed on the edge of the tumor. In this cases, more foreground brushes need to be added. After adding further user edits (nine brushes in Fig. 6e), the result can be better. Considering the characteristic of brain tumor, our method seems more acceptable than GrabCut due to the simplicity of user input.

In our experiments, not all results of the six methods were investigated, but four comparable methods (Fig. 6a∼d) were used to test on the 69 images, for they had the same manual interactive seeds. Here, we also use standard segmentation measure DSC to evaluate the four methods. Figure 7 shows the performance results for these methods.

The statistical results show that our method and grow cut have better results. Grow cut algorithm can give satisfactory segmentation results but also needs the most computing time because of iterative calculating. For one image in our dataset, the computing time of grow cut was about 30 s, which is not desired by the user. Random walker shows good results when the seeds are few or placed near the center of the tumor. It is fast and less affected by the number of the seeds, but is seriously influenced by the seed location, thus leading to a higher std DSC. Traditional graph cuts algorithm needed the least computing time (no more than 1 s) but gave the worst results in our experiment. Compared to graph cuts, the mean DSC of our method is improved from 88.1 to 94.1 %. Our method has the lowest std DSC of the four methods, showing robustness to data changes. Additionally, the computing time for one image is about 4∼8 s (smaller tumor needs less time). This speed can also meet the clinical demand. Considering both the computing efficiency and the result accuracy, our method seems more acceptable by clinical use.