Fully automatic brain segmentation using model-guided level sets and skeleton-based models

. A fully automatic brain segmentation method is presented. First the skull is stripped using a model-based level set on T1-weighted inversion recovery images, then the brain ventricles and basal ganglia are segmented using the same method on T1-weighted images. The central white matter is segmented using a regular level set method but with high curvature regulation. To segment the cortical gray matter, a skeleton-based model is created by extracting the mid-surface of the gray matter from a preliminary segmentation using a threshold-based level set. An implicit model is then built by defining the thickness of the gray matter to be 2.7 mm. This model is incorporated into the level set framework and used to guide a second round more precise segmentation. Preliminary experiments show that the proposed method can provide relatively accurate results compared with the segmentation done by human observers. The processing time is considerably shorter than most conventional automatic brain segmentation methods.


Introduction:
Automated brain segmentation is an important tool for various image-based research of brain structures and diseases, such as Alzheimer's disease and multiple sclerosis (MS).Many brain segmentation and tissue classification methods have been proposed and achieved different level of success in a wide variety of applications; some examples can be found in [1][2][3][4][5].One family of the segmentation methods is the socalled model-based approach, which uses statistical shape priors to guide the evolving contours [6].This approach has been proved to be successful for fully automatically sub-cortical structure analysis [5].However, when it comes to segmentation of the cortical gray matter, such an approach fails to deliver satisfactory results due to the high variation of the brain folding.In [7], we proposed a method to use a skeletonbased implicit model for vessel segmentation where the anatomical variation is also too large to be represented by conventional models.In this paper, we extend the line skeleton to a surface skeleton to generate a personalized gray matter model for each patient and use it to guide the segmentation of the cortical gray matter.The proposed method is implemented in a fully automatic brain segmentation framework where skull stripping, brain ventricle segmentation and basal ganglia segmentation are made using a variation of the conventional statistical model-based level set method.The algorithm was trained on 5 manual segmented brain dataset and tested on 12 datasets provided by the organizer of the MR brain image segmentation challenge [8].
Preliminary results show that the proposed method can achieve reasonable accuracy with much less running time than some conventional methods.

Method:
The proposed brain segmentation method can be divided into four steps.Firstly, the skull is stripped using a model-based level set on T1-weighted inversion recovery images.Secondly, the brain ventricles and basal ganglia are segmented using the same method on T1 weighted images.Thirdly, the central white matter is segmented using a regular level set method but with high curvature regulation.Finally, to segment the cortical gray matter, a skeleton-based model is created by extracting the mid-surface of the gray matter from a preliminary segmentation using a thresholdbased level set.An implicit model is then built by defining the thickness of the gray matter to be 2.7 mm.This model is incorporated into the level set framework and used to guide a second round more precise segmentation.The details of these steps are further described in the following sections.

Skull stripping
In the proposed framework, the statistical model-guided level set method proposed by Leventon et al. [6] is used to strip the skull on the T1-weighted inversion recovery images.In contrast to [6], a threshold-based external speed function is used instead of the gradient-based speed function.The level set equation is summarized in Eq. 1, where C is a clamping function, I is the input image,  is a threshold, M is the statistical model, T is the global transformation and κ is the mean curvature.As described in [6], the statistic model is created by taking the mean of the signed distance functions of each segmented regions ( ) and n prominent variations extracted via Principal Component Analysis (PCA) (M σ1 , M σ2 , … M σn ).Then the model M that fits the current shape can be represented as a weighted combination of these components (Eq.2).The transformation T and the weighing factors  i can be solved by minimizing the squared distance between the model and the level set function, which is also a signed distance map.The optimization of the level set function and the model is usually performed iteratively at the same time, i.e. the model is re-estimated after one or several iterations of the level set evolution.In our implementation, we used the coherent propagation method to evolve the level set function, which is about 10 times faster than the usual sparse field level set algorithm [9,10].It also has the ability to detect convergence of the level set function.This allows us to run the level set evolution until convergence and then performing the model-estimation step and re-running the level-set evolution, therefore it reduces the number of the model-estimation step.To strip the skull,  is set to 2000 in all cases, and α is set to negative so the contour will expand when the voxel's intensity is below the threshold.The statistical model is created from 5 manually segmented brain datasets provided by the challenge organizer.A similar strategy as in the first step is used to segment the brain ventricles and the basal ganglia in T1-weighted MRI images.For the ventricles, the threshold  is set to (P GM +P CSF )/2, where P GM is the intensity of the gray matter peak in the histogram (Fig. 1a) and P CSF is the estimated intensity of cerebrospinal fluid (CSF), which is set to 40 for all cases.For basal ganglia segmentation, two-sided threshold function is used where the lower threshold is set to (P GM +P CSF )/2, the upper threshold is set to V GW , which represents the local minima between the gray matter and white matter peaks in the histogram of T1 images.The statistical models are also extracted from the 5 manually segmented brain datasets.

Central white matter segmentation
To segment the central white matter, the segmented ventricle and basal ganglia regions in the T1-weighted images are set to P WM (representing the intensity of the white matter peak in the histogram of T1 images, cf.Fig. 1a), as the example shown in Fig. 2a.MS lesions are also segmented in the T2-Flair images using a simple thresholding method where the threshold is set to P Flair + 3 Flair .Here P Flair is the intensity of the brain tissue peak in the Flair image, and  Flair is the estimated standard derivation of that peak (Fig. 1b).The corresponding regions of MS lesions in T1 images are also set to P WM (e.g.Fig. 2a).A conventional threshold-based level set method (β in Eq. 1 is set to zero) is used to segment the central part of the white matter in the modified T1 images.Here  is set to V GW , and α is a fixed positive value.A relatively high curvature force (γ = 0.7) is used to avoid holes appearing around the aforementioned segmented structures and lesions.The segmented central white matter region in the T1-weighted images is again set to P WM for the subsequent cortical gray matter segmentation (Fig. 2b).

Cortical gray matter segmentation
As the variation of brain folding is too complicated to be represented by statistical models, we propose to use a surface-skeleton-based model to guide the level set segmentation.Similar to the line-skeleton based model used for vessel segmentation in [7], the implicit gray matter model can be represented using a signed distance transform from the mid-surface of the gray matter.To generate the mid-surface, a threshold-based level set was run twice to segment both the interface between gray matter and white matter and the interface between gray matter and CSF on the modified T1 image from the previous step (Fig. 2b).The threshold is set to VGW and (P GM +P CSF )/2 for those two interfaces respectively, and γ is set to 0.2.The signed distance maps (positive inside) from these two interfaces are denoted by D in and D out , respectively.The mid-surface is thought to be at the zero level of the combined distance map D c using Eq. 3, where  is the thickness of the gray matter (set to 2.7 mm for all cases).
The mid-surface is finally extracted via a third level set segmentation using D c as the external speed function.Note that in the area where the distance between two interfaces is greater than , a platform of zero level will appear in D c (cf. Fig. 3a).This happens often in the tips of a narrow inward/outward folding where the inner or outer surface will stop growing prematurely due to the high curvature force.In such places, the speed function switches to a regular threshold function of the input T1 image (first term in Eq. 1) and  is set to P GM .For this particular step, the input image is smoothed with a Gaussian kernel, so that the gray matter close to white matter will have higher intensity than the one close to the CSF.This strategy helps to push the mid-surface towards the folding direction in most cases (Fig. 3b).After obtaining the mid-surface and its signed distance map D mid , an implicit membrane model can then be built via Eq. 4.
The inner and outer interface segmentation of gray matter is re-run, but this time the model term is added to the level set equation (Eq.1). the generated surface-skeleton (red: white matter; green: gray matter; blue: the mid surface), the arrows point to areas where the white matter segmentation is incorrect and the arrowheads point to areas where the gray matter segmentation is incorrect.Note that the mid-surface is estimated correctly using the proposed method.c.An example of the final brain segmentation.

:
The proposed algorithm was trained on 5 manually segmented brain datasets and tested on 12 datasets provide by the organizer of the MR brain image segmentation challenge [8].As the input images are 3 mm thick slices, they were first resampled to an isotropic resolution of 1 mm using the sinc interpolation method.The output of our segmentation contains three tissue types, CSF (including ventricles), gray matter (including basal ganglia) and white matter.The resulted label images were then down-sampled to the original resolution using the nearest-neighbor method and compared with the manual segmentation results.The Dice coefficients, the 95th percentile of the Hausdorff Distance (HD) and the percentage of Absolute Volume Difference are measured as described in [11] and listed in Table 1.The overall running time of all steps was 2.5-3 minutes on a PC with an Intel i7 CPU and 8G RAM.

Discussions and Conclusion:
The accuracy of the gray matter and white matter segmentation from the proposed method is close to the performance of some widely used brain segmentation software, such as SPM8 and FAST, reported in literature [4,12], but the running time of the proposed method is much shorter than the latter.However, direct comparison will require running that software against the same datasets, which is outside the scope of this paper.
In addition to the tissue classification, the proposed method also outputs the surface-skeleton, which seems to be relatively accurate as judged by preliminary visual examination.It is worth notice that the 95th percentile of HD is around 3 mm, or even less for the gray matter and white matter segmentation, which is close to the image resolution along the z-axis.This, to some extent, confirms that the generated gray matter model is close to the ground truth.This skeleton itself can be useful for analyzing the geometrical features of brain folding.
One limitation of our preliminary validation is that the number of training datasets for generating the statistic model is very small and does not include significant pathological changes like ventriculomegaly.This may limit the algorithm to work only with relatively normal brain data.Large errors in the basal ganglia segmentation were observed in patients with enlarged lateral ventricles.However, the ventricle segmentation seems to work reasonably well in such cases, thanks to the good contrast between CSF and brain tissue.
Resampling is another factor that may affect the evaluation of the proposed method.In theory, the proposed method works better on isotropic datasets.However, the input images used in our experiment were up-sampled from 3 mm scans.Even though high-resolution 1 mm isotropic T1 images were provided in the challenge, they are not well aligned with the 3 mm scans on which the manual segmentation was performed.Applying the proposed algorithm on thin-slice T1 resulted in lower DICE coefficient (about 5% for gray matter) than on the resampled 3 mm T1 images.Some other limitations of the proposed method include using a fixed and unified gray matter thickness in the skeleton-based model and unified thresholds for tissue classification.Introducing some local adaptive strategy and combining image gradient information into the level set function may improve segmentation accuracy.
In conclusion, a fully automatic brain segmentation method has been presented.Preliminary results show that the proposed method can achieve reasonable accuracy with much less running time than some conventional methods.In addition, a surface skeleton of the gray matter which represents the geometric features of brain folding can be obtained.

Fig. 1 .
Fig. 1. a.An example histogram of T1 weighted images, P GM and P WM are the intensity of the gray matter peak and white matter peak.V GW is the local minima between the gray matter and white matter peaks.b.An example histogram of T2 Flair images, P Flair is the intensity of the brain tissue peak.

Fig. 2 .
Fig. 2. a. Modified T1 image before central white matter segmentation.b.Modified T1 image after central white matter segmentation.

Fig. 3 .
Fig. 3. a.An example of the combined distance map D c .The blue lines represent the interfaces, the dash lines indicate the zero level of the map.b.An example of the initial segmentation andthe generated surface-skeleton (red: white matter; green: gray matter; blue: the mid surface), the arrows point to areas where the white matter segmentation is incorrect and the arrowheads point to areas where the gray matter segmentation is incorrect.Note that the mid-surface is estimated correctly using the proposed method.c.An example of the final brain segmentation.

Table 1 .
Automatic segmentation results compared with manual segmentation results