Deep learning-based ensemble model for brain tumor segmentation using multi-parametric MR scans

3.1 Data pre-processing

In this study, we have taken the HGG MRI of the BraTS 2017 dataset for evaluating the results. These MRIs have different artifacts such as motion and intensity inhomogeneity, and are also acquired through different scanners. To validate the algorithm and results, the artifact must be removed through normalization and bias correction.

Initially, to reduce the intensity variation between the surrounding pixels which in turn reduces the noise in an image, a median filter is applied to all 3D MRI modalities of each patient. Then, N4ITK bias field correction [5,25] is applied to all the images of all modalities in order to improve performance by removing unwanted artifacts. Finally, as the intensity values of the images differ from patient to patient, a normalization step is followed which brings the intensities of all the images to a common level. This is done by bringing the mean value of the intensity close to 0 and the standard deviation close to 1. This helps the model to generalize well by removing unwanted biasing. The normalized intensity value I _n of a slice is computed in equation (1) by using the normalization equation, where I is the original intensity value of the slice, µ is the intensity value of I, and σ is the standard deviation of I, given as:

As most of the slices of every patient consisted of completely black masks or masks having very little information regarding the tumor, the model is more inclined to learn background or noise instead of the original tumor part. In order to remove the data imbalance, a threshold value [7] of at least 0.007% tumor informative pixels of the entire image is taken. That means for every patient, the image slices that are having at least 400 pixels out of 57,600 pixels regarding the tumor information are considered for training and testing. Additionally, 1% of the border from all 4 sides was trimmed automatically from every image resulting in decreasing the image slice size from 240 × 240 to 192 × 192. These pre-process not only removed the class imbalance problem but also reduced the training time and memory requirements by a huge margin.

Once the training and testing dataset is created for each of the modalities, it was trained through the U-Net, Encoder–Decoder model, and SegNet models to segment the different tumor regions. In the next section, all the model architecture used in this study is discussed in detail.

3.2 Model architecture

The two-layered ensemble model architecture is shown in Figure 1. This model is based on the three well-known schematic segmentation models which are used widely in medical image analysis. In the first layer, all the models are explored through the multimodal MRI scans, and the maximized feature map is extracted for each individual model. In the second layer, the maximized feature map from each of the model is automatically extracted and it is fused to produce the final segmented output for all the subparts of the tumor.

Figure 1

Proposed two-layered ensemble model for WT (Whole tumor), CT (Core tumor), ED (Edema sub part), and ET (enhancing tumor) segmentation.

The U-Net has been one of the most efficient and well-known architectures [16,17] for medical image segmentation. The model considers the full feature map during up-sampling which makes the network larger in terms of memory. To make the model memory-efficient, the next deep model known as the SegNet encoder–decoder model is explored. This model uses pooling indices during up-sampling due to which it takes less memory but performance degrades. To balance this, we explore another deep model encoder–decoder model which is equivalent to SegNet without pooling indices and uses entire feature maps during up-sampling.

In the first layer, the three models are designed and experimented on all the multimodal scans to get the maximized feature map. In our work, the U-Net model consists of a contraction path, a bottleneck region, and an expansion path. The contraction path contains 4 encoding blocks where each one consists of 2 convolutional layers with a Rectified Linear Unit (ReLU) layer followed by a max-pooling layer having strides 2 × 2 and a kernel size of 2 × 2. After each block of the encoding region, the number of feature channels is doubled. The bottleneck layer contains two convolutional layers which act as an intermediary between the encoding and the decoding region. The decoding region contains a transposed convolutional layer followed by a concatenation layer to concatenate with the corresponding cropped feature map from the encoding region, and two convolutional layers that reduce the number of feature channels to half. In the encoder–decoder model, a typical encoding-decoding architecture [26] is used. The model can be divided into two steps – the encoding steps followed by the decoding steps. The encoding region consists of five encoding blocks where the number of channels is doubled at each block. The first two encoding blocks are similar to the U-Net model, but the next three encoding blocks contain three convolutional layers followed by a max-pooling layer. The decoding region also consists of five decoding blocks where the number of channels is halved at each block. The first 3 decoding blocks contain an up-sampling layer with pool size 2 × 2 followed by 3 convolutional layers. The next 2 decoding blocks contain an up-sampling layer followed by 2 convolutional layers. The implemented SegNet model is a typical SegNet model [27] which is similar to the encoding-decoding model. The encoding and decoding blocks are similar to the encoder–decoder model. Instead of the max-pooling layer, a custom layer termed as Max-Pooling-WithArgMax is used which performs typical max-pooling on the feature maps with a stride of 2 × 2 and kernel size of 2 × 2 but also stores the index position of the max-pooling pixels. Similarly, instead of the up-sampling layer, a custom layer termed Max-UpSampling-WithArgMax is used which uses the corresponding index values stored during pooling to up-sample the feature maps. The final output layer of all the models is a convolutional layer of kernel size 1 × 1 and the number of feature maps are equal to the number of classes and is followed by a sigmoid activation function to get a probability feature map.

Once the best feature map based on MRI scan is computed, in the second layer a feature map fusion is performed based on maximized feature extracted in each of the model. The network connectivity in the ensemble step can be represented as the following equation (2a) and (2b). x(i)[_CT,ET] denotes the intermediate segmented output for the core tumor and enhancing tumor and x(i)[_WT,ED] (equation (2b) denotes the intermediate segmented output for the whole tumor and edema subpart.

and

where H(x) represents a convolutional layer operation followed by a ReLU activation function applied twice on a feature map x. M(x) and U(x) represent the max-pooling layer and the Conv-transpose layer, respectively. Similarly, [] denotes the concatenation layer along with the skip connection. H1(x) and H2(x) represent the operations (2conv + Batch norm. + ReLU) and (3conv + Batch norm. + ReLU) respectively.

3.3 Activation function and loss function

The activation function is one of the most crucial learning parameters for any of the deep learning models. The activation function controls the output, accuracy, and computational efficiency of training a model. In our methodology, the ReLU and sigmoid, two activation functions are used for all the models. The ReLU function (equation (3)) is defined as:

The sigmoid function (equation (4)) is defined as:

The binary cross-entropy loss (BCE) is commonly used in binary classification task (equation (5)). This loss function is computed pixel-wise which compares the actual pixel value y _i (i.e., 0 or 1) with the predicted probability of it being 1 termed as p _i . If y _i is 1, it adds log ( p i ) to the loss and if y _i is 0, it adds log ( 1 − p i ) to the loss.

Another loss function (equation (6)) often used is the dice loss function which measures the overlay between the actual (A) and the predicted mask (B).

For our segmentation task, we took BCE dice loss as the loss function (equation (7)) which is the addition of Binary cross-entropy and dice loss function. This is done to get a faster convergence and better performance.

3.4 Dataset and training details

Different deep learning model performances haves been examined and evaluated using the BraTS 2017 data [9]. Multimodal BraTS has always been focusing on multi-institutional pre-operative MRI scans used for segmentation of Glioma brain tumors which are intrinsically heterogeneous in structural appearance, shape, and histology. There are a total of 210 HGG and 75 LGG a patient subjects available in the BraTS 2017 dataset, whereas this study covers only 210 HGG patient datasets counted by us. Due to the low risk and high survival rate of LGG patients as compared to those with HGG [28], these patients are not included in the evaluation. In BraTS 2017, with the help of different clinical protocols and various scanners, four different multimodal MRI scans known as native (T1), post-contrast T1-weighted (T1ce), T2-weighted (T2), and T2 Fluid Attenuated Inversion Recovery (FLAIR) volumes were acquired for each of the patients from 19 multiple institutes. These four types of modalities are taken as input to each of the models separately and the different tumor part is extracted through training. For each of the patients, the ground truth (mask) is given in the dataset verified by experienced neuron-radiologists. In a manually annotated ground truth mask each tumor sub-part: enhancing tumor (ET), the peritumoral edema (ED), and the necrotic and non-enhancing tumor (NCR/NET), is labeled with label 4, label 2, and label 1, respectively, in the referred data set.

The input image slice is trimmed automatically to 192 × 192 after going through data pre-processing steps. As discussed in the data pre-processing phase, only those image slices are considered for the training which carries at least 0.1% of the information of the tumor in 3D volume. The whole dataset of 210 patients is split into 80% for training and 20% for testing: 170 patients’ data were used for training and 40 patients’ data were used for testing. The training image dataset is again divided into training and validation data in an 8:2 ratio randomly. The details of the number of image slices taken for training a model for each of the sub-tumor parts are given in Table 1.

To avoid over-fitting due to less number of images, data augmentation is applied to the whole training dataset. The augmented dataset is created through 90⁰ rotations, horizontal flip, width and height shifting, and zooming of the images, and a shuffled set of image slices are trained through the model in each step of the epoch. To achieve effective and efficient training, the mini-batch size was set to 32 image slices. For the comparison and fine-tuning of the models, the training was carried out for a fixed number of epochs equal to 100. The training was performed using the Google Collaborator platform with GPU back end in a computer having 4GB RAM.

3.5 Evaluation metrics

Several evaluation metrics were used to evaluate the segmented tumor sub-regions with the segmented ground truth. Dice Score compares the overlays between the predicted and the ground truth. Other metrics include accuracy and F1-Score which are the gold standards for classification tasks used to get deeper insights into the predictions.

Dice Score (equation (8)) can be defined as twice the total overlapping area of the predicted mask and the ground truth divided by the sum of the total number of tumor pixels (i.e., pixel value 1 [foreground]) in both the predicted mask and the ground truth).

Other metrics used require pixel-wise evaluation of the ground truth and the predicted masks. True Positive (TP) can be defined as the total number of positive pixels (belonging to the tumor) in the ground truth which are correctly predicted as a tumor, True Negative (TN) is the total number of negative pixels (belonging to the background) in the ground truth which are correctly predicted negative, False Positive (FP) is the total number of negative pixels which are falsely predicted as positive pixels and finally False Negative (FN) is defined as the total number of positive pixels which are falsely predicted as negative pixels. Accuracy (equation (9)) can be defined as the ratio of the number of pixels correctly classified to the total number of pixels.

Specificity (true negative rate) measures (equation (10)) the actual number of pixels that do not belong to the tumor class that is correctly identified.

Sensitivity (recall) is the ratio of the total number of correctly predicted foregrounds to the total number of actual foregrounds. It shows (equation (11)) the amount of foreground predicted correctly out of all the foreground.

F1-Score (equation (12)) is defined as the harmonic mean of precision and recall and is a widely used metric for the classification of imbalanced classes.

4.1 Effect of network structure

The network structure and the number of trainable and non-trainable parameters of the model are shown in Tables 2 and 3. The Adam algorithm is an adaptive learning rate [29] method that computes individual learning rates for each parameter. In all the models, the Adam optimization [30,31] is used with a default learning rate of 0.001. The model was trained to achieve maximum validation accuracy between the predicted map and the ground truth. Our model is not trained using the mean squared error (MSE) function because it is a non-convex function for binary classification and is not guaranteed to minimize the cost function during training. The hidden convolutional layer uses ReLU function whereas the output layer uses the sigmoid activation function to threshold the tumor and background.

The network parameter is equivalent to the encoder–decoder model and SegNet model as both uses the same number of encoding and decoding layers. In SegNet, a total of 24,957,057 parameters are used out of which 24,942,721 are trainable parameters and 14,336 are non-trainable parameters. In U-Net, the total number of parameters is 28,070,753, all are trainable parameters. While training through SegNet, its Parameters require 2.9 GB memory for a batch size of 32 images, while U-Net (parameters) takes 3.35 GB memory for the same batch size if each of the parameter takes 4 bytes of memory. There is a difference of 12 MB memory requirement to train through these different models for a batch size of one input image.

At the same time, we also measured the time taken to train the entire dataset through each of the models. Each modality took 6.8 h on an average for training the model to do segmentation of four subparts: complete, edema, core, and enhancing using the U-Net model. Similarly, the SegNet model takes 16.5 h on an average. and the encoder–decoder model takes 6.7 h on an average for each of the modalities. The U-Net model stores a full encoding feature map which is used in the decoding level with skip connection. Whereas SegNet stores only the pooling indices of feature maps during max-pooling, which in turn requires less memory than U-Net. But through the experiment, it is observed that it takes higher computation time to reconstruct the feature map from pooling indices in the decoding phase.

4.2 Performance analysis

The performance of the segmentation is evaluated through the most common and useful metrics [32] for image segmentation: sensitivity, specificity, F1 score, dice score, and accuracy measurements. In the primary stage of our experiment, we learned that the SegNet model requires less memory as compared to U-Net but the computational time is very large as compared to the other two methods. A confusion matrix (having four values: true positive, false positive, false-negative, and true negative) is also used to compute the specificity and sensitivity of the segmentation results because it allows visualization of the performance of an algorithm. During training, the confusion matrix for each of the sub-tumor segmentation is computed. Widely used segmentation metrics are sensitivity measures correctly classifying a pixel that belongs to the tumor part, and specificity measures correctly classifying a pixel that does not belong to the tumor part. These two metrics along with the accuracy of the test image dataset is evaluated by using a confusion matrix. The detailed metrics for each of the models is shown in Table 4.

The above metrics calculate the efficiency of the algorithm regarding correctly classifying a pixel if it belongs to the tumor region. Next we have considered the error of prediction i.e., the analysis is also carried out by measuring the mean squared error of the test dataset. In statistics, the MSE of the model measures the average squared difference between the predicted and actual output. The algorithm having a low MSE value will perform better as compared to another model. MSE is measured for each of the predicted masks with respect to the original ground truth. The overall MSE for a model is computed by taking the average of all individual measurements and found that U-Net is having the lowest MSE value of 0.017026377 whereas SegNet is having the maximum MSE value of 0.021142694. The encoder and decoder model is having an average MSE value of 0.018796785, nearly equal to U-Net model.

Dice score is considered as one of the most common metrics to evaluate the performance of the segmentation algorithm. The model is evaluated here by using validation dice scores which are considered to be maximized during training instead of MSE. The dice score is measured for individual images and then computed as the overall dice score for a set of validation data as shown in Table 5 below. It is clearly shown that all the models are having maximum DSC for extracting complete sub-tumor using FLAIR modality, core and enhancing sub-parts using T1ce, and edema region using either FLAIR or T2 weighted modality. But U-Net is performing better as compared to the other two models as per the experimented data.

The research also demonstrates the performance in terms of the Receiver Operating Characteristic (ROC) curve taken between false positive rate vs true positive rate and measures the area under the curve (AUC value) (Figure 2) to analyze the performance of classification of the pixel in detail. The AUC can be used to compare the performance of two or more classifiers in absolute terms: the higher the AUC, the better the segmentation/classification. The details of the AUC value for each of the sub-tumor segmentation through different modalities are shown graphically.

Figure 2

Area under ROC values shown for each of the model for segmenting different sub-tumor parts using FLAIR, T1 weighted, T2 weighted, and T1ce modalities.

There is a similar trend in AUC values for each of the models during segmenting the sub-tumor region using FLAIR, T1 weighted, T2 weighted, and T1ce modalities. For complete and edema region extraction, FLAIR and T2 weighted modalities perform better as compared to T1 weighted and T1ce. Similarly, for core and enhancing, the performance is better in T1ce as compared to the other three modalities. The result also shows the model efficiency in terms of the ROC curve (AUC) for the sub-tumor part with respect to the different modalities. For overall performance of segmentation, from above all the graphs it is clearly shown that the U-Net model is better as compared to the other two models and SegNet (pooling indices) result is satisfactory in comparison to the state of the art but performance is lower than U-Net model.

The segmentation results are demonstrated in Figure 3 on different test images. Figure 3(a–c) depict some semantic segmentation results in visual form using the U-Net model, SegNet model, and encoder–decoder model, respectively, from an axial view. It shows the different sub-tumor parts extracted from FLAIR, T1 weighted, T2 weighted, and T1ce images along with the original ground truth images.

Figure 3

Segmentation results of U-Net, SegNet, encoder–decoder model in (a)–(c). From leftmost, images are: MRI scans, ground truth, tumor-Flair, tumor-T1, tumor-T1ce, and tumor-T2, respectively.

The U-Net model and encoder–decoder model achieve maximum F1 score of 0.94 and 0.93, respectively, for the complete tumor segmentation using FLAIR modality, whereas in SegNet model, maximum of 0.93 F1 score is achieved using T2 weighted modality. Detailed value of F1 score for each of the models is given in Figure 4 represented graphically.

Figure 4

F1 score measurement for all three models from FLAIR, T1 weighted, T2 weighted, and T1ce image.

All the models achieved maximum F1 score for core and enhancing tumor segmentation using T1ce modality (core: 0.95, 0.92, and 0.94 and enhancing: 0.93, 0.92, and 0.91 for U-Net, SegNet, and encoder–decoder models, respectively). For edema region segmentation, U-Net and SegNet achieves 0.87 and 0.84, respectively, using FLAIR modality but encoder–decoder model achieves 0.84 F1 score for both FLAIR and T2 weighted images.

4.3 Statistical analysis

Even if accuracy and F1 score are the most common metrics used in statistics to measure the segmentation performance, sometimes both can be misleading. Both the metrics do not fully consider the size of TP, TN, FP, and FN classes of the confusion matrix in their final computation. To avoid these misleadings and justify our performance comparison, another performance score, the Matthews correlation coefficient (MCC) [35], is considered. It is given in terms of TP, TN, FP, and FN (equation (13)) as follows:

Finally, we conclude our analysis through the MCC metric which is given in Table 6. Along with computed F1 score and accuracy, MCC results also confirm the improvement in the prediction of tumors and show that our segmentation results are very satisfactory with respect to the ground truth.

The statistical analysis is performed through the ANOVA test by using the MCC coefficient value and considering the significance level (α) equal to 0.05. The null hypothesis taken as the mean of each group is equal or there is no significant difference between the performances. From the p values measured for the performance with respect to modalities for U-Net (=0.464793), SegNet (=0.4223184), and encoder-decode model (=0.32607016), it is clearly shown that there is no such evidence to reject the null hypothesis as all the values are greater than 0.05. Hence, there is no significant difference in segmentation results for the different multimodal MRI scans. Similarly, when the test is performed to analyze the segmentation result for different sub-tumor classes, the p-value is computed as 0.019135 for U-Net, 0.03175 for SegNet, and 0.00730 for the encoder–decoder model. All the values are less than 0.05 which indicates there is a significant difference in performance results when segmenting the different sub-parts of the tumor. This analysis is performed between the groups of individual models. Then, this analysis is extended to measure the performance difference between the models as shown in Table 7.

From the ANOVA test, we got a p-value equal to 0.83003 for the performance of three models for segmenting tumors. As it is clearly shown from the different evaluation metrics that there is no huge gap in the segmentation result while comparing the three CNN models, the test also cooperates with the fact. It shows there is no significant difference in segmentation results between the three models. When it comes to the individual sub-tumor parts segmentation or whole tumor segmentation, statistics show that there is a difference between the model performances as the p-value is very much less than 0.05.

In the above section, we compare the three different models through different segmentation metrics to analyze the performance. But no single model is performing better in comparison to the other two for segmenting all the sub-tumor parts. Even if there is no such single MRI modality which will give better accuracy as compared to multi-parametric MRI scans. Each model has its own architecture which will reduce either computational time or memory usage. Similarly, each modality has different characteristics and contrast images which help in segmenting different sub-tumor from different modalities with the highest accuracy.

Hence in ensemble models, the maximized feature maps extracted based on MRI scans and models are fused together to segment the whole tumor of a patient which is time and memory-efficient in some contexts. In the second layer of the proposed model, it selects the model with the corresponding modality for predicting the sub-part of the tumor. The visual output of the model is depicted in Figure 5. The result of the proposed ensemble model is compared with other existing results as shown in Table 8.

Figure 5

Segmentation results of hybrid models: (a)–(g) MRI slices, ground truth, enhancing part, edema region, core tumor, complete part, and predicted whole mask, respectively.