# Using Convolutional Neural Networks for Segmentation of Multiple Sclerosis Lesions in 3D Magnetic Resonance Imaging

Authors: Abdulsattar Abdullah Hamad; Mustafa Musa Jaber; Mohammed Altaf Ahmed; Ghaida Muttashar Abdulsahib; Osamah Ibrahim Khalaf; Zelalem Meraf
Journal: Advances in Materials Science and Engineering
 (2022)

---

## Abstract

Magnetic Resonance Imaging to detect its lesions is used to diagnose multiple sclerosis. Experts usually perform this detection process manually, but there is interest in automating it to speed up the diagnosis and monitoring of this disease. A variety of automatic image segmentation methods have been proposed to quickly detect these lesions. A Gaussian Mixture Model is first constructed to identify outliers in each image. Then, using a set of rules based on expert knowledge of multiple sclerosis lesions, those outliers of the model that do not match the lesions' characteristics are discarded. Furthermore, segmented lesions usually correspond to gray matter-rich brain regions. In some cases, false positives can be detected, but the rules used cannot eliminate all errors without jeopardizing the segmentation’s quality. The second method involves training a convolutional neural network (CNN) that can segment lesions based on a set of training images. This technique can learn a set of filters that, when applied to small sections of an image called “patches,” produce a set of characteristics that can be used to classify each voxel of the image as a lesion or healthy tissue. On the other hand, the results show that the networks are capable of producing results in the worked database comparable to those produced by the algorithms in the literature.


---

## Body

## 1. Introduction

Magnetic resonance imaging (MRI) is a medical imaging technique that allows obtaining images of different parts of the body, such as the brain, heart, or spine, helping diagnose various pathologies that may be affecting a patient [1]. This technique has many advantages over other brain scanning techniques, including the use of nonionizing radiation and the generation of images with a high degree of contrast between the different types of tissues present in the brain (cerebrospinal fluid, white and gray matter), which makes it an excellent support for diagnosing diseases that affect the central nervous system [2], such as epilepsy and dementia.Among these diseases is multiple sclerosis, a chronic disease that causes damage to the central nervous system, characterized by the appearance of regions of inflammation and demyelination of neurons, which result in a progressive loss of the patient’s motor, sensory, and cognitive abilities affected. Through MRI, it is possible to detect the regions affected by this disease and make a diagnosis based on the number and size of the lesions detected [3].This task is generally performed by an expert medical professional or semiautomatically, using computational image segmentation techniques [4], which require parameter adjustment by the expert. This approach tends to present difficulties and disadvantages, such as high variability in the results issued by different experts. Due to this, the medical community’s interest arises in developing automatic segmentation techniques to reduce the time required to identify damaged regions in the brain. There are many challenges in automatic segmentation problems, and their use in practice is quite limited due to the difficulty of correctly validating the predictions made by these tools. In this memory, the effectiveness of different segmentation methods for this type of image will be studied, proposing the use of convolutional neural networks for the segmentation of lesions and comparing them with traditional segmentation methods, analyzing the advantages and disadvantages that each of these methods they bring into practice [5].This work aims to present different techniques of the literature in the area of image segmentation. Then, the implementation of these and their comparison in terms of performance and effectiveness in the segmentation problem of magnetic resonance images will be carried out. To validate results, the set of Boston Children’s Hospital (BCH) training images of the MS Lesion Segmentation Challenge will be used [6], which consists of 10 images in T1, T2, and FLAIR modalities of different patients, accompanied by the segmentation of an expert. Because the segmentation of the test set is not publicly accessible, for supervised learning algorithms, in this case, convolutional neural networks, the set of images will be divided into a training subset, from which the model will be trained of convolutional networks, and a test set, with which the ability of the model to segment lesions in a new image correctly will be tested. This process will be performed five times per technique and an average of the performance measurements will be obtained. The current study aims to experimentally evaluate the performance of convolutional neural networks applied to the problem of automatic segmentation of multiple sclerosis (MS) lesions in multimodal magnetic resonance imaging, comparing it with other literary techniques.


## 2. Methodology

The proposed solution to the problem of segmentation of multiple sclerosis lesions in magnetic resonance images consists of implementing and evaluating two different methods applied in the literature to the problem. The dataset with which the proposal will be validated corresponds to the database of training images used in the multiple sclerosis lesion segmentation challenge organized by the Medical Image Computing and Computer-Assisted Intervention Society [7]. This database is available to the public and comprises magnetic resonance images of twenty different patients in 3 different modalities: T1-Weighted, T2-Weighted, and T2-Fluid Attenuated Inversion Recovery (FLAIR), with their respective masks of segmented lesions. Ten of these images in their 3 modalities, with their respective segmentations, were provided by the University of North Carolina (UNC). The other ten were provided by Children’s Hospital Boston (BCH).The modalities of each image of the data set are coregistered, that is, the voxels of each image correspond to the same spatial coordinate of the brain. However, they lack a greater preprocessing, for which a preprocessing stage is necessary before the analysis application of segmentation algorithms to obtain good results.Once the images have been preprocessed, two literature segmentation algorithms will be implemented. Once state of the art was evaluated, it was decided to implement a technique based on the Gaussian Mixture Models [8] based on the search for outliers in the model and the use of rules to eliminate false positives from the segmentation. On the other hand, the implementation of a convolutional neural network was chosen to segment lesions [9]. This method has become very popular in the literature in recent years.To compare the effectiveness of both methods, the Sensitivity (True Positive Rate) and Precision (Positive Predictive Value) metrics will be used, which are given by(1)TPR=TPTP+FN,PPV=TPTP+FP,where TP is the number of true positives from segmentation (correctly labeled lesions), FP is the number of false positives from segmentation (healthy regions that were labeled as lesions by segmentation), and FN is the number of false negatives from segmentation (lesions that were labeled as healthy tissue). Each of these measurements is quantified at the voxel level. Sensitivity represents the percentage of the correctly segmented lesion, while Precision measures the percentage of hits concerning the segmented region. A high Sensitivity indicates that most of the lesion voxels are correctly segmented. In contrast, a low Sensitivity indicates that lesions are not being detected or only a tiny portion of them are being detected. On the other hand, a high Precision indicates that a tiny percentage of the voxels segmented as lesions correspond to healthy tissue. In contrast, a low Precision indicates that the percentage of healthy tissue segmented as lesions is very high due to a large number of false positives. Finally, the Sorensen Dice similarity (DSC) [10] coefficient will be used, a standard metric used in the literature that allows determining the level of overlap between the obtained segmentation and the segmentation provided by an expert, which is given by the harmonic mean between the two metrics.(2)DSC=21/TPR+1/PPV,=2TP2TP+FP+FN.


## 3. Construction and Implementation of the Solution

### 3.1. Preprocessing

To obtain good results in the segmentation algorithms, it was necessary to perform a preprocessing stage to eliminate image noise and other irregularities present.


#### 3.1.1. Bias Field Correction

MRI images fall victim to a signal that corrupts image quality known as the bias field. This signal is the product of inhomogeneities in the magnetic field used to obtain these images, and its presence results in blurred sectors in the image, adding noise to the intensity values of the image, which negatively affects the quality of the results obtained by algorithms of image processing and segmentation since these disturbances corrupt the distribution of intensities of the different objects and regions that appear in the image. To correct this signal, the N4ITK algorithm was used.The N4ITK algorithm consists of a modification of the N3 algorithm (nonparametric nonuniformity normalization) that seeks to resolve some algorithm’s shortcomings [11]. Given an image, the algorithm models the bias field as a function that multiplies the image without corruption as(3)vx=uxfx+nx,where v (x) is the original image; u (x) is the image without corruption; f (x) is the bias field; n (x) is the noise (assumed to be independent Gaussian noise).Using the notation uˆ = logu, the previous equation applying logarithm is of the form:(4)v^x=u^x+f^x.From this point, the image without corruption is iteratively estimated by(5)u^n=u^n−1−f^rn,=u^n−1−S∗u^n−1−Eu^|u^n−1,where S∗• is an approximation by B-Splines of the bias field residual between the nth and n − 1th iterations. Details of this implementation can be found in [11].


#### 3.1.2. Skull Stripping

To avoid involving the voxels of the images corresponding to the background and the skull, the skull stripping task was performed to obtain a binary mask that only included those voxels corresponding to brain tissues [12]. For this task, BrainSuite software was used. The skull stripping task was performed using the BSE (Brain Surface Extractor) algorithm. Parameter selection was performed manually for each image until an appropriate mask was obtained that did not include skull voxels. The parameters that produced satisfactory results for each of the images in the BCH database were seen in Figure 1 and Table 1.Figure 1
Original FLAIR image, image after applying bias field correction, and image after applying skull stripping.Table 1
Parameters used in the BCH.
ImageEdge constantErosion sizeBCH010.72BCH020.72BCH030.72BCH040.92BCH050.63BCH060.73BCH070.72BCH080.73BCH090.73BCH100.73


### 3.2. Gaussian Mixture Models (GMM)

As in Karimian and Jafari [13] paper, the first method chosen for segmentation is modeling the intensities of the three types of brain tissue using 3-component Gaussian Mixture Models. The segmentation will be done by looking for outliers in the model using the truncated maximum likelihood estimator of the GMM, then applying a set of rules to eliminate false positives.The algorithm proposed by Lee and Song [14] is dependent on several parameters, including the GMM's percentage of outliers and the thresholds used to eliminate them. To choose these parameters, we vary the percentage of outliers in the model h between {0; 0.2; 0.3; 0.4; 0.5}. Then, the Mahalanobis distance, T2 image intensity, and FLAIR image intensity thresholds were varied for each model. The Mahalanobis distance ranged from 0.4 to two. The parameters chosen correspond to the highest DSC combination.A problem that was observed when adjusting the GMMs for the different images in the database is that the intensity distributions for white matter in the T2 and FLAIR images varied a lot between the different images, for which the use of a threshold single fixed for T2 and FLAIR images did not produce good results in segmentation. To solve this problem, an idea similar to that proposed in Jansen [15] of using thresholds based on parameters of the obtained model was taken. In this way, the regions of hyperintensity in the white matter of the T2 and FLAIR images are selected by the thresholds µWM, T2 + kT2 • σWM, T2 and µWM, FLAIR + kFLAIR • σWM, FLAIR, where µWM, i is the mean of the white matter intensities in the image in modality i, σWM, i is the standard deviation of the white matter intensities in the image in modality i, and kT2 and kFLAIR are 2 constants associated with the threshold of their respective image. These last 2 parameters replace the fixed threshold defined by Karimian and Jafari and were made to vary in the range of values {0, 1, 2}. Finally, hypointensity in the T1 image, defined by the negative threshold µWM, T1 − kT1 • σWM, T1, was also added to the set of rules. The parameter kT1 was varied in the same range as kT2 and kFLAIR.Another issue with the method described in the paper is that most predicted lesions were located at the outer edge of the white matter, near the cortical gray matter. The solution proposed by Zhao et al. [16] is to apply a mask obtained from the morphological erosion operation to the white matter mask, allowing the elimination of nearby voxels from the lesion candidates. The internal holes of the white matter mask were filled before applying the erosion mask. The chosen erosion filter was 6 × 6 × 6. A 3 × 3 erosion filter was applied to each slice of their set of images, whose resolution was half that of the images with which they worked, which is why the filter size was initially set to 6 6 6. These voxels were segmented as a lesion and their final size was set to 9 × 9 × 9, see Table 2.Table 2
Parameters of means obtained for the white matter.
Image no.µWM, T2µFLAIR, T2BCH0177.8973.7BCH0237.7662.18BCH0338.0880.85BCH0446.4100.96BCH0545.2663.65BCH0650.0161.64BCH0753.8859.78BCH0843.8860.7BCH0958.3857.84BCH1058.3356.88The initial segmentation to obtain the estimators of the initial parameters of the model was performed on the T1 image using theK-means algorithm.


#### 3.2.1.K-Means Algorithm

TheK-means algorithm is an unsupervised learning algorithm commonly applied in clustering tasks. Given a set of points X = {x1, x2, …, xn} and a parameter k, corresponding to a positive integer, the algorithm is initialized by choosing k random points c1, c2, …, ck within the same domain of X known as centroids. Then, for each point within the set X, the distance between xi and each of the centroids cj is calculated (commonly, the Euclidean distance is used) and the point xi is assigned to the cluster corresponding to the closest centroid. Finally, the centroids are calculated again as(6)Cj=1nj∑xi∈Cj.xi,where Cj is the cluster containing the previously closest points to Cj and nj is the number of points in the cluster Cj. These 2 steps are repeated iteratively until the cluster assignment for each point does not change from one iteration to the next.This algorithm applies to the image segmentation task. Given the set of pixels of a grayscale image, the points to be grouped correspond to the intensities of the pixels. In contrast, the centroids correspond to different intensity values in the range [0, 255]. As a result, theK-means algorithm performs a grouping of pixels with similar intensity values. The K-means algorithm has been used to perform the segmentation of different brain tissues in magnetic resonance images [17].The K-means algorithm was used to perform a preliminary segmentation of the brain in the 3 different types of tissue based on the intensity of the pixels in the T1 image, obtaining, from lower to higher intensity, labels for the pixels corresponding to cerebrospinal fluid, spinal fluid, and gray and white matter.


### 3.3. Convolutional Neural Networks (CNNs)

Second, the first Birenbaum architecture [18] is used to train a convolutional network using a subsampling of FLAIR image patches. Instead of a longitudinal MRI, the paper’s SISTP (Single Image, Single TimePoint) architecture was used. In addition to the proposed architecture, an architecture without the Max pooling layer following the convolution layers was tested to compare performance. Both architectures tested 2 input types: 1 axial image and 1 image for each view type (axial, sagittal, and coronal). Each slice was obtained with a resolution of 0.5 mm and a distance of 1 mm between each slice. The intermediate values were interpolated during the registration process to obtain the final image with a resolution of 0.5 mm. It is then passed to a dense network of 16 neurons followed by the Softmax layer, which concatenates the 48 feature vectors obtained after the first dense layer of each image.The implementation of the network was carried out using Keras [19], a library written in Python for the management and construction of neural networks in a modular way, using the Tensorflow library [20] as a backend.In Figure2, candidate voxels were chosen based on preliminary white matter segmentation from the T1 image using the K-means algorithm. This white matter mask replaced the segmentation obtained from registering an Atlas to a T1 image. WM, FLAIR + kFLAIR • WM, FLAIR was used to define the threshold parameter for the FLAIR image [21]. The parameter was set to match the proportion of candidate voxels in Birenbaum et al. The final number of candidates for the network was 1,242,109 voxels (kFLAIR = 2), which is close to the percentage reported in Birenbaum’s paper, who claims that candidate sampling reduced network computations by 97.5 percent [20]. This subsampling leaves out a percentage of the image voxels corresponding to lesions. In addition to patch sampling, a loss function that penalized false negatives for lesions addressed the class imbalance issue. The weight parameter of this loss function was set at 20% for all trained networks based on the percentage of injuries observed in the training set. We will use leave-one-out cross-validation to evaluate the network's results. The network’s model is then applied to the untrained image, yielding model evaluation metrics for that image. This is done for every image in the training set. Finally, the obtained model’s evaluation metrics are averaged.Figure 2
Architecture for 1 implemented input of the CNN with and without Max pooling layer.


## 3.1. Preprocessing

To obtain good results in the segmentation algorithms, it was necessary to perform a preprocessing stage to eliminate image noise and other irregularities present.


### 3.1.1. Bias Field Correction

MRI images fall victim to a signal that corrupts image quality known as the bias field. This signal is the product of inhomogeneities in the magnetic field used to obtain these images, and its presence results in blurred sectors in the image, adding noise to the intensity values of the image, which negatively affects the quality of the results obtained by algorithms of image processing and segmentation since these disturbances corrupt the distribution of intensities of the different objects and regions that appear in the image. To correct this signal, the N4ITK algorithm was used.The N4ITK algorithm consists of a modification of the N3 algorithm (nonparametric nonuniformity normalization) that seeks to resolve some algorithm’s shortcomings [11]. Given an image, the algorithm models the bias field as a function that multiplies the image without corruption as(3)vx=uxfx+nx,where v (x) is the original image; u (x) is the image without corruption; f (x) is the bias field; n (x) is the noise (assumed to be independent Gaussian noise).Using the notation uˆ = logu, the previous equation applying logarithm is of the form:(4)v^x=u^x+f^x.From this point, the image without corruption is iteratively estimated by(5)u^n=u^n−1−f^rn,=u^n−1−S∗u^n−1−Eu^|u^n−1,where S∗• is an approximation by B-Splines of the bias field residual between the nth and n − 1th iterations. Details of this implementation can be found in [11].


### 3.1.2. Skull Stripping

To avoid involving the voxels of the images corresponding to the background and the skull, the skull stripping task was performed to obtain a binary mask that only included those voxels corresponding to brain tissues [12]. For this task, BrainSuite software was used. The skull stripping task was performed using the BSE (Brain Surface Extractor) algorithm. Parameter selection was performed manually for each image until an appropriate mask was obtained that did not include skull voxels. The parameters that produced satisfactory results for each of the images in the BCH database were seen in Figure 1 and Table 1.Figure 1
Original FLAIR image, image after applying bias field correction, and image after applying skull stripping.Table 1
Parameters used in the BCH.
ImageEdge constantErosion sizeBCH010.72BCH020.72BCH030.72BCH040.92BCH050.63BCH060.73BCH070.72BCH080.73BCH090.73BCH100.73


## 3.1.1. Bias Field Correction

MRI images fall victim to a signal that corrupts image quality known as the bias field. This signal is the product of inhomogeneities in the magnetic field used to obtain these images, and its presence results in blurred sectors in the image, adding noise to the intensity values of the image, which negatively affects the quality of the results obtained by algorithms of image processing and segmentation since these disturbances corrupt the distribution of intensities of the different objects and regions that appear in the image. To correct this signal, the N4ITK algorithm was used.The N4ITK algorithm consists of a modification of the N3 algorithm (nonparametric nonuniformity normalization) that seeks to resolve some algorithm’s shortcomings [11]. Given an image, the algorithm models the bias field as a function that multiplies the image without corruption as(3)vx=uxfx+nx,where v (x) is the original image; u (x) is the image without corruption; f (x) is the bias field; n (x) is the noise (assumed to be independent Gaussian noise).Using the notation uˆ = logu, the previous equation applying logarithm is of the form:(4)v^x=u^x+f^x.From this point, the image without corruption is iteratively estimated by(5)u^n=u^n−1−f^rn,=u^n−1−S∗u^n−1−Eu^|u^n−1,where S∗• is an approximation by B-Splines of the bias field residual between the nth and n − 1th iterations. Details of this implementation can be found in [11].


## 3.1.2. Skull Stripping

To avoid involving the voxels of the images corresponding to the background and the skull, the skull stripping task was performed to obtain a binary mask that only included those voxels corresponding to brain tissues [12]. For this task, BrainSuite software was used. The skull stripping task was performed using the BSE (Brain Surface Extractor) algorithm. Parameter selection was performed manually for each image until an appropriate mask was obtained that did not include skull voxels. The parameters that produced satisfactory results for each of the images in the BCH database were seen in Figure 1 and Table 1.Figure 1
Original FLAIR image, image after applying bias field correction, and image after applying skull stripping.Table 1
Parameters used in the BCH.
ImageEdge constantErosion sizeBCH010.72BCH020.72BCH030.72BCH040.92BCH050.63BCH060.73BCH070.72BCH080.73BCH090.73BCH100.73


## 3.2. Gaussian Mixture Models (GMM)

As in Karimian and Jafari [13] paper, the first method chosen for segmentation is modeling the intensities of the three types of brain tissue using 3-component Gaussian Mixture Models. The segmentation will be done by looking for outliers in the model using the truncated maximum likelihood estimator of the GMM, then applying a set of rules to eliminate false positives.The algorithm proposed by Lee and Song [14] is dependent on several parameters, including the GMM's percentage of outliers and the thresholds used to eliminate them. To choose these parameters, we vary the percentage of outliers in the model h between {0; 0.2; 0.3; 0.4; 0.5}. Then, the Mahalanobis distance, T2 image intensity, and FLAIR image intensity thresholds were varied for each model. The Mahalanobis distance ranged from 0.4 to two. The parameters chosen correspond to the highest DSC combination.A problem that was observed when adjusting the GMMs for the different images in the database is that the intensity distributions for white matter in the T2 and FLAIR images varied a lot between the different images, for which the use of a threshold single fixed for T2 and FLAIR images did not produce good results in segmentation. To solve this problem, an idea similar to that proposed in Jansen [15] of using thresholds based on parameters of the obtained model was taken. In this way, the regions of hyperintensity in the white matter of the T2 and FLAIR images are selected by the thresholds µWM, T2 + kT2 • σWM, T2 and µWM, FLAIR + kFLAIR • σWM, FLAIR, where µWM, i is the mean of the white matter intensities in the image in modality i, σWM, i is the standard deviation of the white matter intensities in the image in modality i, and kT2 and kFLAIR are 2 constants associated with the threshold of their respective image. These last 2 parameters replace the fixed threshold defined by Karimian and Jafari and were made to vary in the range of values {0, 1, 2}. Finally, hypointensity in the T1 image, defined by the negative threshold µWM, T1 − kT1 • σWM, T1, was also added to the set of rules. The parameter kT1 was varied in the same range as kT2 and kFLAIR.Another issue with the method described in the paper is that most predicted lesions were located at the outer edge of the white matter, near the cortical gray matter. The solution proposed by Zhao et al. [16] is to apply a mask obtained from the morphological erosion operation to the white matter mask, allowing the elimination of nearby voxels from the lesion candidates. The internal holes of the white matter mask were filled before applying the erosion mask. The chosen erosion filter was 6 × 6 × 6. A 3 × 3 erosion filter was applied to each slice of their set of images, whose resolution was half that of the images with which they worked, which is why the filter size was initially set to 6 6 6. These voxels were segmented as a lesion and their final size was set to 9 × 9 × 9, see Table 2.Table 2
Parameters of means obtained for the white matter.
Image no.µWM, T2µFLAIR, T2BCH0177.8973.7BCH0237.7662.18BCH0338.0880.85BCH0446.4100.96BCH0545.2663.65BCH0650.0161.64BCH0753.8859.78BCH0843.8860.7BCH0958.3857.84BCH1058.3356.88The initial segmentation to obtain the estimators of the initial parameters of the model was performed on the T1 image using theK-means algorithm.


### 3.2.1.K-Means Algorithm

TheK-means algorithm is an unsupervised learning algorithm commonly applied in clustering tasks. Given a set of points X = {x1, x2, …, xn} and a parameter k, corresponding to a positive integer, the algorithm is initialized by choosing k random points c1, c2, …, ck within the same domain of X known as centroids. Then, for each point within the set X, the distance between xi and each of the centroids cj is calculated (commonly, the Euclidean distance is used) and the point xi is assigned to the cluster corresponding to the closest centroid. Finally, the centroids are calculated again as(6)Cj=1nj∑xi∈Cj.xi,where Cj is the cluster containing the previously closest points to Cj and nj is the number of points in the cluster Cj. These 2 steps are repeated iteratively until the cluster assignment for each point does not change from one iteration to the next.This algorithm applies to the image segmentation task. Given the set of pixels of a grayscale image, the points to be grouped correspond to the intensities of the pixels. In contrast, the centroids correspond to different intensity values in the range [0, 255]. As a result, theK-means algorithm performs a grouping of pixels with similar intensity values. The K-means algorithm has been used to perform the segmentation of different brain tissues in magnetic resonance images [17].The K-means algorithm was used to perform a preliminary segmentation of the brain in the 3 different types of tissue based on the intensity of the pixels in the T1 image, obtaining, from lower to higher intensity, labels for the pixels corresponding to cerebrospinal fluid, spinal fluid, and gray and white matter.


## 3.2.1.K-Means Algorithm

TheK-means algorithm is an unsupervised learning algorithm commonly applied in clustering tasks. Given a set of points X = {x1, x2, …, xn} and a parameter k, corresponding to a positive integer, the algorithm is initialized by choosing k random points c1, c2, …, ck within the same domain of X known as centroids. Then, for each point within the set X, the distance between xi and each of the centroids cj is calculated (commonly, the Euclidean distance is used) and the point xi is assigned to the cluster corresponding to the closest centroid. Finally, the centroids are calculated again as(6)Cj=1nj∑xi∈Cj.xi,where Cj is the cluster containing the previously closest points to Cj and nj is the number of points in the cluster Cj. These 2 steps are repeated iteratively until the cluster assignment for each point does not change from one iteration to the next.This algorithm applies to the image segmentation task. Given the set of pixels of a grayscale image, the points to be grouped correspond to the intensities of the pixels. In contrast, the centroids correspond to different intensity values in the range [0, 255]. As a result, theK-means algorithm performs a grouping of pixels with similar intensity values. The K-means algorithm has been used to perform the segmentation of different brain tissues in magnetic resonance images [17].The K-means algorithm was used to perform a preliminary segmentation of the brain in the 3 different types of tissue based on the intensity of the pixels in the T1 image, obtaining, from lower to higher intensity, labels for the pixels corresponding to cerebrospinal fluid, spinal fluid, and gray and white matter.


## 3.3. Convolutional Neural Networks (CNNs)

Second, the first Birenbaum architecture [18] is used to train a convolutional network using a subsampling of FLAIR image patches. Instead of a longitudinal MRI, the paper’s SISTP (Single Image, Single TimePoint) architecture was used. In addition to the proposed architecture, an architecture without the Max pooling layer following the convolution layers was tested to compare performance. Both architectures tested 2 input types: 1 axial image and 1 image for each view type (axial, sagittal, and coronal). Each slice was obtained with a resolution of 0.5 mm and a distance of 1 mm between each slice. The intermediate values were interpolated during the registration process to obtain the final image with a resolution of 0.5 mm. It is then passed to a dense network of 16 neurons followed by the Softmax layer, which concatenates the 48 feature vectors obtained after the first dense layer of each image.The implementation of the network was carried out using Keras [19], a library written in Python for the management and construction of neural networks in a modular way, using the Tensorflow library [20] as a backend.In Figure2, candidate voxels were chosen based on preliminary white matter segmentation from the T1 image using the K-means algorithm. This white matter mask replaced the segmentation obtained from registering an Atlas to a T1 image. WM, FLAIR + kFLAIR • WM, FLAIR was used to define the threshold parameter for the FLAIR image [21]. The parameter was set to match the proportion of candidate voxels in Birenbaum et al. The final number of candidates for the network was 1,242,109 voxels (kFLAIR = 2), which is close to the percentage reported in Birenbaum’s paper, who claims that candidate sampling reduced network computations by 97.5 percent [20]. This subsampling leaves out a percentage of the image voxels corresponding to lesions. In addition to patch sampling, a loss function that penalized false negatives for lesions addressed the class imbalance issue. The weight parameter of this loss function was set at 20% for all trained networks based on the percentage of injuries observed in the training set. We will use leave-one-out cross-validation to evaluate the network's results. The network’s model is then applied to the untrained image, yielding model evaluation metrics for that image. This is done for every image in the training set. Finally, the obtained model’s evaluation metrics are averaged.Figure 2
Architecture for 1 implemented input of the CNN with and without Max pooling layer.


## 4. Results and Discussion

The results of the selection of parameters for the GMM are presented in the following graph (3), detailing the average DSC obtained for each parameter throughout all the images in the database. The best combination of parameters found, with respect to the DSC obtained, wash = 0.2; dmin = 0.4; kT2 = 0, and kFLAIR = 1.A larger area of the brain is considered a candidate for injury, so the number of true positives increases. There is a segmentation of more false positives, only a fraction of which are eliminated by the false positive reduction rules. Theh parameter also has the issue of losing information in the GMM due to training with fewer voxels, resulting in poor initial segmentation results. The elimination rule for voxels near the gray matter border reduced the number of false positives because many of the outliers detected in the model were in this region and were not filtered by the lesion selection rules. Figure 3 compares the Positive Predictive Value before and after applying this rule; see Figure 4.Figure 3
The results of the selection of parameters for the GMM.Figure 4
Comparison of PPV before and after application of erosion to the white matter mask.Because the original papers from which the implemented techniques came used different databases to validate his experiments, the results obtained by other techniques were investigated in the public database of the MS Lesion Segmentation Challenge. Overall, the results show a high PPV for most of the techniques used, with relatively precise segmentation of multiple sclerosis lesions. The tradeoff is that the TPR obtained was generally low compared to the state-of-the-art techniques. As for the GMM segmentation, this is due to applying a set of rules that filter out false positives and in some cases remove voxels that corresponded to injuries. The low number of true positives in convolutional networks is due to subsampling used to reduce network training times. The main advantage of GMM is the speed of segmentation. The average time required to segment each image was 538.61 15.47 seconds, including reading the file and applying the EM algorithm to the model. Across all test cases, the EM algorithm took 27.96 3.04 seconds per iteration. Convolutional networks performed better on average, with similarity coefficients of 0.33 and 0.29 for 1-view networks and 0.34 and 0.26 for multiview networks, comparable to most results in the state of the art. The network’s average performance was better without the Max pooling layer in terms of results quality. The network underestimates the size of lesions (which explains the low TPR) and simply fails to identify some injuries when the Max pooling layer is added. The number of lesions detected, that is, those for which at least one voxel was classified as a lesion, was as follows.In Figure5, the network with the Max pooling layer was able to detect only a subset of the lesions that the network without it could detect. The only exception to this was image BCH08, in which the network with Max pooling for the axial view managed to find 2 lesions that were not detected by the other network; however, these successes could be classified as noise simply by looking at the results of the segmentation because they appear only in 1 or 2 slices of the image and cover a very small region.Figure 5
Number of segmented lesions.In addition, in these images, it can be seen that the segmented regions corresponding to a lesion are not necessarily contiguous, and multiple segmented lesions may appear that correspond to a single lesion see Table3.Table 3
Number of segmented regions.
ImageCNNCNN (max pooling)Multiview CNNMultiview CNN (MP)No. of regionsSize >6No. of regionsSize >6No. of regionsSize >6No. of regionsSize >6BCH01716173401111536128519115BCH021975714840414103491122BCH0348212043411146511219754BCH04170828020843001374217826162BCH05121823769315090718739992BCH06670165512149599156352115BCH075111084059737086397100BCH0854512946511742010438486BCH09358753869562612416040BCH101593813332309806020It can be seen that the number of regions that are segmented by the CNN is quite large, which is observed in the results as a large number of small-sized regions that appear in the segmentation around larger regions (usually well-defined lesions) classified), regions around the lateral ventricles of the brain and at the outer edges of the white matter, on the border with the gray matter. An alternative, studied later, to mitigate this discrepancy problem consists of performing a postprocessing stage in which morphological operations are applied to the segmentation mask, improving its quality.Another problem observed with the segmentation of networks is the dependence on the initial selection of candidates. Although much of the white matter approximation usingK-means yielded good results, there are some cases in which the hypointense areas in the white matter were classified as another type of tissue.In this case, the GMM segmentation is more accurate because it considers all tissue types when determining injury candidates. Although the GMM, likeK-means, failed to segment white matter, all voxels were considered in the selection of outliers allowing better detection of lesions, not just those predicted to be white matter. One of the benefits of using the Max pooling layer is the reduced training time. The Max pooling layer reduced the feature map size, reducing the network training time by approximately 3. The results of multiview networks were not recorded due to the lengthy training time on the CPU (in hours, minutes, and seconds of execution). The GPU (Graphic Processing Units) reduces the training time. The time required to train the networks, even on the GPU, exceeds the average time required to obtain results from the GMM, which was approximately 9 minutes per image. This time includes file reading, algorithm execution, and segmentation generation. Unlike convolutional networks, which can take days to train on the CPU, the GMM algorithm can be executed on the CPU without compromising execution time; see Figure 6.Figure 6
Example of an image in which the automatic segmentation of WM fails. (a) Preprocessed T1 image. (b) Preprocessed T2 image. (c) Preprocessed FLAIR image. (d) WM obtained fromK-means. (e) Segmentation delivered by the CHB expert. (f) Segmentation obtained from the GMM. (g) Segmentation obtained by CNN (without max pooling). (h) Segmentation obtained by CNN (with Max pooling). (i) CNN multiview segmentation (no Max pooling). (j) CNN multiview segmentation (with Max pooling).


### 4.1. Postprocessing

To improve the quality of the segmentation, morphological operations were applied to the obtained segmentations. The operation performed consists of a closure, which consists of applying a dilation filter followed by applying an erosion filter. The size of this filter was set at 5 × 5 to connect those regions that are close to each other but not enough to join regions that are not related to each other.The results obtained after applying this filter, in addition to the percentage change observed for each metric, are found in Table4.Table 4
Postprocessing results.
Image no.GMMCNNCNN(Max pooling)Multiview CNNMultiview CNN (MP)TPRPPVDSCTPRPPVDSCTPRPPVDSCTPRPPVDSCTPRPPVDSCBCH010.460.80.590.580.570.570.330.60.420.550.590.570.320.640.32BCH020.140.530.230.430.610.50.190.790.30.460.560.50.260.550.35BCH030.20.690.30.390.470.430.350.490.410.40.460.430.250.650.36BCH040.440.260.330.160.110.130.380.160.230.320.180.230.340.290.32BCH050.20.630.30.250.210.230.260.260.260.250.230.240.230.320.26BCH060.120.410.190.30.580.390.270.620.370.290.60.390.240.660.35BCH070.340.620.440.370.690.480.320.730.440.360.730.480.280.770.41BCH080.270.590.370.420.770.550.310.770.440.410.780.550.260.770.39BCH090.040.070.050.240.530.330.30.530.30.30.420.350.110.610.19BCH100.020.120.040.090.720.160.070.610.120.130.610.230.040.790.07Average0.220.470.280.320.530.380.280.560.330.350.520.400.230.610.30The quality of segmentation improved significantly for all methods. That is because using morphological closure helped to alleviate one of the main issues with segmentations, which was an underestimation of the lesions’ size, which was often mistaken for multiple isolated regions. The dilation filter joins adjacent regions by closing. After applying the erosion filter, unconnected regions return to their original state, while connected regions stay connected for the most part.In Figure7, this filter also solves the problem of small holes appearing in some segmentation regions. This problem was most commonly seen in neural networks, where initial candidate selection was left out of the segmentation voxels that were largely surrounded by others chosen as candidates.Figure 7
Changes in the metrics after the application of postprocessing.In Figure8, It should be noted that this filter increases the number of false positives in the segmentation since it allows those regions corresponding to false positives to connect in a similar way to lesions. However, given the improvement observed in the TPR and DSC metrics, it can be said that the advantages of applying the filter outweigh its disadvantages, especially since the filter attacks the weaknesses observed in the techniques, which were able to obtain a segmentation with a low percentage of both false positives and hits to lesions.Figure 8
Example of the benefits of the closing filter.


## 4.1. Postprocessing

To improve the quality of the segmentation, morphological operations were applied to the obtained segmentations. The operation performed consists of a closure, which consists of applying a dilation filter followed by applying an erosion filter. The size of this filter was set at 5 × 5 to connect those regions that are close to each other but not enough to join regions that are not related to each other.The results obtained after applying this filter, in addition to the percentage change observed for each metric, are found in Table4.Table 4
Postprocessing results.
Image no.GMMCNNCNN(Max pooling)Multiview CNNMultiview CNN (MP)TPRPPVDSCTPRPPVDSCTPRPPVDSCTPRPPVDSCTPRPPVDSCBCH010.460.80.590.580.570.570.330.60.420.550.590.570.320.640.32BCH020.140.530.230.430.610.50.190.790.30.460.560.50.260.550.35BCH030.20.690.30.390.470.430.350.490.410.40.460.430.250.650.36BCH040.440.260.330.160.110.130.380.160.230.320.180.230.340.290.32BCH050.20.630.30.250.210.230.260.260.260.250.230.240.230.320.26BCH060.120.410.190.30.580.390.270.620.370.290.60.390.240.660.35BCH070.340.620.440.370.690.480.320.730.440.360.730.480.280.770.41BCH080.270.590.370.420.770.550.310.770.440.410.780.550.260.770.39BCH090.040.070.050.240.530.330.30.530.30.30.420.350.110.610.19BCH100.020.120.040.090.720.160.070.610.120.130.610.230.040.790.07Average0.220.470.280.320.530.380.280.560.330.350.520.400.230.610.30The quality of segmentation improved significantly for all methods. That is because using morphological closure helped to alleviate one of the main issues with segmentations, which was an underestimation of the lesions’ size, which was often mistaken for multiple isolated regions. The dilation filter joins adjacent regions by closing. After applying the erosion filter, unconnected regions return to their original state, while connected regions stay connected for the most part.In Figure7, this filter also solves the problem of small holes appearing in some segmentation regions. This problem was most commonly seen in neural networks, where initial candidate selection was left out of the segmentation voxels that were largely surrounded by others chosen as candidates.Figure 7
Changes in the metrics after the application of postprocessing.In Figure8, It should be noted that this filter increases the number of false positives in the segmentation since it allows those regions corresponding to false positives to connect in a similar way to lesions. However, given the improvement observed in the TPR and DSC metrics, it can be said that the advantages of applying the filter outweigh its disadvantages, especially since the filter attacks the weaknesses observed in the techniques, which were able to obtain a segmentation with a low percentage of both false positives and hits to lesions.Figure 8
Example of the benefits of the closing filter.


## 5. Conclusions

For the MS Lesion Segmentation training database, two alternative methodologies of the state of the art were implemented, along with minor changes. Based on outlier detection using GMM and segmentation using convolutional neural networks, these algorithms had not been tested in this database and had no performance records. Both approaches achieve precise segmentation of MS lesions with good PPV. The fundamental issue with both approaches was the poor average TPR. Convolutional neural networks performed the best, successfully classifying multiple sclerosis lesions from the first set of candidates. It is worth noting that the number of photos used to create these segmentations was rather small, around 1.5% of all possible patches. In this regard, the poor TPR seen due to the small number of patches utilized and the issues observed in automatic matter segmentation may be improved by changing the patch selection for this technique. However, segmentation produced by detecting GMM outliers has several advantages over CNN. CNN's estimate is partially based on the initial injury candidate selection. Using solely white matter voxel-cantered patches presents issues if the white matter segmentation has faults due to picture abnormalities. The GMM could not detect this issue because these abnormalities tend to show as outliers. Due to the minimal number of parameters in the GMM, the EM algorithm executes quickly, allowing for rapid segmentation of lesions. Using the GPU reduces the training time of CNNs to a level comparable to the time necessary to acquire a segmentation using the GMM in the CPU. Finally, applying a closure filter to the segmentation findings improved them significantly, alleviating underestimating the lesion size reported in both techniques.


## 6. Recommendation

Using the suggested methodologies, the segmentation quality can be improved in the state of the art for the database. A contiguous region for a single lesion is not always attainable because the approach does not employ information about the voxel's immediate surroundings, which causes segmentation issues. Instead, the GMM might include voxel-specific properties. A voxel's near surroundings, as well as intensities from each imaging modality, may help segmentation. CNN can be improved greatly. Subsampling candidate voxels removes lesions not used in the model’s training. Using different MRI modalities to identify outliers could improve the candidate voxel sampling step. For generalization, neural networks require a lot of data. The networks’ outcomes may improve if they access more images from varied patients. The neural networks work for MRI at the same resolution they were trained on. Adding more data preprocessing stages would allow input images to match the training images in resolution.



---

*Source: 4905115-2022-04-22.xml*