Brain Tumor Segmentation from Multimodal MR Images Using Rough Sets

Abstract

Automatic segmentation of brain tumors from Magnetic Resonance images is a challenging task due to the wide variation in intensity, size, location of tumors in images. Defining a precise boundary for a tumor is essential for diagnosis and treatment of patients. Rough set theory, an extension of classical set theory, deals with the vagueness of data by determining the boundary region of a set. The aim of this work is to explore the possibility and effectiveness of using a rough set model to represent the tumor regions in MR images accurately, with Quadtree partitioning and simple K-means as precursors to indicate and limit the possible relevant regions. The advantage of using rough sets lie in its ability to represent the impreciseness of set boundaries, which is one of the major challenges faced in tumor segmentation. Experiments are carried out on the BRATS 2013 and 2015 databases and results are comparable to those reported by recent works.

1 Introduction

Glioma is a general term used to refer to brain tumors which arise from the glial cells. Gliomas vary with respect to location, cell characteristics, severity and responsiveness to treatment. The World Health Organization defines a grading system for tumors, where Grade I is the least severe and having the best prognosis, and grade IV is the most severe with worst prognosis of the disease. High-grade gliomas are the ones which are most aggressive in their spread, affording patients an average survival rate of only about two years after discovery and treatment consisting of surgery followed by radiotherapy [1]. Low-grade gliomas offer a much better chance of survival and a longer span of life expectancy. In order to offer patients the best possible chance of recovery, diagnosis of gliomas is thus a time-critical problem.

Brain tumor segmentation is carried out for a wide variety of lesion (area of abnormality in tissue) images. A major challenge in automatic segmentation is the fact that lesions are identified by the difference in their intensity profile with relation to the surrounding normal tissue. This often makes their identification a relative rather than an absolute process, and there are major variations even in case of manual segmentation by experts. Tumor boundaries become ambiguous if the intensity variation is smooth, or if there are artifacts present in the imaging. Artifacts in MRI imaging arise from patient motion (due to breathing, blood flow, heart beat), equipment motion, etc [2]. Multiple modes or types of MR Imaging are available, namely: T1, T2, T1-contrast, FLAIR (Fluid Attenuated Inversion Recovery) etc., each of which provide different anatomical and functional information, which are combined to draw meaningful conclusions about the tumor. Both the tumor and Cerebrospinal fluid (CSF) appear dark in T1 and bright in T2. In terms of image capture, FLAIR is the same as T2, with the added property of reduced CSF brightness. T1-contrast imaging is especially useful in identifying and enhancing active soft tissue lesions by increasing the visual contrast in the image [3].

Intensity normalization [4], registration, skull-stripping [5] are employed as preprocessing steps in MRI. Image intensity and texture are two of the most common features used for segmentation. Positional information is also used in some. Certain algorithms use a combination of these features along with various statistics obtained from them as added features. Early techniques use local or global thresholding [6], region growing [7], morphological edge detection [8], snack methods [9] etc.(Refer to [10] for detailed survey). Threshold-based methods compare ROI intensity with pre-determined intensity thresholds for segmentation. Thresholds can be determined either from prior knowledge or from various image statistics. Region-based methods use the similarity of pixels with their neighbors to combine them and form similar regions. Out of these similar regions, some would be the unhealthy tissue, which can be determined from previous knowledge about the appearance of such tissue. For region growing methods, the initial seed for starting the algorithm is of utmost importance.

Algorithms are based on classification or clustering methods such as Fuzzy C-Means (FCM) [11], Markov Random Fields or Conditional Random Fields, Support Vector Machines [12], Artificial Neural Networks [11], etc., have been developed. Neighborhood, shape or locality constraints can be added to make a more context-aware decision. Support Vector Machine based algorithms have also been widely used in brain tumor segmentation. Good training data, specially from multiple modes of same case, allow the algorithms to learn the expected tumor features and segment accordingly [12]. Atlas-based algorithms [13] use brain atlases constructed by creating a collection of averaged images of a large number of subjects. This allows for use of existing knowledge of normal brain structures in order to segment out the abnormal portions. However, inter-personal variance of the structures can lead to poor performance of these methods.

According to Hirano and Tsumoto [14], rough sets are a good base for tumor segmentation techniques since usage of rough approximations can model the differences between the expected and the observed shape and/or structure of the region of interest. The outer and inner approximations in a rough set can be thought of as the approximations of a region of interest in an image in terms of granules. By using rough sets, one can attempt to represent inconsistency between the existing knowledge of the ROI and the region obtained from the image itself. The authors propose that if there are N types of prior knowledge available about the image. Each of these available knowledge can be used to further refine the rough result obtained by the previous knowledge. Roy and Maji [15] have proposed a post-processing technique for brain tumor segmentation using Rough-Fuzzy clustering and unsupervised feature selection. In their approach, they have used Multiresolution wavelet analysis to extract scale-space feature vectors. The assumption here is that different tissue classes have different texture classes, and all pixels belonging to one class will have similar vectors. Clustering is also carried out using Rough-Fuzzy C-Means method (abbreviated to rRCM). It combines the concept of membership of fuzzy sets and the concept of lower and upper approximations from the theory of rough sets.

The manuscript is organized as follows: Sect. 2 describe the fundamentals of Rough Set theory. Section 3 explains proposed methodology in depth followed by Results over various measures in Sect. 4. Section 5 concludes the manuscript.

2 Rough Set Preliminaries

Conventional sets specify an object can either clearly belonging or clearly not belonging to a set X. Rough sets, first introduced by Pawlak in 1982 [16], is an approximation of sets, which is required where lack of exact knowledge prevent such absolute labeling of an object with respect to X. Rough sets look to define what constitutes the boundary of the set. If the boundary is empty, then the set is considered to be crisp. If boundary is not empty, it indicates that the set could not be defined precisely, hence it is rough.

The lower approximation of a set X with respect to indiscernibility relation $\mathbb {R}$ is defined as the set of all objects, which can be classified as X with certainty. This is the union of all equivalence classes that form a proper subset of X. For all elements x belonging to universe $\mathbb {U}$, this is denoted by: $\underline{\mathbb {R}}X\,=\,\displaystyle \bigcup _{x\in \mathbb {U}}\{\mathbb {R}(x):\mathbb {R}(x) \subseteq X\}$.

The upper approximation of the set X with respect to $\mathbb {R}$ is defined as the set of all objects, which can possibly be classified as X. This formed by all those equivalence classes that have non-empty intersection with X, and is expressed as: $\bar{\mathbb {R}}X\,=\,\displaystyle \bigcup _{x\in \mathbb {U}}\{\mathbb {R}(x):\mathbb {R}(x) \cap X \ne \phi \}$.

Fig. 1.

(a) Original image with grids (b) Showing identified granules

The boundary region of set X can be defined as the set of all objects that cannot be precisely classified as either X or not-X. This is given by the set difference of $\overline{\mathbb {R}}X$ and $\underline{\mathbb {R}}X$, that is, all the elements in the upper approximation that are not present in the lower approximation, i.e. $\mathtt {BND(X): \overline{\mathbb {R}}X - \underline{\mathbb {R}}X}$.

In order to use rough set approximation on images, each image has to be broken up into a collection of granules. The simplest method to do this is to consider the image as the universe $\mathbb {U}$ consisting of pixels. $\mathbb {U}$ is then divided into a number of non-overlapping (or overlapping) blocks, each having multiple pixels. Each such block is considered a granule G. Shown in Fig. 1 is a simple interpretation of segmentation of region of interest using rough sets. In this example, the rough set rules (used to decide which region a granule lies in) are defined on the basis of the number of pixels with intensity 0 in each granule. A granule with all pixels of 0 intensity is said to belong to the Lower Approximation. A granule with some, but not all, black pixels, are part of the Boundary of the set. Any other granule is not part of the set. The rules are mentioned in the Fig. 1 with color coding for easy visualization.

3 Methodolgy

Identification of tumor region and various substructures require information from the different modes to be combined. This section describes first the general architecture, followed by the specifics for identifying different tumor substructures.

Granulation: One way of achieving variable, non-overlapping partitioning of the image space is to use quad-tree decomposition. The image is first divided into four parts (hence the name “Quad” tree). The 256 $\times$ 256 image is initially divided into 16384 2 $\times$ 2 blocks before quad-tree partitioning. These 2 $\times$ 2 blocks are taken as the smallest possible part to be reached via partitioning. Taking 2 $\times$ 2 blocks gives more neighborhood information than is possible from using just one pixel. For partitioning, the correlation vector $v_{i}$ for each block $b_{i}$ is calculated. Each vector has 8 values, $v_{ij}, j=1\,to\,8$, $v_{ij}$ is the correlation between the block $b_{i}$ and its $j^{th}$ neighboring block $b_{ij}$. The final result of calculating the correlation arrays is a 16384 $\times$ 8 matrix of correlation vectors.

The pairwise difference between the correlation vectors of blocks in a granule are calculated. If the discrepancy between the minimum and maximum of these pairwise differences is large, then the granule is partitioned further. This threshold is used to control the granularity that is required to achieve. A higher threshold will allow more dissimilar components to belong to the same granule. In the present discussion, Euclidean distance is used to measure similarity, with greater distance between two vectors implying less similarity between them. Figure 2 shows the varying granularity achieved by varying the partitioning threshold. Thus after all partitioning is complete, a larger granule size means a larger area of comparatively low variation (i.e. a homogeneous area).

Fig. 2.

Size and Number of granules vary as the partitioning threshold takes values (a) 1.75, (b) 1.25, (c) 0.75

Finer granules ensure that the area inside a granule is strictly homogeneous, but it also means a large number of granules will have to be processed. Coarser granules have less homogeneity and are less in number. The amount of granularity has to strike a balance between the accuracy (of homogeneity inside a granule) and the amount of computation involved. Since quad-tree partitioning is performed on a square matrix, all the resultant granules, while unequal, are square matrices too.

After granulation, the next task involves identifying which granules belong to the tumor region (our region of interest). Typically, the region of interest is a small portion of the entire image. For expert judgment, the intensity of the region is the primary indication as to whether or not it is the region of interest, as discussed earlier. Designing a rough set rule based on intensity is difficult since the brightness and contrast of all images are not constant. Hence, it is problematic to assign a single threshold between tumor and non-tumor intensities which works equally well across all cases. The rule has to be obtained from the information afforded by a single case itself.

K Means Clustering: In order to do this, we use simple K-means clustering to infer a crude estimate of the intensity values that are there in the tumor region. We have three classes for clustering: tumor, possibly tumor and not tumor. Each class is initialized with a granule that we are reasonably sure of representing said class. In T2, T1-contrast and FLAIR modes of MRI, the regions of interest (and certain other regions as well) are characterized by high intensity. Calculating the image histogram, the tumor class is initialized by a granule, all of whose pixels are above the 90th percentile in the histogram. The possibly tumor class is initialized by a granule more than half but not all of whose pixels are above the 90th percentile. The granules that are found to belong to these two classes are a smaller subset of all the available granules, and contain the regions of interest.

Rough Sets: From each of the classes obtained from K-Means, we get the mean and standard deviation of intensity ($mean_t$ and $std_t$ for class 1(tumor), $mean_m$ and $std_m$ for class 2(maybe tumor)). Out of the 300 images provided as part of the BRATS database, granules generated by 100 randomly selected images were used to craft the rules which we use to distinguish between Upper and Lower Approximation of rough sets. The Granules with all pixels above $mean_t-f*std_t$ are said to be in the Lower Approximation (i.e. surely of ROI). All the granules of Lower Approximation, plus those granules whose mean intensity is above $mean_m-f*std_m$ is said to be in the Upper Approximation of the ROI. The multiplying parameter f is adjusted according to the percentile of the image intensities that the two means lie at. If there is a small difference between the two, f is kept small, since a smaller range of intensities need to accounted for. If, however, the distance between $mean_t$ and $mean_m$ is large, a larger value of f is chosen to better span the intermediate range.

For identifying the whole tumor region, the FLAIR and T2 mode is used. In order to incorporate the information from both modes, a combined image is produced, having 30% contribution from T2 and 70% from FLAIR. This helps to maintain the hyper-intensity of the tumor region, and at the same time suppress the intensity of other irrelevant areas.

Fig. 3.

(a) T2, (b) FLAIR, (c) T2 and FLAIR combined

Calculating rough sets, as described previously, we obtain the regions where the tumor is located. However, due to various imaging differences (as discussed earlier), there may still be spurious regions detected. If multiple components have similar size and high intensity, then it is probable that they are multiple sites of tumor growth. However, if such multiple sites are highly symmetrical about the central axis, then it likely that it is normal brain structures, like the brain stem or ventricles, since tumors are usually at random locations with uncontrolled growth. Also, spurious bright bands are sometimes present at the edge of the brain image. The following formula, devised experimentally, tries to incorporate these factors as (Fig. 3):

$$\begin{aligned} w = size_c/size_l + 1/(intensity_c-intensity_l) - 0.02*pixels_{center} -0.02*pixels_{edge} \end{aligned}$$

where

• $size_c$ = number of pixels in the component considered

• $size_l$ = number of pixels in the largest connected component

• $intensity_c$ = mean intensities of the component considered

• $intensity_l$ = mean intensities of the largest connected component

• $pixels_{center}$ = pixels of the component very close to the center,

• i.e. are possibly part of the brain-stem and not tumor regions

• $pixels_{edge}$ = pixels of the component on the edge of the brain image

Hence it rewards similarity of size and intensity, while penalizing the areas where stray brightness is usually noticed (i.e. at the center and at the edges). Finally, the components with high score (>0.75, decided by experiments) are retained. Results are shown in Fig. 4.

The enhancing tumor region is the area of active tumor growth with fresh lesions. Since this region has good blood flow, hence usage of a contrasting medium enhances the region in the T1-contrast mode. It lies within the area acquired in the previous step.

Fig. 4.

(a) Results of rough set calculation: Lower Approximation (pink) and Boundary (yellow) (b) tumor Region shown in green. For colours please refer to the pdf. (Color figure online)

For this ROI, the initialization of K-mean centroids and subsequent rough set classification is done as discussed previously. Results are shown in Fig. 5. The sites of dead tissue which previously had tumors, also known as necrotic tissue regions, are usually located inside of adjacent to the enhancing tumor region obtained previously. Other tissue occurring near the active tumor region are termed as the non-enhancing core tissues. These three regions (non-enhancing, enhancing and necrotic) are together termed as the gross tumor core. Process of initialization and rough set calculation is done on the T2 image in a manner similar to the above and the results are shown in Fig. 6.

Fig. 5.

(a) T1-Contrast with tumor area marked (b) Enhancing tumor region: Lower Approximation (pink) and Boundary (yellow) (c) Enhancing tumor region visualized in red over whole tumor region. For colours please refer to the pdf. (Color figure online)

Fig. 6.

(a) T2 with tumor region marked, (b) Gross tumor region: lower approximation (pink) and Boundary (yellow) (c) The necrotic and non-enhancing regions in ochre visualized over previously obtained tissue substructures. For colours please refer to the pdf. (Color figure online)

4 Results and Discussion

The algorithm is tested on the BRATS database [17, 18] which provides a collection of High Grade Gliomas (HGG), and Low Grade Gliomas(LGG) as training data, along with the ground truth segmentations. The database contains fully anonymized images from the Cancer Imaging Archive. The performance measures used here are Dice Coefficient, Positive Predictive Value and Sensitivity. The Dice Coefficient is a similarity measure between two sets A and B, and is given by

$$\begin{aligned} DC = \dfrac{2 \mid A \cap B \mid }{\mid A \mid + \mid B \mid } \end{aligned}$$

Positive Predictive Value is a popular measure used for diagnostic tests that gives the probability of a positive prediction being correct (according to the ground truth), with respect to all positive predictions made by a system. A higher value can be said to correspond to a higher accuracy of the system. It is given by

$$\begin{aligned} PPV&= \dfrac{no.\ of\ true\ positives}{no.\ of\ true\ positives + no.\ of\ false\ positives}\\&= \dfrac{no.\ of\ true\ positives}{no.\ of\ positive\ calls} \end{aligned}$$

Sensitivity is another measure widely used and is also referred to as the hit rate. It gives the probability of a positive case being identified properly by a system. It is the probability of a positive prediction, given that the case is true as per the ground truth. It is given by

$$\begin{aligned} Sensitivity&= \dfrac{no.\ of\ true\ positives}{no.\ of\ true\ positives + no.\ of\ false\ negatives} \end{aligned}$$

The Dice Coefficient, Positive Predictive Values (PPV) and Sensitivity obtained from segmenting the images as compared to the ground truth provided are given in Table 1. As of 2015, the best results reported in the MICCAI-BRATS are in the range of 73 (Hoogi)(for 100 cases) to 88.4 (Bakas)(all cases) mean Dice Coefficient value for the complete tumor region and between 60.6 (Vaidhya)(24 cases) to 82 (Agn)(30 cases) for the gross tumor core region. For the enhancing core region, the average Dice scores lie between 50.9 (Vaidhya) and 75 (Pereira)(all cases) [19]. In the proceedings of the MICCAI Brats 2014 [20], the average Dice Score for the 1st subtask (identifying the whole tumor region), varies between 0.79 (Davy) to 0.87 (Kleesiek and Urban). For the gross tumor core segmentation sub-task, reported average Dice Score has varied from 0.66 (Reza) to 0.79 (Kwon).

Table 1. Results of evaluation on BRATS 2013 and 2015 training databases. Performance measures are given as $Mean \mp Standard Deviation$ and $Median \mp Median Absolute Deviation$ for each scoring method.

Direct comparison between results is difficult since various methods have been tested on different number of cases (and different cases) from the databases, and hence the details of these reported results are also not being included in this work. The algorithm performs at par with the reported methods on the entire dataset. Usage of a combination of K-Means and rough sets help to concentrate the search for the ROI in a much smaller set of possible areas than the whole image.

The following portion shows some of the atypical cases seen during experimentation, which highlight the challenges faced by the algorithm. All the images have four sub-images. The first sub-image shows the mode from which the segmentation has been carried out, and the second shows the truth image for that slice. Third shows the corresponding region as defined by the sub-task and final sub-image shows the obtained region for that particular sub-task. Comparing the last two sub-images will give the reader an idea of the differences (or similarities) between the expected and obtained results.

In Fig. 7, the Dice Score obtained for the whole tumor region is 0.95. Figure 8 shows a case where the brain stem also gets included in the tumor region, since they are very close together and forms one connected component which reduces the Dice score to 0.59. In Fig. 9, there is a region of brightness (within the whole tumor region previously determined) which does not correspond well to the given ground truth, leading to a low Dice score of 0.59. Figure 10 shows a case where the enhancing tumor region is close to the expected region as given in the ground truth, with Dice score 0.89.

Fig. 7.

(a) Flair combined with T2 (b) Ground truth supplied (c) Area to be obtained for whole tumor region (d) Area obtained by algorithm (Dice Score: 0.95)

Fig. 8.

(a) Flair combined with T2 (b) Ground truth supplied (c) Area to be obtained for whole tumor region (d) Area obtained by algorithm (Dice Score: 0.59)

Fig. 9.

(a) Negative T1 combined with T2 (b) Ground truth supplied (c) Area to be obtained for gross tumor region (d) Area obtained by algorithm (Dice Score: 0.59)

Fig. 10.

(a) Flair combined with T2 (b) Ground truth supplied (c) Area to be obtained for enhancing tumor region (d) Area obtained by algorithm (Dice Score: 0.89)

5 Conclusion and Future Scope

In this paper, a Rough Set based approach is presented for Tumor segmentation. The proposed method utilizes the uncertainty handling capabilities of Rough set using approximation spaces. For identification of precise tumor boundaries, K-Means is used as an initial seed algorithm. The novelty of the method lies in the fusion of various modality to exact tumor information available in T2 and Flair images. Although, the present method is its nascent stage with few open parameter. However, result are quite encouraging and the method is competitive enough with state-of-the-art method.