Segmentation of Infant Brain Using Nonnegative Matrix Factorization

Abstract

This study develops an atlas-based automated framework for segmenting infants’ brains from magnetic resonance imaging (MRI). For the accurate segmentation of different structures of an infant’s brain at the isointense age (6–12 months), our framework integrates features of diffusion tensor imaging (DTI) (e.g., the fractional anisotropy (FA)). A brain diffusion tensor (DT) image and its region map are considered samples of a Markov–Gibbs random field (MGRF) that jointly models visual appearance, shape, and spatial homogeneity of a goal structure. The visual appearance is modeled with an empirical distribution of the probability of the DTI features, fused by their nonnegative matrix factorization (NMF) and allocation to data clusters. Projecting an initial high-dimensional feature space onto a low-dimensional space of the significant fused features with the NMF allows for better separation of the goal structure and its background. The cluster centers in the latter space are determined at the training stage by the K-means clustering. In order to adapt to large infant brain inhomogeneities and segment the brain images more accurately, appearance descriptors of both the first-order and second-order are taken into account in the fused NMF feature space. Additionally, a second-order MGRF model is used to describe the appearance based on the voxel intensities and their pairwise spatial dependencies. An adaptive shape prior that is spatially variant is constructed from a training set of co-aligned images, forming an atlas database. Moreover, the spatial homogeneity of the shape is described with a spatially uniform 3D MGRF of the second-order for region labels. In vivo experiments on nine infant datasets showed promising results in terms of the accuracy, which was computed using three metrics: the 95-percentile modified Hausdorff distance (MHD), the Dice similarity coefficient (DSC), and the absolute volume difference (AVD). Both the quantitative and visual assessments confirm that integrating the proposed NMF-fused DTI feature and intensity MGRF models of visual appearance, the adaptive shape prior, and the shape homogeneity MGRF model is promising in segmenting the infant brain DTI.

1. Introduction

The development of computer-assisted diagnostic (CAD) systems for brain disorders requires having reliable segmentation for the various anatomical structures of the brain (e.g., the gray matter (GM), the white matter (WM), and the cerebrospinal fluid (CSF)) [1]. Regarding infants, for instance, analyzing brain magnetic resonance imaging (MRI) scans in order to treat injuries and disorders caused by brain prematurity depends significantly on brain segmentation [2]. Currently, manual segmentation done by experts is the gold standard. However, it is time consuming and can be subjective. To handle these issues, automated and semi-automated segmentation approaches need to be developed. Brain segmentation from MRI can be very challenging due to the similar appearance of neighboring structures inside the brain, which make it difficult to identify boundaries among them, especially in the presence of inhomogeneity, image noise, and partial volume effect. In addition, in diffusion-weighted imaging (DWI), the effects of the distortions arising from the motion of the patients can be significantly emphasized by the diffusion sensitizing gradient [3]. As opposed to segmenting adult brains, infant brain segmentation is harder, as it can be carried out only based on the image intensity, which is the main focus of this study. The tissues of infant brains are not mature enough to give sufficient contrast between the GM and WM [4] along with the presence of much noise [5]. Moreover, distortions from eddy currents and motion of unseated infants add more difficulty to the segmentation problem. Over the past decades, researchers have been working hard to develop techniques for brain segmentation from MRI. Primarily, state-of-the art segmentation techniques deal with adult brains; however, a few of them can fit infant brains. Those segmentation methods can be classified under three categories: statistical [6,7,8,9,10,11,12,13], atlas-based [2,14,15,16,17,18,19,20,21,22], and deformable model-based methods [23,24,25,26,27,28,29,30,31,32]. Statistical methods are the fastest compared to other techniques; however, they cannot handle all real distributions of data as they rely on pre-defined statistics. While atlas-based methods are more accurate than statistical methods, they require heavy computations and are very dependent on the selection of atlases. Deformable-model-based methods outperform other techniques in segmenting connected objects, though these techniques depend on the design of the statistical and geometric guiding forces as well as the model initialization. In this study, we focused on addressing the limitations of the current segmentation techniques that deal with magnetic resonance (MR) infant brains, with a special focus on those which incorporate diffusion tensor imaging (DTI) measurements in their framework. To overcome these limitations, we propose in this study an adaptive, atlas-based technique for segmenting different structures of infant brains automatically from DTI scans (see Figure 1). Details of the proposed segmentation framework will be presented in the following section.

2. Materials and Methods

The segmentation technique we propose in this study includes three main parts: (i) preprocessing, (ii) extracting the brain, and (iii) segmentation using a joint MGRF model. Each part will be discussed in detail in the following subsections below after introducing the used data, and the following notations will be used:

• $\mathbf{R}$$=\left\{\right(x,y,z):0\le x\le X-1,0\le y\le Y-1,0\le z\le Z-1\}$—a three-dimensional raster supporting the MRI data and their associated region maps. For each voxel $s=(x,y,z)$, not located on one of the raster’s edges, we denote its 26-neighborhood by ${\nu }_s$.
• $\mathbf{g}$$:\mathbf{R}\to \mathbf{Q}$—a grayscale image with voxel values in $\mathbf{Q}=\{0,1,\dots ,Q-1\}$.
• $\mathbf{m}:\mathbf{R}\to \mathbf{L}$—a region map (i.e., segmentation) with segment labels in $\mathbf{L}=\{0,\dots ,L\}$.
• $\mathbf{A}=\{a_{i,n}:i=1,\dots ,I;n=1,\dots ,XYZ;a_{i,n}\in {\mathbb{R}}^{+}\}$—a positive matrix comprising I raw image features for each voxel.
• $\mathbf{W}=\{w_{i,j}:i=1,\dots ,I;j=1,\dots ,J;w_{i,j}\in {\mathbb{R}}^{+}\}$—the matrix of basis features in the reduced feature space.
• $\mathbf{H}=\{h_{j,n}:j=1,\dots ,J,n=1,\dots ,XYZ;h_{i,j}\in {\mathbb{R}}^{+}\}$—the matrix of features for each voxel in the reduced feature space.

2.1. Cohort and Imaging Parameters

This study employed MRI and related metadata downloaded from the NIMH Data Archive (NDA) in conformance with its data usage agreement and with approval from the authors’ IRB. Imaging had been originally acquired, and archived on NDA, as part of the Infant Brain Imaging Study (IBIS) [33]. The study subjects are six-month-old infants with high risk of developing autism. Those infants were monitored for further follow-up assessments at one year, and two years, and the final assessment was made at two years of age. Those high-risk subjects were split into two groups: ASD negative (below the ASD cutoff) and ASD-positive (above the cutoff). Dw-MRI brain scans were obtained from 54 individuals using 3-T Siemens TIM Trio scanners (Siemens Medical Solutions, Malvern, PA, USA) with acquisition parameters of: 190 mm FoV, between 75 and 81 slices with 2 mm thickness, 2 mm × 2 mm pixel spacing, TR = 12,800 to 13,300 ms, TE = 102 ms, gradient b-values of 0 and 1000 s/mm${}^2$, and 25 gradient directions.

2.2. Preprocessing

Throughout this study, the DWI quality is controlled through the DTIprep software [34], which is responsible for eliminating different image artifacts, and rectifying the eddy current and motion distortions. Using the detailed report produced by DTIprep and the results obtained, the images with high artifacts were removed. Furthermore, the resulting datasets with a low number of gradient directions that could result in a low SNR were omitted from further processing. The complete DTIprep software details can be found in [34].

Before any DTI features can be extracted, tensor model estimation from the preprocessed DWI datasets is performed. In this study, we selected the weighted linear least square (WLLS) approach because of its ability to provide accurate estimates compared with the traditional linear least square (LLS) methods. The 3D Slicer software [35] was utilized for the DTI to DWI estimation using the WLLS approach. The bad tensors caused by the remainder of acquisition artifacts and noise were handled by shifting the negative eigenvalues.

2.3. Extracting the Brain

The extraction of the brain procedure (AKA skull stripping) includes separating all non-brain regions from brain MR images. Such an approach uses a visual appearance model of the first-order to direct the edge-preserving model for restoring the images, while compensating for any inhomogeneity in the intensities by energy minimization-based optimization.

2.3.1. LCDG Visual Appearance Model

The brain tissue and the non-brain tissue marginal probability distributions can be represented precisely by approximating the empirical gray level distribution for a certain $\mathbf{g}$ by a linear combination of discrete Gaussians (LCDG), whose components can take positive and negative values [36], as the LCDG is capabable of accurately retrieving the transitions between the brain and non-brain tissues, as opposed to using a conventional Gaussian mixture of positive components only. See Appendix A.1 for more implementation details. The final estimation of each class LCDG density is illustrated in Figure 2e.

At the end of this step, a discriminant threshold $\tau$ is calculated in a way that guarantees the maximum separation between the signals of the brain voxels and the non-brain voxels, as illustrated in Figure 2f. This threshold will be used in the next steps to enhance the process of classifying image voxels into either brain or non-brain.

2.3.2. GGMRF Inhomogeneity Reduction Model

The accuracy of the current techniques of infant brain extraction from DWI are significantly hindered by the intensity inhomogeneity. This could be handled by compensating for the intensity inhomogeneity with low frequency. In the proposed brain extraction method, the generalized Gauss–Markov random field (GGMRF) model is used to reduce this inhomogeneity [37], which is applied after DWI data are preprocessed using DTIprep software. In this step, the residual image inconsistencies are removed (smoothed) by compensating for the 3D pair-wise spatial homogeneous interactions along the DWI gray levels. Particularly, the values of the gray levels $q\in \mathbf{Q}$ are treated as samples of a 3D GGMRF model [37] with measurements of the voxel 26-connected neighborhood (neighbors to the pixel from all directions, edges, and corners). The continuity of the q values for each of the brain DWI scans is modified by applying the gradient descent algorithm [38] to search for the closest minimum of the following equation:

$${\hat{q}}_s=arg\underset{{\tilde{q}}_s}{min}\left[{\left|q_s-{\tilde{q}}_s\right|}^{\alpha }+{\rho }^{\alpha }{\lambda }^{\beta }\sum _{r\in {\nu }_s}{\eta }_{s,r}{\left|{\tilde{q}}_s-q_r\right|}^{\beta }\right],$$

(1)

where the original values of the gray levels and their corresponding expectations are denoted respectively by $q_s$ and ${\tilde{q}}_s$ at the observed location $s=(x,y,z)$. The GGMRF potential and the 26-neighborhood voxel set are denoted respectively by ${\eta }_{s,r}$ and ${\nu }_s$. In addition, $\lambda$ and $\rho$ are scaling factors. The level smoothness is controlled by the parameter $\beta \in [1.01,2.0]$. For example, having $\beta =2$ yields to smoothness, while $\beta =1.01$ yields to relatively abrupt edges. The prior distribution is determined by the parameter $\alpha \in \{1,2\}$, where $\alpha =2$ is used for Gaussian and $\alpha =1$ is used for Laplace. See Appendix A.2 for more details about this preprocessing step, yielding a final mask of the brain, which is used to find the final extracted brain (Figure 2).

2.4. The Joint MGRF Segmentation Model

The next step after extracting the infant brain is to parcellate it into its various substructures, which is achieved in this study using a novel 3D joint MGRF. The probability model: $P(\mathbf{g},\mathbf{m})=P\left(\mathbf{g}\right|\mathbf{m}\left)P\right(\mathbf{m})$ is used to jointly describe the input image of the brain $\mathbf{g}$ and its map $\mathbf{m}$, where that input image is co-aligned with the training database. Such probability model is a combination of the images conditional distribution given the map $P\left(\mathbf{g}\right|\mathbf{m})$ and the prior probability distribution of region labels $P\left(\mathbf{m}\right)=P_{\mathrm{sp}}\left(\mathbf{m}\right)P_{\mathbf{V}}\left(\mathbf{m}\right)$, where $P\left(\mathbf{g}\right|\mathbf{m})$ is a NMF-based visual appearance model and $P_{\mathbf{V}}\left(\mathbf{m}\right)$ is the Gibbs probability distribution that describes the spatially homogeneous map $\mathbf{m}$ MGRF model. In addition, $\mathbf{V}$ denotes the potential and $P_{\mathrm{sp}}\left(\mathbf{m}\right)$ is the adaptive shape prior. The components of the model will be explained in details below.

2.4.1. The NMF-Based Visual Appearance Model

After the DTI estimation, we extract various scalar features from the estimated diffusion tensor. These features contain a lot of information about the microstructure of the brain, which is beneficial in significantly distinguishing between different tissues of the brain. Throughout this work, we used the 3D Slicer software [35] to calculate 5 different anisotropy features: fractional anisotropy (FA), mean diffusivity (MD), axial diffusivity (${\lambda }_{\parallel }$), relative anisotropy (RA), and radial diffusivity (${\lambda }_{\perp })$. The two features FA and RA will help to separate WM from other brain tissues, while the other three features (MD, ${\lambda }_{\parallel }$, and ${\lambda }_{\perp })$ can differentiate between GM and CSF tissues. However, the most meaningful information needs to be extracted from that large dimensional space of the DTI.

It has been shown that NMF has strong capabilities to cluster complex data based on extracting features that are capable of learning data classes characteristics [39]. As a result, in the last decade, the use of NMF has been extensively applied to the area of medical image segmentation. A weight matrix $\mathbf{W}$ is determined in NMF-based image segmentation to transform a vector from an input space into a new feature space (H-space) by the factorization of the input matrix $\mathbf{A}$. NMF has been used in a few DTI-based segmentation systems. For instance, Xie et al. [40] utilized NMF for the segmentation of the corpus callosum, spinal cord, and hippocampus from DTI images by k-means clustering [41] in H-space. NMF has been used in some segmentation problems with promising success; however, further investigations are still needed to verify NMF segmentation potential.

In this study, NMF-based feature fusion is used to extract new significant features from the wide-dimensional DTI feature space ($\mathbf{A}$), which consists of one appearance feature (b0) and five anisotropic features (FA, MD, ${\lambda }_{\parallel }$, ${\lambda }_{\perp },\mathrm{RA})$, to more precisely segment various structures of the brain from DTI data. Particularly, our segmentation approach uses NMF to determine the weights for each of the features of the input DTI ($\mathbf{W}$), for the purpose of constructing a new feature space ($\mathbf{H}$) with a decreased number of dimensions, while encoding the information in the training dataset and increasing the separability among different classes. With NMF, two matrices could be factorized out of the input data matrix $\mathbf{A}\in {{\mathbb{R}}^{+}}^{I\times XYZ}$:

$$\mathbf{A}\approx \mathbf{W}\mathbf{H}$$

(2)

where the columns of $\mathbf{W}\in {{\mathbb{R}}^{+}}^{I\times J}$ define the basis vectors, and the encoding of the features in H-space is $\mathbf{H}\in {{\mathbb{R}}^{+}}^{J\times XYZ}$ [39]. The 3D DTI features are processed using NMF by transforming each brain voxel $(x,y,z)$ into vectors of the input data matrix $\mathbf{A}$. To calculate $\mathbf{W}$ and $\mathbf{H}$, the Euclidean distance between $\mathbf{A}$ and $\mathbf{W}\mathbf{H}$ is minimized under the constraint that $\mathbf{W}$ and $\mathbf{H}$ have only non-negative elements; formally,

$$\begin{array}{cccc}\hfill & \underset{\mathbf{W},\mathbf{H}}{\mathrm{minimize}}\hfill & \hfill & \frac{1}{2}{\parallel \mathbf{A}-\mathbf{W}\mathbf{H}\parallel }^2\hfill \\ \hfill & \mathrm{subject}\mathrm{to}\hfill & \hfill & w_{i,j},h_{j,n}\ge 0\mathrm{for}\mathrm{all}i,j,n\hfill \end{array}$$

(3)

Since the introduction of NMF, numerous approaches have been applied for optimizing Equation (3). The most outstanding algorithms applied are multiplicative update, projected gradient descent (PGD), and the alternating least square (ALS) algorithm [42]. In this study, we selected the ALS algorithm, since it is more flexible and has higher speed than other competing algorithms. The ALS algorithm takes advantage from the fact that, although Equation (3) is non-convex in both $\mathbf{H}$ and $\mathbf{W}$, it is convex in either $\mathbf{H}$ or $\mathbf{W}$ when the other is considered as constant. Therefore, a first least square step is required to find one matrix, followed by another least square step in an altering pattern.

In this work, NMF was carried out on an $I\times XYZ$ input matrix of training data with one row per DTI feature and one column per voxel. The achieved $\mathbf{W}$ was considered as the basis vectors for transforming feature vectors into the H-space. Moreover, we used the achieved $\mathbf{H}$ to determine the J-dimensional centroids that correspond to each label of the brain ($C_l$; $l\in \mathbf{L}$). During the testing stage, a new input data matrix of DTI features ($\mathbf{B}$) was presented. These features were transformed into H-space via the pseudo-inverse of $\mathbf{W}$ [43] (see Figure 3), denoted ${\mathbf{W}}^{†}$:

$${\mathbf{H}}_B={\mathbf{W}}^{†}\mathbf{B}.$$

(4)

The segmentation framework that we propose here takes into consideration the new NMF-fused features visual appearance. This visual appearance is modeled using a K-means classifier [44], where the J-dimensional centroids of the classes ($C_l$; $l\in \mathbf{L}$) are computed within the H-space in the training stage. Equation (4) is used to obtain a new J-dimensional test vector $H_{B:x,y,z}$ for any voxel $(x,y,z)\in \mathbf{R}$ in the testing volume $\mathbf{B}$. The Euclidean distance $d_l\left(T_{B:x,y,z}\right)$ between the test vector $H_{B:x,y,z}$ and each of the centroids $C_l$ is computed and used to determine the probabilities of the labels linked to $H_{B:x,y,z}$. For the brain label $l\in \mathbf{L}$, the NMF-based probabilities at the voxel $(x,y,z)\in \mathbf{R}$ are represented by:

$$P_{x,y,z}\left(\mathbf{g}\right|\mathbf{m}=l)=\frac{\frac{1}{d_l\left(H_{B:x,y,z}\right)}}{{\displaystyle \sum _{l=1}^L}\frac{1}{d_l\left(H_{B:x,y,z}\right)}}.$$

(5)

2.4.2. Adaptive Shape Model

To improve the accuracy of the segmentation, the adaptive probabilistic models are utilized to address the challenge of large inhomogeneities in infant brains and to reduce the noise effects. Adaptive probabilistic shape priors are used to constrain the anticipated shapes corresponding to the labels of the brain. To construct the shape atlas, 10 datasets of training images have been acquired from different subjects, including the new NMF-fused features. A twelve-degrees-of-freedom 3D affine transformation is used to co-align the images in a manner that achieves the maximum mutual information (MI) [45]. It is worth mentioning that these degrees of freedom are 3 for the 3D rotation, 3 for the 3D translation, 3 for the 3D shearing, and 3 for the 3D scaling. These shape priors are independent random fields of the co-aligned data region labels that are spatially variant. This is represented by:

$$P_{\mathrm{sp}}\left(\mathbf{m}\right)=\prod _{(x,y,z)\in \mathrm{R}}{p}_{\mathrm{sp}:x,y,z}\left(m_{x,y,z}\right),$$

(6)

where the empirical probabilities of the voxels of each label $l\in \mathbf{L}$ of the brain are denoted by $p_{\mathrm{sp}:x,y,z}\left(l\right)$. We incorporated the IDEA lab [46] probability maps of the infant brain tissues to obtain the atlas labels of the ground truth, where the unified segmentation approach [47] was applied. The implementation was carried out using the statistical parametric mapping (SPM) package [48], which enabled the segmentation of the non-diffusion (b0) scans from the training datasets. The achieved initial segmentation was then refined by an MR expert, in order to generate the final labels of the brain.

The NMF-fused features are used to build the shape prior of each of the input subjects that are to be segmented through an adaptive process that is controlled by the visual appearance. The first step was to choose the most similar subject out of the atlas database, where similarity is defined as the normalized cross-correlation (NCC). The chosen subject is then considered as a reference prototype for co-aligning the input subject by the 3D affine transformation discussed previously. For the test subject, the estimation of each voxel-shape prior probabilities (Figure 4) was performed by following the procedure in Algorithm A2 in Appendix B.1.

2.4.3. 3D Spatial Interaction MGRF Model

To make sure that the segmentation is homogeneous and to account for the noise effects, our proposed approach also included the region labels’ spatial interactions. Such interactions are computed using the well-known Potts model that includes a homogeneous MGRF with the 26 nearest voxel neighborhood and bi-valued Gibbs potentials, as the only consideration is whether neighboring voxels’ labels are the same or different. The used region map second-order 3D MGRF model $\mathbf{m}$ is defined in [49,50,51,52,53,54,55] as:

$$P_{\mathbf{V}}\left(\mathbf{m}\right)=\frac{1}{Z_{{\nu }_s}}\exp {\sum }_{(x,y,z)\in \mathbf{R}}{\sum }_{(\xi ,\eta ,\zeta )\in {\nu }_s}\mathbf{V}(m_{x,y,z},m_{x+\xi , y+\eta , z+\zeta }),$$

(7)

where $Z_{{\nu }_s}$ is the normalization factor which can be approximated by the following equation [56]:

$$Z_{{\nu }_s}\approx \exp {\sum }_{(x,y,z)\in \mathbf{R}}{\sum }_{(\xi ,\eta ,\zeta )\in {\nu }_s}{\sum }_{l\in \mathbf{L}}\mathbf{V}(l,m_{x+\xi , y+\eta , z+\zeta }),$$

(8)

and V denotes the bi-valued Gibbs potential

$$\mathbf{V}(m_{s_1},m_{s_2})=\left\{\begin{array}{cc}{v}_{\mathrm{eq}}\hfill & \mathrm{if}{m}_{s_1}=m_{s_2}\hfill \\ v_{\mathrm{ne}}\hfill & \multicolumn{2}{c}{\mathrm{otherwise}}\end{array},\right.$$

(9)

where $m_{s_i}$ is the label at voxel $s_i$, and $v_{\mathrm{eq}}$ and $v_{\mathrm{ne}}$ are the estimated potentials when labels are equal or unequal, respectively. The initial results of the region map are approximated using these maximum likelihood estimates for the potential values [49]:

$$v_{\mathrm{eq}}=-v_{\mathrm{ne}}\approx 2f_{\mathrm{eq}}\left(\mathbf{m}\right)-1,$$

(10)

where $f_{eq}\left(\mathbf{m}\right)$ refers to the relative frequency of equal labeling of voxel pairs in the 26-neighborhood.

These estimates enable the computation of the probability distribution of the labeling at each voxel $p_s\left(l\right)=P(m_s=l)$. The summary of the overall procedure of the proposed segmentation approach is presented in Appendix B.2.

3. Performance Evaluation Metrics

Throughout this work, three metrics were used for evaluating the segmentation accuracy: (1) the modified Hausdorff distance ($H_{95}$), (2) the Dice similarity coefficient (DSC), and (3) the absolute volume difference (AVD) [57]. These metrics are explained in details below.

3.1. Modified Hausdorff Distance

Distance measures are commonly used for segmentation evaluation. Although the Euclidean distance is usually used, the Hausdorff distance (H) is also a common measure. Considering two subsets of $\mathbf{R}$, $\mathbf{S}$ (segment) and $\mathbf{G}$ (ground truth). The distance from $\mathbf{S}$ to $\mathbf{G}$ is defined as the maximum distance from the nearest point in set $\mathbf{G}$ to $\mathbf{S}$ [57]:

$$H(\mathbf{S},\mathbf{G})=ma{x}_{s\in \mathbf{S}}\left\{mi{n}_{g\in \mathbf{G}}\{d(s,g)\right\}\},$$

(11)

where $d(s,g)$ is the Euclidean distance between the point s in the set $\mathbf{S}$ and the point g in the set $\mathbf{G}$. The bidirectional Hausdorff distance $H_{Bi}(\mathbf{S},\mathbf{G})$ is given by:

$$H_{Bi}(\mathbf{S},\mathbf{G})=max\{H(\mathbf{S},\mathbf{G}),H(\mathbf{G},\mathbf{S})\}.$$

(12)

As the distance decreases, the segmentation is improved. Ideally, a zero-valued bidirectional Hausdorff distance corresponds to perfect segmentation. Throughout this work, to remove the segmentation outlier effect, we used the 95-percentile modified Hausdorff distance ($H_{95}$) in evaluating the accuracy of the segmentation.

3.2. Dice Similarity Coefficient

The Dice similarity coefficient is a measurement of the agreement between two sets. For the sets $\mathbf{S}$ and $\mathbf{G}$, the DSC is the division result between the union size of $\mathbf{S}$ and $\mathbf{G}$; and their average as indicated in the following equation:

$$\mathrm{DSC}(\mathbf{S},\mathbf{G})=\frac{2|\mathbf{S}\cap \mathbf{G}|}{\mathbf{S}\cap \mathbf{G}+\mathbf{S}\cup \mathbf{G}}\times 100\%$$

(13)

While validating a classifier, the DSC is commonly expressed in terms of the number of true positives (TP), false positives (FP), true negatives (TN), and false negatives (FN) [57]:

$$\mathrm{DSC}=\frac{2\mathrm{TP}}{2\mathrm{TP}+\mathrm{FP}+\mathrm{FN}}\times 100\%.$$

(14)

A zero-valued DSC corresponds to no overlapping, while perfect segmentation is achieved when the DSC is equal to 1. Generally, when the value of the DSC increases, the segmentation is improved, as the matching between the ground truth and the segmentation results increases.

3.3. Absolute Volume Difference

In addition to the $H_{95}$ and DSC metrics, the segmentation accuracy could be assessed using a volumetric metric: the percentage absolute volumetric difference, which is obtained by dividing the absolute difference between the segmented volume and original volume by the original volume [57]:

$$\mathrm{AVD}=\frac{|V_s-V_g|}{V_g}\times 100\%,$$

(15)

where $V_s$ and $V_g$ are the segmented volume and the ground truth volume, respectively, and each were calculated by multiplying the number of labeled voxels by the voxel dimensions [57]. The smaller the AVD, the better the segmentation. The proposed techniques were applied on the 54 infants dw-MRI brain datasets to assess the segmentation results, from which nine different datasets having the manually segmented volumes as ground truth that was obtained by an MR expert. The ultimate goal of the proposed segmentation technique is to parcellate the infant brain into WM, GM, CSF, and other (i.e., basal ganglia and other structures). Figure 5 shows an example of the manual segmentation for one of the subjects. The detailed procedure of the proposed segmentation framework is demonstrated in Figure 6.

4. Results

Segmentation Results

Visualization of the segmentation results of the proposed algorithms are presented in Figure 6. As shown in the figure, the proposed algorithm was able to successfully preprocess the raw data (Figure 6a) and ameliorated distortions common in diffusion imaging. Figure 6b shows the preprocessed images after the skull stripping process. Figure 6c presents the results of the initial segmentation using NMF-based visual appearance model and adaptive shape model. Final segmentation results are shown in Figure 6d. The result of the proposed algorithm in 3D is shown in  Figure 7. The complete segmentation results of each subject are shown in Figure 8, Figure 9 and Figure 10. To evaluate the accuracy of the proposed framework, standard segmentation evaluation metrics have been used. These standard metrics are the Dice similarity coefficient, the modified Hausdorff distance, and the absolute volume difference  (Section 3). According to these metrics, the perfect segmentation is obtained if the DSC is close or equal to 1 ($100\%$), the Hausdorff distance decreases (approaches zero), or the the average volume difference is small. As shown in  

Table 1. Segmentation summary results, using DSC (%), H₉₅ (mm), and AVD (%). Results represent $\mu \pm \sigma$ values.

Metric	CSF	GM	WM	Others	Brain
DSC (%)	87.96 ± 3.31	89.92 ± 2.86	95.23 ± 1.18	92.81 ± 5.26	96.64 ± 1.15
H95 (mm)	2.42 ± 0.56	1.98 ± 1.07	1.98 ± 0.01	6.60 ± 2.75	7.17 ± 3.96
AVD (%)	6.10 ± 4.70	9.85 ± 2.34	5.15 ± 2.03	9.74 ± 9.65	2.66 ± 2.80

, where average performance on each tissue class is summarized (expressed as mean ± standard deviation), the proposed segmentation framework recorded promising results compared to ground truth. As seen in  Table 1, the Dice coefficients of similarity are $96.64\pm 1.15\%$ for the whole brain volume, $95.23\pm 1.18\%$ for white matter, $89.92\pm 2.86\%$ for gray matter, and $87.96\pm 3.31\%$ for CSF. Similarly, the modified Hausdorff distance is $7.17\pm 3.96\%$ for the whole brain volume, $1.98\pm 0.01\%$ for white matter, $1.98\pm 1.07\%$ for gray matter, and $2.42\pm 0.56\%$ for CSF. Additionally, the average volume difference recorded $2.66\pm 2.80\%$ for the whole brain volume, $5.15\pm 2.34\%$ for white matter, $9.85\pm 2.34\%$ for gray matter, and $6.10\pm 4.70\%$ for CSF. Moreover, to more highlight the promising performance of our segmentation system, the developed whole-brain segmentation system is compared with other segmentation approaches; namely, slope difference distribution-based clustering [58,59,60], Otsu clustering, K-means clustering, and expectation maximization (EM) clustering; as summarized in  

Table 2. Quantitative evaluation of our brain segmentation system compared to other segmentation approaches using DSC (%), H₉₅ (mm), and AVD (%). Note that SDD and EM stand for slope difference distribution based clustering and expectation maximization clustering, respectively. Results represent $\mu \pm \sigma$ values.

Metric	SDD	Otsu	K-Means	EM	Our System
DSC (%)	$87.82\pm 7.70$	$91.06\pm 3.28$	$72.75\pm 8.64$	$70.81\pm 8.28$	$96.64\pm 1.15$
H95 (mm)	$14.83\pm 8.16$	$12.37\pm 6.67$	$16.88\pm 5.85$	$21.8\pm 8.2$	$7.17\pm 3.96$
AVD (%)	$12.41\pm 10.65$	$14.09\pm 6.49$	$41.92\pm 25.89$	$81.54\pm 31.77$	$2.66\pm 2.80$

. As demonstrated in the table, the EM segmentation system gives the worst accuracy of $70.81\pm 8.28$ DSC, $21.8\pm 8.2$ H₉₅, and $81.54\pm 31.77$ AVD, while the proposed segmentation system outperforms them. These results affirm the promising accuracy results by our segmentation approach. Those experiments demonstrated an accurate identification of the joint MGRF model, enabling successful segmentation of GM and WM from infant DTI images. The current implementation, written in the C++ programming language and running on a typical desktop workstation, took on average about 8.1 s ($\pm 2.53$ s) to process each subject.

5. Discussion

This study presented an adaptive, atlas-based technique for automatically segmenting different structures of infant brains at the isointense age (6–12 months) from DTI scans. The proposed technique integrated NMF-fused DTI feature, intensity MGRF models of visual appearance, adaptive shape prior, and shape homogeneity MGRF model to ensure the accuracy and efficiency of the segmentation process. Moreover, the proposed technique was able to avoid the limitations associated with other segmentation techniques that targeted infant brains. In the next paragraphs, we will discuss different segmentation techniques presented in the literature, along with each technique’s advantages and limitations. Eventually, we will explain how our proposed technique could enhance the accuracy of infant brains segmentation and mitigate limitations of other segmentation methods. As mentioned in the introduction section, brain segmentation from MRI is classified under one of the following categories: statistical techniques, atlas-based techniques, and deformable model-based techniques. While the main goal of most of these techniques is segmenting adult brains, a few of them can fit infant brains. In the subsequent paragraphs, we will overview and address the limitations of the current segmentation techniques that deal with MR infant brains, with a special focus on those which incorporate DTI measurements in their framework. Generally, brain segmentation techniques are categorized as under statistical methods, atlas-based methods, and deformable model-based methods.

Statistical models are one of the most commonly used models in segmentation. For example, prior models are used to characterize the distribution of signals for each structure of the brain. Expectation-Maximization (EM)-based parametric Gaussian density estimation was used by Makropoulos et al. [6], where Markov random field (MRF) regularization accounted for the MRI spatial homogeneity constraints. Automated segmentation of brain structures by Anbeek et al. [7] used probability maps to segment each brain tissue class with a K-nearest neighbor classifier, which employed features such as voxel intensities and coordinates that were manually constructed. Cardoso et al. [8] proposed a posteriori EM multi-modality segmentation algorithm that included the use of a relaxation approach that iteratively adjusts the tissue priors to match the input subject. The segmentation algorithm also included intensity non-uniformity (INU) correction model, an MRF-based homogeneity descriptor, and a correction of partial volume containing voxels. More generalized image priors were introduced by Gui et al. [9], where the creation of these priors was based on a fusion between high-level morphological data such as connectivity and relative spatial location. This morphology information was incorporated in the segmentation process through a neighborhood selection criterion, which guides a region growing algorithm in separating different brain structures. Patrial volume errors are corrected by imposing a condition on WM growing into neighboring extra-ventricular CSF. Wang et al. [10] employed the sparsely represented tissue complementary distribution for segmenting brains from DTI, T1, and T2. The initial segmentation for the various brain structures was performed using patching, along with a ground truth-aligned multi-modal image library. Geometric constraints are then integrated for refinement. Aranguren et al. [11] proposed a multilevel thresholding technique to improve the LSHADE algorithm [12] for the segmentation of magnetic resonance brain images.

The other methodology used for segmentation is using predefined atlases [2,14,15,16,20]. Some of the strongest segmentation approaches employ atlases based on a priori information on the position of the various structures of the brain, so the segmentation could be viewed as a registration problem. An atlas-based segmentation approach from MRI for newborn brains was developed by Prastawa et al. [14]. The approach started with graph robust clustering and estimation of the parameters required to model the initial intensity distributions, then bias correction is performed using the initial intensity distributions as well as the spatial priors. Non-parametric kernel density estimation was used to refine the corrected initial segmentation to produce the final segmentation. An unsupervised deep learning framework for atlas-based registration [20] was introduced for embryonic brain segmentation. In [21], Ruan et al. compared the performance of two different automated brain segmentation in terms of dice coefficient, Hausdorff distance, and the standardized uptake value, to manual segmentation on PET/MR imaging. Weisenfeld et al. [2] proposed an adaptive fully automated classifier fusion technique to segment different tissues of newborn infant brains from MRI, to eliminate bias and variability associated with manual interaction. Tissue class intensity distributions were estimated nonparametrically using samples of the MR intensities of each class representative tissue. The proposed algorithm started with a library of template MRI images, and each template has a large number of manually selected tissue class prototypes. After template–subject registration, different prototype lists lead to nonparametric density estimates that result in different segmentations. These segmentations were then fused to increase the final segmentation accuracy. Shi et al. [15] developed a longitudinal brain segmentation technique for the newborns. Their technique has a bias correction step. It utilizes the probabilistic atlas of the tissues of each subject specifically and builds an atlas longitudinally for each subject from late time scans. To build an atlas for each subject specifically, the late time scans at the ages of 1 and 2 years were segmented using an adaptive fuzzy segmentation algorithm (AFCM) to get the probabilistic maps of tissue classes. The atlas construction is then followed by an iterative process of segmentation and registration carried out jointly, resulting in aligning the atlas. Kim et al. [16] used spatial intensity growth maps computed by linear regression on a voxel basis to compensate for the inhomogeneity present in the infant brain intensity values. Enhanced EM was then used to segment the resulting images. Two registration frameworks that use T2-weighted and DTI images were utilized in [17,18,19] for building the atlas. The registrations within a given subject and between the different modalities were built on a multi-scale using B-spline and affine transformation, where the matching metric was the normalized mutual information (MI). Scans captured at different intervals of time were linearly mapped to the atlas of a 1-year-old subject. Nonlinear transformation mapping was then performed relative to a T2-weighted atlas, where the tensor images were extracted from the aligned DWI. The final DTI atlas was obtained by the log-Euclidean averaging of the tensor images.

Although statistical segmentation is faster than other techniques, it generally does not handle all real distributions of data, as it relies on predefined statistics. Among the factors that affect this are the scanning parameters and the uniqueness of each patient, especially in the case of MRI-scanned brain tissues that have similar intensities, which makes segmentation inaccurate when using intensity alone. While atlas-based segmentation is more accurate than statistical-based segmentation, it is still challenging as it requires heavy computations and appropriate selection of atlases. Moreover, the obtained segmentation accuracy is highly dependable on the involved registration quality. In addition, atlas-based segmentation depends on the feature selection required in linking the test subject and the prior data that were used in constructing the atlas. For example, most of the current techniques use signal intensity to find the correspondence between the prior atlas and the data that need to get segmented, leading to inaccuracy in the results as the signal intensities are variant with scanner, age, and patient.

Segmentation using deformable models is also gaining popularity for its effectiveness [23,24]. Deformable boundaries have been used for the continuous and accurate segmentation of brain structures from MR. Deformable models-based segmentation techniques are based on deforming curves or surfaces (2D image domain) or hyper-surfaces under the effect of both external and internal forces to outline the targeted object boundary. Leroy et al. [23] introduced a surface reconstruction framework for automated segmentation of the boundary between the cerebral cortex and the underlying white matter in infant brains. After the spatial intensity inhomogeneity correction, the initial cortical surface was identified based on two local priors that were not affected by the intensity inhomogeneity. The final reconstruction of the inner cortical surface was obtained by applying deformation on two surfaces located at each side of the cortical region, and the competition between these two surfaces results in a more robust reconstruction. A segmentation framework based on level sets was developed by Wang et al. [24], in which a combination of local intensity information (modeled by the local Gaussian distribution fitting (LGDF) [25]), atlas spatial prior, and a constraint on the thickness of the cortical surface was used to construct the final coupled level sets. An initial segmentation was first obtained by using coupled level set convex optimization based on the statistical information present globally in the images as well as the spatial prior of the altlas. In another study, Wang et al. [26] developed a segmentation technique in which the level sets are guided longitudinally. They used the segmentation result from later time points in guiding the segmentation of neonatal brains. In particular, two-year-old and neonatal brain images are independently segmented using the adaptive fuzzy c-means algorithm [27] and the previously described coupled level sets segmentation framework [24]. Then, the neonatal segmentation was refined by wrapping the older two-year scan to the neonatal space. A 4D multi-modality extension to the last segmentation framework is presented by Wang et al. [28]. They extended the previous work to include 6-month-old scans by integrating both newborn and 1-year-old scans in guiding the segmentation. To handle the associated intensity problems, they used the fractional anisotropy (FA) maps extracted from the DTI. Moreover, temporally consistent segmentation results were obtained by a longitudinally consistent term. Wang et al. [29] proposed another MR level set based neonatal brain segmentation approach that is patch-driven through sparse representation, which is effectively capable of creating an atlas for each subject specifically from an aligned library of images that were manually segmented. The resultant tissue probabilistic maps were then refined by measuring the similarities existing between the patch and its neighbors. The final segmentation was obtained by integrating the refined probabilities into a coupled level set framework. Mesejo et al. [30] introduced a new deformable model that combined region and boundary information with prior shape knowledge introduced using nonlinear registration. Level set parameters were optimized using a genetic algorithm, and scatter search was used to find optimal shape priors. Baghdadi et al. [31] combined both a fuzzy C-means model with a deformable model implemented with the generalized fast marching method (GFMM) to segment brain MR images and compared their technique to the well-known SPM, FSL, and Freesurfer neuroimaging tools. In [32], Wu et al. presented their ROI brain segmentation using deformable shape model guided by an image-based regressor, and compared their performance to multi-atlas-based methods and other learning-based methods on three public brain MR datasets.

Although deformable model-based segmentation is capable of segmenting connected objects with an accuracy higher than other techniques, the achieved accuracy is dependent on the appropriate design of the statistical and geometric guiding forces, as well as the model initialization.

In addition to the limitations associated with each segmentation category, most of the previously mentioned MR segmentation methods for the infant brains were devoted to segment brains either within early infancy (up to 5 months) or early toddler stage (12 months or later), employing T1-weighted MRI, T2-weighted MRI, or both together [2,4,8,9,14,15,16,23,24,26]. Nonetheless, these techniques would come up short for infants within the isointense age (6–12 months), which is the essential concern of this study, since both WM and GM have a high overlap on same levels of intensity in the isointense stage. Hence, a segmentation framework integrating other scan contrasts (e.g., FA maps [10,28]) will significantly distinguish between WM and GM within the isointense age. In any case, these studies have never attempted to consider other DTI scan contrasts in their segmentation frameworks. Moreover, current DTI-based segmentation methods for infant brains endure from some limitations, such as: (i) atlases developed from more than one modality (e.g., T2- and diffusion-weighted) will diminish the accuracy of the segmentation as a result of different contrast levels and inter-slice variability; (ii) utilizing nonlinear warping contrarily influences the shape information, preventing analyzing the shape of the segmented data statistically, which will be advantageous in assessing the relationship between WM and GM morphology [61]; and (iii) most of them depend on the longitudinal information, which is not feasible in most cases.

To overcome these limitations, we proposed in this study an adaptive, atlas-based technique for segmenting different structures of infant brains automatically from DTI scans (see Figure 1). The proposed method involves a novel integration of different models (NMF-based visual appearance model, adaptive shape model, and spatial interactions with MGRF) to account for different information encoded in the brain images, ensure homogeneity, and lessen the effect of noise in the image. This integration allowed an automated segmentation of different brain tissues with comparable output masks to expert’s manual segmentation. A preprocessing step was required to account for the image venetian blind artifacts, eddy current, and noise associated with DWI. Furthermore, the acquisition time of DWI is longer, which results in higher motion artifacts, particularly when the infants are not sedated, in addition to inconsistencies occurring along the slices and the gradient. In addition, the image distortions are amplified while acquiring DWI due to the echo-planar imaging (EPI) fast gradient alteration, which requires DWI screening to remove any bias in the collected data, and hence avoid undesirable consequences in the results [34].

The skull stripping procedure is not a simple task, since the contrast between brain and some other tissues, such as dura mater, can be very slight. Brain extraction is an important preprocessing step in CAD systems that deal with the brain, especially if any form of shape-based analysis is performed as having some parts of the brain missing or failure to expel the dura can underestimate or overestimate the cortical thickness [62]. Several methods for extracting the brain have been presented in the literature. These methods use different techniques, namely, deformable models [63,64,65], atlas-based [66,67], label-fusion based, and hybrid methods [68,69,70]. However, most of the current methods have been developed specifically to work with T1-weighted brain MR images, so they are usually not suitable for extracting the brain directly from DWI scans. Moreover, they are optimized for adult brains and fail to address the additional challenges imposed by MR infant images (e.g., reduced contrast and higher noise).

In order to address these limitations, a new stochastic approach for extracting the brain is proposed in this work.

The results were evaluated using volumetric, overlap, and distance evaluation metrics. Various anatomical structures in the infant brain can be more precisely distinguished during the isointense age range by incorporating diffusivity characteristics of the tissue such as FA. A novel joint Markov–Gibbs random field (MGRF) segmentation model is introduced for modeling the images and the desired region maps of the brain. In particular, the joint MGRF model relies on three image descriptors: (i) a set of diffusivity features derived empirically using nonnegative matrix factorization (NMF), (ii) a shape model that works adaptively, and (iii) a three-dimensional, spatially invariant, second-order MGRF homogeneity descriptor. This segmentation framework utilizes NMF features fusion over the whole space of the DTI features, resulting in fewer dominant features, thus reducing the dimensionality and increasing the separability. The NMF-based visual appearance descriptor is modeled using a K-means classifier with centroids predefined in the training phase. To overcome the atlas bias problem related with most existing techniques, the proposed framework involves the use of an adaptive shape prior. A prior atlas database is used to build our proposed adaptive shape model. Such database is constructed from a subset of the training datasets that are co-aligned and the adaptation is controlled by the first- and second-order visual appearance characteristics of the new NMF-fused image features, which are represented by the voxel-wise intensities of the image and the features of their spatial interaction. Finally, an MGRF spatial interaction model of the second order is added to compensate for the inhomogeneity that exists extensively in infant brains. One limitation of this proposed method is concerning the small sample size. We conducted this study as a pilot study, showing the promising accuracy of our adopted automated segmentation algorithm is comparable to experts manual segmentation. However, we aim to collect more data from the University of Louisville Hospitals to further validate and analyze the performance in the future work. Additionally, dimensionality reduction in the NMF-based feature fusion process will be extensively investigated to find the optimal level of reduction that ensures the best performance of the proposed framework. Moreover, we plan to further investigate the preprocessing step and conduct more enhancement algorithms (e.g., bilateral smoothing [71]), since preprocessing steps play a crucial role in the accuracy of the segmentation process. A future plan is also to integrate advanced DTI features into the proposed framework to convey more information about different structures of the brain and enhance the segmentation accuracy.

6. Conclusions

In this work, we introduced a novel comprehensive, adaptive, atlas-based technique for automated segmentation of different tissue types in the infant brain. The proposed framework provided comprehensive preprocessing steps that accounted for the inherited problems in DTI scans and achieved promising results. The achieved accuracy of $92.5\%$ DSC on average paves the way towards the advancement and development of new research techniques on infants for the detection of various neurological conditions and brain disorders at an early stage.