Anchor-Based Segmentation of Periventricular White Matter for Neonatal HIE

The segmentation and quantitative features of periventricular white matter injury in magnetic resonance imaging are helpful for determining an accurate diagnosis of neonatal hypoxia ischemia encephalopathy (HIE) in clinical practice, but they are tedious tasks. In this paper, we present an interactive segmentation method for clinical magnetic resonance imaging (MRI) using sample estimation at the subvoxel level to establish a quantitative diagnosis. The coarse-grained voxels in clinical MRI are divided into subvoxels, and the intensities of subvoxels are adjusted using gradient adaptive filtering. Because the cerebrospinal fluid is more obvious in grayscale than the white matter and gray matter on T1-weighted images (T1WI), T2-weighted images (T2WI) and fluid-attenuated inversion recovery (FLAIR) images, and its population variance follows a Chi-square distribution, we segment the periventricular region as an anchor by estimating the parameters of the selected data and expand the anchor to obtain periventricular white matter of various widths. After aligning the segmented periventricular white matter in various modalities, we compute the radiomics features of regions of interest for the quantitative diagnosis of HIE. Current clinical practice confirms that our approach is effective and satisfies the clinical diagnostic requirements for neonatal HIE.


I. INTRODUCTION
Due to the lack of quantified medical imaging features, the misdiagnosis rate of neonatal hypoxia ischemia encephalopathy (HIE) is high in clinical practice. Neonatal HIE is closely related to periventricular white matter injury (WMI). WMI refers to the leuko-encephalopathy caused by hypoxia, ischemia or infection in infants born at 24 to 35 weeks. The injury is mainly located in the periventricular region and callosum and includes focal necrosis of deep white matter and a diffuse focus of central white matter. Diffusion-weighted imaging (DWI) is the most sensitive sequence to detect early ischemic lesions [1], and T2-weighted images (T2WI) and fluid-attenuated inversion recovery (FLAIR) images may reflect early WMI. Thus, the segmentation and measurement of periventricular white matter in various magnetic resonance imaging (MRI) modalities is a crucial but challenging task The associate editor coordinating the review of this manuscript and approving it for publication was Wei Wei . in the process of obtaining a diagnosis based on medical images [22]. Image segmentation has long been one of the most important tasks in the processing and analysis of various types of images. Segmentation is also performed as the first step in various clinical imaging diagnoses, but it is a time-consuming task for a radiologist. Although numerous algorithms have been proposed and published, image segmentation of general anatomical structures across a variety of modalities remains a challenging task. Among the two currently used segmentation methods, automatic segmentation is an attractive prospect, but user intervention is inevitable in real applications because of physicians' work habits.
Although automatic segmentation methods have been investigated for many years, they rarely achieve sufficiently accurate and robust results to be useful for clinical applications in medical imaging, mainly due to the poor image quality, inhomogeneous appearances due to pathology, and variability of protocols among clinicians, leading to different definitions of a particularly structural boundary. VOLUME 8, 2020 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ By learning from large amounts of training data, convolutional neural networks (CNNs) have achieved state-of-the-art performance for automatic segmentation [2]- [5]. However, user interaction still must be integrated into CNN frameworks to obtain accurate and robust segmentation of 2D and 3D medical images [6]. However, many automatic approaches are domain-dependent, and obtaining a large-scale training data set of infants with ground truth annotation is not feasible in practice. Since user interactions provide more direct information about the foreground and background than the image itself, user interactions exert a more direct effect on the segmentation result, and the interaction is more consistent with physicians' current working habits. Interactive segmentation methods, which take advantage of the user's knowledge of anatomy and clinical questions to overcome the challenges faced by automatic methods, are widely used and result in increased accuracy and robustness. However, detecting and segmenting WMI on MRI of neonatal brains using the current automatic or interactive segmentation methods is still difficult and is associated with challenges including 1) the lack of a reliable anatomical map or ''atlas'' for the neonatal brain to guide the segmentation process in areas of low contrast and help distinguish tissues of similar intensities; 2) the small size of the neonatal brain, and the duration for which a neonate is scanned is also shorter than adults, resulting in a low contrast-to-noise ratio, signalto-noise ratio and spatial resolution [7]; and 3) the periventricular white matter related to neonatal HIE is difficult to locate.
We present an anchor-based expanding segmentation approach using sample estimation at the subvoxel level (AE3S), which is an atlas-free and interactive approach to segment periventricular WM in infants, to address these problems. First, the coarse-grained voxels are divided into subvoxels, and the intensity of subvoxels is adjusted using gradient adaptive filtering to solve the problem of low spatial resolution in clinical MR images. Second, users from a team of doctors draw regions of interest (ROIs) to guide segmentation and obtain a series of selected sampling points around targets. Since the population of a single tissue exhibits a normal distribution, the population variance follows a Chi-square distribution, and the population mean displays a normal distribution. Thus, we are able to segment each tissue using the sample parameter estimate as a threshold. Third, since cerebrospinal fluid (CSF) values are many times higher than white matter and gray matter values on MRI, and the image of the ventricle is clearer and more stable, we segment the periventricular region as an anchor using sample parameter estimates and obtain the periventricular white matter by extending the anchor along a gradient. By loosely marking the interior regions instead of finely tracing near borders, our interactive segmentation approach has the properties of a higher processing speed and interactive intuitiveness. Moreover, the proposed approach also has other key features, such as the natural handling of N-dimensional images, spontaneous incorporation of user input, and easy implementation.
Fourth, we selected the effective features for a quantitative diagnosis of HIE in clinical practice after computing and comparing the radiomic indexes of the periventricular white matter obtained using various MRI modalities. The aim of this study is to obtain the features of periventricular white matter related to neonatal HIE from T2WI, FLAIR images, DWI and the apparent diffusion coefficient (ADC). Finally, we evaluate the performance of our approach within the context of infant MR data. We compared the findings with manual sketches prepared by physicians in native space as the ground truth. Our evaluation involves 1) an assessment of the ability of the AE3S approach to perform segmentation and obtain parameters such as Dice scores or overlap coefficients, 2) and the determination of the precision or accuracy of radiomics features of each modality in the neonatal HIE diagnosis.
The three main contributions of this paper are described below. 1) Because the cerebrospinal fluid is more obvious in grayscale than white matter and gray matter on MRI, we segment the more stable periventricular region as an anchor and locate the periventricular white matter of various widths by expanding the anchor along a gradient. This approach satisfies the WM location requirements of neonatal HIE. 2) Since the resolution of clinical MR images is low, we divide the voxels into subvoxels and adjust their intensity values using gradient adaptive filtering to satisfy the fine segmentation requirements. 3) We use the seeds of interactive segmentation approaches as sample data and segment ROIs using the sample estimation methods.
Thus, the anchor-based segmentation, higher resolution obtained by the subvoxel decomposition, and ROIs computed using the sample estimation approach are the novel contributions of this paper, which are rarely reported in other studies.
The related literature is reviewed in Section 2. We describe the segmentation and measurement approach for clinical MRI and its implementation in Section 3. The experimental results presented in Section 4 reveal the effectiveness of the proposed method. The conclusions are presented in Section 5.

II. RELATED LITERATURE A. INTERACTIVE IMAGE SEGMENTATION
The general interactive image segmentation problem can be formulated as described below. Given an image and a label propagation algorithm, the user specifies a labeled seed set to produce the expected segmentation. Popular interactive image segmentation algorithms include seeds consisting of line drawings, contours, bounding boxes, or queries [8]. Interactive segmentation methods have been developed, such as graph cuts, level set, and GrowCut.
Graph cut is a popular method that has been employed to solve a variety of image segmentation problems. One of its limitations is that its energy function is able to model only pairwise interactions between pixels or patches. Thus, graph cut algorithms are unable to easily consider global statistics and higher-order geometric properties, such as the curvature [9]. The level set function approximates a trans-formed version of the pixelwise posterior probability of being part of a target object. The evolution of its zero level set is driven by three force terms: region force, edge field force, and curvature force. These forces are based on a probabilistic classifier and an unsigned distance transform of salient edges [10]. By reformulating the GrowCut algorithm as a clustering problem and requiring only local computations of user input, GrowCut propagates these input labels based on principles derived from cellular automata to classify all pixels. Since all pixels must be visited during each iteration, the computational complexity increases quickly as image size increases [11], [12]. The interactive binary image segmentation method uses the Markov Random Field (MRF) framework and the fast bilateral solver (FBS) technique [13]. In interactive image segmentation based on complex networks, nodes representing pixels are connected to their k-nearest neighbors to build a complex network with the small-world property to propagate the labels, and a regular network in a grid format is used to refine the segmentation of the object borders [14]. A generative probabilistic formulation of the model is established in terms of a ''Gaussian Mixture Markov Random Field'', and a pseudolikelihood algorithm is derived that jointly learns the color mixture and coherence parameters for the foreground and background, respectively. Thus, the scheme was developed based on the global histogram to generate the initial probabilities, and the mean of these smoothed probabilities is used as a threshold to classify all pixels into foreground and background pixels with a probabilistic classifier [15]. Superresolution technology provides a constructive direction for solving the low-resolution problem of medical images captured in the clinic [29], [30]. Combined with the gradient information of another high-resolution contrast image, the edge of brain MRI is improved. The gradient relation model between images with different contrast levels was established to recover a high-resolution image from a low-resolution version of the input image. However, it is difficult for users to operate in clinical work. In a unified automatic glioma segmentation algorithm, spatial fuzzy k-means clustering is used to estimate the region of interest in a multimodal MRI image. Based on the new concept of ''affinity'', some seed points are extracted for region growth to improve image segmentation using region merging and an improved distance regularization level set method [31].
However, these segmentation methods do not take full advantage of the characteristics of cerebral MRI for the periventricular WM location of neonatal HIE. Thus, MR values of cerebrospinal fluid are significantly different from the values of gray matter and white matter, and neonatal HIE is commonly related to periventricular white matter [27]; thus, the WM segmentation requirements are quite different from general segmentation methods.

B. DIAGNOSIS OF NEONATAL HIE
The identification of newborns with ''mild HIE'' within 6 h after birth is challenging due to the limitations in the early existing biomarkers of brain injury and the current knowledge gaps in the long-term neuro-developmental outcomes of infants with mild HIE [26]. A review of 102 articles highlights the emergence of MRI and serum-based biomarkers as useful tools to predict long-term HIE outcomes [16]. The combination of MRI and electroencephalography (EEG) is used as a supportive tool for an early and accurate diagnosis and prediction of developmental HIE [17]. The pattern of brain injury depends on the severity and duration of hypoxia and the degree of brain maturation. Mild to moderate HIE results in periventricular leukomalacia in preterm neonates and parasagittal watershed infarcts in full term neonates [18]. A quantitative diffusion tensor imaging (DTI) analysis of the hippocampus in neonates with HIE is a feasible technique to examine neuronal injury [19]. A tract-specific analysis (TSA) is an alternative method for analyzing the WM that is applicable to large-scale studies. TSA produces a skeleton representation of WM tracts and projects the diffusion data of the group onto the skeleton for statistical analysis. The performance of TSA was evaluated in determining the registration accuracy and the data agreement between the native space and template space, and comparing data obtained from preterm infants with the results obtained from native space tractography and tract-based spatial statistics [20]. A fast, unsupervised, atlas-free WMI detection method was used to detect the ventricles as blobs with a fast linear Maximally Stable Extremal Regions algorithm [7]. A reference contour equidistant from the blobs and the brainbackground boundary were used to identify tissues adjacent to the blobs. Assuming a normal distribution of the intensity of this tissue, the outlier intensities in the entire brain region were identified as potential WMI candidates. Thereafter, false positives were discriminated using appropriate heuristics.
The incidence of hypoxic-ischemic injury in preterm infants displays a more complex clinical course with higher rates of adverse neurological outcomes than in term infants [21]. An automatic algorithm using a random forest classifier for detection and quantification of HIE in diffusionweighted MR images was developed and manually trained on the annotations by two independent observers using 1.5 T data obtained from 20 subjects. The method was applied to a larger dataset including 54 additional subjects scanned at both 1.5 T and 3.0 T. The algorithm findings corresponded well with the injury pattern noted by clinicians for both 1.5 T and 3.0 T data and exhibited a strong relationship with the outcome [22]. The Neonatal Research Network (NRN) score of brain injury on MRI and regional DTI are associated with the optimal mean arterial blood pressure (MAP) in neonates with HIE. A higher optimal MAP and lower mean diffusivity may be related because of cytotoxic edema and limited vasodilation during low MAP in the injured brain. Compared to the NRN score, DTI exhibited better detection of injury with an elevated optimal MAP [23]. A previous study [24] confirmed that MRI obtained 3 h after HIE in the CD1 mouse model of perinatal HIE reliably assessed initial brain injury and predicted the histopathological outcome. Another study [32] proposes an innovative processing pipeline for morphological and arterial spin labeling (ASL) data suited to neonates that enables automated segmentation to obtain cerebral blood flow (CBF) values over ROIs. CBF values determined using ASL were higher in the cortical gray matter, basal ganglia and thalamus (BGT), and frontal and parietal lobes in subjects with abnormal MRI results than in subjects with normal MRI results. The change in the CBF pattern can provide prognostic information in addition to the visible structural abnormalities on conventional MRI. However, in the current clinical diagnosis, neonatal encephalopathy is suspected in infants who are depressed at birth and who, in the earliest hours of life, present mainly with disturbed neurological functions, including altered consciousness, seizures, respiration, and depression of tone. The handbook suggests a routine MRI at 7 days of life to better define and assess the extent of injury, which will assist clinicians in determining a likely prognosis [28]. Compared to the traditional practice of treating medical images as pictures intended solely for visual interpretation, radiomics data contain first-, second-, and higher-order statistics. These data are combined with other patient data and mined with sophisticated bioinformatics tools to develop models that may potentially improve the diagnostic, prognostic, and predictive accuracy [25].
In summary, MRI is an important tool for an early, accurate diagnosis and prediction of HIE development, but segmenting the related regions with HIE and computing quantitative indexes of the related white matter are still tedious tasks.

A. SUBJECTS
Permission for this study was granted by the Shandong Provincial Hospital Ethics Committee, and parental consent was acquired prior to imaging.
MR data were collected from 20 preterm subjects who were imaged between August 2018 and January 2020. All images were reviewed by an experienced neuroradiologist, and 7 subjects with major focal lesions were excluded. Twenty subjects (4 female) born at a median gestational age of 33.53 (30.57-36.85) weeks and imaged at a median corrected gestational age of 34.96 (33.57-36.85) weeks were analyzed in this study. These clinical cases are grouped in Table 1 according to the Queensland Clinical Guideline for hypoxic-ischemic encephalopathy [28].

B. DATA ACQUISITION
The axial scan was performed on a 3 T MR system (Siemens) in the imaging unit of Shandong Provincial Hospital. T1-weighted and T2-weighted MRI, FLAIR, and single shot echo planar DWI data were acquired using an 8-channel phase array head coil. The ADC was used to describe the velocity and range of molecular diffusion motion in different directions on DWI. The ADC and DWI image sample sizes are 144 pixels in the row, 144 pixels in the column, and 40 pixels in the slice. Other image sample sizes are 432 * 432 * 20 pixels in three directions. The pulse sequences are shown in Table 2.
All examinations were supervised by a pediatrician experienced in MRI procedures. Infants were sedated with oral chloral hydrate (25-50 mg/kg) prior to scanning, and pulse oximetry, temperature, and electrocardiography data were monitored throughout the procedure. Ear protection was used, comprising earplugs molded from a silicone-based putty placed in the external auditory meatus and neonatal earmuffs.

C. ALIGNMENT AND REGIONS OF INTEREST
Multimodal MR images were visually inspected in orthogonal planes for the presence of motion artifacts. The skulls were removed at a threshold of 0.2, and corrupt diffusion-weighted volumes were excluded before alignment. The signal-to-noise ratio (SNR) was calculated for each subject from the raw MR data to assess the image quality. The SNR was calculated as the mean signal of the brain region divided by the standard deviation on the corner regions of MRI: SNR = 0.66 * MEAN (I brain ) STD(I corners ).
For each subject, a physician's manual sketch and the AE3S algorithm were used to delineate periventricular WM (PWM) as the ROI in any modality, as described in Table 3.
Output: R, f, ctb-three ROIs of periventricular WM, nine radiomics features and the comparison. 1) P(I )-reprocessing the inputted images; 2) S-drawing and calculating the spaces in T1WI/FLAIR; 3) dividing each voxel into 27 subvoxels in the space S and adjusting the intensity values of subvoxels; 4) while continuing to click? 5) sp-selecting seeds in the periventricular region (CSF) of T1WI/FLAIR; 6) th-estimating the confidence interval for 95% of the intensity mean as the threshold; 7) abh-segmenting the periventricular region as an anchor at th; 8) end while 9) R-expanding the anchor to obtain three ROIs (periventricular WM) of different widths on T1WI/FLAIR; 10) R-aligning the ROIs to other modalities; 11) f -computing the radiomics features of ROIs for each modality; 12) ctb-comparing the features of normal and abnormal subjects. 13) Return R, f, ctb.

D. PROCESSING PIPELINE
In this section, we outline the AE3S algorithm and present the iteratively interactive sampling version. The AE3S pipeline is summarized in algorithm 1.
The algorithm shown above is described in detail below. Accordingly, let I ⊂ R D , D ≥ 2 denote the MR image domain. For various types of MR images, the intensity value at point p(x, y, z) is either a scalar V such as in T1WI and FLAIR images, a gradient vector G, or even a tensor T . For p ∈ I , denote N (p) by the Moore neighborhood of p, i.e., 8 directly neighboring subvoxels of p in 2D and 26 subvoxels of p in 3D. The influence of a point q ∈ N (p) on p is defined as ϕ(G q · G p ), a monotonically decreasing function, where G q = ∇I q and G p = ∇I p are gradient vectors at points q and p, respectively.
In algorithm 1, a voxel must be divided, and the intensity value of a subvoxel must be adjusted according to its new neighborhood. First, a voxel is divided into 3 * 3 * 3 = 27 subvoxels in three dimensions. The intensity value of a subvoxel is adjusted by the states of its 26 neighborhoods. We designed a 3D gradient adaptive filter to adjust each subvoxel. Assuming a 3D image I (x, y, z) and its gradient map ∇I = G = (G x , G y , G z ), the filter is represented by the following equation: where ϕ is a function of gradient G. The core idea is that gradient consistency increases the intensity, and gradient variability reduces the intensity. ϕ increases as the gradient difference increases between its neighborhood and itself according to the following equation: where ∇I c and ∇I h are the gradients of a center subvoxel and its neighboring subvoxel, respectively. Thus, we adjust each subvoxel value through the gradient adaptive filter. One advantage of the proposed filtering method is that images must be analyzed in parallel during computation. Anchor-based segmentation at the subvoxel level is a crucial step. The T1WI, T2WI and FLAIR values of cerebrospinal fluid are three times higher than the values of white matter and gray matter; therefore, the image of the ventricle is clearer and more stable. We first segment the ventricle image out as an anchor and locate the periventricular white matter by expanding the anchor contour. This approach meets the WM location requirement for neonatal HIE.
The seeds in the interior of the tissue that are manually specified during the process of segmenting images are actually sampling data. The points selected by users are the sample data of the labeled class from a tissue or lesion population in MR images. In practice, the intensity population of one single tissue on MRI follows a normal distribution. Suppose a subvoxel sample set from a periventricular population is (Ī , S 2 ) and the subvoxel values are independent of each other; then, the population variance displays a Chi-square distribution, and the population mean displays a normal distribution. The population variance and mean intervals are estimated as follows: where the sample size will be sufficient if a team of physicians combine their samples. Here, the sample size must be greater than 30. In the initial iteration, the upper and lower intensity thresholds are set to µ ± 2 * σ to segment the periventricular region. Physicians can repeat the procedure to display the segmented tissue and to sample points. The process from displaying the data to interval estimation continues until the final segmentation is reached. When updating the current segmentation, the new images are based only on the previous result, which significantly reduces the time required for refinement.
Once the anchor of the periventricular region has been identified, we can easily segment the periventricular white matter by expanding its boundary using a combination of VOLUME 8, 2020 operations, such as dilating and closing procedures. Meanwhile, segmentations of other imaging modalities, such as T2WI, DWI and ADC, are obtained by aligning them with T1WI or FLAIR images.
The proposed method was implemented in a MATLAB script, which provides a framework for interactive operations and visualization. The entire implementation is available as a part of an open source package available at [http://graphicslab.sdnu.edu.cn].

E. EVALUATION OF AE3S
We assessed the segmentation aspect of the pipeline using manual sketching by a professional physician as the ground truth to evaluate the AE3S algorithm.
WM values delineated in the T1WI or FLAIR space were warped into the spaces of other modalities using the transformations from cross-entropy-based registration. The segmented area of each subject was converted into binary ROIs, which were compared using Dice scores. The Dice score measures the degree of overlap between two periventricular WM volumes or the numbers of subvoxels.
We computed values in the T1WI, T2WI, FLAIR, DWI and ADC spaces and their radiomics features and constructed histograms and box-and-whisker-plots to determine how accurately the AE3S algorithm projects the MR data from each subject in determining the HIE diagnosis. For the abnormal and normal subject groups, we tested the hypothesis for each combination of features, regions, and modalities. Student's t test was used to compare the differences in paired modalities.
The value of each modality calculated from the normal subject group was considered the ground truth. These values were mainly determined for three periventricular WM regions and nine radiomics features. The statistical values provide a summary of the differences among subjects, whereas the t-value calculation assesses the difference in their deviation from the data within the WM regions. Additionally, we examined the effects of the AE3S algorithm on the metrics of each modality. We expected to observe an increase or a decrease in values based on the presence or absence of HIE, consistent with the predicted pattern observed above

IV. RESULTS AND DISCUSSION
We compared the proposed method of analyzing the MR images of HIE in infants with the regions specified by physicians in the native space as our gold standard. In this test, we focused on three periventricular WM regions and measured their T1WI, T2WI, FLAIR, DWI and ADC values. We divided the clinical FLAIR/T1WI findings into subvoxel data and adjusted the value of each subvoxel. Next, we segmented the periventricular WM regions at the FLAIR/T1WI subvoxel level and obtained the registrations of T2WI, DWI, and ADC images. We achieved a more accurate location by expanding the segmented periventricular region, since the periventricular region is a precise anchor. Finally, we determined the radiomics features and performed a statistical analysis of various modalities in the ROIs.

A. SEGMENTATION EVALUATION
The subvoxel level images had a higher resolution and helped us to refine the next segmentation of the periventricular region. We first segmented the CSF and obtained the periventricular WM on FLAIR images/T1WI and then aligned the images to T2WI, DWI and ADC images obtained using method B in section III to obtain the parameters of the corresponding regions.
The intensity values of the three tissues are a mixture of Gaussian distributions; the CSF values are in the lower range, and the WM and GM values are in the higher ranges. CSF accounts for the vast majority of the area and is approximately 100 times the amount of WM and GM. Therefore, we must roughly identify the ventricular regions through an interaction.
The periventricular CSF with a high contrast on FLAIR images or T1WI is positioned as an anchor, and the periventricular WM is obtained by extending the periventricular boundary along gradient directions. Fig. 1 shows a comparison of the segmentation at the subvoxel and voxel levels. The subvoxel segmentation is smoother and finer.
We compared CSF segmentation using the AE3S algorithm at the voxel and subvoxel levels with the ground truth of manual sketches of experienced physicians at the voxel level, as shown in Fig. 2. The AE3S method is superior in drawing complex contours.
The Dice score, which is described in subsection E of section III, of the voxel segmentation (0.97) on the left was higher than the subvoxel segmentation (0.59) on the right, representing the degree of alignment between the two segmentations. This finding was mainly attributed to the manual sketches used as the ground truth that were performed at the voxel level. However, the subvoxel segmentation method may be more accurate, as shown in Fig. 2. The improvements are likely to be due to both leveraging the full subvoxel information and the clearer contours. Manual methods discard division information, increasing the difficulty of distinguishing neighboring voxels and blurred boundaries.
Our segmentation method explored sampling data that were averaged from the study cohort, which reduces the heterogeneity of the sample data. Our analysis represented a novel approach to increase segmentation accuracy in the neonatal population. Previous studies had focused on using image similarity measures and tissue label overlap scores to assess registration performance.
Notably, our proposed segmentation method has high productivity for the clinical MRI resolution, but a long time is needed to separate voxels into subvoxels, even in a parallel framework.

B. RADIOMICS FEATURES
We calculated and plotted histograms, box-and-whisker plots, and nine radiomics features (mean, std, entropy, uniformity, energy, median, lower quartile, upper quartile and interquartile range) for each modality and for each subject group.
The characteristics of the modality values from three WM regions of the normal group and the abnormal group at the voxel and subvoxel levels were compared using plots and radiomics features. As an example, PWM 1 is shown in Figs. 3 to 6.
Figs. 3-6 show the histograms, box-and-whisker plots, and radiomics features of the five modalities (FLAIR, T2WI, T1WI, DWI, and ADC from left to right) in PWM 1 at the voxel and subvoxel levels. A second group includes the representations of the five modalities in PWM 2, and a third group includes the representations of the five modalities in the deep WM. The latter two groups are not listed here. As shown in the figures, the values of the modalities vary with the width of the periventricular WM region. Figs. 3 and 5 show a comparison of the voxel and subvoxel characteristics of each modality. Although the subvoxel segmentation method achieves accurate contours, the statistical plots and features are similar. Some differences are observed, most noticeably for the outlier points, due to the gradient filtering.
Although the five modality charts show similar trends across the figures, they are very different in each figure.
In each modality chart, the two left columns are selected from the abnormal group, and the two right columns are selected from the normal group. The difference between the abnormal and the normal group is obvious in the FLAIR, T1WI, and T2WI results.
We presented an analysis of the concordance between visual observation and statistical data derived from the AE3S algorithm. The close agreement confirms that the AE3S method is able to more accurately represent neonatal HIE from MRI data. The discrepancy between the visual diagnosis and the AE3S method is most likely due to the (i) low resolution of clinical images, (ii) underdeveloped neonatal brains, and (iii) MRI artifacts. The close agreement of FLAIR and T2WI findings may be due to their ability to clearly reflect abnormal areas.

C. KEY REGIONS AND MODALITIES
Our main purpose was to determine which modality and region show differences between abnormal and normal subjects. Therefore, we needed to determine population differences between the negative and positive groups for various aspects of regions and modalities.
The box-and-whisker plots of the five modalities are shown for the abnormal and normal groups in Fig. 7; subjects from the abnormal group (left panel) exhibited higher values than subjects in the normal group (right panel). The median features of the FLAIR and T1WI modalities in the three regions were potential diagnostic indicators of HIE (p<0.05), as shown in Fig. 8. Other features, modalities, and regions may exert the same effect to some extent, as shown at [http://graphicslab.sdnu.edu.cn]. For example, the energy features of FLAIR and ADC indicate lesions in the three regions. Briefly, these nine radiomics features and the three regions on FLAIR images, T2WI, and T1WI were significantly correlated with HIE lesions. Currently, the indications, which are not explicitly calculated, are empirically applied in the clinical diagnosis of infants with HIE, and all modality trend charts are approximately consistent with the charts obtained from previous clinical diagnoses.
However, target zones consisting of periventricular locations are difficult to identify at birth in term infants. We were able to delineate the PWM using the AE3S algorithm, apart from other segmentation methods, confirming that the AE3S method was successfully applied to neonatal populations.

D. COMPARISON OF VARIOUS DIAGNOSTIC METHODS
Compared to the current clinical visual interpretation of MRI, the AE3S method provides physicians with quantified features of injury due to HIE. It also provides the quantized values of other modalities for the periventricular white matter, for example, CBF and the oxygen extraction fraction (OEF), in addition to the five modalities already mentioned. In the Queensland Clinical HIE Guideline (QCHG), MRI is routinely used for qualitative analysis, and the assessment of neurological status is the major criterion used to determine the probability of HIE in infants [28]. Experienced radiologists are able to distinguish and even score brain injury on MRI, including normal, subnormal and abnormal MRI results. In comparison, our methods may help physicians compute quantitative indexes of various MR modalities and accurately assess the stages of HIE in the next step.
ROIs were derived from a custom anatomical atlas built from ALBERTs data to measure the changes in CBF values in neonates with HIE using ASL. A quantitative CBF map for each subject yielded CBF statistics in the ROIs combined with anatomical segmentation [32]. However, a visual evaluation of the images and segmentation are still required. In the present study, CBF values in the BGT were higher in subjects with abnormal MRI results, and the values were significantly higher than in the cortical gray matter. However, no statistically significant differences were observed between children with normal and abnormal MRI results. Due to the inherent low resolution, the intensity of a single voxel results from a mixture of the measured GM, WM, and CSF signals, and ASL suffers from partial volume effects. Since treatment of the components of CBF effectively limits crosscontamination between brain tissues, more sophisticated methods for partial volume correction should be implemented in the next step. Comparatively, we are able to provide physicians with the radiomics features to distinguish differences in various modalities, ROIs, and subject attributes between normal and abnormal groups. Additionally, considering the prominent impact of CSF, we exclude CSF interference in the ROIs by focusing on the CSF.
The OEF represents an important measurable parameter of brain metabolism and is a hallmark of at-risk tissue in subjects with acute stroke. The OEF within cerebral veins has been measured using advanced quantitative susceptibility mapping (QSM) MRI reconstruction [33], [34]. A significant association was identified between the relative CBF and the absolute OEF in the ischemic hemisphere of patients. T2WIStar imaging also identified an increase in the OEF as the CBF decreases in ischemic tissue. The sensitivity and accuracy of QSM are expected to improve with higher spatial resolution and direct modeling of the vessel shape. Because the OEF values are measured in veins at the hemispheric level and not at the tissue level, vein-based QSM methods are unable to provide voxelwise information about ROIs in the tissue. The most notable limitation regarding the postprocessing of OEF maps is the insufficient spatial resolution. Another obvious limitation is that patients with a previous history of hemorrhage should be excluded due to potential hemosiderin deposition in the brain tissue. We solved the low spatial resolution problem using voxel decomposition methods. We anchored the CSF to segment accurate ROIs and obtain the PWM that is closely related to neonatal HIE. Moreover, we presented the pipeline to compute nine statistical features for each ROI and modality.
More importantly, our approach involved the combination and analysis of modalities and ROIs. Based on these values, physicians will be able to determine which combination of modalities and ROIs is closely related to the HIE stage. A comparison of four methods is shown in Table 4. . Nine radiomics features of the 5 modalities for PWM 1 from four randomly selected subjects are shown at the subvoxel level. The panels from left to right indicate the FLAIR, T2WI, T1WI, DWI, and ADC modalities. The first two subjects are from the abnormal group, and the last two are from the normal group in each panel.

FIGURE 7.
Box-and-whisker plots of the 5 modalities from the abnormal and normal groups in three regions at the subvoxel level. Panels from left to right indicate the FLAIR, T2WI, T1WI, DWI, and ADC modalities. Panels from the top to bottom indicate the PWM 1, PWM 2, and deep WM regions. The abnormal group is shown on the left of each plot, and the normal group is shown on the right.

FIGURE 8.
Results of the t test of the mean features for each combination of the 5 modalities and 3 regions at the subvoxel level. FLAIR, T2WI, T1WI, DWI, and ADC are numbered 1, 2, 3, 4, and 5, respectively. PWM 1, PWM 2, and deep WM are numbered 1, 2, and 3, respectively. If t-value > t-alpha (P<0.05), a significant difference is present in the mean feature for that combination. The mean features of combinations 1, 2, 3, 7, 8, 9, and 15 are possible neonatal HIE indicators.

E. ADVANTAGES AND LIMITATIONS OF THE AE3S METHOD
The AE3S method has several advantages as an analytical tool. It focuses on the stable CSF and uses it as an anchor to segment the ROIs of the PWM. It further divides the clinical images into a fine-grained subvoxel level, approximating research data, and the AE3S method offers improved alignment between different modalities through subvoxel granularity. A major benefit of the AE3S method is that it is anatomically specific and comparable to clinical technology, similar to other interactive methods. A possible application of the AE3S method is to study areas in subjects where the brain atlas lacks data due to a limited number of cases.
A limitation of the AE3S method, as assessed here, is that the subvoxel decomposition requires excessive resources and time, even in a parallel framework. The decomposition processing of each modality required approximately 90 minutes on an Intel i5 1.9 GHz CPU with 8 GB of RAM, and 11 minutes on an Intel hexa-core i7 CPU with 32 GB of RAM and a 4 GB GPU. The subvoxel decomposition is an offline operation. The white matter and the gray matter are intertwined in some periventricular regions. We omitted mixed gray matter in the segmentation process because the AE3S method is ill-suited for these structures.

V. CONCLUSION
In this paper, we presented an effective and efficient image segmentation and analysis method for the clinical diagnosis of HIE in neonates, which was performed by anchor-based expansion, interactive sampling, and parameter estimation. The results revealed the clinical effectiveness and practicability of our proposed method.
In future studies, we plan to more rigorously investigate how closely the proposed method approximates the clinical imaging diagnosis and prognosis of neonatal HIE. Furthermore, we aim to introduce a smoothing procedure into the interactive approach to refine the multilabel segmentation. Regarding the latter point, we are exploring the incorporation of the methodology proposed in popular automated platforms. VOLUME 8, 2020