A fast stochastic framework for automatic MR brain images segmentation.

This paper introduces a new framework for the segmentation of different brain structures (white matter, gray matter, and cerebrospinal fluid) from 3D MR brain images at different life stages. The proposed segmentation framework is based on a shape prior built using a subset of co-aligned training images that is adapted during the segmentation process based on first- and second-order visual appearance characteristics of MR images. These characteristics are described using voxel-wise image intensities and their spatial interaction features. To more accurately model the empirical grey level distribution of the brain signals, we use a linear combination of discrete Gaussians (LCDG) model having positive and negative components. To accurately account for the large inhomogeneity in infant MRIs, a higher-order Markov-Gibbs Random Field (MGRF) spatial interaction model that integrates third- and fourth- order families with a traditional second-order model is proposed. The proposed approach was tested and evaluated on 102 3D MR brain scans using three metrics: the Dice coefficient, the 95-percentile modified Hausdorff distance, and the absolute brain volume difference. Experimental results show better segmentation of MR brain images compared to current open source segmentation tools.

Introduction

Accurate delineation of brain tissues from magnetic resonance (MR) images is an essential step in human brain mapping and neuroscience [1-3]. However, brain MRI segmentation faces challenges stemming from image noise, magnetic field inhomogeneities, artifacts such as partial volume effects, and discontinuities of boundaries due to similar visual appearance of adjacent brain structures. This paper addresses brain segmentation from MR images at different life stages, having the infancy stage the most complicated one due to reduced contrast, higher noise [4], and inverse contrast between the white matter (WM) and gray matter (GM) [5], Fig 1. Segmentation of the brain at later stages might be relatively easier, as the contrast between different types of tissue is much better, and the signal-to-noise ratio (SNR) are improved, Fig 1.

Fig 1

T1-weighted MRI for adult (a) and infant (b) brains.

Lower contrast between infant brain tissue classes stems from the fact that most of the long axons in the WM are not yet myelinated, and WM water content is close to that of GM. Both the WM and GM have the same average T1-weighted intensity at about 9 months of age [6]; hence, it is difficult to classify the infant brain tissues using only the intensity. Furthermore, the unique challenge of imaging the infant patient precludes techniques that might improve contrast or signal-to-noise ratio (SNR) but would also lengthen the time of acquisition. [2]. Many segmentation techniques have been developed for the last two decades in order to address the brain MRI segmentation challenges. These techniques can be roughly classified into the following categories: (i) probabilistic, or statistical methods, (ii) atlas-based methods, (iii) hybrid methods that include (i) and (ii), (iv) deformable-model based methods, and (v) deep-learning based methods.

Probabilistic segmentation

These algorithms involve prior models that describe the signal distributions of each brain structure. Ng et al. [7] segmented MR brain images using the unsupervised K-means clustering of signals and an improved watershed algorithm. A similar approach by Xue et al. [5] employed a parametric Gaussian density estimation with an expectation-maximization (EM) algorithm and constrained spatial homogeneity of the MR images with a Markov random field (MRF) prior. Partial volume averaging effects have been eliminated by predicting the misclassification (e.g., of an “averaged” CSF and GM into an intensity similar to WM). An automated MRI brain segmentation method by Mayer et al. [8] combined spatial and intensity features into a high-dimensional feature space. An adaptive mean-shift classifier extracted a set of convergence modes, i.e. high-density points of a feature space, being good candidates for intensity-based classification. Brain tissues were classified by an intensity-based mode clustering. This approach was very effective with non-convex clustering. Fang et al. [9] developed a tree metrics (TM)-based graph cut algorithm to segment the MRI brain tissues. After a brain MR image is classified using the TM, the goal labeling is inferred by “tree-cutting”. In contrast to most of conventional iterative methods like the EM-based ones, which produce only locally optimal labeling, this algorithm needs no more than one sweep to generate the globally optimal labeling with respect to the TM. An automated segmentation approach by Ortiz et al. [10] classified the brain tissue with no prior information. The segmentation consisted of feature extraction and classification. Extracted first order (pixel/voxel-wise), second order (pair-wise), moment, and scale-invariant features were classified by growing hierarchical self-organizing maps (GHSOM). Li et al. [11] proposed a 3D MGRF model for the segmentation of brain MR images to avoid the shortcomings of the 2D model which is not able to fully capture the spatial information, especially among the slices. An initial segmentation was first obtained by k-means clustering in order to reduce the extensive computations required by the MGRF model. The Iterated Conditional Modes (ICM) algorithm was finally applied to obtain the optimal solution under maximum a posteriori (MAP) criterion. A non-parametric adaptive mean-shift algorithm was proposed by Janney et al. [12] for brain tissue segmentation. The method clustered the joint spatial-intensity feature space, followed by a phase of intensity-based mode clustering into the brain tissue types. Weber et al. [13] segmented the brain tissues using FSL-FAST free software that is based on K-means clustering which provided initial segmentation, followed by the EM algorithm for bias field correction. In order to speed up the process, parallelization to any eligible parts of the software was applied, which needed some adaptation to the algorithms in order to maintain the accuracy obtained by the software package. Mahmood et al. [14] proposed an unsupervised framework for brain tissue segmentation using a combination of Bayesian-based adaptive mean shift that clustered the tissues in the joint spatial-intensity feature space, and fuzzy c-means that is initialized with a priori spatial tissue probability maps to assign the clusters into three tissue types; WM, GM, and CSF.

Infant brain segmentation using statistical-based methods was also addressed in literature. Automated segmentation of brain structures, such as WM, CSF, central GM (CEGM), and cortical GM (COGM) was conducted by Anbeek et al. [15] using T2-weighted and inversion recovery (IR) MRI of the neonatal brains. Probability maps to segment each brain tissue class with a K-nearest neighbor (KNN) classifier using voxel intensities and coordinates as features were constructed manually. A multi-label segmentation process combined the obtained classes. Wang et al. [16] employed a random forest technique to integrate features from different modalities for brain tissue segmentation in infants along with probability maps of GM and CWM. Zhang et al. [17] proposed a deep convolutional neural networks (CNNs) approach for segmenting the neonatal brains from multi-modal MR images, generating the segmentation maps as outputs. The multiple intermediate layers included many operations such as convolution, pooling, and normalization in order to capture the highly nonlinear mappings between inputs and outputs. Moeskops et al. [18] segmented neonatal brains into WM, GM, and CSF using supervised voxel classification.

To recapitulate, statistical-based techniques are generally fast to implement compared to other segmentation methods. However, the fact that actual intensity distributions of brain structures are greatly affected by several factors, such as the unique patient and scanner along with scanning parameters, makes the segmentation hard. Also, due to the similar intensities for the different brain tissue structures of the infant MR brain images, segmentation techniques only based on the intensity remain inaccurate.

Atlas-based segmentation

Atlas-based approaches have emerged as powerful segmentation tools. These approaches are based on a priori knowledge about brain structures, and treat the segmentation problem as a registration one. Ashburner et al. [19] introduced a generative framework that combined image registration, tissue classification, and bias correction. Their framework incorporated a smooth intensity variation and nonlinear registration with tissue probability maps using mixture of Gaussians. Pohl et al. [20] introduced a Bayesian model for simultaneous segmentation and registration. Their framework tried to exploit complementary aspects of registration and segmentation problems. In order to account for different physiological (patient size and weight) and scanning (scanner type and data acquisition protocol) parameters, Han et al. [21] introduced an intensity re-normalization procedure to adjust the prior atlas intensity model to new input data to overcome the problems stemming from using training data acquired from a different scanner that was used for the test data. Artaechevarria et al. [22] proposed a generalized local weighting voting scheme in which the fusion weights were adapted for each voxel based on local estimation of the segmentation performance. The local weighting voting outperformed traditional global strategies that estimate a single value for the segmentation accuracy for the whole image. Sabuncu et al. [23] proposed an automated, label fusion segmentation technique. In order to capture greater inter-subject anatomical variability, each training data set was individually co-registered to the test data set. Then, a nonparametric probabilistic model was employed to fuse the training labels to compute the final segmentation. Morin et al. [24] presented an atlas-based segmentation framework using random walks that combined registration and labeling propagation steps. They used a generative model to provide pixel label probabilities to improve the segmentation for high-confidence labels. To match the target images with atlas images, they used the Affine-Scale Invariant Feature Transform (ASIFT) [25] and Speeded Up Robust Features (SURF) [26] registration techniques. In order to avoid segmentation errors produced by registration imperfection, Lötjönen et al. [27] introduced an optimized pipeline for multi-atlas brain MRI segmentation. They introduced two approaches that combine multi-atlas segmentation and intensity modeling based on using EM and graph cuts for optimization. First, they registered all atlases to the target data and a majority voting was applied to predict the segmentation of the target image. Then, the segmentation was improved using the intensity modeling as a post-processing step. Lijn et al. [28] introduced a segmentation method based on the combination of spatial features and appearance models. They generated a spatial probability map that was obtained from multiple atlas-target image registrations to implement the spatial model. The tissue appearance was modeled by a KNN classifier using Gaussian scale-space features. Then, a Bayesian framework was used to combine both spatial and appearance models and a graph-cut approach [29] was used for optimization. Ledig et al. [30] introduced a framework for labeling whole brain scans by incorporating a global and stationary MRF to ensure consistency of the neighborhood relations between structures with an a priori defined model.

Segmentation of neonatal brains was also conducted in the literature using atlas-based techniques. Segmentation of axial neonatal brain MRI that combined multi-atlas-based segmentation and supervised voxel classification was proposed, [31], in order to segment eight different tissue classes, namely cortical grey matter (CoGM), unmyelinated white matter (UWM), brainstem, cerebellum, ventricles, cerebrospinal fluid in the extra-cerebral space (CSF), basal ganglia, and myelinated white matter (MWM). Some approaches use longitudinal scans at a late-time-point age, where the contrast is much better between different tissue types, from which probabilistic atlases are constructed to guide segmentation of neonatal images [32, 33]. Cherel et al. [34] employed a subject-specific atlas that is based on manually segmented data for brain tissue classes segmentation. The atlases were incorporated with single-atlas expectation maximization (EM) method. Neonatal brains were segmented into 50 regions by [35], where structural hierarchy along with anatomical constraints were employed. Infant brain segmentation using shape priors was also addressed in the literature [36-38].

Atlas-based segmentation techniques show more accuracy compared to statistical-based techniques. Nevertheless, they are still challenged by atlas selection, combination, and the associated heavy computation time. Another major drawback of atlas-based segmentation algorithms is their dependency on the selected features that will be used to link between the test subject and the prior (training) data used in the construction of the atlas. This may lead to inaccurate segmentation, as signals vary due to many factors such as the patient’s age, and the scanning protocol.

Hybrid methods

The literature also shows methods that exploited both probabilistic models along with shape prior (atlas-based) ones. Song et al. [39] proposed a probabilistic neural network (PNN) for segmenting the brain MRI. Probability density functions of the brain tissues were estimated from reference vectors generated by a self-organizing map (SOM). To reduce the partial volume averaging effects, weighting factors were added to the summation layer’s patterns in a weighted probabilistic neural network (WPNN) and soft labeling was performed by a supervised Bayesian classifier. Patenaude et al. [40] proposed a method that used manually labeled image training data, where the principles of both the active shape and appearance models were utilized within a Bayesian framework, allowing probabilistic relationships between shape and intensity to be fully used. Serag et al. [41] employed high-dimensional feature vectors for segmenting brain subjects using a sliding window approach along with a multi-class random forest classifier. Wang et al. [42] segmented T1, T2, and diffusion-weighted brain images using a sparse representation of the complementary tissue distribution. Initially, the brain tissue was segmented into different structures using a patch-based technique with a library of multi-modality images, having been aligned with their ground-truth segmentation maps. Then the segmentation was refined by integrating geometric constraints.

Deformable-model based segmentation

Deformable models have also been employed for brain segmentation in the literature. Angelini et al. [43] introduced a multi-phase level set framework for automated segmentation of brain MRIs. The segmentation of the brain tissues (WM, GM and CSF) was solely based on homogeneity (average grey level) measures. To avoid the need for any prior information and to speed up numerical calculation, a random seed for initialization of the deformable boundaries was used. Colliot et al. [44] proposed a deformable-model based approach that used spatial constraints, represented as fuzzy subsets of the 3D image space, as an external force to control the boundary evolution. To avoid manual selection of the model parameters, a training step was required to estimate the spatial constraints parameters. Miri et al. [45] introduced a topology-preserving deformable model framework for the segmentation of brain MRIs. They employed photometric constraints to guide the deformable model deformations to iteratively reclassify the points located at the evolving boundaries. A deformable model approach for the segmentation of brain regions from MR images was proposed by Liu et al. [46]. The deformable contour was implicitly represented by a set of Wendland’s radial basis functions (RBFs) and was evolved by iterative updates of the locations of the RBFs. Huang et al. [47] introduced an automated, hybrid deformable model framework that integrated both image edge geometry and voxel statistics features to regularize the convergence of the deformable contour. Del Fresno et al. [48] described a hybrid method that combined region growing and deformable models for segmentation of different structures in head MRI and Computed Tomography CT scans. Their approach used a Region-Growing (RG) algorithm to compute an approximation of the objects. This was followed by generating closed and oriented surface meshes to enclose the region of interest. The deformable model method geometry was constructed using the RG-list of boundary voxels generating a hole-free surface mesh. To better detect the structures of interest, the user could select few seeds for RG initial segmentation. Wang et al. [49] proposed a multi-phase level set framework to segment brain MR images with intensity inhomogeneity. They modeled the local image intensities using Gaussian distributions with different means and variances. Then, a variational approach minimized an energy function to compute the means and variances that would guide the contour evolution towards the target boundaries. Bourouis et al. [50] developed a level set framework for segmenting brain tissues. Their framework employed an image registration step and a classification step for the initialization of the deformable boundary. The boundary evolution was controlled by a speed function that accounted for both boundary- and region-based properties. Ciofolo et al. [51] developed an automated framework based on level sets for simultaneous segmentation of multiple structures from brain MRIs. The evolution of each level set was driven by a fuzzy decision system that combined three factors: intensity distribution of the 3D MR volume, the relative position of the evolving contours, and a priori knowledge provided by an anatomical atlas. Wang et al. [52] proposed a multi-layer background subtraction technique with a seed region growing approach which used local texture features represented by local binary patterns (LBP) and photometric invariant color measurements in RGB color space for brain segmentation. Zhao et al. [53] segmented brain tissues using a method adapted from Chan and Vese model, named automatic threshold level set without edges. Thresholds were obtained by fuzzy c-mean algorithm.

Segmentation using deformable models was also exploited in infant brains, [16, 32, 33, 42].

The main advantage of deformable-model based segmentation techniques is the ability to segment connected (non-scattered) objects more accurately than the other segmentation methods. However, the accuracy of this method is based on the accurate design of the guiding forces (statistical, geometric, etc.) as well as the model initialization. A summary for all cited related work is provided in Table 1.

Table 1

Summary of brain segmentation related work.

Group	Category	Methodolody
Probabilistic	Ng et al. [7]	Unsupervised K means clustering.
	Xue et al. [5]	Parametric Gaussian DensityEstimation with Expectation Maximization.
	Mayer et al. [8]	Adaptive mean shiftclassifier and intensity-based mode clustering.
	Fang et al. [9]	Graph cuts algorithm.
	Ortiz et al. [10]	Growing hierarchical self-organisingmap (GHSOM) classifier and probability clustering.
	Li et al [11]	K-means clustering and iterated conditional modes.
	Janney et al [12]	Clustering joint spatial-intensity feature space.
	Weber et al [13]	K-means clustering and Expectation Maximization algorithm.
	Mahmood et al [14]	Bayesian framework and fuzzy c-means.
	Anbeek et al [15]	Probability maps and KNN classifier.
	Wang et al [16]	Probability maps and random forest classifiers.
	Zhang et al [17]	Convolutional neural networks.
	Moeskops et al [18]	Supervised voxel classification.
Atlas-based	Ashburner et al [19]	Registration.
	Pohl et al [20]	Simultaneous registration and segmentation.
	Han et al [21]	Intensity re-normalization for prior shape.
	Artaechevarria et al [22]	Local weighting scheme from fusion.
	Sabuncu et al [23]	Label fusion.
	Morin et al [24]	Labeling propagation steps.
	Lotjonen et al [27]	Atlas with majority voting and intensity models.
	Van et al [28]	Spatial and appearance models.
	Ledig et al [30]	Markov Random Field and prior shape.
	Srhoj et al [31]	Multiple atlases with supervised voxel classification.
	Makropoulos et al [35]	Structural hierarchy and anatomical constraints.
	Cherel et al [34]	Subject-specific atlas with single-atlas expectation maximization.
Hybrid	Song et al. [39]	Probabilistic neural network.
	Patenaude et al [40]	Active shape model within a Bayesian framework.
	Wang et al [42]	Sparse representation of tissue distribution.
	Serag et al [41]	Sliding window and random forest classifier.
Deformable	Angelini et al [43]	Multi-phase level set framework.
	Colliot et al [44]	model with spatial constraints.
	Miri et al [45]	Topology-preserving model with photometric constraints.
	Liu et al [46]	Radial Basis Functions.
	Albert et al [47]	Edge geometry and voxel statistics.
	Del et al [48]	model initialized with region growing.
	Wang et al [49]	Multi-phase level set framework.
	Bourouis et al [50]	Level sets initialized by registration.
	Ciofolo et al [51]	Level sets driven by a fuzzy decision system.
	Wang et al [52]	Region growing with local texture features.
	Zhao et al [53]	Automatic threshold level sets without edges.
DL	Brebisson et al [54]	3D and orthogonal 2D intensity patches.
	Zhang et al [17]	Multi-modality information.
	Chen et al [55]	low-level image appearance features,implicit shape information, and high-level context.

Deep-learning based segmentation

De Brebisson et al. [54] proposed a deep artificial neural network for brain image segmentation from MR scans that assigns each voxel to its corresponding anatomical region. The information to the network is obtained at different scales around the voxel of interest: 3D and orthogonal 2D intensity patches capture a local spatial context. Also, large compressed 2D orthogonal patches and distances to the regional centroids are used so that they would enforce global spatial consistency. Zhang et al. [17] proposed a method that is based on deep convolutional neural networks for multi-modality isointense infant brain image segmentation. Multi-modality information from T1, T2, and fractional anisotropy (FA) images were exploited as inputs to generate the segmentation maps as outputs. The multiple intermediate layers applied convolution, pooling, normalization, and other operations to capture the highly nonlinear mappings between inputs and outputs. Chen et al. [55] proposed a deep voxelwise residual network, referred as VoxResNet, which handles volumetric data. It integrated the low-level image appearance features, implicit shape information, and high-level context together for improving the volumetric segmentation performance.

In summary, brain segmentation work found in the literature suffer from many drawbacks as mentioned above. Moreover, infant brain segmentation techniques depend on either multiple modalities which lengthens the processing time, or on longitudinal studies which are not always available for research purposes. Also the reduction in contrast between CWM and other structures hardens the issue of preserving its edges by most of the available techniques.

To overcome the aforementioned limitations, this paper proposes a novel technique for brain segmentation from MR images. Adaptive probabilistic shape models for the shape and first-order visual appearance of MRI data are employed to initialize the segmentation. This is then combined with a novel higher-order Markov Gibbs random field (MGRF) spatial interaction model (up to fourth order) with analytic estimation of potentials. This joint model guarantees increasing the segmentation accuracy by accounting for the large-scale inhomogeneities and noise in infant brain MRI data. Also, the analytic estimation of potentials generalizes the proposed model to different MRI subjects, unlike using empirical values in most of the work present in literature which would require manual setting for each subject. The strength of the proposed algorithm lies in the fact that it neither depends on multiple modalities for acquiring images nor on longitudinal studies.

Methods

Some brain tissues such as the brain stem and cerebellum are similar to the WM and GM intensities. Therefore, the segmentation technique, Fig 2, makes use of the adaptive shape prior that is co-guided by both a first-order visual appearance descriptor using the estimated LCDG models for each class as well as the 3D spatial relationships between the region labels to segment each label. This forms a 3D joint model that integrates shape, intensity, and spatial information.

Fig 2

Proposed segmentation framework.

Preprocessing

Before segmentation takes place, bias correction is applied to the brain volumes using the nonparametric approach proposed by Tustison at al. [56] to account for nonuniform intensities. This is followed by applying a 3D generalized Gauss-Markov random field (GGMRF) model [57] that reduces noise effects and removes inconsistencies of the scans.

The skull is then removed before segmentation takes place. This work uses the automated approach proposed by our group in [58], which refines the skull stripping results from the brain extraction tool (BET). This is achieved using an additional processing step that is based on the geometric features of the brain. Since the non-brain tissues are brighter than brain tissue, this step exploits the visual appearance features of the MR brain data. Namely, an evolving iso surface-based approach is proposed to remove the non-brain tissues, which is guided by the MR data visual appearance features. First, a set of nested iso-surfaces are generated by the fast marching level sets (FMLS), using the calculated distance map of the extracted brain from the BET step. Then, classification of voxels as brain or non-brain is conducted. Results on both infants and adults using this approach showed the capability of the proposed approach, outperforming four widely used brain extraction tools: BET, BET2, brain surface extractor (BSE), and infant brain extraction and analysis toolbox (iBEAT).

Joint MGRF model of MR brain images

Let R = {(x, y, z) : 0 ≤ x ≤ X − 1, 0 ≤ y ≤ Y − 1, 0 ≤ z ≤ Z − 1}; Q = {0, 1, …, Q − 1}; and L = {0, …, L} denote a finite 3D arithmetic lattice of the size of XYZ supporting grayscale images and their region (segmentation) maps, a finite set of Q integer gray values, and a set of region labels L, respectively. Let g = {g : (x, y, z) ∈ R; g ∈ Q} and m = {m : (x, y, z) ∈ R; m ∈ L} be a grayscale image taking values from Q, i.e., g : R → Q, and a region map taking values from L, i.e., m : R → L, respectively. An input brain image, g, co-aligned to the training data base, and its map, m, are described with a joint probability model: P(g, m) = P(g|m)P(m), which combines a conditional distribution of the images given the map P(g|m), and an unconditional probability distribution of maps P(m) = Psp(m)P(m). Here, Psp(m) denotes a weighted shape prior, and P(m) is a Gibbs probability distribution with potentials V, which specifies a MGRF model of spatially homogeneous maps m. Details of the model’s components are outlined below.

Adaptive shape model Psp(m)

To start the segmentation process, a database is created, where expected shapes of each brain label are constrained with an adaptive probabilistic shape prior. To create the atlases, a training set of images, collected for different subjects (not included as test subjects), are co-aligned by 3D affine transformations with 12 degrees of freedom (3 for the 3D translation, 3 for the 3D rotation, 3 for the 3D scaling, and 3 for the 3D shearing) in a way that maximizes their Mutual Information (MI) [59]. An atlas of 10 subjects containing the three labels to be segmented (WM, GM, and CSF) was constructed for each age group (infants, children, and adults) as described. For each input MR data to be segmented, the shape prior is constructed by an adaptive process guided by the visual appearance features of the input MRI data [60-62]. The shape prior is a spatially variant independent random field of region labels for the co-aligned data: where psp:(l) is the voxel-wise empirical probabilities for each brain label l ∈ L. First, the normalized cross correlation similarity coefficient is used to select the subject from the shape database that has the best match with the input subject (i.e., highest similarity). The selected subject is then used as a reference prototype to co-align the input subject using the 3D affine transformation described above. In order to estimate the shape prior probabilities for each voxel in the test subject, the steps summarized in Algorithm 1 are followed. Fig 3 show the calculated probabilistic maps for each structure using the proposed adaptive shape prior.

Fig 3

The calculated shape probability for the CSF(a), GM(b), and WM(c).

Algorithm 1 Steps of the shape prior segmentation.

1. Align the test subject with the shape database to get the 3D affine transformation matrix T.

2. For each slice i, i = 1 to n

I. For each voxel v in slice i

(a) Transform v to the atlas domain using the transformation matrix T.

(b) Initialize a 3D cube, C, of size N1 × N2 × N3 centered around the mapped voxel (vmapped).

(c) Search C for voxels with corresponding grey level in all training sets with equalized intensities that fall within a predefined tolerance ±τ in w.

(d) If no voxels are found using Step (c), increase size of C and repeat step (c) until correspondences are found or the maximum size allowed for C is reached.

(e) Calculate the shape probability for each structure at location r based on the found voxels and their labels.

End for

End for

3. Return the constructed 4D shape probabilities.

First-order intensity model P(g|m)

The first-order visual appearance of each brain label is modeled by separating a mixed distribution of voxel intensities of the brain scans into individual components associated with the dominant modes of the mixture. The latter is precisely approximated with a Linear Combinations of Discrete Gaussians (LCDG) [63] with positive and negative components, which is based on a modified version of the classical Expectation-Maximization (EM) algorithm.

Let Ψ = (ψ(q|θ) : q ∈ Q) denote a discrete Gaussian (DG) with parameters θ = (μ, σ), integrating a continuous 1D Gaussian density with mean μ and variance σ2 over successive gray level intervals. The LCDG with four dominant positive DGs and Cp ≥ 4 positive and Cn ≥ 0 negative subordinate DGs is [63]: where all the weights w = [wp:, wn:] are non-negative and meet an obvious constraint . All LCDG parameters, including Cp and Cn, are estimated from the mixed empirical distribution to be modeled using the modified EM algorithm [63]. Fig 4 shows an example of the empirical density for different brain tissues using the LCDG model for an infant subject, Fig 4(a), and an adult one Fig 4(b).

Fig 4

Normalized empirical density using the LCDG model for an infant subject (a), and an adult one (b).

Note that dashed = empirical, red = CSF component, green = GM component, blue = WM component.

MGRF model with second- and higher–order cliques PV(m)

In addition to the first-order visual appearance model, the spatial interactions between the brain voxels are also taken into account. Using spatial models that are only second-order-clique based, (e.g., [64]), will not enable accounting for the spatial inhomogeneity of brain MR images, especially for infants. Therefore, in this paper we propose a higher-order Markov-Gibbs Random Field (MGRF) spatial interaction model that adds the families of the triple and quad cliques to the pairwise cliques (Fig 5(b) and 5(c)), along with analytical estimation of the potentials. The proposed approach accounts for the spatial inhomogeneity of the brain scans, especially for those of infants, thus, reducing noise effects and increasing segmentation accuracy. Details of the proposed higher-order MGRF model are described below.

Fig 5

Samples of the second- (a), third- (b), and fourth-order (c) cliques for the 26-neighborhood (graph cliques are shown in different colors for visualization purpose).

Let C denote a family of s-order cliques of an interaction graph with nodes in the 3D lattice sites (x, y, z) and edges connecting the interacting, or interdependent, sites (see Fig 5). To account for the scan inhomogeneities, especially with infant MRI, the label interactions are modeled by a spatially homogeneous MGRF with up to fourth-order interactions over the nearest 26-neighborhoods of voxels: where A clique families describe the geometry of the graph interactions, V = [V : {0, …, L} → (−∞, ∞) : a = 1, …, A] is a collection of Gibbs potential functions V for the families C, and the partition function Z normalizes the probabilities over the parent population of all maps. An initial region map m, obtained by the voxel-wise classification, allows for analytically approximating the maximum likelihood estimates of the potentials and computing the voxel-wise probabilities of the region labels. For symmetry sake, only equality or inequality of the labels in clique c is taken into account. The second- third- and forth–order potentials are given by Eqs (4), (5) and (6), respectively: where , and F(m°) = [ρF(μ1, …, μ|m°) : (μ1, …, μ) ∈ {0, …, L}; a = 1, …, A] is the collection of scaled relative frequencies of co-occurrences of configurations (μ1, …, μ) of the labels in the cliques of each family C over a given training map m°. where where and where m is the region map label at the voxel p = (x, y, z). The proposed analytical approximation of the Gibbs potentials from a given map m extends earlier second-order MGRFs (e.g., [64]) to the higher-order models. The complete proof for the higher-order MGRF model is provided in the supplement.

Finally, the region map m is improved using Iterative Conditional Mode (ICM) algorithm [65] that maximizes the probabilities of the 3D joint model. The complete steps of our segmentation framework are summarized in Algorithm 2. Also Fig 2 illustrates the whole proposed framework.

Experimental results

The proposed segmentation framework was tested on different databases at different ages to show its generality and robustness. 42 subjects from the Kennedy Krieger Institute (KKI)(8–12.8 years), 20 subjects from the university of California (UCLA) (8.4–17.9 years), and 20 subjects from the NYU Langone Medical Center (6.5–39.1 years) [66] were used to validate the segmentation approach. Moreover, 20 infants from the NDAR/IBIS database (aged 6 months) [67] were segmented using the proposed framework. For all subjects used throughout this study, there was no available information that would reveal the identity of the individual participants during or after data collection.

The IBIS database comprises T1-weighted images, and were acquired on a 3 tesla scanner with TR = 2400 millisecond (ms), TE = 3.16 ms TI = 1200 ms, and flip angle = 8. 160 sagittal slices were acquired at 1 millimeter (mm) thickness, with each slice being 224 × 256 pixels with 1 mm resolution.

The UCLA database includes participants between 8.4 and 17.9 years of age. T1-weighted images were acquired using MPRAGE with TR = 2300 ms, TE = 2.84 ms, and flip angle = 9. Sagittal slices were acquired at 1.2 mm thickness. The pixels of each 256 × 256 slice were 1 mm each side.

Algorithm 2 Steps of the Proposed segmentation framework.

○

(a) Use the automated approach in [58] to remove the skull from the MR images.

(b) Construct the shape database for each age group through a co-alignment of the biased-corrected training volumes (both grey scale and their ground truth).

○

(a) For each Atlas:

i. Estimate the adaptive shape prior probability (P(m)) using Algorithm 1

ii. Approximate P(g) using an LCDG with four dominant modes.

iii. Form region map m using marginal estimated density and prior shape.

iv. Find the Gibbs potentials for the MGRF model from the initial map m.

v. Improve m using the iterative conditional mode (ICM) algorithm [65].

(b) Apply majority voting to fuse the segmentation results of the three atlases.

Finally, the NYU data includes participants between 6.5 and 39.1 years of age. T1-weighted images were acquired on a 3 tesla Allegra with TR = 2530 ms, TE = 3.25 ms, and flip angle = 7. Sagittal slices were acquired at 1.33 mm thickness. The pixels of each 256 × 256 slice were 1.3 mm, and 1 mm. A summary of all the databases used in this work is provided in Table 2.

Table 2

Summary of databases used to validate the proposed method.

Database	Subjects	Age	Scan Parameters
IBIS	20	6 months	B = 3 T, TR = 2400 ms, TE = 3.16 msflip angle = 8°, 1 mm thickness.
KKI	42	8–13 years	B = 3 T, TR = 8 ms, TE = 3.7 ms,flip angle = 8°, 1.33 mm thickness.
UCLA	20	8.4–17.9 years	TR = 2300 ms, TE = 2.84 ms,flip angle = 9°, 1.2 mm thickness.
NYU	20	6.5–39.1 years	B = 3 T, TR = 2530 ms, TE = 3.25 ms,flip angle = 7°, 1.33 mm thickness.

The proposed segmentation approach was evaluated on all subjects above using their manually segmented ground truth created by an MR expert. Special care was given to the infant subjects, since these were the hardest to delineate due to the unmyelinated nature of WM at this early age. Results in the upcoming pages show that the joint model combining the intensity, spatial, and shape information is in general a better performer than having only one or two of three models. The intensity information alone would often fail to differentiate between different tissue types of the infant MRI scans. This is enhanced after using higher-order MGRF, where edges of different tissue types are better retained. Comparing the visual segmentation results of using the intensity model alone for the IBIS database, with those for other databases, it can be inferred that major parts of the CWM bundles were missing for infant scans, which was not the case for the other databases. This result is expected, since the contrast is extremely poor between the different tissue types in infant MRI scans, whereas it gets better with older ages represented by the other databases. Results also show the merit of using the higher-order MGRF model with infant scans, which was compared against the second-order model and achieved better results. The performance of the two MGRF models with other age groups didn’t show a notable difference, which is expected since those scans are more homogeneous and less noisy than infant ones.

The segmentation performance was evaluated using 3 metrics: (i) the Dice similarity coefficient (DSC) [71], (ii) the 95-percentile modified Hausdorff distance (MHD) [72], and (iii) the absolute brain volume difference (ABVD), by comparing to the ground truth segmentation. Tables 3, 4, and 5 summarize the accuracy results obtained using the three metrics for the WM, GM, and CSF of different databases. As these tables show, the DSC metric was used to evaluate the proposed approach using both the second- and the higher-order MGRF. The reported accuracies show the advantages of using the higher-order MGRF, especially with the IBIS database, Table 3. Accuracies of 89.5% and 90.9% were achieved using the second-order model for WM and CGM segmentation respectively, whereas they enhanced to 94.7% and 95.2% with the higher-order model. The accuracies were also increased using the higher-order model with the other databases (Tables 4 and 5), yet not with the same rate it did with the IBIS database. This is acknowledged to the fact that the infant MR volumes suffer from more noise and image inhomogeneities which could be accounted for using the higher-order MGRF model. Segmentation results samples from each database are shown in Figs 6, 7, 8, 9 and 10, where the three extracted labels (WM, GM, CSF) using the proposed method and iBEAT [68], FSL [69], and FreeSurfer [70] are displayed along with the ground truth segmentation.

Table 3

Accuracy of our segmentation approach using Dice Similarity Coefficient (DSC)(%), the modified Hausdorff Distance (MHD)(mm), and Absolute Brain Volume Difference (ABVD) (%) for the WM, GM, and CSF of the IBIS database.

Metrics are represented as Mean±Standard Deviation. Results for the proposed approach are shown using both the second- and higher-order MGRF model. Age of this group is 6 months.

		Segmentation Method
Struct.	Metric	Proposed Method		iBEAT	FSL	FS
		2nd-order MGRF	Higher-order MGRF
WM	DSC	89.5±2.43	94.7±1.53	73.3±1.27	80.4±1.57	85.6±2.3
	MHD	10.5±0.7	7.3±1.23	18.27±1.53	13.6±1.28	11.2±1.0
	ABVD	6.5±2.54	3.17±1.73	37.94±0.61	15.8±0.6	10.2±3.7
GM	DSC	90.9±1.56	95.2±0.13	81.6±3.5	89.5±0.65	86.5±1.3
	MHD	6.7±1.5	3.5±0.24	23.3±0.52	15.7±0.86	10.6±3.7
	ABVD	5.23±3.1	1.62±1.24	34.46±0.18	14.5±1.17	12.1±2.5
CSF	DSC	93.2±1.2	94.58±0.44	79.65±1.38	82.4±1.5	88.5±1.4
	MHD	5.6±2	4.35±1.1	27.23±1.43	9.7±0.3	7.8±1
	ABVD	2.5±0.73	1.9±0.11	21.07±0.98	7.5±2	5.3±1.2

Table 4

Accuracy of our segmentation approach using Dice Similarity Coefficient (DSC)(%), the modified Hausdorff Distance (MHD)(mm), and Absolute Brain Volume Difference (ABVD) (%) for the WM, GM, and CSF of the KKI database.

Metrics are represented as Mean±Standard Deviation. Results for the proposed approach are shown using both the second- and higher-order MGRF model. Age range of this group is 8–13 years.

		Segmentation Method
Struct.	Metric	Proposed Method		FSL	FS
		2nd-order MGRF	Higher-order MGRF
WM	DSC	94.1±0.8	95.8±1.5	88.4±2.5	91.5±2.3
	MHD	5.8±0.9	5.2±1.5	9±1.2	7.5±1.5
	ABVD	3.2±0.8	2.5±1.6	9.8±1.6	9.1±1.8
GM	DSC	94.9±0.77	96.7±1.2	92.3±1.85	92.2±1.2
	MHD	4±1.5	3.1±1.7	11.8±1.2	8.2±1.25
	ABVD	1.8±0.4	1.12±1.1	8.3±2.7	5.2±1.4
CSF	DSC	95.5±1.25	96.5±1	92.1±2.7	93.2±3.7
	MHD	5±1.3	3.7±2	7.25±2.1	7±1.5
	ABVD	3±1.25	2.5±1	9±2.0	8.25±2.1

Table 5

Accuracy of our segmentation approach using Dice Similarity Coefficient (DSC)(%), the modified Hausdorff Distance (MHD)(mm), and Absolute Brain Volume Difference (ABVD) (%) for the WM, GM, and CSF of the NYU and UCLA databases.

Metrics are represented as Mean±Standard Deviation. Results for the proposed approach are shown using both the second- and higher-order MGRF model. Age range of this group is 6.5–39.1 years.

		Segmentation Method
Struct.	Metric	Proposed Method		FSL	FS
		2nd-order MGRF	Higher-order MGRF
WM	DSC	95.7±1.1	96.1±1.5	92.4±1.5	93.5±0.2
	MHD	2.7±0.5	2.3±1.2	11.8±1.8	7.5±2.0
	ABVD	1.9±0.72	1.5±1.7	10.7±1.6	7.8±2
GM	DSC	97.1±0.6	97.8±0.13	93.2±1.3	94.4±1.8
	MHD	1.3±0.44	0.9±0.2	9.1±1.6	5±0.8
	ABVD	1.88±0.76	1.2±1.4	10.6±1.73	7±2.3
CSF	DSC	96.25±1	97.5±0.7	91.5±2	93±1.5
	MHD	2±0.65	1.5±0.8	8.5±0.6	7±1.2
	ABVD	1.5±0.2	1±0.3	9.5±3	8.2±0.9

Fig 6

Segmentation results; for CSF (a), GM (b), and WM (c); projected onto axial plane for a six-month-old subject from the IBIS database for infants: Segmentation using our proposed method (first row); using the iBEAT method (second row); and Ground truth (third row).

Fig 7

Segmentation results; for CSF (a), GM (b), and WM (c); projected onto axial plane for a additional six-month-old subject from the IBIS database for infants: Segmentation using our proposed method (first row); using the iBEAT method (second row); and Ground truth (third row).

Fig 8

Segmentation results; for CSF (a), GM (b), and WM (c); projected onto axial plane for a nine-year-old subject from the KKI database: Segmentation using our proposed method (first row); using the FSL method (second row); using the FreeSurfer method (third row); and Ground truth (fourth row).

Fig 9

Segmentation results; for CSF (a), GM (b), and WM (c); projected onto axial plane for a 16-year-old subject from the UCLA database: Segmentation using our proposed method (first row); using the FSL method (second row); using the FreeSurfer method (third row); and Ground truth (fourth row).

Fig 10

Segmentation results; for CSF (a), GM (b), and WM (c); projected onto axial plane for a sample from the NYU database: Segmentation using our proposed method (first row); using the FSL method (second row); using the FreeSurfer method (third row); and Ground truth (fourth row).

The performance of the proposed segmentation approach is highlighted by comparing it against the software package (iBEAT) [68], that performs bias correction followed by brain segmentation. Moreover, segmentation was done using the FSL package [69], and also using FreeSurfer software [70]. Segmentation accuracies for the iBEAT, FSL, and FreeSurfer results are also summarized in Tables 3, 4, and 5, where accuracies are reported for CWM, GM, and CSF. These results emphasize the efficiency of the proposed approach that is required for possible subsequent processes such as shape analysis. Table 6 breaks down the timing expended by each approach to perform segmentation for a single subject. The proposed approach turned out to be the fastest, which adds up to its advantages.

Table 6

Summary of the time required by the proposed approach and other approaches for segmenting a brain subject.

Approach	Required execution time
iBEAT [68]	40 minutes
FSL [69]	22 minutes
FreeSurfer [70]	18 hours
Proposed approach	92 seconds for skull stripping,370 seconds for alignment,65 seconds for segmentation,Total is about 9 minutes.

It is worth mentioning that in some cases where the anatomical changes in the brain structures are irregular/rapid, there will be a lack of support from the 3D relative neighbors which will affect the 3D spatial dependencies between each voxels and its neighbors. This smoothness assumption by the MGRF model may results in an underestimation of thin white matter strands. Fig 6 shows an example of such effect for an axial cross section from one subject from the IBIS database, while Fig 7 shows another axial cross section of enhanced results.

Conclusions and future work

This paper proposed a 3D automated approach for brain segmentation from MR images in subjects spanning different ages. The segmentation method integrates intensity, shape, and spatial information in a hybrid model. The MGRF spatial model, the higher-order one in particular, accounts for scan inhomogeneities and noise that are drastically affecting infant scans. The novelty of the proposed algorithm lies at using the adaptive shape model along with the higher-order MGRF model. The work has been tested and validated on 102 MR scans and compared against state-of-the-art approaches. The metrics used to evaluate the segmentation showed that the proposed approach is a better performer in terms of accuracy and time.