Multiple Sclerosis Lesion Segmentation Using Statistical and Topological Atlases

This paper presents a new fully automatic method for segmentation of brain images that possess multiple sclerosis (MS) lesions. Multichannel magnetic resonance images are used to delineate multiple sclerosis lesions while segmenting the brain into its major structures. The method is an atlas based segmentation technique employing a topological atlas as well as a statistical atlas. An advantage of this approach is that all segmented structures are topologically constrained, thereby allowing subsequent processing with cortical unfolding or diffeomorphic shape analysis techniques. Validation on data from two studies demonstrates that the method has an accuracy comparable with other MS lesion segmentation methods, while simultaneously segmenting the whole brain

Multiple Sclerosis (MS) is a demyelinating disease of the central nervous system that commonly leads to inflammatory and atrophic pathology, often causing cognitive impairment [1,2].It is primarily expressed as focal lesions in the white matter (WM) of the brain, but the state and progression of the disease is also correlated with cerebral atrophy [1,3].Because of its superior contrast, magnetic resonance (MR) imaging is the modality of choice for clinical evaluations of MS.Quantitative analysis of MR images to measure and monitor the lesion load and tissue volumes has also become invaluable for patient follow-up and evaluation of therapies.Manual delineation of MS lesions, however, is both challenging and time-consuming since three-dimensional information from several MR contrasts must be integrated.
Several techniques have been proposed for automated MS lesion segmentation.Most of these techniques rely on multichannel MR acquisitions and then employing pattern recognition techniques to detect pixels that are either outliers to the standard intensity profiles of healthy brain tissue [4,5], or are similar to lesion intensity profiles derived from a training set [6,7].Since the intensities of lesions often overlap with the intensities of other structures in the brain, additional contextual processing is often utilized to minimize false positives.The problem of segmenting MS lesions is closely related to the detection of white matter signal abnormalities that frequently occur in Alzheimer's disease and older populations [8,9,10].
Most of these techniques focus solely on lesion segmentation, even though they may recover a tissue classification as well [5,6,7,11].The extent and location of brain atrophy, important for monitoring the progression of the disease, is either not computed or is not subject to validation.Standard processing techniques for measuring these quantities are often not applicable to data that possess lesions.Furthermore, none of these methods segment the sub-cortical structures of the brain.Advanced analysis of cortical structure [12,13] is also not readily applicable as the topology of the brain is changed by the lesions.Another disadvantage of many lesion segmentation algorithms, particularly those that employ training data to model lesion intensity profiles, is that they are dependent on a specific acquisition pulse sequence.These approaches must be modified or re-trained to process data acquired using alternative pulse sequences.
In this paper, we propose a new technique for segmenting white matter lesions in MS that provides a detailed and topologically consistent segmentation of the brain into its main cortical and sub-cortical components.
The method incorporates both spatial and intensity information to segment multichannel MR images without post-processing.Moreover, this method enforces topological constraints in such a way that the segmented images reach topological equivalence with those of healthy subjects, allowing the direct use of techniques for subsequently performing cortical reconstruction and cortical unfolding [13], as well as diffeomorphic shape analysis [14].Although it utilizes multichannel acquisitions, training data is not required to model the intensity distributions, allowing the algorithm to be flexible enough to be applied to data originating from a variety of pulse sequences.
The method extends our previous work on topology-preserving anatomical segmentation [15] by introducing several modifications that effectively model the lesions while maintaining the constraints provided by topological and statistical atlases.The key observation from which the new approach is formulated is that topological outliers, such as lesions, can be addressed in a topology-preserving framework when they are grouped together with the appropriate classes.We validate this approach using manual delineations from two studies of MS patients, demonstrating good performance.

Statistical and Topological Atlas-based Segmentation in Healthy Anatomy
Our algorithm segments the brain into its major structures (cerebral gray and white matter, cerebellar gray and white matter, basal ganglia, ventricles, and brainstem) while delineating MS white matter lesions.It is based on our previous work that incorporated information from statistical and topological atlases into an intensity-based classification technique [15].This method is capable of segmenting brains while preserving the global topology of these structures.In this section, we briefly explain the different components of this method, before introducing the specific improvements required for handling lesions.
The segmentation method utilizes a statistical atlas as well as a topological atlas.The statistical atlas is built from a set of 18 manual delineations of the structures of interest, based on the IBSR dataset [16].The boundary of each structure was blurred to make a smooth probability map and to account for anatomical variations beyond those present within the training set.The topological atlas is a parcellation of the brain edited to encode a specific topology for each structure and group of structures, based on statistical atlases and anatomy textbooks.The topological atlas is used for preserving topology as well as lowering the influence of competing intensity clusters in regions that are spatially disconnected, while the statistical atlas affects the segmentation of adjacent structures having similar intensity.Prior to the segmentation, the atlas is rigidly registered to the studied MR image, and the registration is further updated at each iteration of the segmentation algorithm.
Given a multichannel MR image {I i }, where i denotes the channel index, the segmentation algorithm performs two interleaved processes consisting of (1) performing a fuzzy segmentation, (2) defining topologically consistent regions using fast marching.The fuzzy segmentation is obtained by minimizing the following energy function: (1) with respect to fuzzy membership functions u jk for each pixel j and structure k, a gain field g i j that models intensity inhomogeneities for each channel i, and the intensity centroids c i k for each structure.The first term in (1) is data driven, the second term enforces smoothness on the memberships [17], and the third term controls the influence of the statistical atlas.The parameters β and γ control the relative weighting of each term and are set empircally.The exponent q is a parameter controlling the "hardness" of the membership functions and is typically set to 2. The gain field is a smoothly varying function modeled as a low degree polynomial (see [18] for details).The variables p jk are probabilities derived from the statistical atlas that represent the prior probability of pixel j being inside structure k.The variables w km and r jk are weighting the impact of the atlas and the structure relationships as described in the following.
The atlas weight w km between two classes is a function of distance between their centroids: where s w is a parameter and typically set to 0.1.The atlas weight is close to one when c k c m but goes to zero when c k = c m , so that the priors influence the segmentation only where the intensity contrast between structures is low.
The relationship weights r jk take into account global and local relationships between the structures, and are defined by The relationship weights penalize against membership configurations that are inconsistent with the topology atlas.
The energy function (1) is used to compute membership functions for each structure in a coordinate descent fashion, similar to FANTASM (Fuzzy And Noise Tolerant Adaptive Segmentation Method) [17].In addition, we compute a "hard" segmentation that is derived from the memberships but is constrained to be homeomorphic to the topology template.This hard segmentation is performed using two successive iterations of a fast marching front propagation technique, where the speed function is modulated by the memberships and the topology is preserved.The first step thins the structures into a skeleton-like object and then the second step grows the structures back to find the optimal boundary, using a minimal path strategy [19,15].
The complete algorithm follows these steps: 1. Align atlases to image and set initial segmentation to the topological atlas.
2. Compute r jk from current hard segmentation 3. Compute the memberships u jk , centroids c k and the inhomogeneity field g j .
4. Thin structures using the fast marching algorithm.
5. Grow back the structures and update the segmentation.
6. Refine the alignment of the atlases.
The convergence criterion is the relative amount of change in the energy E with each iteration.We choose a threshold of 10 −4 or maximum of 20 iterations.

Lesions in Topology-preserving Segmentation
The segmentation method above is based on anatomical priors, which seems to conflict with the globally distributed nature of MS lesions.Since lesions can occur anywhere in the WM, we cannot associate to them a specific topological or statistical model, unlike other structures.Similarly, the topology of the WM is modified arbitrarily by the appearance of lesions.However, if we make the observation that the MS lesions must appear inside the WM region, we can then assume that the structures made of WM grouped with the lesions have the same shape and topology as healthy WM.Thus, in an anatomical sense, WM and lesions are treated as a single structure.This allows our previously generated topological and statistical atlases to be directly applicable to the segmentation of brain images possessing lesions, even though they were created from images of healthy subjects.
To address lesions in the algorithm, we first add an additional lesion class to our model, with a corresponding lesion membership function and intensity centroid.We also need to modify the definitions of the atlas weights w km and the relationship weights r jk to reflect the new model, as described later in this section.The values of the statistical atlas prior p jm for lesions are set to the WM value.The speed function for evolving the hard segmentation of the grouped WM and lesions uses the sum of membership functions for WM and lesions.The lesion and WM are then separated by selecting, inside the grouped region, whichever has the higher membership value.With these alterations, we preserve the topology of brain structures while computing a competitive process between WM, lesions, and neighboring structures.
The atlas weights should be large when the centroids of two nearby structures are close to each other.For the lesion atlas weight, we use the lesion centroid to compute w km in (2).When comparing lesions and WM, however, the atlas priors should have no influence, so we set w lesion,wm = 0. Similarly, the relationship function r j,lesion can be replaced by the function of the underlying region, i.e. r j,lesion = r j,wm .However, in some regions we further modify this function to reduce the potential for false postives.
Classifying lesions using a tissue classification or clustering technique often suffers from a large amount of false positives.This is due to the fact that lesions can have intensity profiles close to those of other structures in the brain.For example in T1 images, MS lesions appear as hypointense WM voxels whose intensity is close to GM intensity (see Fig. 1).Because we compute a hard segmentation of brain structures at each step of our algorithm, we can use our knowledge about areas where false positives commonly occur to define an appropriate relationship function.In many cases, it is unlikely that MS lesions appear right at the boundary between cortical or sub-cortical GM and WM.However, because the intensity centroid of the lesion class lies between centroids of these structures, their boundaries are common areas for false positives.Similarly, the boundary between WM and ventricles may appear bright in T2, PD and FLAIR (FLuid Attenuated Inversion Recovery) images and yet not contain any lesion.
One area where false positives commonly occur is in the WM between ventricles.This occurs particularly because this area tends to enhance on FLAIR even when lesions are not present.By modifying the statistical atlas to treat this region as a separate class, we are able to reduce these false positives.(see Fig 2) Another area where false positives often occur is along the GM/WM boundary.To address this problem, we update the relationship function between lesion and this structures as a function of the distance from the boundary of GM , ventricles and inter-ventricular WM: )r j,wm d j,GM ≤ d max,GM and d j,V EN > d max,V EN , r j,wm otherwise. ( where d j,V EN and d j,GM are the distance to ventricles and GM structures.
In addition to modifying the lesion relationship function with respect to distance from GM and ventricles, we also modify it with respect to distance from inter ventricular WM: r j,lesion otherwise. ( where d j,W M int is the distance from boundary of inter-ventricular WM.
In addition, the relationship function for the cortical GM is also modified in the region where d j,V EN < d max,V EN , to preserve the boundary between WM and CSF : The relationship function for sub-cortical structures are unchanged, as they legitimately share a boundary with the ventricles.The distances d max,V EN and d max,GM used must be small, as lesions may still appear in the vicinity of the boundaries.In our experiments, we set d 2 max,V EN = d 2 max,GM = 3 voxels 2 .However we chose the distance from inter-ventricular WM bigger because lesions can not appear too close to the boundary of this region, we chose d max,W M = 5 voxels.
Finally, an initialization for the intensity centroids of the different structures is necessary.We employ a simplistic intensity template of the expected centroids for each modality.When processing a multichannel image, we first estimate the robust minimum and maximum of the intensity in each channel (the intensity values at 5% and 95% of the histogram, respectively), and normalize the profiles such that 0 corresponds to the minimum and 1 to the maximum.We initialize the centroids to empirically determined values stemming from the expected intensities of the three major tissue classes and lesions.After the centroids are initialized, however, they are allowed to evolve freely to minimize the energy of Eq 1.In some data sets where the the initial segmentations are particularly poor and intensity inhomogeneities are strong, the lesion centroid computation can become unstable.To address this issue, we employ a weighted update scheme to compute the centroids for only the lesion class as follows: λ is chosen with respect to the amount of inhomogeneity in data set.We set λ = 0.8 in this work.

Experiments
We applied our method to the two different studies of MS patients from the MICCAI Challenge.Study 1 includes 10 patients while Study 2 consists of 14 cases.In both studies T1, T2 and FLAIR acquisitions with 0.5mm isotropic resolution are available.To reduce the computational burden, we subsampled the data to 1mm cubic volumes.Results were upsampled to the original resolution after processing.Two data sets were manually segmented by two human experts and the automated segmentation was compared with the segmentation from these two raters.Four metrics were computed to evaluate the performance of automated methods: Volume Difference, Average Distance, True Positive and False Positive.The volume difference captures the absolute percent volume difference to the expert rater segmentation.The average distance captures the symmetric average surface distance to the expert rater segmentation.True positive rate is the number of lesions in the automated segmentation that overlaps with a lesion in the expert segmentation divided by the number of overall lesions in the expert segmentation.False positive rate is the number of lesions in the automated segmentation that do not overlap with any lesion in the expert segmentation divided by the number of overall lesions in the automated segmentation.Moreover, using all manual segmentations as well as the high quality segmentations provided by several automated methods, a consensus segmentation has been computed via the STAPLE methodology.The STAPLE values of sensitivity, specificity and predictive value (posterior probability value) were also computed.Positive predictive value (PPV) is the ratio of true positives to the sum of true positives and false positives.Finally all metrics are scored in relation to how expert raters compare against each other.A score of 90 for any of the metric would equal performance akin an expert rater.An overall score of 80 was achieved by our method.Table1 shows the complete results of our method.The algorithm seemed to perform poorly on 4 cases.On Case CHB03, the volume difference score was lower than expected.In this case, our method found some false positives on the WM/GM boundary.This problem was likely due to the low quality of the T1 channel which makes the down-weighting of the GM boundary less accurate.Moreover, we observed that the FLAIR channel suffers from a large amount of inhomogeneity that caused some of the voxels on the GM/WM boundary to have an intensity close to that of lesions.
On cases CHB05, CHB11 and UNC09, the method achieved low scores for True and False Positives.These three cases posses a very low amount of lesions.In these cases our method also segmented very few voxels as lesions which causes the centroid computation to be somewhat unstable.The centroid for the lesion class can easily move far away from its optimal value given a few misclassified pixels.The large inhomogeneity in the FLAIR channel also contributed to this problem.We expect a more robust estimation procedure for the centroids and perhaps a more sophisticated inhomogeneity correction scheme would be needed.In our future work, we will also investigate the automated identification of such cases in order to select the most adequate parameter settings from the data.
Additional validation results on computational phantoms and real data are described in [20].In this work, a previous version of the algorithm achieved a Dice coefficient of 0.720 for lesions and Dice coefficients of 0.900, 0.925, 0.899, 0.720, 0.871, 0.879 and 0.774 for 7 other structures (cortical csf, cerebral WM, cerebral GM, brainstem grouped with cerebellar WM, cerebellar GM, ventricles, sub-cortical structures including  caudate, putamen and thalamus, respectively) when applied to the Brainweb MS phantom [21] under default noise and inhomogeneity settings.Moreover, an intraclass correlation of 0.772 was achieved on a data set of 10 real MR images acquired from MS patients.Overall, the algorithm was shown to be comparable to other state of the art algorithms while yielding a topologically consistent segmentation of multiple structures.

Conclusion
In this paper, we presented a new fully automatic segmentation technique for detecting MS lesions.The performance of the method was evaluated on data sets from two MS patient studies, and exhibited good performance on most of the cases.The proposed technique not only gives us a reasonable segmentation of lesions, but also gives us a topologically correct segmentation of the main brain regions, unlike previous lesion segmentation algorithms.Preserving topology has been of central interest in neuroimaging for applications ranging from computational anatomy to cortical reconstruction and mapping.Until now, these methods have focused on brains devoid of lesions and other structures that would alter the overall topology.In Multiple Sclerosis, brain atrophy has been observed in addition to the appearance of lesions.The proposed algorithm will enable volumetric analysis, cortical thickness analysis, and diffeomorphic shape analysis in MS populations, potentially leading to new insights on the disease.

Figure 1 :
Figure 1: MS lesions contrasts in MR images: hypointense WM in T1 looks similar to GM, hyperintense WM on T2 and PD looks similar to CSF.Lesions are brighter than other tissues on FLAIR, but the boundary of the ventricles is also hyperintense (see areas pointed by red arrow).

FLAIRFigure 2 :
Figure 2: Example of removing false positives using an interventricular-WM class.White arrows show the common false positives (caused by hyper-intensities between ventricles) which have been removed by adding the extra class.

Fig 3 and
Fig 4  shows an actual result computed using our method on selected cases from the Study 1 and Study 2, respectively.

Figure 3 :
Figure 3: Example of simultaneous tissue classification and lesion segmentation from Study 1, showing the classification of major structures, as well as segmented lesions and inter-ventricular WM.

Table 1 :
Result of comparison of automated lesion segmentation with manual segmentation of two human experts