The neonatal databaseWe evaluated the performances of our proposed method using the first release of the developing human connectome project (dHCP) database which is publicly available. The image acquisition was done using 3 T Philips Achieva including 32 channel phased array head coil24. By anatomical imaging of 40 non-sedated term neonates (gestational age of 37–44 weeks), the T2-weighted MR images were acquired in sagittal and axial slice stacks within plane resolution 0.8 × 0.8 mm2 and 1.6 mm slices overlapped by 0.8 mm. Other parameters are: TR = 12 s; TE = 156 ms; SENSE factor: axial = 2.11, sagittal = 2.58. Furthermore, we applied the proposed method to segment brain tissues of 40 neonates with the same age range at acquisition time from dHCP data release four (dHCPR4)25. This release includes images of 783 neonates with age range 23–44 weeks (gestational age) which have been acquired with the same imaging parameters of the first release.The motion corrected, reconstructed T2-weighted images were bias corrected using the N4 algorithm26. Manual atlases, annotated by a neuroanatomist expert, were registered to the T2-weighted images, and their labels were fused in the subject space to derive structure priors. Subsequently, segmentation was performed employing an EM scheme27 to extract 87 regions, which were further consolidated into 8 tissues, including cortical gray matter (cGM), deep gray matter (dGM), WM, CSF, ventricle, cerebellum, brainstem, and hippocampus. The GM was created combining the cGM and dGM tissues. GM, WM, and CSF masks were then reviewed by neuroanatomists and any mis-segmentation was subsequently manually corrected in all scans and utilized as the “ground truth”. These pre-processed T2-weighted images and their corresponding “ground truth” images were then used to assess the accuracy of the proposed methods proposed in this study.Proposed methodThe proposed neonatal brain tissue segmentation method consisted of three stages: 1) atlas selection, 2) subject-specific atlas creation using RF classifier, and 3) brain tissue segmentation using the EM-MRF method.Stage 1: atlas selectionIn Stage 1, the uniform atlas selection (UAS) algorithm28 was used to uniformly select a subset of \(K\) (\(K \le M\)) images from the entire library as atlases to reduce the size of the training library.The affine registration (Implemented in SPM12 toolbox29) was used to transform each image \(I_{{tr_{i} }}\) and its corresponding label image \(L_{{tr_{i} }}\) into the common coordinate space as defined by neonatal brain atlas template \(T_{avg}\)30. Each label image \(L_{{tr_{i} }}\) consisted of \(c = \{ 1,2,…,C\}\) brain tissues, i.e. GM, WM, and CSF. The library set \(D = \{ D_{i} \left| {i = 1,…,M\} } \right.\) was then created using \(D_{i}\) containing the affine aligned image \(I_{{tr_{i} }}^{A}\) and its corresponding label image \(L_{{tr_{i} }}^{A}\). As explained in Noorizadeh et al.30, the atlas set \(S = \{ S_{k} \left| {k = 1,…,K;S_{k} \in D} \right.\}\) was selected using the UAS algorithm. Each atlas \(S_{k}\) consisted of the \(I_{{trs_{k} }}^{A}\) and label image \(L_{{trs_{k} }}^{A}\).Stage 2: subject specific atlas creation using RF classifierThe goal of the second stage was to make a subject-specific atlas for each target image \(I_{ta}\). To achieve this goal as illustrated in Fig. 1, \(I_{ta}\) was first transformed to the template space using the affine registration. Then, the subject-specific atlas of the aligned image (\(I_{ta}^{A}\)) was created using developed window-based segmentation method as follows:Figure 1Schematic diagram of Stage 3 used to create a subject-specific atlas for each target image. \(T_{avg}\), neonatal template; \(I_{ta}\), target image; \(I_{ta}^{A}\), affine aligned target image; (\(I_{{trs_{k} }}^{A}\), \(L_{{tr_{i} }}^{A}\)), affine aligned atlas image and its corresponding label image; (\(I_{{trs_{k} }}^{N}\), \(L_{{tr_{i} }}^{N}\)), nonlinearly aligned atlas image and its corresponding label image; \(sw_{j}\), \(j^{th}\) non-overlapping window; \(j\), index of windows; \(R\), number of windows; \(f_{k}^{{sw_{j} }} (.)\), feature vector of \(k^{th}\) atlas in \(sw_{j}\); \(F_{tr}^{{sw_{j} }}\), atlas feature matrix; \(l_{k}^{{sw_{j} }} (.)\), label of \(k^{th}\) atlas in \(sw_{j}\); \(L_{tr}^{{sw_{j} }}\), atlas label vector; \(f_{ta}^{{sw_{j} }} (.)\), target feature vector in \(sw_{j}\); \(F_{ta}^{{sw_{j} }}\), atlas feature matrix; \(SSA_{ta}\), subject-specific atlas.Step 1: Transform the atlas images into the target image space: In order to reduce the computational time, only the images in the atlas set \(S\) were used to create the subject-specific atlases. The affine aligned images \(I_{{trs_{k} }}^{A}\) and their corresponding label images \(L_{{trs_{k} }}^{A}\) were non-linearly transformed to the \(I_{ta}^{A}\) space using the method proposed by Ashburner et al. and implemented in SPM12 toolbox29 to generate the transformed images (\(I_{{trs_{k} }}^{N}\)) and label images (\(L_{{trs_{k} }}^{N}\)).Step 2: Non-overlapping windowing: In neonatal MR images, WM contains high water content in some regions (e.g. near to ventricles) and therefore its signal intensity is similar to that of CSF in these regions31. The classification methods based on classifiers trained using global information have shown to perform poor on tissues with similar intensities (like CSF or WM) due to the ambiguous tissue intensity distribution14,15,32. To overcome this weakness, local classification methods have been introduced33. With a relatively high computational load, the local strategy has been shown to provide a higher segmentation accuracy by training independent classifiers in a voxel-wise manner. We implemented a segmentation method adopted from Serag et al.13 based on the trade-off between the performance improvement and computational cost. In this method, the images domain \(\Omega\) were split into \(w \times w \times w\) non-overlapping windows \(sw_{j}\) (\(j = 1,…,R;sw_{j} \in \Omega\)). A local classifier was then trained on each window \(sw_{j}\) to reduce misclassification caused by local intensity changes. The label assignment was done simultaneously by each trained classifier for all the \(w^{3}\) voxels in each window \(sw_{j}\) to reduce the computational cost.Step 3: Feature extraction from atlas images: We used different gradient- and intensity-based features to improve the segmentation accuracy by using attributes providing higher order descriptions of images especially around boundaries as well as the image characteristics at the voxel-wise level.For each voxel \(v{}_{l}\)(\(l = 1,…,w^{3} ;v_{l} \in sw_{j}\)) in the window \(sw_{j}\) of the atlas images \(I_{{trs_{k} }}^{N}\), a feature vector \(f_{k}^{{sw_{j} }} (v_{l} )\) containing 10 gradient- and intensity-based features (gray scale intensity, first and second order derivatives in \((x,y,z)\) direction, gradient magnitude, azimuth, and zenith angles) was extracted (see supplementary Sect. 1)13.An atlas feature matrix \(F_{tr}^{{sw_{j} }} = \{ f_{k}^{{sw_{j} }} (v_{l} )\left| {k = 1,…,K;l = 1,…,w^{3} } \right.\}\) was then constructed with columns representing the feature vectors \(f_{k}^{{sw_{j} }} (v_{l} )\) of all \(w^{3}\) voxels belonging to the window \(sw_{j}\) of all the atlas images. The corresponding label, extracted from the label image \(L_{{trs_{k} }}^{N}\) for each voxel \(v_{l}\)(\(l_{k}^{{sw_{j} }} (v_{l} )\)) in the window \(sw_{j}\), was also included in the atlas label vector \(L_{tr}^{{sw_{j} }}\) as:$$L_{tr}^{{sw_{j} }} = \{ l_{k}^{{sw_{j} }} (v_{l} ) = L_{{trs_{k} }}^{N} (v_{l} )\left| {k = 1,…,K;l = 1,…,w^{3} } \right.\}$$
(1)
The atlas feature matrix \(F_{tr}^{{sw_{j} }}\) and the atlas label vector \(L_{tr}^{{sw_{j} }}\) were then used to train the local classifier for the \(j^{th}\) window \(sw_{j}\).Step 4: Training RF classifiers: An RF classifier with an ensemble of \(tre = [1,TRE]\) decision trees was used to create a subject-specific atlas for each target image. In this classifier, each decision tree contains two kinds of nodes: non-leaf (internal) and leaf nodes. By using the stored split function at each internal node, the incoming input data was sent to the left or right child node. The final prediction of each decision tree was also stored in each leaf node. During training, a weak class predictor \(p_{tre} (L_{tr}^{{sw_{j} }} \left| {F_{tr}^{{sw_{j} }} } \right.)\) was trained in each decision tree using \(F_{tr}^{{sw_{j} }}\) and \(L_{tr}^{{sw_{j} }}\). To this end, a subset of all possible features was first randomly sampled. The selected features were then compared with the split function in each internal node to assign them to the left or right child node. Therefore, the purpose of the training was to optimize both the selected features and split functions by maximizing the information transfer at each internal node. By splitting other internal nodes, the tree continued growing. The tree stopped growing at a specific depth (\(D^{TRE}\)) or when the leaf node contained less than a certain number of training samples (\(S_{\min }^{TRE}\)). By repeating this procedure for all other decision trees, the RF classifier was trained on \(j^{th}\) window.Step 5: Feature extraction and label assignment in target image: For each voxel \(v{}_{l}\)(\(l = 1,…,w^{3} ;v_{l} \in sw_{j}\)) in the window \(sw_{j}\) of the target image \(I_{ta}^{A}\), a feature vector \(f_{ta}^{{sw_{j} }} (v_{l} )\) containing the same features as those extracted from the atlas images (see supplementary Sect. 1) was used to construct the target feature matrix \(F_{ta}^{{sw_{j} }} = \{ f_{ta}^{{sw_{j} }} (v_{l} )\left| {l = 1,…,w^{3} } \right.\}\) Each voxel \(v{}_{l}\) was then pushed through each trained decision tree (\(tre = [1,TRE]\)). Once it passed all internal nodes and arrived at a leaf node, the leaf node distribution was used to determine the tissue probability (\(p_{ta}^{{sw_{j} ,tre}} (L\left| {v_{l} } \right.)\)) of voxel \(v{}_{l}\) in the \(tre^{th}\) tree. By averaging \(p_{ta}^{{sw_{j} ,tre}} (L\left| {v_{l} } \right.)\) over all decision trees, the final tissue probability of each voxel \(v{}_{l}\) in the window \(sw_{j}\) was computed as:$$p_{ta}^{{sw_{j} }} (L\left| {v_{l} } \right.) = \frac{1}{TRE}\sum\limits_{tre = 1}^{TRE} {p_{ta}^{{sw_{j} ,tre}} (L\left| {v_{l} } \right.)}$$
(2)
By repeating steps 3–5 for other windows, the initial estimation of GM, WM, and CSF was generated and used as the subject-specific atlas (\(SSA_{ta}\)) for the target image.Stage 3: tissue segmentation using EM-MRFIn the final stage of the proposed method, we performed the brain tissue segmentation (\(L_{ta} = \{ l_{c} |c = 1,…,C\}\)) in the target image using the unified segmentation algorithm based on the EM-MRF algorithm34 as implemented in SPM12 toolbox29. In this algorithm, the spatial normalization, segmentation, and bias correction were done using the generative model, which is based on a mixture of Gaussians (For more details see supplementary Sect. 2).Parameter selectionThe performance of the RF classifiers used in the subject-specific atlas creation stage depended on four parameters: number of atlases (\(K\)), tree depth (\(D^{TRE}\)), minimum number of training samples in leaf nodes (\(S_{\min }^{TRE}\)), and number of trees (\(TRE\)). Since the atlas selection strategy used in this study for brain tissue segmentation was similar to the method described in Noorizadeh et al. for brain extraction30, in order to have a trade-off between the accuracy and the computational complexity, the number of atlases (\(K\)) was set to 10. The RF parameters such as the tree depth (\(D^{TRE}\)), and the minimum number of training samples in leaf nodes (\(S_{\min }^{TRE}\)) were set to 50 and 8, respectively, according to Wang et al.14. However, the effect of the number of trees (\(TRE\)) was evaluated by varying \(TRE\) in the range of [5–30] with an increment of five. The accuracy of the proposed tissue segmentation depended also on the windows size (\(w \times w \times w\)). We changed \(w\) in the range of [3–9] with an increment of two to investigate the impact of \(w\) on the performance improvement and computational time.Evaluation measuresWe used the \(DSC\)35, Jaccard coefficient (\(JC\))36, and Conformity coefficient (\(CC\))37 to measure the similarity between automatically segmented GM, WM, and CSF and manually delineated tissues as follows:$$DSC = \frac{2TP}{{2TP + FP + FN}}$$
(3)
$$JC = \frac{TP}{{TP + FP + FN}}$$
(4)
$$CC = 1 – \frac{FP + FN}{{TP}}$$
(5)
where \(TP\), \(FP\), and \(FN\) denote true positive, false positive, and false negative, respectively. In addition, the sensitivity (\(SE\)) and specificity (\(SP\)) are computed to assess the amount of correctly identified \(TP\) and \(TN\) as:$$SE = \frac{TP}{{TP + FN}}$$
(6)
$$SP = \frac{TN}{{TN + FP}}$$
(7)
where \(TN\) denotes true negatives. In addition to the similarity metrics, we computed the overall accuracy (\(OAC\)) across all tissues38:$$OAC = \frac{TP + TN}{{TP + TN + FP + FN}}$$
(8)
The “leave-one-out cross-validation” (LOOCV) strategy was used to evaluate the performance of the proposed tissue segmentation method on T2-weighted MR images from the first and fourth release of dHCP database. In each round, one image was selected as the test image while the remaining 39 images were considered as training images. By averaging the evaluation measures over the 40 runs, the performance of the proposed method was computed. All programs were developed in MATLAB R2019b software and run on a standard single-threaded PC comprising an Intel ® Core i7 (6700HQ CPU @ 2.60 GHz) with 12 GB RAM.The accuracy of the segmented tissues by the proposed method was compared with SEGMA, which is a multi-atlas segmentation method based on multi-class RF classifiers13 as well as with the EM-MRF algorithm29,34 with three different atlases used as prior knowledge:
AtlasI: the probabilistic atlas created by Murgasova et al.39 based on MR images of 153 neonates born at 40 weeks of gestational age (GA).
AtlasII: the probabilistic atlas created by averaging the aligned label images of 40 neonates from the first release of dHCP database (Fig. 2a).
AtlasIII: the probabilistic atlas created by averaging the aligned label images of \(K\) (= 10) neonates, selected as atlases by the UAS algorithm28 for each target image (Fig. 2b).
Figure 2Illustration of the probabilistic atlases used as prior knowledge in the EM-MRF algorithm. (a) AtlasII, and (b) AtlasIII created for a selected subject in axial, coronal, and sagittal views using dHCP first release.To achieve an optimal trade-off between the segmentation performance and computational complexity, we set the SEGMA parameters, the number of atlases (\(K\)), tree depth (\(D^{TRE}\)), minimum number of training samples in leaf nodes (\(S_{\min }^{TRE}\)), and number of trees (\(TRE\)) to 10, 50, 8, and 20, respectively.We further investigated significant differences between the performance of the proposed method and that of SEGMA, EM-MRF (AtlasI), EM-MRF (AtlasII), and EM-MRF (AtlasIII) using two sample t-tests (p < 0.001).