Novel solutions toward high accuracy automatic brain tissue classification in young children

This paper addresses a need for the development of automatic and accurate brain tissue classification techniques and validation of classification in MR brain images obtained during early childhood development. Namely, we consider the National Institutes of Health (NIH) pediatric database that represents the largest demographically diverse U.S. pediatric cohort with respect to ethnicity, gender, race and income level (Almli et al.,, 2007, Evans,, 2006). Within the larger cohort from birth to adolescence, the ”Objective-2” (O2) MR subset consists of 69 typically developing subjects in the age range from 10 days to 4.5 years scanned longitudinally.
Almli et al., (2007) and Fonov et al., (2011) have reported that the dynamic changes in intensity and brain shape with age in the O2 data pose a challenge to using general anatomical image processing pipelines such as FSL (Smith et al.,, 2004), CIVET (MacDonald et al.,, 2000) and statistical analysis tools provided by SPM (Ashburner and Friston,, 1997).
There is a high variability of MR signal across times and within each tissue class brought about by rapid changes in signal intensity of grey matter and white matter (Figure 1.1). In particular, in the course of time white matter reverses in intensity as axons become myelinated from occipital to frontal lobes, and later from central to subcortical white matter.
Secondly, the contrast-to-noise ratio between grey matter and white matter in the pediatric brain MRI can be as low as half of that in adult MRI. Automated state-of-the-art classification algorithms fail to delineate accurately tissue class clusters in brain subdomains with significantly overlapping intensity histograms.
Thirdly, due to the risk of motion in non-sedated infants during brain scanning, a trade-off between spatial accuracy and motionless acquisition time was achieved at a resolution of $1\times 1\times 3\>mm^{3}$. The MRI resolution is critical to the quality of segmentation. Thicker axial slices accentuate the partial volume artifact, a voxel that contains signal from two or three different tissues.
There is a need for new quantitative analysis techniques for the O2 data that will accommodate their unique characteristics such as high tissue intensity variation and changing grey/white matter contrast. As these changes occur rapidly over time during brain development, analytical methods for early pediatric MRI must employ age-specific anatomical templates (Figure 1.1).

In this paper, we implement a novel KFDA-based classification method (Portman et al.,, 2014) to automatically generate GM, WM and CSF labels (Figure 1.2) in age-dependent pediatric O2 templates. We started with the youngest brain template for ages 4.5 to 8.5 years from NIH Objective-1 (O1) pediatric database (Fonov et al.,, 2011) that represents the average anatomy closest in age to the oldest brain template for ages 44 to 60 months from the O2 database. For the initial estimation of GM, WM and CSF spatial patterns, brain tissue probability maps were transferred from the youngest O1 T1w template to the oldest O2 T1w template using mni$\_$autoreg (Collins et al.,, 1994). Then, hard GM, WM and CSF labels were obtained from the transferred tissue probability maps by choosing the tissue class with a highest probability at a given brain voxel.

Having combined the best guess on GM, WM and CSF location with multi-channel, T1-, T2- and PD-weighted template volumes, tissue labels for the oldest O2 template were generated using the KFDA algorithm. The process can then be repeated to map these labels to the second oldest template for the age range 33-44 months and so on in a pairwise comparison to the lowest age group.
Figure 1.1 shows drastic intensity changes in WM as it reverses in intensity due to the myelination process and complicates the pairwise template diffeomorphic registration. The subcortical WM of the frontal and temporal lobes is not myelinated for the age range 5-8 months and the T1w signal of WM in these regions is isointense with GM. For ages 2 to 5 months the myelinated area spreads from the internal capsule to the occipital lobe along cortico-spinal tracts, and for an earlier age it occupies the internal capsule only. For these early childhood templates, an extension of the proposed algorithm to identification of myelinated and unmyelinated WM in addition to GM and CSF is needed, and is beyond the scope of this paper.

The development of segmentation techniques in newborn and infant brain MRI over the past decade has shown that atlas-based segmentation methodology for adult brain MR images can be applicable to child brain MRI when tuned to the specific properties of the dataset under study. There are three major child MRI segmentation frameworks, Expectation-Maximization, Registration-based method and Adaptive Label Fusion:
Expectation-Maximization (EM) framework. This widely used framework assumes a global Gaussian mixture model of the distribution of major tissue class intensities (Van Leemput et al.,, 1999). The EM method relies on spatial priors from probabilistic atlases that describe anatomical variability of the brain population and constrains classification during parameter estimation step.
Since application of adult prior probability maps to pediatric brain atlas construction generates errors due to anatomical differences, studies have focused on generation of a probabilistic atlas that would appropriately capture brain anatomy in early childhood. In particular, age-specific atlases of GM, WM and CSF tissue classes for O1 database of children aged 4.5 to 18.5 years have been created (Fonov et al.,, 2011).
Various techniques for creation of population-specific brain atlases for MRI study of early brain development have been proposed (Evans et al.,, 2012, Murgasova et al.,, 2007, Kuklisova-Murgasova et al.,, 2011, Prastawa et al.,, 2005, Xue et al.,, 2007). Most of these atlases are derived from small cohorts that do not capture the full range of normal anatomical variability. Furthermore, they do not cover the entire O2 developmental epoch from birth to 4.5 years of age. The creation of such standard atlas still remains a significant challenge since it requires semi-automated segmentation of large datasets.
The successful EM classification of child brain MRI is achieved at a huge computational cost associated with building an appropriate brain atlas and overcoming restrictions of the global Gaussian mixture modelling (Prastawa et al.,, 2005, Van Leemput et al.,, 1999, Wells et al.,, 1996, Xue et al.,, 2007, Tohka et al.,, 2004, Choi et al.,, 1991). In reality, voxel intensity distribution within infant brain MRI differs from a Gaussian due to partial volume effects and biological intensity variation.
Registration-based framework. Registration-based brain segmentation methods are seemingly attractive in applications to the population of young children since tissue intensity variations do not affect the segmentation results and there is no need for tissue intensity probabilistic models. A powerful technique known as ANIMAl+INSECT (developed by (Collins et al.,, 1999)) non-rigidly registers an adult brain atlas to an individual subject via ANIMAL while simultaneously correcting for intensity inhomogeneity in the subject. The ANIMAL procedure iteratively estimates a 3D non-linear deformation field iteratively in a multiscale hierarchy by optimal matching of both Gaussian blurred volumes and subsequent refining of the resulting displacement field. The estimated transformation is then used for tissue label transfer from the atlas to the subject.
(Murgasova et al.,, 2007) examined the suitability of a registration-based technique for segmentation of child brain MRI at ages 1 and 2 years. The method outperformed EM (Van Leemput et al.,, 1999) in segmentation of subcortical brain structures (such as the thalamus) that exhibit substantial intensity variation. However, EM was more successful in cortical areas due to the difficulty of registering complex cortical folds in atlas and subject brains.
As with EM, a barrier to implementation of registration methods is the dependency on a pediatric brain atlas with accurate measures of spatial boundary uncertainty for tissue classes that the O2 dataset does not possess.
Adaptive label fusion. Adaptive label fusion is a machine learning approach that estimates subject-specific tissue intensity distributions non-parametrically and therefore, captures complex distributions of voxel intensities in young child MRI (Weisenfeld and Warfield,, 2009). Given a library of newborn template MRI each containing manually selected tissue class prototypes, the tissue prototype labels of each template are transferred to the subject through the registration process. Then the conditional (on the subject intensities) probability density of the tissue labels is non-parametrically estimated based on each template prototype labels leading to candidate segmentations of the subject. The candidate segmentations are then fused using the maximum-likelihood method that computes the label weights based on how reliable the estimated candidate segmentation is.
The label fusion technique uses global intensity statistics that requires for each tissue type to produce similar intensities in different parts of the image. To manage the tissue intensity variability and partial volume effects a probabilistic atlas constructed from 15 newborn brain MRI is incorporated into the fused segmentation model (Warfield et al.,, 2000).
The label fusion approach is flexible as it may or may not be atlas-dependent. However, it requires a library of infant brain MRI templates of the appropriate developmental age with manually chosen tissue prototypes that is not yet available for the pediatric O2 dataset.

Portman (Portman et al.,, 2014) introduced a local atlas-independent framework using modern trends in Computer Vision such as Kernel Fisher Discriminant Analysis (KFDA) for pattern recognition and structural similarity index (SSIM; (Wang et al.,, 2004)) for perceptual image quality evaluation. In this framework, within-class tissue intensity variations are handled by optimal partitioning of the brain into overlapping subdomains having different average intensities.
KFDA-based methods explore complex non-linear structures of the tissue clusters in the input or intensity space to find optimal separating surfaces between the class clusters. KFDA attempts to make the input data more separable by non-linearly mapping them from the input intensity space to an abstract, high-dimensional feature space where the classification is performed.
Furthermore, the proposed method automatically identifies voxel categories such as PVE voxels and tissue class prototypes, and predicts the most likely membership of PVE voxels from the set of prototypes.
We aim to classify age-specific pediatric O2 templates into two major tissue classes (GM, WM) and the CSF. The classified representative templates can then be used for the classification of a subject brain MRI of a developmental age similar to the age range of the template.
An important question that arises with classification is validation. In the absence of a “ground truth”, we assess the quality of classification via SSIM that predicts image quality as perceived by the Human Visual System (HVS) (Wang and Bovik,, 2009) which is regarded as an optimal information extractor that seeks to identify objects in the image. As such, the HVS must automatically identify structural distortions and compensate for the nonstructural distortions (e.g., image corruption by noise). The SSIM is an objective image quality metric that simulates this functionality and computes the degree of structural similarity between reference and distorted images. It has been shown that the SSIM is well-matched to the subject ratings of image databases, and therefore, it is a good approximation to the perceived image quality (Wang et al.,, 2003).
Given the objectivity and effectiveness of the SSIM we apply it for comparison of Partial Volume Estimation (Tohka et al.,, 2004) and KFDA-based classification algorithms, as well as for automatic monitoring of the quality of classification. That is, we compute the structural closeness of classified 3D brain images containing GM, WM and CSF mean intensity values with their counterparts seen in an MR image type of a highest contrast and regarded as references. For the age range from 8 to 60 months, we rely on T1w as the most informative of all MR image types. However, for an earlier age group T1w template serves as a reference only for the CSF pattern (Figure 1.1) due to a poor GM/WM contrast. For GM and WM, T2w data can be used as a reference since they exhibit higher GM/WM contrast.

In the following we will first illustrate the pipeline structure of the proposed segmentation algorithm, and discuss a need for local segmentation. We will then describe the algorithm step by step, namely, optimal brain partitioning, subsequent delineation of the CSF and WM clusters, and image stitching. Finally, we will report on classification results for O2 brain templates for ages from 44 to 60 months and for ages 8 to 11 months, compare the performance of PVE and KFDA methods and discuss possible improvements and further applications of our method.

The problem of classification into GM, WM and the CSF is partitioned into two binary classification problems.

First, the multichannel (T1w, T2w, pdw) image data are separated into the CSF and G+WM using a kernel function that best describes the non-linear behaviour of the CSF cluster in the 3D intensity space. We postulate that in this way we can delineate the CSF with a higher accuracy. As reported in Bouix et al., (2007) the major drawback of the existing automatic classification algorithms is incorrect classification of the CSF.
Second, the G+WM data is separated into GM and WM. These are major steps of the proposed KFDA-based method whose pipeline structure is illustrated in Figure 2.1.

Healthy children aged 10 days to 4.5 years (n=106) were enrolled for the NIH-funded O2 MRI study of early normal brain development. In this project, T1w,T2w and pdw data were acquired approximately at quarterly time intervals (Almli et al.,, 2007) at a 1.5 T Siemens Sonata scanner with a $1\times 1\times 3$ mm spatial resolution. The data were resampled to $1$ mm³ grid using tri-cubic interpolation. The brain image quality control was applied to the MR data that reduced the original sample to 69 subjects. T1w, T2w and PDw average atlases have been created for these age ranges (Fonov et al.,, 2011). All images were N3 non-uniformity corrected (Sled et al.,, 1998), linearly normalized to have the same intensity range as the ICBM152 template by a linear histogram scaling (Nyùl Làszlò et al.,, 2000) and registered to the ICBM152 stereotaxic space using MINC mni_autoreg non-rigid registration tool (Collins et al.,, 1994). Brain masks were created using BET from the FSL package (Smith,, 2002). Brain masking was applied to all templates to remove non-brain tissue prior to classification.

We examined infant MR image data and made a simple but important observation that led to the optimal brain partitioning algorithm presented below. For the O2 T1-weighted template for ages 8 months to 11 months with GM, WM and CSF probability maps obtained by global PVE we performed hard labelling into GM, WM and CSF by choosing the highest probability class at each brain voxel. We took samples of about 10,000 voxels contained in disjoint parallelepiped-like brain subdomains and plotted the corresponding T1w intensity profiles along with the tissue and subdomain intensity means as seen in the low panel of Figure 2.2.
For simplicity, we illustrate a variable behaviour of these brain data samples using examples of two brain subvolumes denoted by A and B, respectively, and shown in the upper panel of Figure 2.2. Notice by comparison of diagrams 2.2.c and 2.2.d that the tissue means are approximately the same in their values as a result of global estimation of the probability distribution of tissue classes across the entire brain. However, subdomains A and B differ in their overall mean intensities.

The observations of various brain data samples suggest that there is more uncertainty in identification of tissue boundaries in subvolumes with lower overall means. Such samples usually contain lower intensities of both GM and WM classes.
Figure 2.3.a shows the classified image subdomain A with underestimated WM and the CSF. In this subdomain, WM has a range of intensities appearing bright in the central part of the brain and fading along WM fibers. A low intensity signal of myelinated WM can be seen in the temporal lobe. As expected, it is not detected by the global PVE due to low GM/WM contrast in the temporal brain region. Similarly, a low GM/CSF contrast causes misclassification of the CSF into GM. We conclude that a lower subdomain intensity mean is an indicator of a “problematic” subvolume with a lower contrast. Conversely, the global classification yields quite accurate WM, GM and CSF patterns in subdomain B with a higher mean intensity as observed by comparison to GM, WM and CSF that we can visually identify in T1w data of subdomain B given in Figure 2.2.b (Figure 2.3.b).

We then applied a KFDA-based algorithm (Portman and Evans,, 2013) locally to classify low intensity image subdomain A into GM, WM and the CSF. We used the PVE-based classification given in Figure 2.3.a to initialize our algorithm.

Figure 2.4.a shows a significant improvement in tissue identification over the global PVE. The local KFDA approach reveals more complex WM structure with WM streaks dispersing from the central part of the brain and the presence of more myelinated WM in the temporal lobe. In subvolumes with high intensity means (consisting of mostly WM) the WM pattern appears no different from the one detected by the global PVE as seen in Figure 2.4.b.
Furthermore, we determined experimentally that compared to the global KFDA (Portman and Evans,, 2013) the local KFDA classification is capable to capture spatial tissue patterns with a higher accuracy. A comparison of tissue classifications obtained by global and local KFDA methods in the O2 brain template from 44 to 60 months of age is shown in Figure 2.5. Observe that WM fibers in the posterior brain are detected by the local KFDA approach and omitted by the global one.

Therefore, we tackle intra-class MR signal variability by splitting the brain into parallelepiped-like subdomains having different average T1w intensity values. In this way, we are able to control the quality of classification locally, namely, in subdomains with low GM/WM or/and GM/CSF contrast.

The algorithm for 3D brain image partitioning used here is based on maximization of mutual information (MI) between the histogram bins of the T1w image and the subdomains of the partitioned image. The algorithm progressively subdivides the 3D image into subdomains using binary space partitioning (Rigau et al.,, 2004), that is, each subdomain is subdivided into two subdomains with coronal, sagittal or axial planes according to the maximum MI gain. At each partitioning step, the MI optimality criterion finds a planar cut separating each parallelepiped-like volume into two subvolumes whose intensities best match the histogram bins (splitting the data into two non-overlapping intensity ranges). Thus, the obtained subvolumes will have different average intensity values.
To maintain the continuity of the classified image subdomains across their boundaries, we added two slices to each planar boundary of each subdomain thus creating an overlap of 4 slices between adjacent brain subdomains.
Examples of the optimally partitioned brain templates for age groups 8 to 11 months and 44 to 60 months are shown in Figures 2.6.a-b. Since the partitioning algorithm is intensity-based, different brain regions and their total number can be obtained for different individual brain scans.

The optimal brain partitioning algorithm allows to identify lower intensity regions. The computed GM/WM contrast-to-noise ratio (CNR) in image subdomains of the template for ages 8 to 11 months show that GM/WM CNR takes lower values in the posterior brain and higher values in the central part of the brain (Figure 2.7).

The fundamental idea of KFDA is to exploit naturally occurring class structures in the input space to find an accurate non-linear separating boundary. We performed KFDA in two steps.
Step1: Classification into the CSF and G+WM. Through extensive experimentation we determined that the sigmoidal kernel function $K(\vec{I}_{m},\vec{I})=\tanh(a({\vec{I}_{m}}^{T}\cdot\vec{I})+b)$ (with $a$ and $b$ being user-specified parameters) is the best choice for delineation of the CSF. Since KFDA computes an optimal decision surface between the CSF and G+WM it easily identifies PV voxels that lie near or on the boundary between the classes.

A visualization of the overlapping set (coloured in black and green) in the stereotaxic space given in Figure 2.8 shows that the overlapping voxels are located within the boundary regions between G+WM and the CSF. We recognize a particularly problematic brain area around the sulcus where the CSF is usually poorly detected due to the presence of PVE. Therefore, it is reasonable to assume that the overlapping voxels contain a mix of both tissue classes.

Step2: Classification into GM and WM.
Having delineated the CSF we classify G+WM into GM and WM. Experiments with various kernel functions showed that the Gaussian radial basis function $K(\vec{I}_{i},\vec{I})=\exp\left(-\frac{(\vec{I}_{i}-\vec{I})^{T}(\vec{I}_{i}-\vec{I})}{2\sigma^{2}}\right)$ is the best choice to model the non-linear structure of WM and GM clusters.

We applied a Simulated Annealing technique to stitch the local classifications together into a cohesive global picture of the classified brain. First of all, for each brain slice we collected the constituent subimages. Each pair of overlapping subimages contained a joint image region $I_{s}$ of size $4\times ncols$ or $nrows\times 4$ that is to be optimally estimated from two sets of observations $I_{l}$, $I_{r}$ or $I_{up}$, $I_{low}$ as illustrated in Figures 2.9.a-b. We implemented the iterative Simulated Annealing (SA) algorithm (Grenander,, 1996) to find the most probable joint region for each pair of overlapping classified subimages in a brain slice.
An example of SA application is shown in Figure 2.9.c. The leftmost columns of $I_{r}$ and $I_{l}$ appear slightly different in presence of the CSF and GM.The rightmost columns of $I_{l}$ and $I_{r}$ only differ in the value of a single central pixel. The optimal labelling $I$ preserves the label configuration of its first and last columns as they appear in their respective overlapping regions $I_{l}$ and $I_{r}$.

We started with the upper left subimage and stitched its neigbouring subimages vertically and horizontally. Then for each of the neigbouring subimages we identified overlapping ones in both horizontal and vertical directions, and glued them together using the SA algorithm. By progressive stitching of constituent image parts we assembled the entire brain slice as illustrated in Figure 2.9.d. The proposed SA algorithm yields seamless estimates of the joint regions.

We applied the proposed segmentation pipeline to the O2 3D brain template for ages 44 to 60 months with initialization obtained by a label transfer from the brain template for ages 4.5 to 8.5 years. Figure 3.1 illustrates the performance of KFDA using classified brain slices. The significant improvement in WM pattern detection was observed in the cerebellum (Figure 3.1.c) and validated using MSSIM (Figure 3.1.b). There is a remarkable closeness of the stripe-like WM pattern identified by KFDA to its poorly visible counterpart in a low contrast T1w brain slice.
The comparison of KFDA and initial classifications of the central axial slice demonstrates the adaptability of KFDA to individual brain tissue intensity distributions (Figures 3.1.d-f). The initialization derives from the labelling of the older brain closest in age to the brain template under study. With KFDA, we obtained more accurate CSF detection, in particular in areas between the insula and precentral gyrus where the initial CSF is underestimated. The anterior horns of the lateral ventricles and the third ventricle that are not present initially in the CSF are observed in the T1w input and captured by KFDA. Also, left and right frontal lobes appear to be separated by the CSF and are detected by KFDA.
Contours of the WM in the frontal lobe are more complex, having more visible cusps in the T1w image Figure 3.1.f). However, in the left frontal lobe, the WM is overestimated as it spreads into the middle frontal gyrus. WM fibers passing between the caudate nucleus and putamen appear connected in both the reference (T1w) and KFDA-classified images.

There are WM fibers separated from the main fiber bundle in the anteromedial temporal lobe with KFDA classification that are in agreement with the fading T1w signal intensity in this area of the brain. However, in the initialization generated by a label transfer from an older brain template they appear to be a part of the whole WM. Overall, KFDA tends to accurately reproduce the segmented brain anatomy according to its natural appearance in the input MR images.

We applied the local KFDA-based method to segment the O2 brain template for ages 8 to 11 months initially classified into GM, WM and the CSF using PVE. The proposed approach leads to a significant improvement in the CSF detection throughout the brain almost doubling the initial CSF volume. The initial CSF cluster determined by PVE consisted of 26717 voxels and has increased to 53681 voxels with application of KFDA.
The comparison of KFDA and PVE classified results clearly demonstrates that CSF is underestimated in PVE classification (Figures 3.2 and 3.3). KFDA reveals complex CSF patterns in the occipital lobe and on the border between cerebellum and temporal lobes as seen in Figures 3.3.c and 3.3.f. However, it does not capture the CSF completely on the border between the cerebellum and the left temporal lobe (Figures 3.3.c and 3.3.f) due to the less pronounced intensity differences between GM and the CSF in the left hemisphere.
Figure 3.2.e shows that in contrast to PVE, KFDA detects most WM seen in low-contrast temporal lobes and lower intensity WM streaks that stretch out from the central part of the cerebellar WM (Figure 3.2.a). This example demonstrates the capability of the proposed method to identify tissue patterns structurally similar to their counterparts that are visible in low contrast subdomains of T1w images.