Probabilistic Multiple Sclerosis Lesion Detection using Superpixels and Markov Random Fields Detección Probabilística de Lesiones de Esclerosis Múltiple usando Superpixeles

Multiple Sclerosis (MS) is the most common neurodegenerative disease among young adults. Diagnosis and mo-nitoring of MS is performed with T2-weighted or T2 FLAIR magnetic resonance imaging, where MS lesions appear as hyperintense spots in the white matter. In recent years, multiple algorithms have been proposed to detect these lesions with varying success rates, which greatly depend on the amount of a priori information required by each algorithm, such as the use of an atlas or the involvement of an expert to guide the segmentation process. In this work, a fully automatic method that does not rely on a priori anatomical information is proposed and evaluated. The proposed algorithm is based on an over-segmentation in superpixels and their classification by means of Gauss-Markov Measure Fields (GMMF). The main advantage of the over-segmentation is that it preserves the borders between tissues, while the GMMF classifier is robust to noise and computationally efficient. The proposed segmentation is then applied in two stages: first to segment the brain region and then to detect hyperintense spots within the brain. The proposed method is evaluated with synthetic images from BrainWeb, as well as real images from MS patients. The proposed method produces competitive results with respect to other algorithms in the state of the art, without requiring user assistance nor anatomical prior information.


INTRODUCTION
Neurodegenerative diseases are one of the most critical issues for the health sector. Not only elderly people are the most affected by neurodegenerative diseases, but also young people can suffer from them. Multiple sclerosis (MS) is a neurodegenerative disease that mainly affects people between 20 and 40 years old, with high incidence in the general population. In fact, it is the second in incidence, epilepsy being the first [1] . The cause of MS is controversial but seems to depend on genetic and environmental factors and may also have a strong auto-immune component. The diagnosis and prognosis are well established nowadays by neurologists. The symptoms are described by the patient, and evidence is found through physical examination. In clinics, the neurologist verifies the symptoms of the patient (for example weakness, blurred vision, ataxia, etc.), and then requests a brain imaging study. Magnetic resonance imaging (MRI) is highly recommended in this case, with the most common protocols for this purpose being: T1-w, T2-w, Proton Density (PD) and Fluid Attenuated Inversion Recovery (FLAIR). In brain images, the neurologist manually annotates MS lesions, which in FLAIR images are shown as high-intensity spots on the white matter. It is important to mention that manual annotation and counting of hyper-intense spots is often an extensive and tedious process because the clinician needs to check dozens of images and may find several spots for a single patient. Hence, there is a large interest in designing algorithms that can automatically detect MS lesions or assist the expert during the process [2] [3] . In the past decades, various methods to segment MS lesions in MRI images have been published; however, some of these methods suffer from low accuracy or have such a large number of parameters to tune that they are not user-friendly. Other approaches rely heavily on the user's participation, or need atlas databases, requiring additional preprocessing time (for instance, to align the atlas with the input images). For these reasons, it is interesting to develop fully automatic and user-friendly methods to aid in the detection of MS lesions.
There are a few reviews of methods for MS lesion segmentation in the literature [4] [5] . Some of those methods are based on probabilistic approaches [6] [7] , support vector machines (SVM) [8] , region growing [9] , K-Nearest Neighbors [10] or neural networks [11] , while some methods may also use additional prior information such as atlases [12] or clinical information [13] . Many of these methods use pre-processing steps to prepare the input images for MS lesion detection, such as image denoising and non-uniformity correction. Also, since some non-brain tissues such as scalp and optic nerve are also shown with high intensity on T2-w images, a skull stripping step is often required; to this end, several methods use the Brain Extraction Tool (BET) [3] .
Among the algorithms for MS lesion detection, several methods are based on Expectation-Maximization (EM) to segment MS lesions [14] [15] due to its good accuracy and easy implementation. In the EM-based method proposed by Garcia-Lorenzo et al. [15] , the algorithm is divided in three steps. In the first and second steps, there is a non-uniformity correction of the input image, followed by a skull stripping process. Finally, in the last step, the MS lesion detection is applied using clinical rules to select potential regions with good results (reported Specificity of 0.9954). Another interesting method, which uses Markov Random Fields (MRF), is proposed by Khayati et al. [6] [7] . They developed an MS lesion detector by estimating a conditional probability density function for each class, which was trained using the adaptive mixtures method (AMM).
For validating their results, in [6] they use a cross-validation approach where the first MS reader was used as the gold standard, leading to very accurate results since they achieved an average Dice Similarity Coefficient (DSC) of 0.75 [16] [17] . Other proposal was developed by Lao et al. [18] , where they first perform affine registration of T1-w, T2-w PD and FLAIR images by maximization of the mutual information [19] and skull stripping based on affine registration using the BET algorithm [3] . In the training process, the proposed method combines T1-w, T2-w, PD and FLAIR in an attribute vector (AV) for each voxel, along with the information of the neighboring voxels. A set of manually segmented scans is used to train a support vector machine (SVM) using the AdaBoost method [8] [20] ; once trained, the SVM outputs, for each voxel, a scalar measure of abnormality that is binarized by applying a tuned threshold to discriminate lesions from normal tissue. Although similarity values (e.g., Dice index) were not reported, specificity and sensitivity look very promising which are also complemented by good visual results. Finally, in [21] the proposed method is based on Artificial Neural Networks (ANN), and the implemented software is freely available online at https://med.inria.fr/the-app/downloads. After selecting the segmentation option in this interface, the user can upload data and interactively click on the MS lesion to enhance them. One disadvantage of this software is that sometimes the algorithm segments the entire white matter region, particularly when MS lesions have blurred borders. Despite the diversity of algorithms for MS lesion detection, it is difficult to find one that fulfills the requirements of public health institutions, particularly when images are acquired with low resolution, few slices, and using a single modality to reduce costs. These requirements include border preservation, robustness to noise and blurred borders, good accuracy with single-modality (T2-w or T2-FLAIR) images, a reduced number of tuning parameters, and an implementation that does not rely on atlases or databases.
In this paper, a new method for automatic MS lesion segmentation is proposed. The method works with single-modality low resolution images, it does not need prior information related to anatomical structures or user annotations, and it only requires a few parameters to be tuned. In order to achieve good border preservation and robustness, an over-segmentation in superpixels (SPs) is performed as the first step [22] , followed by a post-processing stage to achieve connectedness and eliminate spurious SPs produced by noise. Each SP is then classified by means of a Gauss-Markov Measure Field (GMMF) model [23] . This segmentation approach is applied twice: first to isolate the brain region, and a second time to segment MS lesions in the brain. This paper is organized as follows: Section 1 contains the introduction and a description of the goals of this work. Section 2 presents the details of the proposed segmentation algorithm, which is called SP-GMMF, and the methodology for MS lesion segmentation.
Results from applying the proposed method to the analysis of synthetic and real MR images are shown in Section 3, with a comparison against a state-of-the-art algorithm based on Expectation Maximization (EM).
Finally, the conclusions of this work are presented in Section 4.

MATERIALS AND METHODS
In this proposal, segmentation of MS lesions is achieved in two stages: the first stage is to isolate the region of interest (in this case, the brain) from the rest of the image (skull, meninges, or bone cavities); then, in the second stage, MS lesions within the brain region are detected. Each stage employs a segmentation algorithm. The segmentation method proposed here is a combination of an over-segmentation of the image in superpixels using the Simple Linear Iterative Clustering (SLIC) [22] and a probabilistic labeling of the superpixels by means of the GMMF [23] . Besides the main steps of the segmentation process, there are other pre and post-processing steps that are important for the efficiency of the proposed algorithm. The details of the complete algorithm are described below.
The SLIC method is a clustering algorithm that combines spatial location and intensity information to subdivide an image in a relatively small number of groups of pixels that have similar color and are spatially coherent, commonly called superpixels [22] [24] [25] .
This algorithm is a variant of the k-means algorithm that uses a reduced search space to form each SP, and whose main advantages are high computational speed and border-preserving SPs. The SLIC method works as follows: let us define l( ) as the image over the lattice L (that is, l(r) ∈ L), and each superpixel is defined as S k = [C k , D k ] where C k and D k are the average intensity and the geometric center of the superpixel, respectively.
Indexes of superpixels are denoted by k, that is k = 1, 2, …, K, with K as the number of superpixels over L. The initialization of each D k is given by a regular hexagonal grid in L, and C k as the intensity at pixel site C k . For each superpixel S k , a square neighborhood of size 2M × 2M is defined with center at D k and M = √|L|/K', where K' is the desired number of SPs given by the user; as a rule of thumb, this parameter can be defined as the total number of pixels in the image divided by the area (in pixels) of the smallest lesion that is observed; for instance, we are using K' = 3000 SPs for all images in this study. The actual number of superpixels K will vary across all the steps in the algorithm. Only those pixels that belong to the neighborhood of D k can be assigned to S k according to the distance measure: where the first term measures the intensity distance, and it is normalized by the dynamic range of the data m, which can be computed as m = max (l ( )) -min (l ( )).
The second term is the spatial distance term and it is normalized by the size of the neighborhood. Finally, γ is a hyperparameter that weighs the importance between both terms. Once every pixel is assigned to some S k , C k and D k are updated with the average intensity and spatial position of the pixels that belong to the k-th SP, and δ k is computed again; this process is iterated until convergence. In our experience, the algorithm converges in 5 to 10 iterations. Although SPs adhere well to borders, for an adequate choice of γ, the SLIC algorithm can sometimes produce fragmented superpixels, which may lead to unconventional neighborhoods for the Markovian Fields process [23] . For that reason, a relabeling process by means of connected component algorithm is applied [26] , so that each connected fragment of a fragmented SP is considered as a new, individual superpixel. Once the relabeled field has been obtained, small SPs which are often due to noise, are fused to the most similar (in terms of average intensity) neighboring superpixel. In our case, a superpixel is considered too small when its area is less that 3% of the SP average area M 2 ; the 3% threshold was found experimentally as the percentage for which the number of SPs after fusion approximates better the number of desired SPs K'. With this fusion, not only the number of variables, but also the noise was reduced. Once the SPs are obtained, the next step for the segmentation process is to classify them with respect to their intensity. For that purpose, we consider that each class is defined by a Gaussian distribution of pixel intensities with a given mean and variance. Note that, although we aim for a binary classification at each stage (brain vs non-brain in the first stage, lesion vs non-lesion in the second stage), there are several types of tissues represented in the images; for this reason, multiple classes must be considered using the GMMF model [23] , the goal is to estimate, for each SP S k , the probability p j (S k ) that it belongs to class j, for each j = 1, …, C, where C is the number of classes. Under this model, one can obtain the probability field for each class j by minimizing the following energy function U given by: where the first term is a data term which enforces p j (k) to be similar to the normalized likelihood between the k-th superpixel and the j-th class, which is given by where v j (k) is a likelihood function which measures how well the k-th SP fits in class j.
Assuming that classes follow a Gaussian distribution, the likelihood can be obtained as: where μ j and σ j are the mean and variance of the j-th class, respectively. Notice that the μ can be automatically initialized by the k-means method. Additionally, there is a hyper-parameter κ which controls the overall variance of all classes. The second term is a regularization term that promotes the similarity in the neighborhood. The neighborhood N k for the k-th superpixel can be obtained, first, inspecting the borders of the image of superpixel labels and then inspecting the image of labels; i.e., when a border is detected, the labels in that point (SP) are added to its SP neighborhood. Finally, λ is a hyper-parameter that weighs the importance between terms. To solve (2), one can calculate its derivative and equal it to zero to obtain a linear equation system, which can be iteratively solved by the Gauss-Seidel method where the solution of p j (k) is given by: After each Gauss-Seidel iteration, the means and variances of all classes are updated with a forgetting factor where t is the current iteration). This method generally stops until convergence is obtained, but in our experience, it converges from 5 to 10 iterations. Once convergence is achieved, each superpixel S k is assigned to the class given by: Figure 1 shows the flow chart of the SP-GMMF Segmentation proposal.

MS Lesion Classification
The SP-GMMF segmentation proposal takes advantages from SLIC and GMMF methods to implement a general-purpose segmentation method. In the proposed two-stage algorithm for MS lesion detection, which is illustrated in Figure 2, the first stage is ori- where the hemispheres are joined into a single region, the next largest region will be another structure whose area will be significantly smaller than the largest one (with a ratio less than 0.7), and thus the algorithm will consider it as part of the brain. Once the brain area mask has been obtained, a hole-filling algorithm is applied to it in order to recover any dark structure within the brain that might have been discarded by the previous operation. This mask, applied to the MRI image, automatically isolates the brain for the rest of the image (Figure 2e).
On the other hand, MS lesions appear in T2 images as hyperintense spots; for this reason, an intermediate step is to apply a contrast-enhancing intensity adjustment to the isolated brain image to saturate the hyperintensities. This operation scales the voxel intensities by a factor such that the lowest 3% percentile will be saturated to zero, and the highest 3% percentile will be saturated to 255 (in an 8-bit grayscale image), which will facilitate the MS lesion detection for the next segmentation procedure.
In the second stage, the SP-GMMF segmentation is applied to the intensity adjusted brain image with a finer resolution of the SPs; that is, a larger desired number of SPs K' (Figures 2f and 2g). To obtain the binary mask that contains the MS lesions, the algo-  algorithm [4] . The main reason for choosing EM is due to a good balance between simplicity, popularity and good results in this task. In general, this approach assumes that intensities belonging to the structures on MRI images follow a Rician distribution that can be fairly approximated by a Gaussian distribution [13] [14] . Additionally, another proposal for MS lesion segmentation based on EM is presented here. Specifically, the two-step algorithm can be implemented, first using SP-GMMF to segment the brain, and then using EM at pixel level classification to segment MS lesions. For simplicity, this algorithm will be called EM*. Notice that, both EM* and EM methods follow the same MS lesion classification of the proposed SP-GMMF method.

EM* Classification
In other words, EM* and EM methods, in the first step of the algorithm obtain a mask of the brain, then an intensity adjustment operation is applied to the masked region, and finally in the second step, the EM algorithm is applied to the intensity-adjusted brain, and the same discrimination criteria is applied to obtain the MS lesions.

RESULTS AND DISCUSSION
For comparison purposes, both SP-GMMF and EM* were tested and compared against EM using both synthetic and real brain MRI images. A brief summary of the details for each algorithm is shown in Table 1.
Experimental results were obtained using the BrainWeb dataset [27] and several MRI images (axial and sagittal sequences) from two subjects positively  The parameters for the first step (brain peeling) with SP-GMMF and EM* methods were: K' = 1000 desired

Results from synthetic images
Data for these experiments consisted on simulated 181×217×181 T2-weighted MRI volumes from the MS Lesion Brain Database in Brain Web [27] Table 2.    Table 3. Figure 5 shows boxplots of the DSC and SEN results.

Tabla 4
Method DSC for BrainWeb Images

Requires multiple modality
Requires prior data TOADS [12] 0.79 0.63 X X Graph Cuts [28] ~0.7 0.63 X AWEM [30] 0.7 X CGMM+CE [31] 0.78 X X Khayati [ results are summarized in Table 2, shows that, in general, EM has a significantly lower performance than EM* and SP-GMMF, possibly due to errors in the brain peeling stage. On the other hand, no significant differences between EM* and SP-GMMF were found.
In the literature, some methods are applied to the BrainWeb data [24] [28] [29] [30] [31] with a high performance in DSC. However, their good results are not clearly explained, for example, some of them do not mention the slice (or slices) used for the tests. In [29] , experiments using T1 and T2 images with 3% and 5% of noise are shown, obtaining a DSC of 0.782 in the best case. These experiments illustrate the main advantage of using prior information from an atlas; in general, using an atlas is a good option, however, algorithms without an atlas could reach similar performance; for instance, in [29] they report a DSC = 0.74 with an atlas and DSC = 0.75 without atlas. Moreover, they do not indicate which atlas is used or how the subject brain is registered to the atlas, which suggests they might have used the healthy BrainWeb data as atlas in order to avoid the registration process; in that case, applying the method proposed in [29] to real images would require further pre-processing steps which could introduce additional errors. In [31] , the authors show a methodology based in Gauss-Markov model followed by curve evolution. Experiments were achieved using 61 slices (from the 60th to the 120th) using the T1, T2 and PD images. They obtained very good results since DSC is over 0.78 for 3% to 9% of noise. The results obtained by Garcia-Lorenzo et al. [28] present an interesting behavior for INU = 0%. The method yields a lower DSC = 0.24 for low noise levels (1%), increases its performance for 3% and 5% of noise (DSC of 0.65 and 0.79 respectively), and then decreases at 7% and 9% with DSC = 0.6 and DSC = 0.25, respectively. These results, along with our own results from Figure 4, suggest that some algorithms perform better when there is a mild amount of noise in the images. For instance, the EM and EM* algorithms estimate the mean and standard deviation of gray intensities for each class; in the absence of noise, the standard deviation will approach zero, which may cause numerical instabilities; on the other hand, the algorithm in [28] is gradient-based and may also be affected by untextured, noiseless regions.  [30] , they report results from three different patients for which they can obtain DSC values of 82%, 56%, and 52%, respectively. In [28] , experiments with 10 patients yield DSC values between 40% and 75%. Considering that the authors of these works use multi-modal images (combining T1-w, T2-w, FLAIR and PD) with high resolution, we consider our results (DSC between 45% and 69%) to be competitive for low-resolution single modality imaging.
Finally, we discuss the differences in the false positive lesions reported by the proposed methods, as shown in

CONCLUSIONS
An automatic algorithm to detect and segment Multiple Sclerosis lesions from T2 MR images was presented. The proposed method is based on a segmentation process where the image is subdivided into connected clusters (superpixels) which are then labeled according to their average intensity using Gauss-Markov Measure Field model. The segmentation pro-cess is applied two-fold: first to isolate the brain region, and then to detect hyperintense spots within the brain region. Finally, some of the false positives are discarded depending on their area and eccentric- ity. An experimental test using synthetic images from the BrainWeb database shows that the proposed segmentation has strong advantages against the popular EM method, even under noisy conditions. While SP-GMMF is fairly robust to noise, it is very sensitive to non-uniformity, so additional pre-processing might be required for some images to deal with the spatial intensity variations. In the case of real MRI images, SP-GMMF maintains its advantages against EM, although using EM for detecting hyperintense spots (once the brain region has been isolated) has shown benefits, such as a more accurate lesion count. In brief,