2.1 Data Set

To show how the proposed method performs we used the longitudinal multiple sclerosis lesion segmentation data set [³].^fn The training set contains 4-5 volumes of multi-modal MR images of five patients.

The available modalities were acquired on a 3T MRI scanner and corresponds to T1-weighted, T2-weighted, PD-weighted and FLAIR.

Successive volumes were acquired in time intervals that allow to reflect disease progression. For each image, the corresponding manual lesion segmentation is also provided.

2.2 Preprocessing Pipeline

The preprocessing steps that were applied correspond to the standard preprocessing pipeline typically performed to prepare MR images, i.e. correction of the MR intensity inhomogeneities, skull-stripping, dura mater stripping, and rigid-body registration to a 1 mm isotropic MNI template [³].

In addition to the previously described standard MRI preprocessing pipeline, a series of additional operations were performed with the purpose of generating high quality image patches to train the CNN model.

These operations are summarized below. The image intensities were truncated to the quantiles within the range [0.01,0.9995], and then scaled to the range [0,1].

The truncation quantiles were determined experimentally in such a way that they allow to achieve better correspondence between MR image intensities, especially in high intensity voxels in T2 and FLAIR images which are potentially part of MS-related lesions.

We found that this process allows to achieve results even better than other approaches used to maximize the contrast between image intensities (e.g., histogram matching).

The quantiles were determined experimentally in this work. We found values relatively close to quantiles found by other researchers [²]. To train the CNN model, around of 650.000 patches were generated from the available images.

This was done by selecting image patches that contain voxels that satisfy two conditions: (i) the patch contains hyper-intense voxels on FLAIR images which have the potential of being a MS-related lesion [⁵], and (ii) the patch contains white matter.

To satisfy the condition (i) a set of supra-threshold voxels were determined from the FLAIR images, and to satisfy the condition (ii) the subset of voxels belonging to white matter were determined.

This was done by registering the FLAIR images to the ICBM452 probabilistic atlas using a rigid-body registration algorithm implemented in FSL [⁹].

2.3 Network Architecture and Training

The CNN used has 7 convolutional layers. The input of the network consists of three patches of size 33×33. The selected CNN architecture was initially inspired by the architecture presented in [²], but the number of hidden layers was increased.

This increment of the model’s hidden layers was motivated by our empirical findings in terms of the Dice similarity criterion. The inputs of the model correspond to patches generated from the T1-weighted, T2-weighted, and FLAIR images.

The PD-weighted images were discarded as they do not contribute to the CNN model’s performance according to our empirical study. Each convolutional layer (except for the output layer) uses Leaky ReLU activation function with negative slope coefficient β=0.3.

The output layer corresponds to a 1-dimensional sigmoidal activation function that retrieves the probability of a voxel belonging to the MS-related lesion class.

By evaluating a whole image, a probability map can be generated with the purpose of solving the segmentation problem as a voxel-to-voxel binary classification task. The network contains 2 layers of MaxPooling of 2×2 with stride 2.

Dropout after each convolutional layer was used as regularization method. In order to increase the depth of the network without generating adverse effects such as the problem of the vanishing gradient, a single ResNet block was included into the model.

We empirically found that by increasing the number of ResNet blocks the overall model’s performance is not improved. Figure 1 summarizes the network architecture. The image sampling performed to generate patches is highlighted with green squares.

Fig. 1 Red blocks: 2×2 maxpooling layers, Blue blocks: convolutional layers, Yellow blocks: Leaky-ReLU activation functions, except the classification layer use sigmoidal activation. One residual connection present in the network 

Considering that our CNN model has a one dimensional output, the loss functions used for training was the binary cross entropy and focal loss with parameter values γ=1.0, γ=2.0, γ=3.0 and γ=5.0.

In addition, a mechanism to improve the CNN model performance in difficult voxels near the edges of the lesions was implemented. This mechanism consisted in balancing the number of patches of the non-lesion class and lesion class with the ratio 1:1.

A stratified sampling process was performed considering a proportion p for voxels located at a distance not greater than 4 mm to the edges of the MS lesions and a proportion 1−p for voxels randomly sampled from voxels candidates that are distant more than 4 mm to the MS lesions. The training of the CNN model was made by using the stochastic gradient descent strategy.

The learning rate selected was 0.01 with decay learning rate of 1×10−6, the momentum parameter was 0.9 with Nesterov momentum and batch size of 64. Dropout with probability 0.1 was used.

The implementation and training of the CNN was done using the Keras Deep Neural Network Library [⁴] with TensorFlow as backend numerical engine. The model was trained using NVIDIA GeForce GTX 1080TI Graphic Card.

To compare the segmentation performance, we implemented a group five-fold cross validation with five identical models. For each fold, the training and testing sets were determined is such a way no patches of the same patient belongs to the testing and training sets simultaneously.

The details are presented in Table 1. The number of epochs was determined in such a way that the Dice similarity coefficient reaches a stable value in the validation set (a total of 100 epochs was applied for training, see Figure 2).

Fig. 2 Evolution of the Dice similarity coefficient during training using a validation example 

2.4 Focal Loss Function

Our model was trained using the focal loss function that was proposed in [¹⁰], with the purpose of addressing the class imbalance problem, specifically in CNN models trained in one-stage.

The focal loss function is based on the binary cross entropy loss function, which is given by:

where y∈{0, 1} is the ground truth, y=0 and y=1 correspond to the non-MS and MS-lesion classes, respectively, and y^∈[0, 1] is the probability of the class with label y=1 estimated by the model. For notation convenience, y^t is defined as:

Thus, for the t-th input instance, the BCE can be rewritten as: BCE(y^t)=−log⁡(y^t).

In the focal loss function, a weighting term given by (1−y^t)γ is added to the binary cross entropy loss function. According to this, the focal loss function can be expressed as:

where γ is a parameter that modulates the weighting effect provided to this loss function.

The focal loss weights dynamically the training examples losses, down-weighting easy-to-classify training examples and assigning higher losses to examples with low prediction accuracy, i.e. hard-to-classify training examples.

3.1 Evaluation

The accuracy of the CNN model was computed by comparing the segmentation mask predicted by the model S and the gold standard G manually segmented by the expert. The evaluation metrics used in this work aim to highlight the spatial-overlap and spatial-distance and are defined below [¹¹].

Spatial-overlap metrics:

where TP, FN and FP correspond to the amount of true positive, false negative and false positive values.

Spatial-distance metrics:

The Hausdorff distance is usually adversely affected by the presence of outliers since in its construction it includes the max statistic, which is non-robust. Due to this fact, the Hausdorff distance and the percentile 95% are usually used in the literature.

where d(x,∂S)=min⁡y∈∂G‖x−y‖2 and d(y,∂G)=min⁡x∈∂S‖x−y‖2, |∂G| and |∂S| denote the number of points on the surfaces of G and S, respectively.

3.2 Results and Discussion

Table 2 presents the results obtained by our proposed CNN network (Figure 1) using stratified sampling together with the binary cross entropy loss functions (BCE) and focal loss (FL).

The parameter s corresponds to the proportion of voxels sampled near to the lesions edges respect to the total number of voxels of the non-lesion class and γ is the focal loss parameter (see Eq. (3)).

The metrics used correspond to, Hausdorff distance (Eq. ((7)), Average symmetric surface distance (Eq. (9)), Sensivity (Eq. (5)), Precision (Eq. (6)) and Dice (Eq. (4)).

Regarding the results obtained by the spatial-distance metrics, it can be noted from Table 2 that the HD and ASSD metrics are sensitive to the presence of outliers, due to the fact that the mere presence of a lesion not detected by CNN or a lesion detected erroneously can generate a large number of voxels greatly distant from the lesions in the gold standard.

That is why this metric is often used in evaluating the image segmentation of big and high-contrast organs. For metrics based on spatial-overlap, the best result about the sensitivity metric corresponded to 0.7143(0.0127) using s=0.2 and γ=1.0.

For the precision metric, the best result was 0.7593(0.0146) obtained using the parameters s=0.2 and γ=3.0. These results show that this methodology reduces false negative and false positive errors.

In the particular case of false positive errors they are reduced thanks to the stratified sampling strategy that increases the representation of difficult voxels close to the MS lesions.

The best results in terms of the Dice similarity coefficient corresponded to 0.7150(0.0094), obtained using s=0.2 and γ=3.0, which indicates that using focal loss, the increased representation of difficult-to-classify voxels by the stratified sampling strategy will affect positively the model’s performance compared to uniform sampling, which has s=0.026 (in average).

The best result in the Dice similarity coefficient was compared with all other combinations of parameters by means of a t-test with the assumption of different and unknown variances (Welch’s test), where the null hypothesis corresponds to the assumption that there are no differences of the Dice means with the different combinations.

The alternative hypothesis is that this result has a higher Dice mean. It must be noted in Table 2 that the p-values were less than 0.01 when using a uniform sample (s=−), i.e. a non-stratified sample, independent of γ parameter of the focal loss.

A low p-value is also observed in stratified samples s=0.2 and s=0.4 only when binary cross entropy loss (BCE) since using focal loss increases the value of the Dice similarity coefficient, thus decreasing the power of the test (p-value>0.1), except for s=0.4 where FL has a marginal incremental effect.

We also observed that the use of focal loss drastically influences the decrease of the binarization threshold (bt) necessary to maximize the results of the spatial-distance metrics and the Dice similarity coefficient, especially when stratified sampling was used.

This increases the proportion of non-lesion voxels close to the edges MS-lesions in the training set, which supports the hypothesis that these areas contain most difficult-to-classify voxels.

Some segmentation results are shown in Figure 3. Axial FLAIR MRI views of a patient with multiple sclerosis are presented with the segmentation mask obtained with our proposed algorithm overlayed. True positive are represented in red, false negative errors are presented in green and false positive errors in yellow.

Fig. 3 Automatic segmentation: (a)-(d) pacient 1 (e)-(f) pacient 5 

Figure 3 (a)-(d) corresponds to images from patient 1 and (e)-(h) corresponds to images from patient 5. In Figure 3 (a) and (e) the parameters are s=0.0 and γ=0.0 and in Figure 3 (b) and (f) the parameters are s=0.0 and γ=3.0. It can be noted that in Figure 3 (b) there was a slight decrease in false negative errors which is attributable to the application of focal loss.

In Figure 3 (c) and (g) the results of automatic segmentation are presented using the parameters s=0.2 and γ=0.0 and in Figure 3 (d) and (h) s=0.2 and γ=3.0. In Figure 3 (d) and (h) an increase in true positive and a decrease in false negative errors and false positive errors with respect to all the previous combinations of parameters. This is attributable to the use of the stratified sampling mechanism and the proper of the γ parameter.