This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TMI.2014.2322280, IEEE Transactions on Medical Imaging 1

Automatic Whole Brain MRI Segmentation of the Developing Neonatal Brain Antonios Makropoulos, Ioannis S. Gousias, Christian Ledig, Paul Aljabar, Ahmed Serag, Joseph V. Hajnal, A. David Edwards, Serena J. Counsell* and Daniel Rueckert

Abstract—Magnetic resonance (MR) imaging is increasingly being used to assess brain growth and development in infants. Such studies are often based on quantitative analysis of anatomical segmentations of brain MR images. However, the large changes in brain shape and appearance associated with development, the lower signal to noise ratio and partial volume effects in the neonatal brain present challenges for automatic segmentation of neonatal MR imaging data. In this study, we propose a framework for accurate intensity-based segmentation of the developing neonatal brain, from the early preterm period to term-equivalent age, into 50 brain regions. We present a novel segmentation algorithm that models the intensities across the whole brain by introducing a structural hierarchy and anatomical constraints. The proposed method is compared to standard atlas-based techniques and improves label overlaps with respect to manual reference segmentations. We demonstrate that the proposed technique achieves highly accurate results and is very robust across a wide range of gestational ages, from 24 weeks gestational age to term-equivalent age. Index Terms—neonatal brain MRI, image segmentation, Expectation-Maximization, partial volume, hierarchical modelling, model averaging

A. Challenges in neonatal brain MRI Automated segmentation of infant MR imaging data, and especially those of preterm infants, is considerably more challenging than adult brain segmentation. MR images of the neonatal brain have a much lower contrast-to-noise ratio (CNR), frequently have lower signal-to-noise ratio due to the small size of the neonatal brain and vary enormously in terms of brain shape and appearance as a result of rapid brain development during this period. In addition, the partial volume effect present due to the inverted signal intensity in white matter (WM) presents an obstacle for tissue classification [8]. The partial volume (PV) of gray matter (GM) and cerebrospinal fluid (CSF) in the cortical gray matter (CGM)CSF boundary leads to intensities similar to those of WM, which is predominantly unmyelinated in the neonatal brain. This PV effect leads to voxels mislabelled as WM at the CGMCSF interface [8]. The partial volume correction in this study addresses specifically this partial volume effect. B. Related work on neonatal brain segmentation

I. I NTRODUCTION EGMENTATION of MR brain images is increasingly being used to assess brain growth and development in the neonatal brain [1]–[10]. However, manual segmentation of magnetic resonance (MR) images is extremely time consuming and thus an expensive process. Furthermore, manual labelling is subject to inter and intra-observer variability. These limitations of manual approaches present an obstacle in labelling large cohorts of subjects that are frequently required for neuroimaging studies. There is therefore a need for an accurate automatic technique to parcellate the brain into multiple structures.

S

Copyright (c) 2010 IEEE. Personal use of this material is permitted. However, permission to use this material for any other purposes must be obtained from the IEEE by sending a request to [email protected]. A. Makropoulos, C. Ledig and D. Rueckert are with the Biomedical Image Analysis Group, Department of Computing, Imperial College London, London SW7 2AZ, United Kingdom. A. Makropoulos, P. Aljabar, J. V. Hajnal, A. D. Edwards and S. J. Counsell are with the Centre for the Developing Brain, Division of Imaging Sciences and Biomedical Engineering, King’s College London, Kings Health Partners, St. Thomas’ Hospital, London SE1 7EH, United Kingdom. I. S. Gousias is with the Centre for the Developing Brain, Imperial College London and MRC Clinical Sciences Centre, Hammersmith Hospital, London W12 0NN, United Kingdom. A. Serag is with the Advanced Pediatric Brain Imaging Research Laboratory, Childrens National Medical Center, Washington, DC, USA. Corresponding author e-mail: [email protected] (Serena J. Counsell)

I. Tissue segmentation Previous work in the neonatal segmentation field is primarily focused on the segmentation of tissues (CSF, CGM, WM) and large subcortical structures (subcortical GM, brainstem, cerebellum, ventricles). Prastawa et al. [5] introduced a method for tissue segmentation based on non-parametric kernel density estimates. Novelties of the work included differentiation between the myelinated and unmyelinated WM class according to a graph based clustering technique (minimum spanning tree) and removal of outliers with the use of the Minimum Covariance Determinant estimator. Spatial tissue priors were propagated from a probabilistic atlas. A drawback is that the atlas was created by averaging semi-automatic segmentations from three subjects and therefore can not capture the large differences in the neonatal population. In previous work by our group, [8], we implemented an EM scheme and introduced a Markov Random Field (MRF) scheme similar to [11] to enforce a smooth labelling. One novelty of this work was the introduction of a knowledge-based method that adjusts the EM priors to correct for connected components of mislabelled CSF voxels as WM in the CGM-CSF interface. Further, subjectspecific tissue priors were estimated with k-means clustering eliminating the use of an atlas for the spatial priors. However, Xue et al. only segmented the CSF, CGM and WM tissues and masked out the deep grey matter structures. Weisenfeld et al. [7] used an iterative sample editing process, estimating

0278-0062 (c) 2013 IEEE. Translations and content mining are permitted for academic research only. Personal use is also permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.

This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TMI.2014.2322280, IEEE Transactions on Medical Imaging 2

Fig. 1. Different slices from a label-based encephalic ROI template (ALBERT) [9] (the manually segmented atlases used in this study) in the sagittal and axial plane (top and bottom row respectively) with the manually segmented labels overlaid.

the segmentation with the use of STAPLE [12]. Partial volume correction was performed using a method similar to [8] and spatial homogeneity modelled with a MRF. The tissue segmentation of Weisenfeld et al. depends on the initial alignment of atlases, which is a non-trivial problem due to the large developmental differences in the neonatal brain at different ages. Other approaches for neonatal brain MRI segmentation includes the use of graph-cuts [6] where the intensity likelihood term was estimated with the use of the Parzen-window non-parametric density estimation technique. In their work [6], Song et al. adopted extensive preprocessing of the images. Song et al. only segmented the grey and white matter having excluded the skull and CSF with manual delineation. Cardoso et al. [1], [2] proposed an EM-MRF scheme that adapts the atlas priors similar to [13]. The atlas priors provided by atlases were modelled as samples drawn from a Dirichlet distribution (in essence a hyperprior) and were adapted according to the posteriors of each EM iteration. The partial volume was modelled as mixed distributions among the different tissues. Moreover, Cardoso et al. deviated from the classic Gaussian modelling by introducing a semi-conjugate Gaussian prior over the tissue Gaussian means initialised by sampling manually selected patches. Cardoso et al. avoided the dependence on the atlas alignment with a prior relaxation scheme. However, when segmenting structures with very similar intensity profiles it is unclear how the relaxation scheme will adapt the propagated priors, for example in deep grey matter structures. II. Structural segmentation In order to assess regional brain volume changes with time, more detailed segmentations of a large number of brain regions are required. Recently, a first approach to achieve this has been developed in [10] which shows that neonatal T1-weighted MR images can be automatically and reliably segmented using multi-atlas approaches.

In atlas-based segmentation methods, the MR images of manually segmented atlases are registered, typically nonrigidly, to the subject’s MR image and their labels are propagated to it, based on the calculated transformations [14]. The propagated labels of multiple atlases can then be fused to yield probabilistic spatial priors or to provide the final segmentation result. A commonly used technique, (majority voting) label fusion [3], [10], [14], [15], labels each voxel with the structure that is favoured by most of the atlases. More sophisticated fusion techniques weight the vote of each atlas locally according to the similarity of the atlas MR image to the image to be segmented [16], [17]. Label fusion (with majority or local-weighted voting) is among the techniques that provide best results when segmenting multiple regions [10], [17], [18]. We previously manually segmented neonatal T1-weighted MR brain images to create label-based encephalic ROI templates (ALBERTs), a set of labelled brain atlases [9]. The ALBERTs were used to automatically segment the termequivalent neonatal brain T1-weighted MR images into 50 regions. Label fusion of the transformed labels of the ALBERTs and alignment of a maximum probability atlas were investigated [10] and demonstrated accurate results.

C. Contributions In the present work, we propose a novel multi-structure Expectation-Maximization (EM) based segmentation technique for the subdivision of the whole brain. We combine the transformed labels of the ALBERTs with subject-specific tissue priors obtained as in [8]. Regions whose intensity distribution needs to be described with more than one Gaussian distribution are split into separate classes. The spatial priors are combined with an intensity model to segment the structures in a fashion similar to [19], [20]. We propose a number of improvements to the standard EM algorithm employed for

0278-0062 (c) 2013 IEEE. Translations and content mining are permitted for academic research only. Personal use is also permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.

This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TMI.2014.2322280, IEEE Transactions on Medical Imaging 3

brain segmentation ( [11], [19], [20]) to accurately segment the neonatal brain: • The manual protocol of the atlases used in this study [9] subdivides the combined WM and CGM space into 32 structures. Each one of the 32 structures contains both the WM and CGM part of the respective region of the brain (see Figure 1 for an example of this manual parcellation). The boundaries between the WM/CGM structures are based on the morphology of the cortex and are not separated according to intensity. To overcome this issue we introduce a hierarchy within anatomical structures. We model the whole CGM and WM tissues as ’superlabels’ (aggregations) of their subdivisions. This hierarchical modelling allows the accurate sampling of the tissues’ intensities. Furthermore, the boundaries between the WM/CGM structures are defined according to the spatial priors and are not modified according to the intensity. We only allow for a small amount of smoothing to alter these boundaries. In contrast to this, the boundaries of the WM/CGM structures to the other structures (e.g. CSF, subcortical GM) are adapted with the EM segmentation. • Segmentation of the subcortical grey matter into different regions based on the intensity is difficult for some parts of the brain. Subcortical GM structures have very similar intensity profiles and parts of their boundaries are often hard to distinguish. We tackle this challenge by defining the predictive (posterior) distribution as a model averaging of label fusion and the Gaussian Mixture Model (GMM). The segmentation decision in homogeneous regions in the deep grey matter is primarily based on the label fusion vote rather than the GMM since the intensity is less reliable in these areas.

18 subcortical regions [9]. MR imaging was performed on a Philips 3T system using an 8 channel head coil. The MPRAGE MR images were acquired using the following parameters: repetition time (TR): 17 ms; echo time (TE): 4.6 ms; flip angle 13◦ ; slice thickness: 0.8 mm; field of view: 210 mm; matrix: 256 × 256 (voxel size: 0.82 × 0.82 × 0.8 mm). T2weighted MRI were acquired using the following parameters: TR: 8670 ms; TE: 160 ms; flip angle 90◦ ; slice thickness 1 mm; field of view: 220 mm; matrix: 256 × 256 (voxel size: 0.86 × 0.86 × 1 mm). The T2-weighted MRI of the atlases were co-registered to their corresponding T1 images [9]. This allows the segmentation to utilize any of the two modalities for the labelling of unseen images. The T2-weighted images were used for the segmentation in this study which provides better contrast for the differentiation of tissues in neonates. 2) Subjects: T2-weighted MR images were acquired for 198 premature infants born at a median (range) GA of 29+2 (24+3 - 35+2 ) weeks at birth and imaged at 39+1 (26+5 44+2 ) weeks PMA at scan and 36 term controls born at a median (range) GA of 39+2 (36 - 41+5 ) weeks at birth and imaged at 40+6 (36+3 - 44+4 ) weeks PMA at scan. The median (range) weight at birth of the subjects was 1.18 (0.6 - 3.7) kg for the premature infants and 3.3 (2.2 - 4.3) kg for the term controls. The head circumference of the infants at gestational age was 27.15 (22 - 36.8) cm and at scan age 34 (23 - 39.6) cm . Ear protection was provided for each infant with earplugs moulded from a silicone-based putty (President Putty, Coltene Whaledent, Mahwah, NJ) and neonatal ear protectors (MiniMuffs, Natus Medical, San Carlos, CA). All examinations were supervised by a pediatrician experienced in MRI procedures and pulse oximetry, temperature and electrocardiography data were monitored throughout the procedure. B. Segmentation pipeline

The contributions of this paper are three-fold: • We present a novel multi-structural segmentation algorithm for parcellation of the developing brain into 50 regions. • The proposed model improves the accuracy of our previous study [10] and outperforms current state-of-the-art techniques used in the segmentation field. • Our algorithm allows for the first time accurate regional segmentation of the neonatal brain from 24 weeks gestational age to term-equivalent age. The algorithm is successfully applied to 234 T2-weighted images of different gestational ages. II. M ETHODS

Figure 2 presents an overview of the proposed segmentation pipeline. The unlabelled MRI is initially brain-extracted with the Brain Extraction Tool [21] and corrected for field inhomogeneity with the N4 algorithm [22]. Afterwards, the ALBERTs’ MR images are registered to the target image and the atlas labels are propagated to the image. The propagated labels are averaged in a locally-weighted scheme and subdivided with the use of subject-specific tissue priors obtained with k-means clustering. The proposed EM scheme is used for the estimation of the resulting labelling. Finally, a postprocessing step results in a labelling similar to that of the atlases. Section II.C. and II.D. describe the atlas propagation and the probabilistic modelling of the whole brain and Section II.E. to Section II.J the EM framework implemented in this work.

A. Data acquisition 1) Atlases: The ALBERTs [9] were obtained by manual delination by an expert on MPRAGE MR brain images in 20 preterm infants born at a median (range) gestational age (GA) (determined from the date of the last menstrual period and confirmed on early ultrasound) of 29+2 (26 - 38) and imaged at 40+6 (36+4 - 44+6 ) weeks (+days ) post-menstrual age (PMA) at scan (GA plus number of weeks (days) after birth), dividing the brain into 32 WM/CGM structures and

C. Atlas propagation The ALBERTs’ MR images are rigidly, affinely and nonrigidly registered to each subject’s image in the T2 space. The non-rigid registration is carried out using free-form deformations with control point grid spacings of 20mm, 10mm, 5mm and 2.5mm and normalized mutual information (NMI) as similarity measure [23] (registration parameters can be found in [24]). The ALBERTs’ labels are then propagated to the

0278-0062 (c) 2013 IEEE. Translations and content mining are permitted for academic research only. Personal use is also permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.

This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TMI.2014.2322280, IEEE Transactions on Medical Imaging 4

Atlases

Subject

Brain Extraction Registration

N4 Bias Field Correction

Image Preprocessing

k-means + simple EM

Tissue Priors

Brain Subdivision Modelling

Atlas Priors

Atlas weighting

Transformation

Label Priors Preprocessing Expectation step Repeat until convergence

Correction step

Expectation-Maximization

Maximization step

Merge substructures

Label Postprocessing

Fig. 2. Segmentation pipeline. After the preprocessing of the subject’s MR image, the ALBERTs’ MR images are registered to it and the atlas labels are propagated to the image. The atlas priors are combined with subject-specific tissue priors. The resulting spatial priors of the structures are introduced to the proposed EM algorithm that provides the final segmentation, after a label postprocessing step.

image space and are averaged to form a probabilistic spatial prior for each structure (for each voxel i): PA a wia γik πik = PK a=1 PA a j=1 a=1 wia γij where a = 1..A denotes the transformed labels of each a ALBERT and k ∈ [1..K] denotes the different structures, γik is the vote for structure k produced by atlas a calculated as ( 1, if voxel i belongs to structure k in atlas a a γik = 0, else and wia denotes the weight of each ALBERT. wia may be set to a uniform value of 1 globally for all voxels or may be spatially varying in a local weighting scheme. In this work we consider both cases. The local weighting is based on the

sum of squared differences (SSD) of the intensity-normalized images of the atlas and the unseen image [17], in the local neighbourhood of each voxel (in patches of 3 × 3 × 3 size). D. Brain subdivision modelling EM segmentation involves the accurate modelling of the intensities of the whole space covered by the different regions of the brain. In the standard EM segmentation, every region is assumed to follow a Gaussian distribution over the intensities of the voxels. In this section we describe the subdivision of regions that need to be described with more than one Gaussian distribution. The subdivision is performed in the space of the image to be segmented by splitting the spatial priors propagated from the ALBERTs. As the manual segmentation protocol [9] did not divide the WM/CGM structures of the brain into CGM and WM,

0278-0062 (c) 2013 IEEE. Translations and content mining are permitted for academic research only. Personal use is also permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.

This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TMI.2014.2322280, IEEE Transactions on Medical Imaging 5

these regions contain both tissues. Therefore, in order to describe each structure with a single Gaussian intensity model, we divided each WM/CGM structure into its corresponding tissue parts (WM, CGM) using a soft segmentation. The prior probability of each WM/CGM structure in the individual subject space is multiplied with a WM and a CGM probability map to define the corresponding WM and CGM parts. The required subject-specific tissue probabilities are obtained as the posteriors of a simple EM scheme with 4 classes (CSF, GM, WM and extra-cranial space) in the intensity domain. In this EM scheme, the means and standard deviations are initialised using k-means clustering [25]. The k-means clustering is performed in the intensity space with k=4 clusters. An alternative approach would be to segment the WM/CGM structures of the ALBERTs into WM and CGM and use the modified atlases for the propagation of the CGM,WM components of these structures. However, the registration between different neonatal brains does not align accurately the cortical ribbon (Figure 10.A. demonstrates this effect). The term-equivalent atlases have very different cortical complexity from early neonatal brains and the registration can not capture the large deformations occuring as a result of the rapid development. Therefore, we have chosen to differentiate between the tissue types in the space of the target image with intensity clustering. The background of the manual segmentations contains the extra-cranial space, CSF, as well as parts of the subcortical GM space as can be seen in Figure 1. The background prior probability is subdivided accordingly into 3 parts. The CSF and extra-cranial space are defined from the tissue probability maps. We further exclude the ventricles as separate regions from the CSF. The subcortical background mainly represents the internal capsule and is segmented by masking the background with a deep grey matter binary map (obtained by transforming and thresholding the subcortical GM map of [26]). The ventrolateral nuclei of the thalamus deviates significantly in terms of intensity from the rest of the thalamus [27]. To account for this, we subdivided the right and left thalamus into the low intensity (in T2 space, inversely for the T1 space) part of the thalamus that represents the ventrolateral nuclei and the rest of the thalamus, with another simple EM scheme initialised with a 2 cluster k-means technique. The remaining subcortical structures retain their propagated definitions from the ALBERTs. It should be noted that the subcortical GM background and the thalamus parts are segmented using a hard segmentation in the atlas space and then propagated to the subject image, as these regions are difficult to differentiate in early preterm brains. The initial labels and the ones resulting from the subdivision are presented in Figure 3. The 50 atlas labels are subdivided into 87 labels. Once the EM estimation presented in the next section has converged, we merge the subdivided parts of the structures. E. Expectation-Maximization formulation In this section we present the standard EM technique. Our notation is similar to [11]. The image intensity likelihood is approximated by a Gaussian Mixture Model of K = 87

Atlas labels (50 + background)

Labels after brain subdivision (87)

WM/CGM structure 1 (WM part) WM/CGM structure 1 (CGM part) WM/CGM structure 2 (WM part) WM/CGM structure 2 (CGM part)

WM/CGM structure 1 WM/CGM structure 2

WM/CGM structure 32

WM/CGM structure 32 (WM part) WM/CGM structure 32 (CGM part)

Subcortical structure 1 Subcortical structure 2

Subcortical structure 1 Subcortical structure 2

Subcortical structure 16

Subcortical structure 16 Left thalamus without the ventrolateral nuclei Left ventrolateral nuclei Right thalamus without the ventrolateral nuclei Right ventrolateral nuclei

Left thalamus Right thalamus

Extra-cranial space CSF Subcortical background

Background

Fig. 3. Initial labels of the ALBERTs and labels after the subdivision.

structures. The label of each voxel i is denoted by a 1-of-K indicator variable zi . The prior distribution of zi ’s is denoted with the K × 1 vector parameter πi , i.e. P (zi = ek ) = πik

(1)

Assuming the observed intensities yi are independent, the problem can be formalised as a Maximum a Posteriori (MAP) estimation of the means µk and standard deviations σk of the structures’ Gaussians, Φy = {µ1 , .., µk , σ1 , ..., σk } [11]. The parameters Φy are estimated by the EM algorithm, at each iteration m: Expectation step: [m+1]

pik

= P (zi = ek | yi , Φy ) [m]

P (yi | zi = ek , Φy )P (zi = ek ) = PK [m] j=1 P (yi | zi = ej , Φy )P (zi = ej ) where P (yi | zi = ek , Φ[m] y ) = G(yi | µk , σk )

(2) (3)

and G(. | µk , σk ) represents a Gaussian distribution with parameters µk , σk . Maximization step: PN [m+1] pik yi [m+1] (4) µk = Pi=1 [m+1] N i=1 pik [m+1] 2 (σk )

PN =

i=1

[m+1]

[m+1] 2

pik (yi − µk PN [m+1] i=1 pik

)

.

(5)

F. Markov Random Field regularization In practice, the label of a voxel is dependent on the label assignment of its surrounding voxels. To account for this we include a spatial regularisation term that enforces a smooth labelling among neighbouring voxels [11]. Our implementation follows the work of [1] adopting the mean field approximation,

0278-0062 (c) 2013 IEEE. Translations and content mining are permitted for academic research only. Personal use is also permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.

This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TMI.2014.2322280, IEEE Transactions on Medical Imaging 6

where the Markov Random Field (MRF) is expressed with a multiclass extension of the Potts model:

WM/CGM structure 1 (WM part) WM/CGM structure 2 (WM part)

WM

−UM RF (ek |Ni )

πik e P (zi = ek ) = PK −UM RF (ej |Ni ) j=1 πij e

}

SWM

}

SCGM

WM/CGM structure 32 (WM part)

(6)

WM/CGM structure 1 (CGM part) WM/CGM structure 2 (CGM part)

CGM

with Ni the first-order neighbours of voxel i. The energy function UM RF is defined as:

WM/CGM structure 32 (CGM part) CSF

UM RF (ek | Ni ) = K X

Akj (

j=1

X l∈Nix

sx plj +

Extra-cranial space

X l∈Niy

sy plj +

X

sz plj )

(7)

l∈Niz

where the variable s = {sx , sy , sz } , models the anisotropic voxel spacing [1]. Akj is an a priori defined connectivity strength matrix that penalises the interaction of structures according to their spatial proximity. The spatial proximity is automatically defined from the ALBERTs. A pair of structures is defined as ”neighboring” if the structures have neighboring voxels in at least one set of atlas labels. Since the atlases did not differentiate the WM/CGM structures into WM, CGM, these parts were automatically subdivided prior to the definition of the neighboring structures. In this framework, Akj is defined as:  0, if structure k is the same as j      a, if structure k is neighbouring j,      and both k and j belong to WM/CGM Akj = b, if structure k is neighbouring j,      and either k or j do not belong to WM/CGM      c, if structure k is distant from j with a < b < c. Larger values attribute a stronger effect of the neighboring voxels resulting in a smoother labelling. To preserve the propagated anatomical information at the boundaries between the WM/CGM structures while removing isolated voxels, we only allow weak smoothing to alter the boundaries between these structures. The connectivity strength parameters are set to a = 0.3, b = 1, c = 1.7 in our experiments. G. Superlabels The WM/CGM structures on the ALBERTs were parcellated according to the position of the sulci and gyri. Therefore, the boundaries between these regions are determined by anatomical landmarks and not according to intensity. In the standard EM formulation the different structures that belong to the same tissue (WM, CGM) compete with each other based on their prior and intensity likelihoods. As a result, the WM part of a WM/CGM structure can expand into the space of neighboring WM parts (and the CGM part of a WM/CGM structure can expand into the space of neighboring CGM parts) driven by its intensity likelihood. To prevent this from happening we constrain the expansion amongst structures of the same tissue using only the prior probability, by defining a single intensity likelihood for all the divisions of the same tissue.

Subcortical background

MR Image

Subcortical structure 1 Subcortical structure 2 Subcortical structure 16 Left thalamus without the ventrolateral nuclei Right thalamus without the ventrolateral nuclei Left ventrolateral nuclei Right ventrolateral nuclei

Fig. 4. Label hierarchy.

In this study, we defined superlabels for the WM and CGM, as aggregations of their subdivisions. Figure 4 presents the hierarchical tree of the structures implemented in this study. Let the WM regions belong to the WM superclass SW M and the CGM regions to the CGM superclass SCGM . The WM and CGM prior and posterior probabilities are defined as the summation of the prior and posterior probabilities of their subdivisions (at the Expectation step (2) ): X

πiCGM =

[m+1]

πik

, piCGM =

k∈SCGM

X

πiW M =

k∈SW M

[m+1]

X

pik

k∈SCGM [m+1]

πik

, piW M =

X

[m+1]

pik

k∈SW M

The intensity of each superlabel is described with a single Gaussian (GCGM , GW M ) which is updated in the standard Maximization step in equations (4)-(5) as an additional class. The superlabel’s Gaussian is thus sampled from the whole space of its subdivisions. This allows for a more accurate estimation of the tissue parameters. In order to consider only the prior probabilities for the separation of the WM subdivisions and the CGM subdivisions, we define the same intensity likelihood for all the subdivisions of each superlabel as: k ∈SCGM : [m] P (yi | zi = ek , Φ[m] y ) = P (yi | zi = eCGM , Φy )

k ∈SW M : [m] P (yi | zi = ek , Φ[m] y ) = P (yi | zi = eW M , Φy )

0278-0062 (c) 2013 IEEE. Translations and content mining are permitted for academic research only. Personal use is also permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.

This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TMI.2014.2322280, IEEE Transactions on Medical Imaging 7

The segmentation problem is thus reduced to the estimation of Φh = {µH(1) , .., µH(k) , σH(1) , ..., σH(k) }, where |Φh | < |Φy |. H is a hierarchy function similar to [28] and indicates the parent of the label in the hierarchical tree of the structures:  CGM, if k ∈ SCGM   W M, if k ∈ SW M H(k) =   k, else The estimation of the Gaussian parameters in the proposed hierarchical framework is similar to the hippocampal subfield segmentation in [28]. Van Leemput et al. [28] define a global WM, CGM and CSF tissue type with a shared mean and variance over the tissue subparts as we propose in this study. In our approach the hierarchical model is defined for the whole brain instead of a region of interest [28]. The proposed hierarchical modelling is further similar to [29] where Pohl et al. define a whole-brain hierarchical model. In their work, the Gaussian parameters are estimated for all the substructures. Then the intensity distribution for structures with children is defined as the weighted sum of the intensity distribution of its sub-structures, where the weight is equal to the atlas prior. Our framework deviates from the hierarchical modelling in [29] in the sense that the Gaussian parameters are estimated from top to bottom in the hierarchy rather than starting from the substructures. As shown in [28], the EM update equations still hold by setting: [m] P (yi | zi = ek , Φ[m] y ) = P (yi | zi = eH(k) , Φy )

(8)

H. Partial Volume correction The spatial regularisation introduced by the MRF allows the removal of isolated voxels that may be mislabelled due to partial volume signal averaging and noise in the image. However, there is a large number of adjacent voxels that can be misclassified when labelled based on the intensity in the neonatal brain. Since most of these voxels are neighboring each other, they favor their neighbors through the MRF energy. As a result, the MRF alone is not sufficient to remove them. A way to overcome this problem as shown in [8] is to introduce a knowledge-based approach designed specifically to lower the probability of PV in unwanted regions. Table I summarizes the knowledge-based rules to detect partial volume voxels and the tissues they should belong to instead. The connected components (ccs) are computed with connected component labelling. The partial volume correction in [8] was performed by reducing the prior probability of the misclassified tissue in favor of the appropriate tissue(s). Let r represent the misclassified tissue class and cj the appropriate tissue classes. The prior probabilities of the detected voxels are adjusted as following: 0 = πicj + (1 − λ)wicj πir πic j 0 πir

where

= λπir

πic wicj = P j n πicn

and λ is set to 0.5 in all the experiments. Essentially, the PV correction redistributes the unwanted tissue probability mass to the appropriate tissue(s) after each EM iteration. The PV correction in this work is implemented to follow the top to bottom hierarchy described in the previous section. The correction for labels that have subdivisions is initially estimated for the superlabels and then passed to their subdivisions with respect to their prior probability: πik 0 0 πit t ∈ [CGM, W M ], k ∈ St : πik =P j∈St πij We propose a modified EM algorithm and implement the partial volume correction in a Correction step (C-step) after the Expectation step. I. Model averaging A challenge in the GMM modelling for multi-class whole brain labelling arises when segmenting the deep grey matter into subcortical structures. These structures are very similar in terms of intensity and the differentiation between them is very difficult. The Gaussian distributions of deep grey matter structures overlap significantly (see Figure 7). The problem becomes more significant for images of early preterm infants where the contrast is limited. To account for this, we model the predictive distribution of the data as model averaging of an ensemble of two models, label fusion and GMM. We incorporate a weighting scheme that reduces the update based on the GMM model in voxels belonging to homogeneous regions (with low gradient with respect to their neighborhood) of the deep GM. In these regions the segmentation is primarily based on label fusion of the ALBERTs. This modelling prevents the expansion of structures in homogeneous areas due to intensity and favors the anatomical definition based on the a-priori information of the propagated labels of the atlases. The predictive (posterior) distribution pC ik of the ensemble of models, indexed by M , with overall parameters Φ = {ΦM } can be written as X pC P (M )P (zi = ek | yi , M, ΦM ) ik = P (zi = ek | yi , Φ) = M

= P (MGM M )P (zi = ek | yi , MGM M , ΦMGM M ) + [1 − P (MGM M )]P (zi = ek | yi , Mf us , ΦMf us ) where MGM M denotes the GMM model with parameters ΦMGM M and Mf us denotes the fusion model with parameters ΦMf us . The predictive distribution of label fusion is independent of the intensity modelling of the image. It is effectively equal to the spatial prior formed as the averaged propagated label sets from the atlases: P (zi = ek | yi , Mf us , ΦMf us ) = πik The weighting of the GMM is spatially defined as P (MGM M ) = τi . In this study τi is set to 1 to allow for the standard EM update with the GMM, except for the subcortical regions where τ is defined as τ = 1 − (1 + L(∇I, σe ))−1

(9)

0278-0062 (c) 2013 IEEE. Translations and content mining are permitted for academic research only. Personal use is also permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.

This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TMI.2014.2322280, IEEE Transactions on Medical Imaging 8

A

C

B

Fig. 5. Example of segmented WM/CGM structures (A) and the corresponding superlabels for CGM (B) and WM (C) in the axial plane of a T2-weighted MRI of a 31 weeks GA brain (acquired at 39 weeks age at scan). Knowledge-based rules Misclassified tissue Appropriate tissues WM ccs, WM voxels in the CSF - extra-cranial space boundary WM CSF WM voxels in the CGM - CSF, CGM - extra-cranial space boundaries WM CSF, CGM CSF ccs mainly surrounded by WM CSF WM CGM voxels in the CSF - extra-cranial space boundary CGM CSF, extra-cranial space TABLE I KNOWLEDGE - BASED RULES FOR DETECTING MISCLASSIFIED VOXELS AND THEIR APPROPRIATE TISSUES ACCORDING TO TISSUES ’ CONNECTED COMPONENTS ( CCS )

where ∇I is the gradient magnitude of the image. The Lfunction represents the robust Lorentzian error norm in [30] with σe being a spatially-variant scale calculated in the local neighborhood of the voxels (an example of the weighting factor can be seen in Figure 6). Modelling the segmentation as averaging of label fusion and GMM alters the predictive distribution of the original EM model as pik = pC (10) ik = (1 − τi )πik + τi pik The posteriors with the new modelling are updated at each iteration of the modified EM algorithm in the Correction step. J. Modified Expectation-Maximization algorithm Algorithm 1 presents the proposed EM algorithm. The EM steps (Expectation, Correction, Maximization) are repeated until convergence.

Algorithm 1 The proposed EM algorithm Repeat Expectation step πik e−UM RF (ek |Ni ) G(yi | µH(k) , σH(k) ) pik = PK −UM RF (ej |Ni ) G(y | µ i H(j) , σH(j) ) j=1 πij e P t ∈ [CGM, W M ] : pit = k∈St pik Correction step pik = (1 − τi )πik + τi pik Detect partial volume voxels i and cj , r ∈ [CGM, W M, CSF ] with connected component labelling and update: πic πicj = πicj + (1 − λ) P j πir n πicn πir = λπir t ∈ [CGM, W M ], k ∈ St : πik = P

j∈St

III. E VALUATION S TRATEGY Leave one out evaluation using the atlases was carried out by specifying in turn one of the ALBERTs as the target and using the remaining ALBERTs to segment it. The atlases were segmented using the proposed EM algorithm and label fusion with both atlas weighting schemes (global and local). The atlas segmentations were carried out on the T1 images as the manual delineations of the ALBERTs are defined in the T1 space. The T2-weighted MR images present a different intensity contrast from the T1-weighted images in the tissue boundaries. Therefore, performing the segmentations in the T1-weighted images allows a more accurate comparison to

πik πij

πit

Maximization step Calculate parameters Φh : PN piH(k) yi µH(k) = Pi=1 N p PNi=1 iH(k) 2 p iH(k) (yi − µH(k) ) 2 σH(k) = i=1 PN i=1 piH(k) Until convergence

0278-0062 (c) 2013 IEEE. Translations and content mining are permitted for academic research only. Personal use is also permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.

This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TMI.2014.2322280, IEEE Transactions on Medical Imaging 9

Fig. 6. Axial slice of a T2-weighted MRI (A), magnified region of the deep grey matter (B) and the weighting of EM in the model ensemble (C). In red areas the segmentation decision is primarily based on label fusion while in blue areas the result is dominated by the GMM model.

4.5

Hippocampus Amygdala Cerebellum Brainstem Caudate nucleus Thalamus Subthalamic nucleus Lentiform Nucleus Corpus Callosum Lateral Ventricle CSF CGM WM

4 3.5

likelihood

3 2.5 2 1.5 1 0.5 0 4.5

5

5.5

6 6.5 intensity (log domain)

7

7.5

Fig. 7. Gaussian intensity likelihood of the subcortical structures and the tissues in the log domain (in the T2 space)

the gold standard which follows the T1 intensity profile for the structure delineation. The methods were evaluated with the Dice coefficient to compare overlaps of the automatically computed labels with the manually parcellated gold standard. A primary analysis explores the segmentation accuracy of the two techniques in labelling the whole brain, obtained as the union of the segmented areas with each method. Further exploratory analysis assesses the improvement in regional accuracy of the proposed technique over label fusion. In order to explore the structures that have the most significant improvement with the introduction of the proposed technique we further present the regional Dice scores and include a conservative analysis which accounts for multiple testing. The results are presented in Section IV.A. The ALBERTs are all scans around term-equivalent age and registrations between pairs of atlas images are more likely to succeed given that the morphology of the subjects corresponds to very similar stages of development. However, early preterm brains vary largely in the cortical appearance and subcortical structures compared to the term-equivalent brains. To test performance of the methods across a range of gestational ages, the proposed segmentation scheme was applied to 198

preterm subjects with age ranging from 24+3 to 35+2 weeks GA and 36 term controls with GA ranging from 36 - 41+5 . The segmentations were performed in the T2 images since the T2 intensity profile provides more contrast for the neonatal images between the subcortical GM, unmyelinated WM and cortical GM. The T2 images allow a better intensity characterization and therefore more reliable intensity segmentation. All the 234 segmentations were rated by an observer blinded to the processing method with a scale from 1 to 5 (1=segmentation failure, 2=extensive region inaccuracies/alignment problems, 3=moderate region inaccuracies/alignment problems, 4=regional boundaries inaccuracies, 5=accurate segmentation) with both label fusion and the proposed EM algorithm in the locally-weighted scheme based on 7 equally spaced axial slices. The rating performance of the 2 techniques is presented in Section IV.B. 15 cases were further selected to qualitatively evaluate the method at different scan ages (26+5 , 27+6 , 28+2 , 29, 30, 31+1 , 32, 33, 34, 35+2 , 36, 37+2 , 39+6 , 42, 44+2 weeks) in all the structures and all the image slices. The qualitative evaluation was based on the region definitions in the protocol used for the manual delineation of the ALBERTs [9]. In order to evaluate quantitatively the performance of the proposed

0278-0062 (c) 2013 IEEE. Translations and content mining are permitted for academic research only. Personal use is also permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.

This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TMI.2014.2322280, IEEE Transactions on Medical Imaging 10

Structure hippocampus amygdala cerebellum brainstem caudate nucleus thalamus subthalamic nucleus lentiform nucleus corpus callosum lateral ventricle frontal lobe occipital lobe WM/CGM structures subcortical structures overall

MV-Fusion 0.758 0.818 0.922 0.920 0.840 0.902 0.741 0.872 0.706 0.790 0.918 0.850 0.790 0.828 0.804

MV-EM 0.772 0.822 0.925 0.922 0.844 0.900 0.743 0.874 0.710 0.829 0.927 0.860 0.801 0.836 0.814

LW-Fusion 0.780 0.820 0.925 0.917 0.843 0.901 0.741 0.873 0.717 0.813 0.927 0.857 0.802 0.835 0.814

LW-EM 0.790 0.826 0.929 0.922 0.846 0.901 0.744 0.876 0.727 0.838 0.932 0.860 0.808 0.842 0.820

TABLE II D ICE COEFFICIENT OF LEAVE ONE OUT CROSS - VALIDATION WITH THE ALBERT S . T HE ATLAS SEGMENTATIONS ARE CARRIED OUT IN THE T1- WEIGHTED IMAGES . F USION AND EM ARE COMPARED WITH BOTH MAJORITY VOTE (MV) AND LOCAL WEIGHTING (LW) OF THE TRAINING ATLASES ( BOLD = SIGNIFICANTLY BETTER AT P < 0.05).

algorithm at different gestational ages we manually segmented the ventricles, frontal and occipital lobe for each of the 15 selected cases. An initial version of the labels was semiautomatically obtained with label fusion and automatic tissue segmentation [8]. These labels were then refined with manual editing to result in the final labels. The manually segmented cases were compared against label fusion and EM with both weighting schemes. Section IV.C. presents the qualitative and quantitative evaluation of the 15 cases. The rater that performed the qualitative assessment and manual segmentations of the data was different from the expert who performed the manual segmentation of the ALBERTs. All statistical tests for significance in Section IV are reported according to a two-sided paired t-test. The results are conservatively adjusted for multiple comparisons using Bonferroni’s correction where stated. IV. R ESULTS A. Atlas cross-validation EM significantly (p < 10−4 ) improves the whole-brain segmentation accuracy with a mean Dice score of 0.961 over 0.957 obtained with label fusion in the local weighting of the ALBERTs (0.958 over 0.949 with majority vote weighting, p < 10−5 ). The results of the leave-one-out cross-validation of the ALBERTs for all the 50 structures are shown in Appendix A, Table IV with both weighting schemes. To allow for an easier comparison we averaged the results of the WM/CGM structures and the subcortical regions over the 2 hemispheres and present the evaluation of major structures in Table II. Local weighting of the atlases performs significantly better in both label fusion and EM. LW-EM performs equally or better in all the structures, with significantly better results (p < 0.05) for 36 out of the 50 structures with majority vote weighting (24 after Bonferroni correction), and 39 out of the 50 structures with locally-weighted voting (25 after Bonferroni correction). The structures with stronger intensity differentiation from their surrounding regions benefit most by using the EM technique, such as the lateral ventricles, hippocampus and corpus callosum and the WM/CGM structures. The boundaries of these structures can be more clearly

separated in terms of intensity and are improved with the presented EM technique through the GMM modelling. Label fusion relies on accurate registration of the atlases and does not model the intensity of the structures. However, accurate registration of neonatal brain images is challenging, even in areas of large intensity contrast, due to the large developmental changes. As a result, label fusion is sensitive to misalignments as it cannot remove the intensity outliers in the boundaries of the structures. Areas with difficult boundaries such as the thalamus, caudate nucleus and subthalamic nucleus are not significantly different to label fusion. In these areas the weighting of label fusion in the model ensemble is increased through the Lorentzian error norm. This avoids expansion of structures into areas without sufficient intensity differentiation. To test the overall significance of improvement we compared the average Dice scores of each structure obtained with EM against label fusion. EM is significantly better than label fusion with p < 10−8 for global and p < 10−10 for local weighting of the atlases. Another experiment was performed to compare the modelling of the predictive distribution as a model averaging as opposed to the Gaussian Mixture Model formulation alone (essentially the effect of Equation (10)). The results are summarized in Appendix A, Table V. The proposed modelling with the model ensemble results in significantly better overlap for 31 out of the 50 structures and worse for one structure with majority vote weighting (better for 22 structures after Bonferroni correction), and significantly better score in 20 out of the 50 structures with locally-weighted voting (10 after Bonferroni correction). All of the subcortical structures except for the hippocampus and corpus callosum have significantly better (p < 0.05) overlap with respect to the gold standard. Figure 8 allows a visual comparison of the results of the atlas cross-validation presented in this section. B. Evaluation across a large cohort Results of the blinded rating (1-5) of all 234 segmentations of locally-weighted label fusion and the proposed technique are presented in Figure 9. The segmentations produced with the presented method are greatly improved with a mean rating of 4.9 across ages compared to 3.8 judged for label fusion. The majority of the segmentations with the proposed EM implementation are rated with 5 (accurate segmentation) at all scan ages (92% of the cases), in contrast to label fusion mostly voted with 4 (regional boundaries inaccuracies) for brains around term (81.9% of the cases) and declines to 3 (moderate region inaccuracies/alignment problems) for early preterm brains (16% of the cases). The increase is consistent and significant (p < 0.05) for all the ages with a minimum increase in the rating of 0.9 (for the age of 40 weeks) up to 1.8 (for the age of 30 weeks). Example segmentations of neonates at 28 weeks, 31 weeks and 32 weeks GA with the two techniques are shown in Figure 10. Registration between term-equivalent brains to early preterm brains is challenging, especially in the cortical areas due to different degree of cortical complexity. Utilizing the intensity information, the tissue boundaries are refined (see Figure 10.A).

0278-0062 (c) 2013 IEEE. Translations and content mining are permitted for academic research only. Personal use is also permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.

This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TMI.2014.2322280, IEEE Transactions on Medical Imaging 11

5 4 3 2 1

60

50 Number of images

0.9

40

30

20 0.85

10

0

0.8

28

30

32

34 36 38 40 Age at scan (weeks)

42

44

Fig. 9. Rating (5=best) of label fusion (left columns) and the proposed EM segmentation method (right columns) with local weighting across the scan ages of the cohort. Each age at the x axis covers the span [age−1 week, age+ 1 week). Structure lateral ventricle frontal lobe occipital lobe

0.75

0.7

M M V− V EM−Fu −G sio n LW LWMV MM −E −F−EM M us − io LWGM n −EM M

Hippocampus r Hippocampus l Amygdala r Amygdala l Cerebellum r Cerebellum l Brainstem Caudate nucleus l Caudate nucleus r Thalamus l Thalamus r Subthalamic nucleus l Subthalamic nucleus r Lentiform Nucleus l Lentiform Nucleus r Corpus Callosum Lateral Ventricle r Lateral Ventricle l Frontal lobe l Frontal lobe r Parietal lobe l Parietal lobe r Occipital lobe l Occipital lobe r Anterior temporal lobe M r Anterior temporal lobe M l Anterior temporal lobe L r Anterior temporal lobe L l Insula l Insula r Cingulate g A l Cingulate g A r Cingulate g P l Cingulate g P r Superior temporal g middle r Superior temporal g middle l Superior temporal g P l Superior temporal g P r Medial and inferior temporal g A r Medial and inferior temporal g A l Medial and inferior temporal g P l Medial and inferior temporal g P r Gyri parahippocampalis A r Gyri parahippocampalis A l Gyri parahippocampalis P l Gyri parahippocampalis P r Lateral occipitotemporal g A r Lateral occipitotemporal g A l Lateral occipitotemporal g P l Lateral occipitotemporal g P r

Fig. 8. Dice coefficient of leave one out cross-validation with the ALBERTs. Fusion, EM with the GMM modelling and the proposed EM (as model averaging) are compared with both majority vote (MV) and local weighting (LW) of the training atlases.

C. Manually segmented cases Label fusion and EM were compared in the 15 subjects selected at different scan ages. The method produced segmentations that were judged more accurate than locally-weighted label fusion in an average of 21 out of the 26 structures (structures that are defined separately for each hemisphere were merged) in the 15 cases, and equally accurate in the rest of the structures. Appendix A, Table VI summarizes the qualitative results. The volumetric comparison in terms of the Dice score of the manually segmented cases against label fusion and EM is presented in Table III. The proposed technique significantly

MV-Fusion MV-EM LW-Fusion LW-EM 0.690 0.821 0.712 0.806 0.905 0.952 0.934 0.963 0.805 0.870 0.903 0.945 TABLE III D ICE COEFFICIENT OF THE VENTRICLES , FRONTAL AND OCCIPITAL LOBE IN THE 15 MANUALLY SEGMENTED CASES . F USION AND EM ARE COMPARED WITH BOTH MAJORITY VOTE (MV) AND LOCAL WEIGHTING (LW) OF THE TRAINING ATLASES ( BOLD = SIGNIFICANTLY BETTER AT P < 0.05).

outperforms label fusion in the manually segmented regions in both employed weighting schemes. V. D ISCUSSION AND C ONCLUSIONS With improvements in neonatal intensive care, increasing numbers of infants who are born preterm now survive [31], [32]. However, approximately 10% of infants who are born preterm will develop cerebral palsy [33] and up to 50% will develop cognitive and/or behavioural problems in childhood [34], [35]. Quantitative neuroimaging studies are increasingly being used to assess brain growth and development in this vulnerable population [9], [36], [37]. These studies have identified disturbances in the growth of cortical grey matter [38], [39], white matter microstructure [37], [40]–[42] and regional brain growth [9], [36], [37] compared to infants born at term. These quantitative MR measures, obtained at term equivalent age or older, are associated with neurodevelopmental performance [37], [39], [42]. However, approaches to assess regional brain growth prior to term equivalent age are lacking. In this study we have presented an automatic segmentation approach that allows for the first time detailed parcellation of the neonatal brain from the early preterm period to term equivalent age. In contrast to existing methods in the neonatal field, our approach performs whole brain intensity-based segmentation into 50 structures, enabling the detailed assessment of regional brain growth and development in these infants.

0278-0062 (c) 2013 IEEE. Translations and content mining are permitted for academic research only. Personal use is also permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.

This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TMI.2014.2322280, IEEE Transactions on Medical Imaging 12

Fig. 10. Axial slices of T2 weighted MRIs of 28(A), 31(B) and 32(C) weeks GA brains at different scan ages (28, 31, 42 weeks) with labels overlaid. The labels are estimated with label fusion (left images) and the proposed EM scheme (right images). The arrows show areas of improvement with the proposed technique.

In this work we adapt an EM scheme to segment a large number of structures simultaneously for the neonatal brain. The proposed method improves the segmentation utilizing the image intensity and provides more accurate results than state of the art techniques used in the segmentation field (majority-vote, local-weighted fusion). The novelties of the EM implementation compared to previous approaches ( [11], [19], [20]) include the hierarchical intensity modelling of the whole brain and constraints allowing for the restriction of the EM update to regions that are difficult to be differentiated from

the signal intensity. A careful modelling was adopted to subdivide the subject’s brain MRI into distinct intensity regions, that allowed each structure to be modelled with a single Gaussian. Further, a hierarchical approach is introduced to model the WM subdivisions and CGM subdivisions to avoid competition between them in terms of the intensity likelihood. The hierarchical modelling adopted in this study is essential for the segmentation of different structures that belong to the same tissue. Modelling the predictive distribution as a model averaging between label fusion and the GMM prevents subcortical labels from expanding in homogeneous regions. This is essential as in homogeneous areas the intensity distributions of the subcortical structures overlap significantly and this would introduce errors in the segmentations. The model averaging is necessary due to the low contrast in the subcortical space in the neonatal brains, especially in early preterm MR images. Our method incorporates bias field correction, spatial regularization, partial volume correction similar to previous studies in the field. Despite the hierarchical modelling of the proposed model, our algorithm is flexible in the modification of the atlas used for the prior information. Further work will investigate the benefit of integrating different modalities in a multi-modal segmentation of the neonatal brain similar to [43], [44]. To evaluate our technique we compared it with label fusion, one of the most commonly used techniques, using leave-oneout cross-validation in a global and local weighting of the atlases. A qualitative and quantitative evaluation is further performed in a wide range of gestational ages, from early preterm to term-equivalent, as the ALBERTs cover a spectrum of ages around term equivalent age. The propagated labels within our framework are clearly improved in terms of intensity and regions problematic for registration. It should be noted here that automatic quantitative evaluation methods would facilitate a time-efficient and large-scale evaluation of segmentation techniques. Due to the large developmental changes in the developing neonatal brain, evaluation needs to be performed at different ages to validate the robust performance of a segmentation technique. In summary, the proposed framework allows the segmentation of neonatal MR brain images with gestational ages ranging from the early preterm period to term-equivalent age, a time when the developing brain is changing considerably in terms of brain volume, cortical folding and signal intensity on MR images. Furthermore, we have extended an EM scheme previously used to assess CSF, CGM and WM to enable the brain to be divided into multiple brain structures. The proposed algorithm outperforms current approaches. This improvement becomes even more significant in the brain of early preterm infants that vary considerably from the brain at term-equivalent age. The method has been successfully applied to a large cohort of subjects and will enable to assess the impact of aberrant brain growth on subsequent neurodevelopment following preterm birth.

0278-0062 (c) 2013 IEEE. Translations and content mining are permitted for academic research only. Personal use is also permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.

This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TMI.2014.2322280, IEEE Transactions on Medical Imaging 13

A PPENDIX A Region Hippocampus right Hippocampus left Amygdala right Amygdala left Cerebellum right Cerebellum left Brainstem, spans the midline Caudate nucleus left Caudate nucleus right Thalamus left Thalamus right Subthalamic nucleus left Subthalamic nucleus right Lentiform Nucleus left Lentiform Nucleus right Corpus Callosum Lateral Ventricle right Lateral Ventricle left Frontal lobe left Frontal lobe right Parietal lobe left Parietal lobe right Occipital lobe left Occipital lobe right Anterior temporal lobe, medial part right Anterior temporal lobe, medial part left Anterior temporal lobe, lateral part right Anterior temporal lobe, lateral part left Insula left Insula right Cingulate gyrus, anterior part left Cingulate gyrus, anterior part right Cingulate gyrus, posterior part left Cingulate gyrus, posterior part right Superior temporal gyrus, middle part right Superior temporal gyrus, middle part left Superior temporal gyrus, posterior part left Superior temporal gyrus, posterior part right Medial and inferior temporal gyri anterior part right Medial and inferior temporal gyri anterior part left Medial and inferior temporal gyri posterior part left Medial and inferior temporal gyri posterior part right Gyri parahippocampalis et ambiens anterior part right Gyri parahippocampalis et ambiens anterior part left Gyri parahippocampalis et ambiens posterior part left Gyri parahippocampalis et ambiens posterior part right Lateral occipitotemporal gyrus, gyrus fusiformis anterior part right Lateral occipitotemporal gyrus, gyrus fusiformis anterior part left Lateral occipitotemporal gyrus, gyrus fusiformis posterior part left Lateral occipitotemporal gyrus, gyrus fusiformis posterior part right

MV-Fusion 0.767 0.749 0.813 0.822 0.921 0.922 0.920 0.837 0.843 0.895 0.908 0.738 0.744 0.871 0.873 0.706 0.785 0.795 0.917 0.918 0.870 0.866 0.854 0.846 0.815 0.806 0.783 0.788 0.856 0.843 0.766 0.816 0.804 0.799 0.816 0.816 0.682 0.667 0.821 0.814 0.795 0.777 0.782 0.786 0.689 0.701 0.718 0.721 0.686 0.672

MV-EM 0.781 0.764 0.817 0.828 0.925 0.925 0.922 0.842 0.846 0.894 0.907 0.740 0.746 0.873 0.876 0.710 0.822* 0.836* 0.925* 0.929* 0.883 0.877* 0.864* 0.856* 0.837* 0.836* 0.814* 0.815* 0.860 0.848* 0.765 0.819 0.806 0.802 0.831* 0.825* 0.681 0.667 0.838* 0.827* 0.799 0.784* 0.796* 0.801* 0.702* 0.714* 0.734* 0.739* 0.696* 0.679

LW-Fusion 0.788 0.772 0.816 0.824 0.924 0.926 0.917 0.841 0.844 0.895 0.907 0.737 0.746 0.872 0.874 0.717 0.803 0.822 0.926 0.928 0.884 0.876 0.862 0.852 0.826 0.819 0.802 0.797 0.862 0.852 0.770 0.820 0.803 0.798 0.833 0.828 0.686 0.671 0.840 0.830 0.799 0.778 0.797 0.805 0.712 0.721 0.745 0.740 0.702 0.688

LW-EM 0.797 0.783 0.821 0.831 0.929* 0.930* 0.922* 0.845 0.848 0.895 0.908 0.740 0.748 0.874 0.877 0.727* 0.829* 0.846* 0.930* 0.934* 0.887 0.881 0.865 0.856* 0.842* 0.834 0.821* 0.811* 0.866* 0.856* 0.774 0.826* 0.808* 0.803* 0.841* 0.832 0.685 0.671 0.847* 0.836 0.800 0.781 0.810* 0.814* 0.718 0.729* 0.752* 0.749* 0.705 0.695*

TABLE IV D ICE COEFFICIENT OF LEAVE ONE OUT CROSS - VALIDATION WITH THE ALBERT S . F USION AND EM ARE COMPARED WITH BOTH MAJORITY VOTE (MV) AND LOCAL WEIGHTING (LW) OF THE TRAINING ATLASES ( BOLD = SIGNIFICANTLY BETTER AT P < 0.05, * = SIGNIFICANTLY BETTER AFTER B ONFERRONI CORRECTION ).

0278-0062 (c) 2013 IEEE. Translations and content mining are permitted for academic research only. Personal use is also permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.

This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TMI.2014.2322280, IEEE Transactions on Medical Imaging 14 Region Hippocampus right Hippocampus left Amygdala right Amygdala left Cerebellum right Cerebellum left Brainstem, spans the midline Caudate nucleus left Caudate nucleus right Thalamus left Thalamus right Subthalamic nucleus left Subthalamic nucleus right Lentiform Nucleus left Lentiform Nucleus right Corpus Callosum Lateral Ventricle right Lateral Ventricle left Frontal lobe left Frontal lobe right Parietal lobe left Parietal lobe right Occipital lobe left Occipital lobe right Anterior temporal lobe, medial part right Anterior temporal lobe, medial part left Anterior temporal lobe, lateral part right Anterior temporal lobe, lateral part left Insula left Insula right Cingulate gyrus, anterior part left Cingulate gyrus, anterior part right Cingulate gyrus, posterior part left Cingulate gyrus, posterior part right Superior temporal gyrus, middle part right Superior temporal gyrus, middle part left Superior temporal gyrus, posterior part left Superior temporal gyrus, posterior part right Medial and inferior temporal gyri anterior part right Medial and inferior temporal gyri anterior part left Medial and inferior temporal gyri posterior part left Medial and inferior temporal gyri posterior part right Gyri parahippocampalis et ambiens anterior part right Gyri parahippocampalis et ambiens anterior part left Gyri parahippocampalis et ambiens posterior part left Gyri parahippocampalis et ambiens posterior part right Lateral occipitotemporal gyrus, gyrus fusiformis anterior part right Lateral occipitotemporal gyrus, gyrus fusiformis anterior part left Lateral occipitotemporal gyrus, gyrus fusiformis posterior part left Lateral occipitotemporal gyrus, gyrus fusiformis posterior part right

MV-EM-GMM 0.776 0.763 0.800 0.812 0.916 0.915 0.914 0.822 0.824 0.886 0.900 0.715 0.723 0.859 0.865 0.711 0.815 0.828 0.926 0.929 0.881 0.874 0.862 0.853 0.835 0.836 0.810 0.819 0.855 0.844 0.764 0.815 0.803 0.799 0.829 0.824 0.681 0.664 0.831 0.824 0.797 0.782 0.783 0.791 0.699 0.709 0.726 0.732 0.690 0.672

MV-EM 0.781 0.764 0.817* 0.828* 0.925* 0.925* 0.922* 0.842* 0.846* 0.894* 0.907* 0.740* 0.746* 0.873* 0.876* 0.710 0.822* 0.836* 0.925 0.929 0.883 0.877* 0.864 0.856 0.837 0.836 0.814 0.815 0.860* 0.848* 0.765 0.819 0.806 0.802 0.831 0.825 0.681 0.667 0.838* 0.827 0.799 0.784 0.796* 0.801 0.702 0.714 0.734* 0.739 0.696 0.679*

LW-EM-GMM 0.794 0.788 0.814 0.823 0.927 0.929 0.919 0.832 0.834 0.890 0.901 0.720 0.729 0.866 0.871 0.727 0.826 0.844 0.930 0.934 0.887 0.881 0.864 0.856 0.842 0.833 0.821 0.811 0.864 0.852 0.773 0.825 0.807 0.803 0.841 0.832 0.684 0.670 0.847 0.836 0.799 0.781 0.807 0.813 0.717 0.728 0.751 0.749 0.705 0.694

LW-EM 0.797 0.783 0.821 0.831 0.929 0.930 0.922 0.845* 0.848* 0.895* 0.908* 0.740* 0.748* 0.874* 0.877* 0.727 0.829 0.846 0.930 0.934 0.887 0.881 0.865 0.856 0.842 0.834 0.821 0.811 0.866* 0.856* 0.774 0.826 0.808 0.803 0.841 0.832 0.685 0.671 0.847 0.836 0.800 0.781 0.810 0.814 0.718 0.729 0.752 0.749 0.705 0.695

TABLE V D ICE COEFFICIENT OF LEAVE ONE OUT CROSS - VALIDATION WITH THE ALBERT S . T HE PROPOSED EM AS AN ENSEMBLE OF MODELS IS COMPARED WITH THE EM WITH GMM MODELLING ALONE WITH BOTH MAJORITY VOTE (MV) AND LOCAL WEIGHTING (LW) OF THE TRAINING ATLASES ( BOLD = SIGNIFICANTLY BETTER AT P < 0.05, * = SIGNIFICANTLY BETTER AFTER B ONFERRONI CORRECTION ).

Region Hippocampus Amygdala Cerebellum Brainstem, spans the midline Caudate nucleus Thalamus Subthalamic nucleus Lentiform Nucleus Corpus Callosum Lateral Ventricle Frontal lobe Parietal lobe Occipital lobe Insula Medial and inferior temporal gyri, anterior part Medial and inferior temporal gyri posterior part Cingulate gyrus, anterior part Cingulate gyrus, posterior part Gyri parahippocampalis et ambiens anterior part Gyri parahippocampalis et ambiens posterior part Lateral occipitotemporal gyrus, gyrus fusiformis anterior part Lateral occipitotemporal gyrus, gyrus fusiformis posterior part Superior temporal gyrus, middle part Superior temporal gyrus, posterior part Anterior temporal lobe, medial part Anterior temporal lobe, lateral part

26+5 o EM EM EM EM EM o o o EM EM EM EM EM EM EM o EM EM o EM EM EM EM EM EM

27+6 EM o EM EM EM EM o EM EM EM EM EM EM o EM EM EM EM EM EM EM EM EM EM EM EM

28+2 EM EM EM EM o o EM EM EM EM EM EM EM EM EM EM o EM o EM EM o EM EM EM EM

29 EM EM EM EM o o o EM EM EM EM EM EM EM EM EM o o EM EM EM EM EM EM EM EM

30 o EM EM EM EM EM EM EM EM EM EM EM EM EM EM EM EM EM EM EM EM EM EM EM EM EM

31+1 o o EM EM EM EM o EM EM EM EM EM EM o EM EM EM EM EM EM EM EM EM EM EM EM

32 EM EM EM EM EM EM o EM EM EM EM EM EM EM EM EM EM EM EM EM EM EM EM EM EM EM

33 EM o EM EM EM o o o o EM EM EM EM o EM EM o o o EM EM EM EM EM EM EM

34 o o EM EM EM EM o EM EM EM EM EM EM EM EM EM EM EM EM EM EM EM EM EM EM EM

35+2 o o EM EM o o o o o EM EM EM EM EM EM EM EM EM EM EM EM EM EM EM EM EM

36 EM EM EM EM o EM o EM EM EM EM EM EM EM EM EM EM EM EM EM EM EM EM EM EM EM

37+2 o o EM EM EM o o EM EM EM EM EM EM EM EM EM EM EM o o o EM EM EM EM EM

39+6 o o EM EM o EM o o o EM EM EM EM EM EM EM EM EM o o EM EM EM EM EM EM

42 EM EM EM EM EM EM o EM EM EM EM EM EM EM EM EM EM EM EM EM EM EM EM EM EM EM

44+2 o o EM EM o o o EM o EM EM EM EM EM o o EM EM EM EM EM EM EM EM EM EM

TABLE VI Q UALITITATIVE COMPARISON OF FUSION AND EM WITH LOCAL WEIGHTING OF THE TRAINING ATLASES . T HE BEST SEGMENTATION RESULT FOR EACH LABEL IS PRESENTED ACCORDING TO THE OBSERVER . L ABELS WHERE THE TWO SEGMENTATION TECHNIQUES ARE EQUALLY ACCURATE ARE DENOTED WITH ’ O ’.

0278-0062 (c) 2013 IEEE. Translations and content mining are permitted for academic research only. Personal use is also permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.

This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TMI.2014.2322280, IEEE Transactions on Medical Imaging 15

ACKNOWLEDGMENT The authors acknowledge financial support from: the Department of Health via the National Institute for Health Research (NIHR) Comprehensive Biomedical Research Centre award to Guys & St Thomas NHS Foundation Trust in partnership with Kings College London and Kings College Hospital NHS Foundation Trust; Imperial College Healthcare Comprehensive Biomedical Research Centre Funding Scheme; the Medical Research Council Foundation; the Medical Research Council (UK); Action Medical Research. The work was supported by the Department of Perinatal Imaging & Health at Kings College London. R EFERENCES [1] M. J. Cardoso, A. Melbourne, G. S. Kendall, M. Modat, C. F. Hagmann, N. J. Robertson, N. Marlow, and S. Ourselin, “Adaptive neonate brain segmentation,” Medical image computing and computer-assisted intervention (MICCAI), vol. 14, no. Pt 3, pp. 378–386, 2011, PMID: 22003722. [2] M. J. Cardoso, A. Melbourne, G. S. Kendall, M. Modat, N. J. Robertson, N. Marlow, and S. Ourselin, “AdaPT: an adaptive preterm segmentation algorithm for neonatal brain MRI,” NeuroImage, Aug. 2012, PMID: 22906793. [3] I. S. Gousias, D. Rueckert, R. A. Heckemann, L. E. Dyet, J. P. Boardman, A. D. Edwards, and A. Hammers, “Automatic segmentation of brain MRIs of 2-year-olds into 83 regions of interest,” NeuroImage, vol. 40, no. 2, pp. 672–684, Apr. 2008, PMID: 18234511. [4] L. Gui, R. Lisowski, T. Faundez, P. S. Hppi, F. Lazeyras, and M. Kocher, “Morphology-driven automatic segmentation of MR images of the neonatal brain,” Medical image analysis, Jul. 2012, PMID: 22921305. [5] M. Prastawa, J. H. Gilmore, W. Lin, and G. Gerig, “Automatic segmentation of MR images of the developing newborn brain,” Medical image analysis, vol. 9, no. 5, pp. 457–466, Oct. 2005, PMID: 16019252. [6] Z. Song, S. P. Awate, D. J. Licht, and J. C. Gee, “Clinical neonatal brain MRI segmentation using adaptive nonparametric data models and intensity-based markov priors,” Medical image computing and computer-assisted intervention (MICCAI), vol. 10, no. Pt 1, pp. 883–890, 2007, PMID: 18051142. [7] N. I. Weisenfeld and S. K. Warfield, “Automatic segmentation of newborn brain MRI,” NeuroImage, vol. 47, no. 2, pp. 564–572, Aug. 2009, PMID: 19409502. [8] H. Xue, L. Srinivasan, S. Jiang, M. Rutherford, A. D. Edwards, D. Rueckert, and J. V. Hajnal, “Automatic segmentation and reconstruction of the cortex from neonatal MRI,” NeuroImage, vol. 38, no. 3, pp. 461–477, Nov. 2007, PMID: 17888685. [9] I. S. Gousias, A. D. Edwards, M. A. Rutherford, S. J. Counsell, J. V. Hajnal, D. Rueckert, and A. Hammers, “Magnetic resonance imaging of the newborn brain: Manual segmentation of labelled atlases in term-born and preterm infants,” NeuroImage, vol. 62, no. 3, pp. 1499–1509, Sep. 2012. [10] I. S. Gousias, A. Hammers, S. J. Counsell, L. Srinivasan, M. A. Rutherford, R. A. Heckemann, J. V. Hajnal, D. Rueckert, and A. D. Edwards, “Magnetic resonance imaging of the newborn brain: Automatic segmentation of brain images into 50 anatomical regions,” PLoS ONE, vol. 8, no. 4, p. e59990, Apr. 2013. [11] K. Van Leemput, F. Maes, D. Vandermeulen, and P. Suetens, “Automated model-based tissue classification of MR images of the brain,” IEEE transactions on medical imaging, vol. 18, no. 10, pp. 897–908, Oct. 1999, PMID: 10628949. [12] S. K. Warfield, K. H. Zou, and W. M. Wells, “Simultaneous truth and performance level estimation (STAPLE): an algorithm for the validation of image segmentation,” IEEE transactions on medical imaging, vol. 23, no. 7, pp. 903–921, Jul. 2004. [13] N. Shiee, P.-L. Bazin, J. L. Cuzzocreo, A. Blitz, and D. L. Pham, “Segmentation of brain images using adaptive atlases with application to ventriculomegaly,” Information Processing in Medical Imaging, vol. 22, pp. 1–12, 2011, PMID: 21761641. [14] R. A. Heckemann, J. V. Hajnal, P. Aljabar, D. Rueckert, and A. Hammers, “Automatic anatomical brain MRI segmentation combining label propagation and decision fusion,” NeuroImage, vol. 33, no. 1, pp. 115–126, Oct. 2006, PMID: 16860573. [15] T. Rohlfing, R. Brandt, R. Menzel, and J. Maurer, Calvin R, “Evaluation of atlas selection strategies for atlas-based image segmentation with application to confocal microscopy images of bee brains,” NeuroImage, vol. 21, no. 4, pp. 1428–1442, Apr. 2004, PMID: 15050568. [16] P. Aljabar, R. A. Heckemann, A. Hammers, J. V. Hajnal, and D. Rueckert, “Multi-atlas based segmentation of brain images: atlas selection and its effect on accuracy,” NeuroImage, vol. 46, no. 3, pp. 726–738, Jul. 2009, PMID: 19245840. [17] X. Artaechevarria, A. Munoz-Barrutia, and C. Ortiz-de Solorzano, “Combination strategies in multi-atlas image segmentation: application to brain MR data,” IEEE transactions on medical imaging, vol. 28, no. 8, pp. 1266–1277, Aug. 2009, PMID: 19228554. [18] K. O. Babalola, B. Patenaude, P. Aljabar, J. Schnabel, D. Kennedy, W. Crum, S. Smith, T. Cootes, M. Jenkinson, and D. Rueckert, “An evaluation of four automatic methods of segmenting the subcortical structures in the brain,” NeuroImage, vol. 47, no. 4, pp. 1435–1447, Oct. 2009, PMID: 19463960. [19] C. Ledig, R. Wolz, P. Aljabar, J. Lotjonen, R. Heckemann, A. Hammers, and D. Rueckert, “Multi-class brain segmentation using atlas propagation and EM-based refinement,” in 2012 IEEE International Symposium on Biomedical Imaging (ISBI), May 2012, pp. 896 –899. [20] J. M. L¨otj¨onen, R. Wolz, J. R. Koikkalainen, L. Thurfjell, G. Waldemar, H. Soininen, and D. Rueckert, “Fast and robust multi-atlas segmentation of brain magnetic resonance images,” NeuroImage, vol. 49, no. 3, pp. 2352–2365, Feb. 2010, PMID: 19857578. [21] S. M. Smith, “Fast robust automated brain extraction,” Human brain mapping, vol. 17, no. 3, pp. 143–155, Nov. 2002, PMID: 12391568. [22] N. J. Tustison, B. B. Avants, P. A. Cook, Y. Zheng, A. Egan, P. A. Yushkevich, and J. C. Gee, “N4ITK: improved n3 bias correction,” IEEE transactions on medical imaging, vol. 29, no. 6, pp. 1310–1320, Jun. 2010, PMID: 20378467. [23] D. Rueckert, L. I. Sonoda, C. Hayes, D. L. Hill, M. O. Leach, and D. J. Hawkes, “Nonrigid registration using free-form deformations: application to breast MR images,” IEEE transactions on medical imaging, vol. 18, no. 8, pp. 712–721, Aug. 1999, PMID: 10534053. [24] I. Gousias, A. Hammers, R. Heckemann, S. Counsell, L. Dyet, J. Boardman, A. Edwards, and D. Rueckert, “Atlas selection strategy for automatic segmentation of pediatric brain MRIs into 83 ROIs,” in 2010 IEEE International Conference on Imaging Systems and Techniques (IST), 2010, pp. 290–293. [25] J. Macqueen, “Some methods of classification and analysis of multivariate observations,” in Proceedings of the Fifth Berkeley Symposium on Mathematical Statistics and Probability, 1967, pp. 281–297. [26] A. Serag, P. Aljabar, G. Ball, S. J. Counsell, J. P. Boardman, M. A. Rutherford, A. D. Edwards, J. V. Hajnal, and D. Rueckert, “Construction of a consistent high-definition spatio-temporal atlas of the developing brain using adaptive kernel regression,” NeuroImage, vol. 59, no. 3, pp. 2255–2265, Feb. 2012, PMID: 21985910. [27] S. J. Counsell and M. A. Rutherford, “Magnetic resonance imaging of the newborn brain,” Current Paediatrics, vol. 12, no. 5, pp. 401–413, Oct. 2002. [28] K. Van Leemput, A. Bakkour, T. Benner, G. Wiggins, L. L. Wald, J. Augustinack, B. C. Dickerson, P. Golland, and B. Fischl, “Automated segmentation of hippocampal subfields from ultra-high resolution in vivo MRI,” Hippocampus, vol. 19, no. 6, pp. 549–557, Jun. 2009, PMID: 19405131. [29] K. M. Pohl, S. Bouix, R. Kikinis, and W. E. L. Grimson, “Anatomical guided segmentation with non-stationary tissue class distributions in an expectationmaximization framework,” in 2004 IEEE International Symposium on Biomedical Imaging (ISBI), 2004, p. 8184. [30] M. J. Black, G. Sapiro, D. Marimont, and D. Heeger, Robust Anisotropic Diffusion, 1998. [31] R. M. Schwartz, A. M. Luby, J. W. Scanlon, and R. J. Kellogg, “Effect of surfactant on morbidity, mortality, and resource use in newborn infants weighing 500 to 1500 g,” New England Journal of Medicine, vol. 330, no. 21, pp. 1476–1480, May 1994, PMID: 8164699.

0278-0062 (c) 2013 IEEE. Translations and content mining are permitted for academic research only. Personal use is also permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.

This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TMI.2014.2322280, IEEE Transactions on Medical Imaging 16

[32] D. Wilson-Costello, H. Friedman, N. Minich, A. A. Fanaroff, and M. Hack, “Improved survival rates with increased neurodevelopmental disability for extremely low birth weight infants in the 1990s,” PEDIATRICS, vol. 115, no. 4, pp. 997–1003, Apr. 2005. [33] M. Hack and A. A. Fanaroff, “Outcomes of children of extremely low birthweight and gestational age in the 1990s,” Seminars in Neonatology, vol. 5, no. 2, pp. 89–106, May 2000. [34] M. Delobel-Ayoub, C. Arnaud, M. White-Koning, C. Casper, V. Pierrat, M. Garel, A. Burguet, J.-C. Roze, J. Matis, J.-C. Picaud, M. Kaminski, and B. Larroque, “Behavioral problems and cognitive performance at 5 years of age after very preterm birth: The EPIPAGE study,” Pediatrics, vol. 123, no. 6, pp. 1485–1492, Jun. 2009. [35] N. Marlow, D. Wolke, M. A. Bracewell, and M. Samara, “Neurologic and developmental disability at six years of age after extremely preterm birth,” New England Journal of Medicine, vol. 352, no. 1, pp. 9–19, Jan. 2005, PMID: 15635108. [36] G. Ball, J. P. Boardman, D. Rueckert, P. Aljabar, T. Arichi, N. Merchant, I. S. Gousias, A. D. Edwards, and S. J. Counsell, “The effect of preterm birth on thalamic and cortical development,” Cerebral cortex (New York, N.Y.: 1991), vol. 22, no. 5, pp. 1016–1024, May 2012, PMID: 21772018. [37] J. P. Boardman, C. Craven, S. Valappil, S. J. Counsell, L. E. Dyet, D. Rueckert, P. Aljabar, M. A. Rutherford, A. T. M. Chew, J. M. Allsop, F. Cowan, and A. D. Edwards, “A common neonatal image phenotype predicts adverse neurodevelopmental outcome in children born preterm,” NeuroImage, vol. 52, no. 2, pp. 409–414, Aug. 2010, PMID: 20451627. [38] O. Kapellou, S. J. Counsell, N. Kennea, L. Dyet, N. Saeed, J. Stark, E. Maalouf, P. Duggan, M. Ajayi-Obe, J. Hajnal, J. M. Allsop, J. Boardman, M. A. Rutherford, F. Cowan, and A. D. Edwards, “Abnormal cortical development after premature birth shown by altered allometric scaling of brain growth,” PLoS Medicine, vol. 3, no. 8, p. e265, Aug. 2006, PMID: 16866579. [39] R. Rathbone, S. J. Counsell, O. Kapellou, L. Dyet, N. Kennea, J. Hajnal, J. M. Allsop, F. Cowan, and A. D. Edwards, “Perinatal cortical growth and childhood neurocognitive abilities,” Neurology, vol. 77, no. 16, pp. 1510–1517, Oct. 2011, PMID: 21998316. [40] M. Anjari, L. Srinivasan, J. M. Allsop, J. V. Hajnal, M. A. Rutherford, A. D. Edwards, and S. J. Counsell, “Diffusion tensor imaging with tract-based spatial statistics reveals local white matter abnormalities in preterm infants,” NeuroImage, vol. 35, no. 3, pp. 1021–1027, Apr. 2007, PMID: 17344066. [41] G. Ball, S. J. Counsell, M. Anjari, N. Merchant, T. Arichi, V. Doria, M. A. Rutherford, A. D. Edwards, D. Rueckert, and J. P. Boardman, “An optimised tract-based spatial statistics protocol for neonates: applications to prematurity and chronic lung disease,” NeuroImage, vol. 53, no. 1, pp. 94–102, Oct. 2010, PMID: 20510375. [42] B. J. M. van Kooij, L. S. de Vries, G. Ball, I. C. van Haastert, M. J. N. L. Benders, F. Groenendaal, and S. J. Counsell, “Neonatal tract-based spatial statistics findings and outcome in preterm infants,” American Journal of Neuroradiology, vol. 33, no. 1, pp. 188–194, Jan. 2012, PMID: 21998101. [43] B. B. Avants, N. J. Tustison, J. Wu, P. A. Cook, and J. C. Gee, “An open source multivariate framework for n-tissue segmentation with evaluation on public data,” Neuroinformatics, vol. 9, no. 4, pp. 381–400, Dec. 2011. [44] L. Wang, F. Shi, P.-T. Yap, J. H. Gilmore, W. Lin, and D. Shen, “4D multi-modality tissue segmentation of serial infant images,” PLoS ONE, vol. 7, no. 9, p. e44596, 2012.

0278-0062 (c) 2013 IEEE. Translations and content mining are permitted for academic research only. Personal use is also permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.

Automatic whole brain MRI segmentation of the developing neonatal brain.

Magnetic resonance (MR) imaging is increasingly being used to assess brain growth and development in infants. Such studies are often based on quantita...
7MB Sizes 24 Downloads 4 Views