Optimizing Automated Brain Extraction for Moderate to Severe Traumatic Brain Injury Patients: The Role of Intensity Normalization and Bias-Field Correction

Abstract

Accurate brain extraction is crucial for the validity of MRI analyses, particularly in the context of traumatic brain injury (TBI), where conventional automated methods frequently fall short. This study investigates the interplay between intensity normalization, bias-field correction (also called intensity inhomogeneity correction), and automated brain extraction in MRIs of individuals with TBI. We analyzed 125 T1-weighted Magnetization-Prepared Rapid Gradient-Echo (T1-MPRAGE) and 72 T2-weighted Fluid-Attenuated Inversion Recovery (T2-FLAIR) MRI sequences from a cohort of 143 patients with moderate to severe TBI. Our study combined 14 different intensity processing procedures, each using a configuration of N3 inhomogeneity correction, Z-score normalization, KDE-based normalization, or WhiteStripe intensity normalization, with 10 different configurations of the Brain Extraction Tool (BET) and the Optimized Brain Extraction Tool (optiBET). Our results demonstrate that optiBET with N3 inhomogeneity correction produces the most accurate brain extractions, specifically with one iteration of N3 for T1-MPRAGE and four iterations for T2-FLAIR, and pipelines incorporating N3 inhomogeneity correction significantly improved the accuracy of BET as well. Conversely, intensity normalization demonstrated a complex relationship with brain extraction, with effects varying by the normalization algorithm and BET parameter configuration combination. This study elucidates the interactions between intensity processing and the accuracy of brain extraction. Understanding these relationships is essential to the effective and efficient preprocessing of TBI MRI data, laying the groundwork for the development of robust preprocessing pipelines optimized for multi-site TBI MRI data.

1. Introduction

Traumatic brain injury (TBI) is a significant public health challenge, hospitalizing approximately 230,000 people annually in the United States alone and leading to a wide range of long-term cognitive and physiological impairments [1,2,3]. Among the most common and most debilitating sequelae of TBI is the development of post-traumatic epilepsy (PTE), a condition characterized by recurrent seizures persisting for more than one week after a traumatic insult to the brain [4]. Epidemiological studies have shown that TBI is a major risk factor for epilepsy, with the likelihood of seizure development increasing with the severity of the injury [5]. The pathophysiological mechanisms underlying PTE are complex and involve changes in brain architecture, neuroinflammation, and alterations in neuronal network excitability [6]. Research in this area is ongoing, with studies employing advanced neuroimaging techniques, animal models, and clinical trials to unravel the intricacies of PTE and improve health outcomes for patients who suffer from it [4]. Previous work has shown that the development of PTE typically follows a period of latency, extending from several months to years after the initial insult [7]. As of now, there are no established preventive measures available to counteract the emergence of seizures during this critical period, largely because substantive biomarkers allowing for the early detection of patients at risk for PTE have not yet been identified. Pioneering the search for these biomarkers, the Epilepsy Bioinformatics Study for Antiepileptogenic Therapy (EpiBioS4Rx) is an international, multidisciplinary study dedicated to the detection of the biomarkers of PTE to identify at-risk patients before seizures emerge, thus providing the prerequisite foundation necessary to launch large-scale clinical trials aimed at the development of antiepileptogenic drugs capable of preventing PTE [8]. Through EpiBioS4Rx, we have acquired substantial amounts of neuroimaging (MRI), electrophysiological (EEG), and molecular data. The focus of this current investigation is MRI, which is a key data modality for PTE biomarker identification because it provides the granular spatial resolution necessary for the precise quantitative characterization of the heterogeneous brain lesions resulting from TBI.

Heterogeneity, from patient presentations to the pathophysiology of brain lesions, has been acknowledged as a hallmark of TBI [9]. However, the same heterogeneity that grants MRI relevance in the search for PTE biomarkers also introduces significant limitations and technical challenges to its analysis [9,10]. TBI may present with extensive lesions, deformations, edema, and skull fractures [11,12,13]. TBI can also induce the displacement of brain ventricles and complicate the intensity distributions of MRIs as well [14]. Each case of TBI manifests a unique set of morphometric abnormalities, precipitating substantive alterations in brain architecture that can obscure the anatomical landmarks typically prominent in healthy brains. The complexity and variability of traumatic brain injuries is a major issue for automated methods of processing TBI MRI data [11,15]. Bennett et al. [16] demonstrate the necessity for modified automated lesion segmentation methods within the EpiBioS4Rx framework. Their work highlights that existing segmentation algorithms struggle with deformities introduced by TBI, thus necessitating manual segmentation in the absence of reliable automated TBI lesion segmentation methods.

Challenges to Preprocessing of Multi-Site TBI MRI Data

Automated brain extraction is a critical step in MRI preprocessing where TBI-related challenges manifest. Brain extraction, or skull strip**, refers to the separation of brain tissue from non-brain components, such as the skull, scalp, and dura mater, and is often one of the first and most critical preprocessing steps in MRI analysis. While manual brain extraction is often more accurate than automated brain extraction for TBI MRI data due to the limitations of contemporary automated methods, manual extractions are susceptible to human error and are plagued with inter-rater reliability issues. Moreover, manual methods are far more logistically complex and labor-intensive, taking from weeks to months to accomplish what an automated method can in seconds. A pivotal development in the domain of brain extraction has been the automated Brain Extraction Tool (BET) from the FMRIB Software Library (FSL) [17,18]. BET utilizes a deformable model guided by Bayesian priors to generate an outline of the brain’s surface. In the healthy brain, this strategy is overwhelmingly effective and has made the tool very popular, with over 8000+ studies citing it to date [19].

While BET performs extractions effectively in the healthy brain, its reliance on Bayesian priors severely reduces its accuracy in brains with pathology (see Figure 1c). This is because Bayesian priors rely on predefined assumptions of what a “typical” brain should look like. In brains with TBI, these assumptions are often broken, ultimately resulting in inaccurate brain extractions [20]. Inaccurate brain extraction is a critical misstep in MRI analyses. The precision of brain extraction directly influences the accuracy of all volume-based brain analyses. However, because brain extraction is one of the first steps in MRI preprocessing pipelines, inaccurate brain extraction is also detrimental to the accuracy of subsequent preprocessing steps, such as tissue segmentation and registration. Thus, poor-quality brain extraction prompts a cascade of errors that fundamentally compromise the validity of findings with these data, obscuring true insights into brain anatomy and function and potentially leading to misinterpretations of neural substrates underlying cognitive processes and disease mechanisms.

Several methods that improve brain extraction for patient data have been proposed. These techniques aim to overcome the shortcomings of conventional automated methods in extracting the brain. ROBEX from Iglesias et al. [21], SPECTRE from Carass et al. [22], BEaST from Eskildsen et al. [23], optiBET from Lutkenhoff et al. [20], and MONSTR from [24] have shown improved performance in TBI data compared to conventional techniques like BET. Of these, optiBET is particularly notable for its outstanding performance in 74 human MRIs with TBI—outperforming machine learning-based algorithms such as ROBEX [20]. OptiBET is a brain extraction script optimized for brains with pathology that, in addition to its robust testing on TBI data, also stands out from other extraction tools in being entirely composed of freely available and near-universally employed neuroimaging software. These strengths have made OptiBET a widely adopted tool for TBI brain extraction. To perform brain extraction, the optiBET script first employs BET to make the initial brain outline estimate, then it performs linear and non-linear registration to the MNI template space, and finally, an MNI brain mask is back-projected to the original subject space to mask out non-brain tissue from the original image [20]. While optiBET demonstrates significantly improved performance in TBI brain extractions, it is unreliable, requiring extensive manual correction in some cases—shown in Figure 1d. OptiBET especially shows a decline in performance on T2-FLAIR MRIs, which is a specialized MRI sequence crucial for the detailed visualization and segmentation of TBI-induced brain lesions [16,24,25].

Intensity normalization and intensity inhomogeneity correction are also two essential steps in the preprocessing of multi-site TBI MRI data. While intensity inhomogeneity correction and intensity normalization represent two distinct preprocessing techniques, for simplicity, they are often collectively referred to as “intensity processing” in later sections. Intensity inhomogeneity correction mitigates spatial variations in signal intensity that are not inherent to the anatomical structure but rather result from the variations in the MRI scanner’s magnetic field, called the bias field [26,27]. Conversely, intensity normalization aims to standardize the distribution of intensity values. Normalizing the intensity distribution across different images allows for inter-scan and inter-scanner comparisons [28,29]. In raw MRI data, variations in intensity can be attributed both to tissue properties and to the MRI scanner settings used during acquisition. Due to this dependency on scanner settings, comparing intensities across MRIs is problematic, especially when different settings or scanners are involved [28]. Variability persists among all scans, even when using identical scan and scanner settings, making intensity normalization essential for meaningful comparisons across MRI scans [28,29]. Intensity normalization and intensity inhomogeneity correction are critical for the accuracy and consistency of MRI data analysis and, thus, are essential components of MRI preprocessing pipelines [26,27].

Considering the limitations of contemporary automated techniques in performing brain extractions in the presence of TBI, the fundamental importance of accurate brain extraction to the validity of MRI analyses, and the necessity for intensity processing for multi-site MRI data sets, there is an urgent need for the development of more reliable and efficient automated brain extraction methods and a greater understanding of the interplay of different preprocessing steps in general to create an optimized MRI preprocessing pipeline for TBI data. To address these challenges, the current investigation examines the utility of N3 bias-field correction from Sled et al. [30] and three different intensity normalization algorithms from Reinhold et al. [31] in refining the precision and reliability of T1-weighted Magnetization-Prepared Rapid Gradient-Echo (T1-MPRAGE) and T2-weighted Fluid-Attenuated Inversion Recovery (T2-FLAIR) TBI MRI brain extractions with optiBET and BET.

This study offers several critical contributions to the field of neuroimaging preprocessing, specifically addressing challenges in the automated extraction of patient brains. Our primary objective is to develop and validate optimized preprocessing techniques for automated brain extraction in patients with moderate to severe TBI, addressing a critical need in the analysis of large-scale, multi-site TBI MRI data. We provide a comprehensive evaluation of the effects of intensity normalization and bias-field correction on the performance of automated brain extraction tools (BET and optiBET) using a large, multi-site dataset of patients with TBI. This systematic analysis across 154 preprocessing pipelines and two MRI sequence types (T1-MPRAGE and T2-FLAIR) allows for robust comparisons and far-reaching insights. We identify specific parameter configurations of N3 bias-field correction that significantly improve brain extraction accuracy for TBI data, with notable differences between T1-MPRAGE and T2-FLAIR images. These findings can directly guide the optimization of preprocessing pipelines for TBI MRI analysis. We find that alternative BET configurations (BET|R and BET|B,f0.2,g0.3) outperform the widely used BET|B,f0.1 setting for TBI data. This highlights the need for TBI-specific validation of preprocessing techniques and suggests avenues for improving tools like optiBET. We demonstrate the nuanced effects of multiple intensity normalization methods (Z-score, WhiteStripe, KDE-based) on brain extraction performance, providing a basis for further research. Collectively, our results delineate the complex interplay between intensity processing and brain extraction in the presence of TBI, offering valuable insights for the development of robust, optimized preprocessing pipelines tailored to the challenges of TBI MRI data. The detailed analysis of technique interactions and TBI-specific effects informs both the refinement of existing methods and the creation of new, specialized tools for TBI neuroimaging research.

2. Materials and Methods

In this section, we describe the dataset used in the experiment, the different automated brain extraction pipelines we compared, the inhomogeneity correction, normalization, and brain extraction configurations that form the basis of these pipelines, and the method of extraction quality assessment.

2.1. Participants

This research received approval from the Institutional Review Board of the University of California, Los Angeles, and the local ethics committees of all participating institutions within the EpiBioS4Rx Study Group. EpiBioS4Rx successfully enrolled 250 participants from 13 international locations with advanced clinical TBI programs. The demographic data and enrollment criteria of the participants in the EpiBioS4Rx study have been thoroughly described in previous work [16]. In brief, the full EpiBioS4Rx cohort has an average age of 44.28 years (SD = 21.2) and consists of 192 males and 58 females, with an average arrival Glasgow Coma Scale (GCS) score of 7.89 (SD = 3.91), and the MRI post-injury was taken, on average, 12.16 days (SD = 12.54) after the injury. Regarding the enrollment criteria, subjects were admitted within 72 h of experiencing a moderate to severe TBI. The inclusion criteria specified patients between 6 and 100 years old with acute TBI evident from intracranial, cortical, and/or subcortical bleeding on CT imaging, a GCS score ranging from 3 to 13 upon emergency department arrival, and the capability of undergoing continuous EEG monitoring for seven days post-injury, as well as MRI scanning within 18 days of the injury. The exclusion criteria ruled out individuals with diffuse axonal injury without hemorrhagic contusions; known cases of HIV/AIDS or Hepatitis B or C; those who are pregnant; those with pre-existing neurological diseases, CNS malignancies, epilepsy/seizure disorders, or dementia; cases of isolated anoxic brain injury, severe cervical spine injuries, or brain death; those in situations of current or impending incarceration; and those with positive tests for COVID-19.

Here, we report analyses performed on T1-MPRAGE and T2-FLAIR MRIs from EpiBioS4Rx subjects who had gold-standard manually corrected brain extractions available. A total of 131 subjects had T1-MPRAGE, and 83 subjects had T2-FLAIR gold-standard extractions available. Six T1-MPRAGE and eleven T2-FLAIR sequences were determined to be unusable. Therefore, the total number of subjects with a manually corrected, gold-standard brain extraction available was 125 T1-MPRAGE and 72 T2-FLAIR images. The cohort has an average age of 42.52 years (SD = 20.67) and consists of 77.62% males and 22.38% females, with an average arrival Glasgow Coma Scale (GCS) score of 7.58 (SD = 3.79), and the MRI post-injury was taken, on average, 11.99 days (SD = 7.31) after the injury.

2.2. MRI Acquisition

Structural MRI scans were obtained using 3-Tesla or 1.5-Tesla scanners following the EpiBioS4Rx protocol described in [8]. Approximately two weeks post-injury, 3D T2-FLAIR images were captured (average = 11.81, standard deviation = 7.74), with settings as follows: slice thickness of 1 mm, a field of view of 256 mm, a frequency of 256 Hz, a flip angle of 90° or greater than 120°, a repetition time (TR) exceeding 5000 ms, an echo time (TE) between 80 and 140 ms, an inversion time (TI) of 2000 to 2500 ms, no inter-slice gap, an isotropic resolution of 1 mm, and an excitation number (NEX) of at least 1, employing up to 2x parallel imaging. For the acquisition of 3D T1-weighted images, a Magnetization-Prepared Rapid Gradient-Echo (MPRAGE) sequence was utilized. The imaging parameters included a slice thickness of 1 mm, a field of view of 256 mm, a frequency of 256 Hz, flip angles ranging from 8° to 15°, a TR of 1500 to 2500 ms, the smallest possible TE, an inversion time of 1100 to 1500 ms, no inter-slice gap, an isotropic resolution of 1 mm, a minimum NEX of 1, and the use of up to 2x parallel imaging.

2.3. Gold-Standard Brain Extraction Protocol

Manual brain extractions are needed to validate and assess the accuracy of automated brain extraction methods because manual extraction is currently the most accurate method of brain extraction for TBI cohorts. These images serve as a definitive reference for the accuracy of the automated extraction, enabling the assessment of pipeline performance. The protocol used to create our “gold-standard” brain extractions is as follows. First, reorientation was carried out with FSL’s fslreorient2std, and, following this, two different brain extraction techniques were employed: BET with options B and f = 0.4 (Figure 1c) and optiBET (Figure 1d) [32]. Two trained researchers were assigned to independently review the extraction quality for the same subject. The extractions were visually examined using graphical user interfaces (GUIs) such as MRIcron, FSLEYES, or ITK-SNAP [33,34,35]. Both researchers reviewed the BET and optiBET images and selected the superior brain extraction for further refinement. For any subject where there were conflicting ratings, a staff researcher reviewed the images and made the final decision.

Then, several techniques were attempted to refine brain extraction quality. Initially, various methods were used to try to improve extraction quality, such as adjustments to BET option “f” and dilation and erosion with fslmaths. If the extraction was still of unacceptable quality, manual correction was performed as a last resort to fill in or remove voxels as necessary. Manual corrections were performed by a group of trained researchers using fsleyes [35]. The change in extraction quality from Figure 1c,d to Figure 1b exemplifies the degree of corrections made during this part of the procedure. After completing the manual corrections, a final dual-review process was carried out to ensure the quality of the corrected brain extractions. A staff researcher and a trained volunteer researcher with three years of neuroimaging experience independently reviewed each manually corrected image. In instances where the quality of the corrected brain extraction was uncertain or poor, a final evaluation was performed by a medical doctor with extensive neuroimaging experience, who made final decisions regarding extraction quality and the need for further corrections.

2.4. Brain Extraction Pipelines

We ran 154 brain extraction pipelines on each of the 125 T1-MPRAGE and 72 T2-FLAIR images for a total of 30,338 extractions. A brain extraction pipeline comprises one of the extraction configurations (listed in

Table 1. A description of the BET and optiBET brain extraction configurations. The Primary Option column contains the brain extraction tool used (BET or optiBET) and the primary BET option, “B” or “R”, if applicable. The Option “f” and “g” Value columns contain the fractional intensity threshold (“f”) and the vertical gradient in the fractional intensity threshold (“g”) values if they differed from their respective defaults of f = 0.5 and g = 0. The Label column contains the short-hand notation used when referencing the extraction configuration. The Description column provides a brief description of how the configuration’s specified options affect the brain extraction.

Primary Option	Option “f” and “g” Values		Label	Description
BET	-	-	BET	Default BET configuration (f = 0.5 and g = 0).
	f 0.2	-	BET|f0.2	Generates larger brain outline estimates.
	f 0.8	-	BET|f0.8	Generates smaller brain outline estimates.
	-	g -0.3	BET|g-0.3	Generates smaller brain outline estimates at the bottom of the image and a larger outline at the top.
	-	g 0.3	BET|g0.3	Generates a larger outline at the bottom of the image and a smaller outline at the top.
	f 0.2	g 0.3	BET|f0.2,g0.3	Generates a larger brain outline estimate, with the bottom of the image having a larger estimate relative to the top.
BET|B	-	-	BET|B	Performs FAST bias-field correction and standard-space masking for image bias-field reduction and neck voxel cleanup.
	f 0.1	-	BET|B,f0.1	Performs bias-field correction and generates a larger brain outline estimate.
	f 0.2	g 0.3	BET|B,f0.2,g0.3	Performs bias-field correction and generates a larger brain outline estimate, with the bottom of the image having a larger estimate relative to the top.
BET|R	-	-	BET|R	Runs BET iteratively for robust brain center estimation.
optiBET	f 0.1	-	optiBET	Performs BET|B,f0.1 for the initial extraction and then FLIRT and FNIRT to generate the final extraction by masking the input image with a back-projected standard brain mask.

) and one of the intensity processing procedures (listed in

Table 2. A description of the intensity processing procedures. IP procedures were executed prior to brain extraction with one of the configurations detailed in Table 1. The IP (intensity processing) column contains the intensity processing technique used. The Iter. (iterations) and Mask columns contain the number of iterations that a technique was run for and the brain extraction configuration used to generate the brain mask that was applied during intensity processing, if applicable. The Label column contains the short-hand notation used when referencing the IP procedure. The Description column provides a short description of each IP procedure.

IP	Iter.	Mask	Label	Description
N3	1	-	N3 iter = 1	1 iteration of N3 inhomogeneity correction without a mask.
	4	-	N3 iter = 4	4 iterations of N3 without a mask (FreeSurfer default).
ZS	1	-	ZS iter = 1	1 iteration of Z-score intensity normalization without a mask.
	4	-	ZS iter = 4	4 iterations of Z-score int. norm. without a mask.
	1	BETf0.1	ZS m = BETf0.1	1 iteration of Z-score with a BETf0.1-extracted mask.
	1	BETf0.5	ZS m = BETf0.5	1 iteration of Z-score with a BETf0.5-extracted mask.
	1	BETf0.8	ZS m = BETf0.8	1 iteration of Z-score with a BETf0.8-extracted mask.
KDE	1	BETf0.1	KDE m = BETf0.1	1 iteration of KDE-based intensity normalization with a BETf0.1-extracted mask.
	1	BETf0.5	KDE m = BETf0.5	1 iteration of KDE with a BETf0.5-extracted mask.
	1	BETf0.8	KDE m = BETf0.8	1 iteration of KDE with a BETf0.8-extracted mask.
WS	1	BETf0.1	WS m = BETf0.1	1 iteration of WhiteStripe intensity normalization with a BETf0.1-extracted mask.
	1	BETf0.5	WS m = BETf0.5	1 iteration of WhiteStripe with a BETf0.5-extracted mask.
	1	BETf0.8	WS m = BETf0.8	1 iteration of WhiteStripe with a BETf0.8-extracted mask.
-	-	-	-	No intensity processing procedure (control).

). The intensity processing procedure refers to the first step of each pipeline and is executed prior to the actual brain extraction. The extraction configuration refers to the second step of each pipeline, which executes the actual brain extraction via BET or optiBET. Of the 154 pipelines tested, 140 performed brain extraction with BET, and 14 used optiBET. These two steps together make up the full brain extraction pipeline. For each extraction configuration, one pipeline had no intensity processing procedure, which is referred to as the extraction configuration’s control pipeline.

Extraction Configurations

One hundred forty pipelines performed brain extraction with one of the ten BET configurations described in Table 1. These ten BET configurations differ in their specific parameters or “options” set in the BET command. BET option “B” incorporates intensity inhomogeneity correction, also called bias-field correction. BET option “R” robustly locates the brain’s center of gravity more accurately by performing multiple iterations of BET to refine the center-of-gravity estimation. Options “B” and “R” are mutually exclusive and optional. BET option “f” adjusts the fractional intensity threshold, which adjusts the size of the brain outline estimation, with values closer to 0 yielding larger outlines and values closer to 1 yielding smaller outlines. BET option “g” adjusts the vertical gradient in the fractional intensity threshold, modifying brain outline estimations along the vertical axis such that positive values increase the outline size at the brain’s bottom and decrease it at the top, while negative values do the opposite. Unless different options are specified, BET runs with options f = 0.5 and g = 0 by default. The ten BET configurations we tested utilize one or none of the mutually exclusive primary options, “B” and “R”, along with a combination of “f” and/or “g” option values.

Fourteen pipelines performed brain extraction with the default configuration of optiBET from Lutkenhoff et al. [20]. The difference between optiBET and BET is that optiBET is a shell script designed to perform brain extractions even in brains with significant pathology. The first step of the optiBET script is performing a preliminary brain outline estimate using BET with options “B” and “f = 0.1” (BET|B,f0.1 in Table 1). This specific configuration of BET was found by Popescu et al. [36] to generate the most accurate BET extractions in MRIs with multiple sclerosis. BET|B,f0.1 was further shown to improve brain extractions over default BET in T1-MPRAGE MRIs with TBI by Lutkenhoff et al. [20]. Following the initial BET extraction, optiBET uses both linear and non-linear registration to map the extraction onto the MNI template space. Third, it inverts the projection of a standard MNI brain-only mask from the template space back to the subject’s original space. Finally, this brain mask is applied to the input image to mask out non-brain tissue, generating the final brain. The configuration of optiBET we tested was the default, which uses an MNI152 T1 1 mm brain-only mask for masking and performs the initial extraction with FSL using the described BET configuration [20].

2.5. Intensity Processing Procedures

For each of the 10 BET configurations, plus optiBET, 14 different intensity processing procedures (see Table 2) were implemented prior to extraction with the BET configuration. One inhomogeneity correction technique, Nonparametric Nonuniform Intensity Normalization (N3), and three intensity normalization techniques, Z-score, KDE-based, and WhiteStripe normalization, were compared. KDE-based and WhiteStripe intensity normalization both work exclusively with images of the brain and thus require an initial mask prior to the final brain extraction. For Z-score, KDE-based, and WhiteStripe normalization, we varied the pre-normalization mask by the fractional intensity threshold at the values f = 0.1, f = 0.5 (default), and f = 0.8. Additionally, we tested N3 inhomogeneity correction and Z-score normalization at one and four iterations. FreeSurfer performs N3 normalization through the MINC tool mri_nu_correct, which runs four iterations of N3 by default. These configurations are described in Table 2. The intensity-normalization package from Reinhold et al. [31] was used to execute the intensity normalization techniques. For N3, we utilized the MINC tool mri_nu_correct from Sled et al. [30].

2.5.1. N3 Inhomogeneity Correction

Nonparametric Nonuniform Intensity Normalization (N3), from Sled et al. [30], is a widely used technique for addressing intensity inhomogeneity in MRI images. N3 has been extensively employed in various neuroimaging studies and is often considered a standard preprocessing step in many MRI analysis pipelines. For instance, N3 is the default inhomogeneity correction method in popular neuroimaging software packages such as FreeSurfer [37] and is commonly used in large-scale, multi-site neuroimaging studies like the Alzheimer’s Disease Neuroimaging Initiative (ADNI) [38]. Moreover, N3 has been shown to effectively reduce intensity inhomogeneity in MRI images of various patient populations, including those with multiple sclerosis [39], Alzheimer’s disease [40], and brain tumors [41].

N3 addresses intensity inhomogeneity using a nonparametric approach without presupposing specific tissue distribution models or necessitating initial segmentation. Central to this method is the principle that a uniform tissue’s intensity histogram, unaffected by inhomogeneity, should ideally manifest as unimodal. Through iterative refinement, N3 enhances the sharpness of this histogram, targeting the elimination of the bias field distorting the true signal. Initially, N3 operates under the presumption that deviations from the expected unimodal distribution stem primarily from an underlying, smoothly varying bias field, attributing intensity variations to this inhomogeneity. The correction process unfolds iteratively, alternating between estimating the true tissue intensity distribution, assuming that the bias field is known, and estimating the bias field, assuming that the true tissue intensity distribution is known. Key to N3 is the deconvolution step, which separates the true signal from the bias field by conceptualizing the observed image intensities as the multiplication of the true signal by the bias field and then solving for the bias field in the frequency domain, where convolution simplifies to multiplication. This process is accompanied by the histogram sharpening of the image, aiming to counteract the bias field’s blurring effect by iteratively adjusting the image to achieve a more pronounced unimodal histogram. Notably, the N3 method is designed to be independent of the MRI pulse sequence and is robust to pathological data that might otherwise compromise model-based assumptions. The iterative process is designed to continue until subsequent adjustments yield negligible changes in histogram sharpness, indicating successful compensation for the bias field and the restoration of the true tissue intensities. This automated, minimal-assumption approach renders N3 a versatile and effective tool for enhancing tissue appearance uniformity across an MRI image.

2.5.2. Z-Score Intensity Normalization

Z-score (ZS) normalization is a technique designed to standardize the intensity values of MR images to a set range, making images from different subjects or scanning sessions comparable. By calculating the mean ${\mu }_{\mathrm{zs}}$ and standard deviation ${\sigma }_{\mathrm{zs}}$ of the intensities within the image, Z-score normalization adjusts each voxel value $I\left(x\right)$ as follows:

$$I_{\mathrm{zs}}\left(x\right)=\frac{I\left(x\right)-{\mu }_{\mathrm{zs}}}{{\sigma }_{\mathrm{zs}}}$$

(1)

This formula ensures that the resulting image $I_{\mathrm{zs}}$ has a mean intensity and a standard deviation set to one [31]. However, because Z-score normalization can generate an intensity distribution with negative intensity values, which violates image requirements for some brain extraction procedures, the distribution was shifted with fslmaths to have a minimum intensity value of zero after normalization [18].

2.5.3. KDE Intensity Normalization

Kernel Density Estimate (KDE)-based normalization advances the concept of intensity normalization by estimating the probability density function of the intensities within the brain mask using kernel density estimation with a Gaussian kernel. This method smooths the intensity histogram, allowing for the precise identification of peaks corresponding to different tissue types. The process is delineated as follows:

$$I_{\mathrm{kde}}\left(x\right)=c_3·\frac{I\left(x\right)}{\rho },$$

(2)

where $\rho$ represents the intensity peak associated with white matter, and $c_3=1000$ sets this peak to a standardized value post-normalization [31]. The KDE method not only smooths out fluctuations in the histogram but also enables the identification of tissue-specific intensity peaks. This capacity for detailed peak identification makes KDE-based normalization particularly effective in images with complex intensity distributions.

2.5.4. WhiteStripe Intensity Normalization

WhiteStripe (WS) normalization specifically targets the normalization of intensities of normal-appearing white matter (NAWM). This method isolates a band of intensities (the “white stripe”) centered around the peak intensity of NAWM, determined by analyzing the smoothed histogram of the image [39]. The normalization is applied as follows:

$$I_{\mathrm{ws}}\left(x\right)=\frac{I\left(x\right)-{\mu }_{\mathrm{ws}}}{{\sigma }_{\mathrm{ws}}},$$

(3)

where ${\mu }_{\mathrm{ws}}$ is the mean intensity of the white stripe, and ${\sigma }_{\mathrm{ws}}$ is its standard deviation [31]. Like Z-score normalization, WhiteStripe normalization can generate an intensity distribution with negative intensity values, so the distribution was shifted with fslmaths to have a minimum intensity value of zero after normalization as well.

2.6. Quality Assessment

To evaluate the efficacy of the automatic brain extraction pipelines, the brain-extracted image generated by these techniques can be compared to the manually corrected brain extraction of the same image. We employed the dice similarity coefficient (DSC) as our primary metric for comparison. The DSC is a statistical tool that can quantify the degree of similarity or overlap between two MRIs. A DSC score close to 1 indicates a high overlap between two images, while a score closer to 0 indicates very little overlap. Thus, the automated extraction method that produces the highest DSC on average with the gold-standard image is the most similar to it and, therefore, is considered the more accurate automated brain extraction technique. The dice similarity coefficient is calculated as follows:

$$DSC=\frac{2|X\cap Y|}{\left|X\right|+\left|Y\right|},$$

(4)

where X and Y are two sets representing a binary brain-extracted MRI. Fslmaths was used to execute the calculations of the dice similarity coefficients [18].

3. Results

Table 3. A ranked table of all 154 brain extraction pipelines by mean DSC scores for T1-MPRAGE and T2-FLAIR MRIs. Pipelines consist of an IP procedure from Table 2 and an Ex. Config. (extraction configuration) from Table 1. Pipelines at the top of the table are closer to DSC = 1, meaning they exhibit higher overlap with the gold standard and are considered more accurate. The SD columns report the standard deviation in DSC scores for each pipeline. Additionally, Figure 1 can be used to help visualize and compare brain extraction quality represented by different DSC scores.

T1-MPRAGE				T2-FLAIR
Brain Extraction Pipeline		DSC		Brain Extraction Pipeline		DSC
Ex. Config.	IP Procedure	Mean	SD	Ex. Config.	IP Procedure	Mean	SD
optiBET	N3 iter = 1	0.9780	0.0259	optiBET	N3 iter = 4	0.9557	0.0544
optiBET	-	0.9779	0.0455	optiBET	N3 iter = 1	0.9525	0.0676
optiBET	WS m = BET|f0.8	0.9769	0.0488	optiBET	ZS iter = 4	0.9475	0.0845
optiBET	ZS iter = 1	0.9766	0.0447	optiBET	-	0.9471	0.0882
optiBET	ZS m = BET|f0.8	0.9766	0.0471	optiBET	ZS iter = 1	0.9470	0.0874
optiBET	WS m = BET|f0.5	0.9766	0.0454	optiBET	WS m = BET|f0.1	0.9438	0.0926
optiBET	ZS iter = 4	0.9763	0.0462	optiBET	ZS m = BET|f0.8	0.9427	0.0944
optiBET	ZS m = BET|f0.1	0.9762	0.0450	optiBET	WS m = BET|f0.5	0.9410	0.1025
optiBET	KDE m = BET|f0.5	0.9761	0.0459	optiBET	KDE m = BET|f0.5	0.9409	0.1049
optiBET	N3 iter = 4	0.9759	0.0265	optiBET	WS m = BET|f0.8	0.9409	0.1030
optiBET	KDE m = BET|f0.1	0.9759	0.0460	optiBET	KDE m = BET|f0.8	0.9408	0.1052
optiBET	ZS m = BET|f0.5	0.9739	0.0545	optiBET	ZS m = BET|f0.1	0.9408	0.1025
optiBET	KDE m = BET|f0.8	0.9711	0.0668	optiBET	KDE m = BET|f0.1	0.9407	0.1057
optiBET	WS m = BET|f0.1	0.9690	0.0988	optiBET	ZS m = BET|f0.5	0.9404	0.1040
BET|B	N3 iter = 1	0.9391	0.0299	BET|B,f0.2,g0.3	N3 iter = 4	0.9384	0.0249
BET|R	N3 iter = 4	0.9383	0.0466	BET|R	N3 iter = 1	0.9343	0.0493
BET|B	N3 iter = 4	0.9382	0.0314	BET|B,f0.1	N3 iter = 4	0.9342	0.0293
BET|B	WS m = BET|f0.1	0.9368	0.0363	BET|B,f0.2,g0.3	N3 iter = 1	0.9333	0.0474
BET|R	N3 iter = 1	0.9352	0.0480	BET|R	KDE m = BET|f0.1	0.9324	0.0458
BET|B	WS m = BET|f0.8	0.9325	0.0652	BET|R	KDE m = BET|f0.5	0.9322	0.0460
BET|B	ZS m = BET|f0.8	0.9317	0.0389	BET|R	KDE m = BET|f0.8	0.9320	0.0464
BET|B	ZS m = BET|f0.1	0.9311	0.0381	BET|R	ZS iter = 1	0.9318	0.0463
BET|B	-	0.9306	0.0403	BET|R	ZS iter = 4	0.9318	0.0458
BET|B	ZS iter = 4	0.9294	0.0409	BET|R	ZS m = BET|f0.8	0.9317	0.0466
BET|B	ZS m = BET|f0.5	0.9292	0.0419	BET|R	ZS m = BET|f0.1	0.9316	0.0458
BET|B	ZS iter = 1	0.9286	0.0405	BET|R	ZS m = BET|f0.5	0.9315	0.0466
BET|B,f0.2,g0.3	N3 iter = 1	0.9262	0.0234	BET|R	WS m = BET|f0.5	0.9315	0.0464
BET|B,f0.2,g0.3	N3 iter = 4	0.9249	0.0256	BET|R	WS m = BET|f0.1	0.9314	0.0467
BET|B,f0.2,g0.3	-	0.9226	0.0330	BET|R	-	0.9312	0.0467
BET|B,f0.2,g0.3	ZS m = BET|f0.8	0.9223	0.0343	BET|R	WS m = BET|f0.8	0.9309	0.0461
BET|B,f0.2,g0.3	WS m = BET|f0.1	0.9223	0.0335	BET|B,f0.1	N3 iter = 1	0.9294	0.0491
BET|B,f0.2,g0.3	ZS m = BET|f0.5	0.9222	0.0341	BET|B,f0.2,g0.3	ZS iter = 4	0.9280	0.0510
BET|B,f0.2,g0.3	ZS m = BET|f0.1	0.9221	0.0342	BET|B,f0.2,g0.3	ZS iter = 1	0.9274	0.0518
BET|B,f0.2,g0.3	ZS iter = 4	0.9220	0.0340	BET|R	N3 iter = 4	0.9273	0.0616
BET|B,f0.2,g0.3	ZS iter = 1	0.9218	0.0343	BET|B,f0.1	ZS iter = 4	0.9255	0.0523
BET|B,f0.2,g0.3	KDE m = BET|f0.1	0.9207	0.0382	BET|B,f0.1	ZS iter = 1	0.9249	0.0527
BET|B,f0.2,g0.3	KDE m = BET|f0.5	0.9206	0.0378	BET|B,f0.2,g0.3	-	0.9247	0.0615
BET|B	WS m = BET|f0.5	0.9205	0.0810	BET|B,f0.2,g0.3	WS m = BET|f0.1	0.9244	0.0615
BET|B	KDE m = BET|f0.1	0.9193	0.0513	BET|B,f0.2,g0.3	ZS m = BET|f0.8	0.9235	0.0614
BET|B	KDE m = BET|f0.5	0.9182	0.0519	BET|B,f0.1	-	0.9227	0.0621
BET|B,f0.2,g0.3	WS m = BET|f0.8	0.9181	0.0622	BET|B,f0.2,g0.3	WS m = BET|f0.5	0.9220	0.0693
BET|R	WS m = BET|f0.5	0.9179	0.0694	BET|B,f0.1	WS m = BET|f0.1	0.9213	0.0620
BET|B,f0.2,g0.3	KDE m = BET|f0.8	0.9174	0.0502	BET|B,f0.2,g0.3	ZS m = BET|f0.5	0.9211	0.0694
BET|R	ZS iter = 4	0.9158	0.0729	BET|B,f0.2,g0.3	ZS m = BET|f0.1	0.9207	0.0694
BET|R	KDE m = BET|f0.1	0.9157	0.0731	BET|B,f0.1	ZS m = BET|f0.8	0.9206	0.0621
BET|R	ZS m = BET|f0.5	0.9156	0.0731	BET|B,f0.2,g0.3	KDE m = BET|f0.1	0.9203	0.0676
BET|R	ZS m = BET|f0.8	0.9155	0.0732	BET|B,f0.2,g0.3	KDE m = BET|f0.8	0.9202	0.0673
BET|R	ZS m = BET|f0.1	0.9155	0.0731	BET|B,f0.1	KDE m = BET|f0.1	0.9202	0.0676
BET|B	KDE m = BET|f0.8	0.9154	0.0625	BET|B,f0.2,g0.3	KDE m = BET|f0.5	0.9199	0.0672
BET|R	ZS iter = 1	0.9154	0.0732	BET	ZS m = BET|f0.8	0.9197	0.0518
BET|B,f0.2,g0.3	WS m = BET|f0.5	0.9152	0.0570	BET	KDE m = BET|f0.8	0.9196	0.0523
BET|R	KDE m = BET|f0.8	0.9150	0.0734	BET|B,f0.2,g0.3	WS m = BET|f0.8	0.9195	0.0709
BET|B,f0.1	N3 iter = 1	0.9142	0.0280	BET|B,f0.1	KDE m = BET|f0.8	0.9195	0.0673
BET|R	-	0.9133	0.0768	BET|B,f0.1	KDE m = BET|f0.5	0.9195	0.0673
BET|R	KDE m = BET|f0.5	0.9133	0.0769	BET	ZS iter = 4	0.9193	0.0525
BET|R	WS m = BET|f0.8	0.9132	0.0793	BET	KDE m = BET|f0.1	0.9193	0.0527
BET|R	WS m = BET|f0.1	0.9131	0.0771	BET	-	0.9193	0.0531
BET|B,f0.1	N3 iter = 4	0.9122	0.0301	BET	KDE m = BET|f0.5	0.9193	0.0534
BET|B,f0.1	ZS m = BET|f0.8	0.9121	0.0339	BET	ZS iter = 1	0.9191	0.0524
BET|B,f0.1	-	0.9121	0.0339	BET	WS m = BET|f0.8	0.9190	0.0538
BET|B,f0.1	KDE m = BET|f0.5	0.9121	0.0369	BET	ZS m = BET|f0.1	0.9190	0.0522
BET|B,f0.1	ZS m = BET|f0.5	0.9120	0.0344	BET|B,f0.1	WS m = BET|f0.5	0.9188	0.0696
BET|B,f0.1	KDE m = BET|f0.1	0.9120	0.0372	BET	WS m = BET|f0.5	0.9188	0.0536
BET|B,f0.1	ZS iter = 4	0.9119	0.0339	BET	WS m = BET|f0.1	0.9187	0.0535
BET|B,f0.1	ZS iter = 1	0.9118	0.0348	BET|B,f0.1	ZS m = BET|f0.1	0.9186	0.0700
BET|B,f0.1	ZS m = BET|f0.1	0.9117	0.0350	BET|B,f0.1	ZS m = BET|f0.5	0.9184	0.0700
BET|B,f0.1	WS m = BET|f0.1	0.9113	0.0348	BET	ZS m = BET|f0.5	0.9184	0.0537
BET	N3 iter = 4	0.9097	0.0639	BET	N3 iter = 1	0.9181	0.0549
BET|B,f0.1	KDE m = BET|f0.8	0.9090	0.0485	BET|B,f0.1	WS m = BET|f0.8	0.9159	0.0710
BET	N3 iter = 1	0.9075	0.0653	BET|g-0.3	KDE m = BET|f0.1	0.9051	0.0672
BET|B,f0.1	WS m = BET|f0.8	0.9073	0.0620	BET|g-0.3	-	0.9050	0.0672
BET|B,f0.1	WS m = BET|f0.5	0.9052	0.0533	BET|g-0.3	KDE m = BET|f0.5	0.9048	0.0678
BET|g-0.3	N3 iter = 4	0.9029	0.0537	BET|g-0.3	ZS iter = 4	0.9042	0.0680
BET|g-0.3	N3 iter = 1	0.9025	0.0560	BET|g-0.3	ZS m = BET|f0.1	0.9042	0.0683
BET	WS m = BET|f0.5	0.9000	0.0711	BET|g-0.3	WS m = BET|f0.5	0.9040	0.0682
BET|g-0.3	KDE m = BET|f0.5	0.8991	0.0527	BET|g-0.3	ZS m = BET|f0.8	0.9040	0.0685
BET|g-0.3	ZS m = BET|f0.8	0.8990	0.0527	BET|g-0.3	ZS iter = 1	0.9039	0.0683
BET|g-0.3	ZS iter = 1	0.8990	0.0530	BET|g-0.3	ZS m = BET|f0.5	0.9038	0.0678
BET|g-0.3	WS m = BET|f0.8	0.8989	0.0530	BET|g-0.3	WS m = BET|f0.1	0.9031	0.0698
BET|g-0.3	WS m = BET|f0.1	0.8989	0.0528	BET|g-0.3	N3 iter = 1	0.9020	0.0654
BET|g-0.3	-	0.8988	0.0538	BET|g-0.3	WS m = BET|f0.8	0.9011	0.0691
BET|g-0.3	ZS iter = 4	0.8987	0.0536	BET	N3 iter = 4	0.8980	0.0793
BET|g-0.3	ZS m = BET|f0.1	0.8986	0.0534	BET|B	ZS iter = 1	0.8963	0.0823
BET	WS m = BET|f0.8	0.8984	0.0759	BET|g-0.3	KDE m = BET|f0.8	0.8917	0.1263
BET|g-0.3	KDE m = BET|f0.1	0.8984	0.0531	BET|B	WS m = BET|f0.1	0.8917	0.1156
BET	ZS m = BET|f0.8	0.8983	0.0766	BET|B	WS m = BET|f0.5	0.8916	0.1048
BET	WS m = BET|f0.1	0.8979	0.0763	BET|B	WS m = BET|f0.8	0.8904	0.1098
BET|g-0.3	ZS m = BET|f0.5	0.8979	0.0536	BET|B	ZS iter = 4	0.8895	0.1030
BET	ZS iter = 1	0.8978	0.0769	BET|g-0.3	N3 iter = 4	0.8883	0.0736
BET	-	0.8978	0.0762	BET|B	N3 iter = 4	0.8870	0.1142
BET	ZS m = BET|f0.5	0.8977	0.0758	BET|B	ZS m = BET|f0.8	0.8857	0.1128
BET	KDE m = BET|f0.5	0.8977	0.0769	BET|B	ZS m = BET|f0.5	0.8844	0.1073
BET|g-0.3	KDE m = BET|f0.8	0.8976	0.0544	BET|B	N3 iter = 1	0.8815	0.1148
BET	ZS iter = 4	0.8976	0.0767	BET|B	ZS m = BET|f0.1	0.8811	0.1076
BET	ZS m = BET|f0.1	0.8975	0.0755	BET|B	-	0.8724	0.1211
BET	KDE m = BET|f0.1	0.8971	0.0766	BET|B	KDE m = BET|f0.5	0.8609	0.1182
BET|g-0.3	WS m = BET|f0.5	0.8970	0.0533	BET|B	KDE m = BET|f0.1	0.8609	0.1207
BET	KDE m = BET|f0.8	0.8964	0.0769	BET|B	KDE m = BET|f0.8	0.8604	0.1174
BET|f0.2	WS m = BET|f0.5	0.7776	0.1711	BET|f0.2	ZS iter = 4	0.8361	0.1394
BET|f0.2	WS m = BET|f0.8	0.7732	0.1702	BET|f0.2	KDE m = BET|f0.1	0.8361	0.1397
BET|f0.2	KDE m = BET|f0.5	0.7732	0.1704	BET|f0.2	-	0.8361	0.1398
BET|f0.2	WS m = BET|f0.1	0.7730	0.1704	BET|f0.2	KDE m = BET|f0.8	0.8359	0.1398
BET|f0.2	ZS iter = 4	0.7730	0.1703	BET|f0.2	WS m = BET|f0.1	0.8353	0.1400
BET|f0.2	ZS iter = 1	0.7729	0.1705	BET|f0.2	WS m = BET|f0.8	0.8352	0.1398
BET|f0.2	-	0.7727	0.1709	BET|f0.2	ZS m = BET|f0.8	0.8348	0.1406
BET|f0.2	KDE m = BET|f0.1	0.7726	0.1707	BET|f0.2	KDE m = BET|f0.5	0.8346	0.1412
BET|f0.2	ZS m = BET|f0.8	0.7724	0.1709	BET|f0.2	ZS m = BET|f0.1	0.8345	0.1410
BET|f0.2	ZS m = BET|f0.5	0.7717	0.1715	BET|f0.2	ZS m = BET|f0.5	0.8344	0.1420
BET|f0.2	KDE m = BET|f0.8	0.7715	0.1696	BET|f0.2	ZS iter = 1	0.8341	0.1421
BET|f0.2	ZS m = BET|f0.1	0.7704	0.1719	BET|f0.2	WS m = BET|f0.5	0.8341	0.1411
BET|f0.2	N3 iter = 1	0.7699	0.1724	BET|f0.8	KDE m = BET|f0.8	0.8333	0.0999
BET|f0.2	N3 iter = 4	0.7618	0.1700	BET|f0.8	WS m = BET|f0.1	0.8333	0.0995
BET|f0.8	WS m = BET|f0.8	0.7536	0.1135	BET|f0.8	KDE m = BET|f0.1	0.8332	0.1001
BET|f0.8	ZS m = BET|f0.1	0.7536	0.1136	BET|f0.8	-	0.8332	0.1001
BET|f0.8	ZS iter = 1	0.7532	0.1138	BET|f0.8	ZS m = BET|f0.5	0.8331	0.0997
BET|f0.8	-	0.7531	0.1138	BET|f0.8	KDE m = BET|f0.5	0.8331	0.1000
BET|f0.8	ZS iter = 4	0.7531	0.1142	BET|f0.8	ZS m = BET|f0.8	0.8330	0.1002
BET|f0.8	KDE m = BET|f0.5	0.7529	0.1146	BET|f0.8	WS m = BET|f0.8	0.8330	0.1002
BET|f0.8	KDE m = BET|f0.1	0.7529	0.1146	BET|f0.8	ZS m = BET|f0.1	0.8329	0.1002
BET|f0.8	WS m = BET|f0.1	0.7529	0.1145	BET|f0.8	WS m = BET|f0.5	0.8329	0.1003
BET|f0.8	ZS m = BET|f0.8	0.7528	0.1144	BET|f0.8	ZS iter = 1	0.8327	0.1003
BET|f0.8	ZS m = BET|f0.5	0.7528	0.1146	BET|f0.8	ZS iter = 4	0.8327	0.1002
BET|f0.8	KDE m = BET|f0.8	0.7495	0.1181	BET|g0.3	KDE m = BET|f0.8	0.8237	0.1029
BET|g0.3	WS m = BET|f0.5	0.7446	0.1794	BET|g0.3	KDE m = BET|f0.5	0.8236	0.1028
BET|g0.3	ZS m = BET|f0.1	0.7432	0.1819	BET|g0.3	KDE m = BET|f0.1	0.8235	0.1032
BET|g0.3	-	0.7430	0.1822	BET|g0.3	ZS m = BET|f0.1	0.8234	0.1032
BET|g0.3	ZS iter = 1	0.7429	0.1821	BET|g0.3	-	0.8233	0.1032
BET|g0.3	KDE m = BET|f0.1	0.7429	0.1823	BET|g0.3	ZS iter = 4	0.8232	0.1033
BET|g0.3	ZS m = BET|f0.8	0.7427	0.1824	BET|g0.3	WS m = BET|f0.1	0.8231	0.1032
BET|g0.3	WS m = BET|f0.1	0.7425	0.1822	BET|g0.3	ZS iter = 1	0.8231	0.1032
BET|g0.3	WS m = BET|f0.8	0.7425	0.1824	BET|g0.3	ZS m = BET|f0.5	0.8230	0.1034
BET|g0.3	ZS m = BET|f0.5	0.7425	0.1826	BET|g0.3	WS m = BET|f0.5	0.8230	0.1048
BET|g0.3	ZS iter = 4	0.7421	0.1829	BET|g0.3	ZS m = BET|f0.8	0.8227	0.1048
BET|g0.3	KDE m = BET|f0.5	0.7419	0.1831	BET|g0.3	WS m = BET|f0.8	0.8213	0.1046
BET|g0.3	KDE m = BET|f0.8	0.7414	0.1818	BET|f0.2	N3 iter = 1	0.8208	0.1508
BET|f0.8	WS m = BET|f0.5	0.7401	0.1229	BET|f0.8	N3 iter = 1	0.8194	0.1049
BET|g0.3	N3 iter = 1	0.7395	0.1810	BET|g0.3	N3 iter = 1	0.8094	0.1158
BET|f0.8	N3 iter = 1	0.7375	0.1203	BET|f0.8	N3 iter = 4	0.7888	0.1093
BET|f0.8	N3 iter = 4	0.7283	0.1145	BET|f0.2	N3 iter = 4	0.7794	0.1773
BET|g0.3	N3 iter = 4	0.7265	0.1782	BET|f0.2,g0.3	WS m = BET|f0.8	0.7672	0.1780
BET|f0.2,g0.3	WS m = BET|f0.5	0.7128	0.1612	BET|f0.2,g0.3	WS m = BET|f0.5	0.7670	0.1761
BET|f0.2,g0.3	KDE m = BET|f0.8	0.7077	0.1545	BET|f0.2,g0.3	-	0.7669	0.1770
BET|f0.2,g0.3	ZS iter = 4	0.7074	0.1552	BET|f0.2,g0.3	KDE m = BET|f0.8	0.7663	0.1765
BET|f0.2,g0.3	KDE m = BET|f0.5	0.7059	0.1565	BET|f0.2,g0.3	KDE m = BET|f0.5	0.7660	0.1773
BET|f0.2,g0.3	ZS m = BET|f0.1	0.7057	0.1562	BET|f0.2,g0.3	ZS iter = 4	0.7658	0.1773
BET|f0.2,g0.3	WS m = BET|f0.1	0.7057	0.1560	BET|f0.2,g0.3	KDE m = BET|f0.1	0.7657	0.1768
BET|f0.2,g0.3	ZS m = BET|f0.8	0.7054	0.1562	BET|f0.2,g0.3	ZS m = BET|f0.1	0.7654	0.1772
BET|f0.2,g0.3	ZS iter = 1	0.7045	0.1560	BET|f0.2,g0.3	WS m = BET|f0.1	0.7653	0.1777
BET|f0.2,g0.3	WS m = BET|f0.8	0.7043	0.1567	BET|f0.2,g0.3	ZS m = BET|f0.8	0.7651	0.1775
BET|f0.2,g0.3	-	0.7042	0.1560	BET|f0.2,g0.3	ZS iter = 1	0.7645	0.1774
BET|f0.2,g0.3	ZS m = BET|f0.5	0.7042	0.1561	BET|f0.2,g0.3	ZS m = BET|f0.5	0.7632	0.1784
BET|f0.2,g0.3	KDE m = BET|f0.1	0.7034	0.1562	BET|g0.3	N3 iter = 4	0.7631	0.1424
BET|f0.2,g0.3	N3 iter = 1	0.7027	0.1523	BET|f0.2,g0.3	N3 iter = 1	0.7421	0.1943
BET|f0.2,g0.3	N3 iter = 4	0.6948	0.1470	BET|f0.2,g0.3	N3 iter = 4	0.6920	0.2104

ranks and summarizes the performance of the 14 optiBET and 140 BET brain extraction pipelines with T1-MPRAGE and T2-FLAIR MRIs in terms of the mean DSC. Referencing the DSC scores of the brain extractions in Figure 1 can be helpful for imagining the degree of inaccuracies in a brain extraction based on the mean DSC scores reported in Table 3. There was significant variability in performance across extraction configurations, intensity processing procedures, and MRI sequences. Pipelines using optiBET generated more accurate extractions than pipelines extracting with a configuration of BET for all intensity processing procedures and both MRI types. Of these top pipelines, the highest-performing pipelines for both image sequences were those employing optiBET with either one iteration of N3 inhomogeneity correction (denoted as N3 iter = 1) or four iterations (denoted as N3 iter = 4). Specifically, optiBET–N3 iter = 1 generated the most accurate extractions for T1-MPRAGE MRIs, with a mean DSC score of 0.9780. Similarly, for the T2-FLAIR images, optiBET–N3 iter = 4 was the highest-performing pipeline (mean DSC = 0.9557). These results are reflected in the first row of Table 3. It should further be noted that for T1-MPRAGE, the optiBET control pipeline was the second-highest-performing pipeline, barely falling short of optiBET–N3 iter = 1 (mean DSC = 0.9779), while, in T2-FLAIRs, optiBET–N3 iter = 4 showed a more substantial increase from the performance of its control pipeline (mean DSC = 0.9471).

While Table 3 generally exhibits a smooth, continuous decline in mean DSC between pipelines, there are notable points of discontinuity where a pronounced decrease in mean DSC is observed between two neighboring pipelines in the table. For the T1-MPRAGE MRI, the first discontinuity is observed between the poorest-performing optiBET pipeline, optiBET – WS m = BET|f0.1, which produced a mean DSC score of 0.9690, and the highest-performing BET pipeline, BET|B – N3 iter = 1, which produced a mean DSC score of 0.9391 (p < 0.001). In T2-FLAIR images, this gap between optiBET and BET performance was not observed. The lowest-performing optiBET pipeline in T2-FLAIR images was OptiBET – ZS m = BET|f0.5 with a mean DSC of 0.9404, and the highest-performing BET pipeline was BET|B,f0.2,g0.3 – N3 iter = 4 with a mean DSC of 0.9384. Figure 2 and Figure 3 visualize the results summarized in Table 3, displaying pipeline performance based on DSC scores by percentile rank for T1-MPRAGE and T2-FLAIR MRIs, respectively, with an emphasis on the lower end of the DSC score distribution. Subfigures (a) through (k) are ranked from lowest to highest by the mean DSC score of all pipelines with the same extraction configuration.

3.1. Effect of BET Parameter Configurations

The best-performing BET pipeline in terms of the mean DSC was BET|B – N3 iter = 1 for T1-MPRAGE and BET|B,f0.2,g0.3 – N3 iter = 4 for T2-FLAIR images. While previous investigations into BET parameter optimization have found that BET|B,f0.1 has outstanding performance in the presence of other kinds of pathology, as shown in Figure 2, we found that BET|B,f0.1 performs worse than BET|B, BET|R, and BET|B,f0.2,g0.3 for T1-MPRAGE MRIs, regardless of the intensity processing procedure [36]. The top BET|B pipeline (BET|B – N3 iter = 1) resulted in a mean DSC score of 0.9391, the mean DSC score of the top BET|R pipeline (BET|R – N3 iter = 4) was 0.9383, and the top BET|B,f0.2,g0.3 pipeline (BET|B,f0.2,g0.3 – N3 iter = 1) resulted in a mean DSC score of 0.9262, while the top BET|B,f0.1 pipeline (BET|B,f0.1 – N3 iter = 1) had a mean DSC of 0.9142 (differences were p < 0.001). In contrast, BET|B,f0.1 performed markedly better in T2-FLAIR. The top BET|B,f0.1 pipeline, BET|B,f0.1 – N3 iter = 4, produced the third highest DSC score among all BET pipelines for T2-FLAIR images, with a mean DSC = 0.9342. However, BET|B,f0.2,g0.3 – N3 iter = 4 and BET|R – N3 iter = 1 both still perform marginally better, producing DSC scores of 0.9384 and 0.9343, respectively.

Another notable difference in pipeline performance between T1MPRAGE and T2FLAIR is seen in BET|B, which performs significantly worse with T2-FLAIR images. As Figure 3 shows, BET|B is the only BET configuration with a primary option (“B” or “R”) to perform worse than pipelines without a primary option. Specifically, the top default BET pipeline (BET – ZS m = BET|f0.8) and the top BET|g−0.3 pipeline (BET|g−0.3 – KDE m = BET|f0.1) had mean DSC scores of 0.9197 and 0.9051, respectively, while the highest-performing BET|B pipeline (BET|B – ZS iter = 1) had a mean DSC score of only 0.8963.

For both MRI types, default BET and BET|g−0.3 had distinctly better performance relative to other BET configurations with no primary option—namely, BET|g0.3, BET|f0.8, and BET|f0.2,g0.3. The highest-performing pipelines using default BET were BET – N3 iter = 4 for T1-MPRAGE images, with a mean DSC of 0.9097, and BET – ZS iter = 4 for T2-FLAIR images, with a mean DSC of 0.9193. For BET|g−0.3, the highest-performing pipelines were BET|g−0.3 – N3 iter = 4 for T1-MPRAGE images and BET|g−0.3 – KDE m = BET|f0.1 for T2-FLAIR images, which had mean DSC scores of 0.9029 and 0.9051, respectively. On the other hand, BET|g0.3, BET|f0.8, and BET|f0.2,g0.3 exhibited worse performance. In Table 3, the top BET|f0.2 marks the second point of discontinuity for T1-MPRAGE images and the first major point of discontinuity for T2-FLAIR MRIs, representing a drastic decrease in pipeline performance from the preceding pipelines for both MRI sequences. For T1-MPRAGE, the poorest-performing default BET pipeline, BET – KDE m=BET|f0.8, has a mean DSC score of 0.8964 and precedes the top BET|f0.2 pipeline, BET|f0.2 – ZS iter = 4, in which the mean DSC plummets to 0.7776. This discontinuity can also be observed in the difference between the slopes of the pipelines for each extraction configuration in Figure 2.

An analogous but less drastic drop in mean DSC scores occurs at this point for T2-FLAIR images as well, where the lowest-performing BET|B pipeline (BET|B – KDE m = BET|f0.8) precedes the top BET|f0.2 pipeline (BET|f0.2 – ZS iter = 4), with the mean DSC decreasing from 0.8604 to 0.8361. In contrast to the T1-MPRAGE results, this was not the largest point of discontinuity for the T2-FLAIR images. Rather, the largest singular differential for T2-FLAIR images was observed between the two lowest-scoring pipelines in Table 3, drop** from 0.7421 to 0.6920 between BET|f0.2,g0.3 – N3 iter = 1 and BET|f0.2,g0.3 – N3 iter = 4. These two pipelines also have the lowest performance for the T1-MPRAGE images; however, the mean DSC only decreased from 0.7027 to 0.6948 between them.

3.2. Effect of Intensity Processing Procedures

While there was generally greater variability between extraction configurations, intensity processing procedures had varying effects between extraction configurations and other IP procedures, especially at the lower end of the DSC score distribution, as shown in Figure 2 and Figure 3. In this regard, N3 intensity inhomogeneity correction had a markedly different impact on brain extraction than KDE-based, WhiteStripe, and Z-score intensity normalization, which had more idiosyncratic and less significant effects. Figure 4 and Figure 5 show the impact of the intensity inhomogeneity correction and intensity normalization techniques on the mean DSC score of each extraction configuration relative to the control pipeline. For T1-MPRAGE images, BET|B,f0.2,g0.3, BET|B,f0.1, BET|g−0.3, BET|f0.2, BET|f0.8, and BET|f0.2,g0.3 cluster around similar mean DSC values regardless of intensity processing, while other BET configurations, like BET|B, BET|R, and BET, show larger changes in mean DSC depending on the intensity processing procedure.

3.2.1. Intensity Inhomogeneity Correction

N3 intensity homogeneity correction had a notable effect on the variance of the pipelines into which it was incorporated. While N3 offered virtually no improvement in the mean DSC of optiBET in T1-MPRAGE images, it did reduce the standard deviation by roughly half, decreasing it from 0.0455 to 0.0258 for N3 iter = 1 and 0.0265 for N3 iter = 4. This effect is also represented in Figure 2, which shows that pipelines with N3 have higher DSC scores in the bottom 10th percentile of DSC scores relative to control pipelines and other IP procedures for optiBET and several BET extraction configurations. Specifically, this effect is observed in the BET, BET|R, BET|B, BET|B,f0.2,g0.3, and BET|B,f0.1 extraction configurations (Figure 2g–k). For these configurations, N3 iter = 1 showed a robust decrease in the standard deviation relative to the control pipeline, while N3 iter = 4 showed a smaller reduction in the standard deviation (Table 3). However, N3 does not have this effect in the five lowest-performing BET configurations in Figure 2.

In addition to the differences in variance shown in Table 3 and Figure 2, pipelines incorporating N3 homogeneity correction had the largest effect on the mean DSC from their respective control pipelines for T1-MPRAGE MRIs. As Figure 4 shows, N3 iter = 1 and N3 iter = 4 both significantly increased DSC relative to the control pipeline across the default BET, BET|R, and BET|B configurations (p < 0.001). Notably, N3 iter = 4 had a significantly larger positive effect than N3 iter = 1 on default BET (p < 0.02) and BET|R (p < 0.001). For BET|g(-0.3), both N3 iterations showed significant improvements, with p-values < 0.02 for N3 iter = 1 and p < 0.05 for N3 iter = 4. Similarly, both iterations increased the DSC for BET|B,f0.1, with N3 iter = 1 outperforming N3 iter = 4 (p < 0.01). In the case of BET|B,f0.2,g0.3, N3 iter = 1 had a more substantial positive effect compared to N3 iter = 4 (p < 0.05). Conversely, for optiBET, N3 iter = 1 exhibited a significant positive effect, whereas N3 iter = 4 had a slight negative impact (p < 0.001). Similarly, for BET|f0.2, BET|f0.2,g0.3, and BET|f0.8, N3 iter = 4 demonstrated a significant decrease in DSC relative to the control pipelines and showed a larger decrease in DSC compared to N3 iter = 1 (p < 0.001 for BET|f0.2 and BET|f0.2,g0.3; p < 0.001 for BET|f0.8). Lastly, for the BET|g0.3 configuration, N3 iter = 1 had a significant negative effect, and so did N3 iter = 4 (p < 0.001), with N3 iter = 4 having a significantly larger negative effect than N3 iter = 1 (p < 0.001) (Figure 4).

For T2-FLAIR images, N3 shows a pattern similar to but distinct from its effects on T1-MPRAGE images (represented in Figure 5). The decreasing effect of N3 on the standard deviation is similarly observed with T2-FLAIR images (Figure 5), with N3 reducing the standard deviations of BET|B, BET|B,f0.1, BET|B,f0.2,g0.3, and optiBET by over 50 percent in some cases. However, opposite to the results for T1-MPRAGE images, reductions in the SDs of these configurations were larger with N3 iter = 4 than with N3 iter = 1. Moreover, the addition of N3|iter = 4 to default BET and BET|R (Figure 5f,j) resulted in increases in SDs from 0.0531 and 0.0467 to 0.0793 and 0.0616, respectively, rather than decreases, as was the case for the T1-MPRAGE MRIs.

For BET|B and BET|B,f0.2,g0.3, both iterations significantly increased the DSC relative to the control pipeline (p < 0.05). However, for BET|B,f0.1 and optiBET, the increases produced by both iterations were not statistically significant. Interestingly, in the BET|R configuration, N3 iter = 1 demonstrated a non-significant increase in DSC, while N3 iter = 4 showed a non-significant decrease; these were the highest- and lowest-performing pipelines, respectively, for BET|R. For BET configurations without a primary option, N3 iter = 4 consistently showed a larger negative effect than N3 iter = 1. Default BET, BET|f0.2, BET|f0.2,g0.3, BET|f0.8, BET|g(-0.3), and BET|g0.3 resulted in significant decreases in DSC relative to the control pipeline. For default BET and BET|g(-0.3), only N3 iter = 4 led to significant declines (p < 0.01), while N3 iter = 1 showed non-significant decreases. Conversely, in BET|f0.2, BET|f0.2,g0.3, BET|f0.8, and BET|g0.3, both iterations resulted in significant decreases, with each configuration having p-values less than 0.01 for N3 iter = 1 and less than 0.001 for N3 iter = 4. Across these configurations, the negative impact of N3 iter = 4 on the DSC was notably larger than that of N3 iter = 1 (p < 0.001 for each).

3.2.2. Intensity Normalization

In comparison to N3 inhomogeneity correction, Z-score, KDE-based, and WhiteStripe intensity normalization demonstrated more idiosyncratic effects on the pipelines for T1-MPRAGE MRIs. However, significant positive differences were also observed. KDE-based normalization with m = BET|f0.5 increased the DSC for BET|f0.2,g0.3, BET|f0.2, and BET|g-0.3 configurations (p < 0.05) (Figure 4). Additionally, KDE m = BET|f0.8 increased the DSC for BET|f0.2,g0.3 (p < 0.05). WhiteStripe normalization with m = BET|f0.1 increased the DSC for BET|B and BET|f0.2,g0.3 (p < 0.001 and p < 0.05, respectively). WS m = BET|f0.5 increased the DSC for BET|f0.2,g0.3, BET|f0.2, BET|R, and BET, with significance ranging from p < 0.05 to p < 0.01. WS m = BET|f0.8 increased the DSC for BET|f0.2 and BET|f0.8 (p < 0.05). Z-score normalization with ZS iter = 1 increased the DSC for BET|g-0.3 (p < 0.05). ZS iter = 4 increased the DSC for BET|f0.2,g0.3, p < 0.05. ZS m = BET|f0.1 increased the DSC for BET|f0.8, BET|B, and BET|f0.2,g0.3, with significance from p < 0.05 to p < 0.01. Lastly, ZS m = BET|f0.8 increased the DSC for BET|B and BET|f0.2,g0.3, with p < 0.01 and p < 0.05, respectively (Figure 4).

Negative effects were also observed with different normalization procedures in T1-MPRAGE MRIs. KDE m = BET|f0.1 resulted in a decrease in DSC for BET|B (p < 0.001), BET|B,f0.2,g0.3 (p < 0.01), and optiBET (p < 0.05). Additionally, KDE m = BET|f0.5 led to a decrease in DSC for BET|B (p < 0.01) and optiBET (p < 0.02), while KDE m = BET|f0.8 caused a decrease in DSC for BET (p < 0.001), BET|B (p < 0.01), and optiBET (p < 0.02) (Figure 4). WhiteStripe normalization also showed detrimental effects; WS m = BET|f0.1 resulted in a decrease in DSC for optiBET (p < 0.001) and BET|B (p < 0.001). WS m = BET|f0.5 led to a decrease in DSC for BET|B (p < 0.001), optiBET (p < 0.001), BET|B,f0.1 (p < 0.01), and BET|g0.3 (p < 0.05). Z-score normalization was not exempt from negative impacts either; ZS iter = 1 resulted in a decrease in DSC for optiBET (p < 0.001) and BET|B,f0.2,g0.3 (p < 0.01). ZS iter = 4 led to a decrease in DSC for optiBET (p < 0.001) and BET|g0.3 (p < 0.05). Additionally, ZS m = BET|f0.1 caused a decrease in DSC for optiBET (p < 0.001), BET|B,f0.2,g0.3 (p < 0.01), and BET (p < 0.05), while ZS m = BET|f0.5 and ZS m = BET|f0.8 each led to a decrease in DSC for BET|g0.3 (p < 0.05) (Figure 4).

For T2-FLAIR MRIs, there were fewer significant positive effects. KDE m = BET|f0.8 demonstrated an increase in DSC with a p-value of less than 0.05. For the BET|g-0.3 configuration, an increase was also observed with KDE m = BET|f0.1, showing a significant improvement (p < 0.001) (Figure 5). In the BET|B configuration, several normalization techniques resulted in notable increases in DSC: WS m = BET|f0.1 and WS m = BET|f0.5 both achieved increases with p < 0.001 and WS m = BET|f0.8 with p < 0.01. Additionally, ZS iter = 1 and ZS iter = 4 in the BET|B configuration led to significant increases in DSC, both p < 0.001. Further enhancements in the same configuration were seen with ZS m = BET|f0.8, ZS m = BET|f0.5, and ZS m = BET|f0.1 (all p < 0.01) (Figure 5).

Negative effects were also observed across different procedures for the T2-FLAIR images. In the BET|B configuration, KDE m = BET|f0.8 and KDE m = BET|f0.5 both led to a decrease in DSC, each with (p < 0.05) (Figure 5). Similarly, in the BET|f0.2,g0.3 configuration, KDE m = BET|f0.1 and KDE m = BET|f0.5 resulted in decreases in DSC with p < 0.01 and p < 0.05, respectively. Moreover, KDE m = BET|f0.5 in the BET|g-0.3 configuration showed a significant decrease in DSC (p < 0.001). WhiteStripe normalization also demonstrated detrimental effects; in the BET|B,f0.2,g0.3 configuration, WS m = BET|f0.1 and WS m = BET|f0.5 both resulted in decreases in DSC, with p < 0.01 and p < 0.05, respectively. Additionally, WS m = BET|f0.8 in the BET|f0.2 configuration and WS m = BET|f0.5 in the BET|f0.8 configuration each led to decreases in DSC with p < 0.05. WS m = BET|f0.1 and WS m = BET|f0.5 in the BET|f0.2,g0.3 configuration and WS m = BET|f0.5 in the BET|g-0.3 configuration also showed significant decreases (both p < 0.001). Z-score normalization exhibited similar negative trends; in the BET|B,f0.2,g0.3 configuration, ZS m = BET|f0.1 and ZS m = BET|f0.5 led to decreases in DSC with p < 0.02 and p < 0.05, respectively. Further, ZS procedures with the BET|f0.2 and BET|f0.2,g0.3 configurations showed decreases in mean DSC, with p-values ranging from p < 0.05 to p < 0.02. In the BET|g-0.3 configuration, ZS iter = 4, ZS m = BET|f0.1, ZS m = BET|f0.5, and ZS iter = 1 each resulted in decreases in DSC (all p < 0.01) (Figure 5).

4. Discussion

This study investigated the impact of different intensity processing procedures and brain extraction configurations on the efficacy of automated brain extraction from MRI images of human subjects suffering from traumatic brain injury. In our analyses, we primarily use the mean dice similarity coefficient to evaluate the accuracy of brain extraction pipelines. While other studies evaluating automated brain extraction methods generally report the median, we deliberately report the mean. The frequency distribution of DSC scores across a set of images typically exhibits a left skew, particularly in high-performing brain extraction pipelines. The performance of these pipelines does not decrease linearly but rather tends to drop sharply. Thus, a portion of the data consist of outliers, where the extraction results are exceptionally poor. These very poor extractions have an outsized impact on preprocessing efficiency, as they require extensive manual corrections. Extensive manual corrections, even if only for a few dozen images, necessitate the recruitment of significant auxiliary resources. Large datasets, such as EpiBios4Rx’s, can require the recruitment and training of volunteer researchers. Thus, the poorest extractions introduce the most significant logistical complexities and, in turn, are the greatest burden on the efficiency of the TBI preprocessing pipeline. Reporting the median would obscure these critical outliers, providing a less complete picture of the pipeline’s overall performance. On the other hand, using the mean ensures that these critical cases are represented in our results.

Our comparison of brain extraction pipelines with optiBET and BET configurations on T1-MPRAGE and T2-FLAIR MRIs produces several critical findings with significant implications for the analysis of patient MRI data, particularly for the precise and efficient preprocessing of subjects with traumatic brain injury. As expected, we found that optiBET performed superior to BET in all pipelines for both MRI types. We found notable differences among the ten different BET configurations, with several results contradicting and significantly adding to previous findings. For both sequence types, we found that BET configurations with a primary option (“B” or “R”) almost always outperformed configurations without a primary option. For T1-MPRAGE specifically, we found that BET|B was the most accurate BET pipeline, followed by BET|R and BET|B,f0.2,g0.3. This contrasts with Popescu et al. [36], who found that BET|B,f0.1 was the highest-performing pipeline for 3D T1 images in multiple sclerosis. Possible reasons for this discrepancy include differences in pathology, the MRI acquisition protocol, procedures for generating gold-standard extractions, and evaluation with the mean vs. median (however, using the median would not change the relative performance of the pipelines for our data). Additionally, because Popescu et al. [36] did not evaluate BET configurations with alterations to option “g”, BET|B,f0.2,g0.3 would not be accounted for. The divergence in optimal pipeline performance underscores the potential need for tailored preprocessing approaches for TBI data. These findings further prompt a reevaluation of assumptions on the efficacy of preprocessing techniques, which are often based on their performance on datasets with different pathologies. This suggests that a “one-size-fits-all” approach to preprocessing patient MRI data may be inappropriate, particularly in complex cases such as TBI, where a technique’s effectiveness should be rigorously validated against the specific challenges presented by the pathology. Alternatively, generalizations about technique performance on patient data may be viable, but these assumptions can only be made when they are sourced from investigations on pathologies of equal or greater complexity. Given that TBI is regarded as one of the most complex brain pathologies, the development of TBI-optimized MRI preprocessing methods could be key to highly efficient and effective generalized preprocessing pipelines for patient MRI data on the whole.

Concerning the effect of the intensity processing techniques, N3 inhomogeneity correction stands out as a significantly beneficial technique for improving the accuracy of brain extractions—especially in the more accurate configurations incorporating options “B” or “R”. This is evidenced by the higher mean DSCs and reduced standard deviations observed in pipelines using N3 relative to their controls. OptiBET combined with N3 inhomogeneity correction achieved the highest DSC scores for the BET and optiBET pipelines, specifically optiBET – N3 iter = 1 for T1-MPRAGE and optiBET – N3 iter = 4 for T2-FLAIR. While the optiBET pipeline with one iteration of N3 showed only a slight advantage over its control, for T2-FLAIR, optiBET with four iterations of N3 significantly outperformed the control. Similarly, the BET pipelines with option “B” (BET|B, BET|B,f0.1, and BET|B,f0.2,g0.3) showed greater accuracy with one iteration of N3 relative to four iterations for the T1-MPRAGE MRIs but showed greater accuracy with four iterations of N3 for T2-FLAIR images. This suggests that small, precise adjustments to the bias field can significantly influence the accuracy of brain extractions, potentially enhancing the reliability of TBI-related morphometric analyses if applied appropriately. The relevance of these findings could extend to various aspects of TBI processing and analysis research. A detailed understanding of how different inhomogeneity correction iterations affect extraction outcomes can guide researchers in selecting the most appropriate settings for their specific needs. This is particularly valuable in TBI cases, where brain morphology can be drastically altered, and conventional methods are less effective.

4.1. Limitations

A potential limitation of this investigation’s findings is the manually corrected, gold-standard brain extractions. As previously described, these brain extractions are a benchmark used to evaluate the accuracy of automated brain extraction pipelines. Our gold-standard extractions were generated from an initial extraction by optiBET or BET|B,f0.4 and then manually corrected using various techniques until they were of sufficient quality. Therefore, DSC scores could be skewed in favor of pipelines that are more similar to the optiBET and BET|B,f0.4 control pipelines. However, because subsequent manual corrections were made for all images, the impact of this bias on our results should be negligible. Regardless, we would like to acknowledge the potential for this bias and address its two consequences on the interpretation of this study’s findings. First, as evidenced by the inter-rater reliability issues associated with manual tracings of the brain, specific brain structures, and brain lesions from TBI, there is no single “correct” brain extraction [16,42]. Rather, there is a range of acceptable brain extractions that can all be considered correct. Consequently, any change in the extracted brain resulting from the addition of an intensity processing procedure could be interpreted as a decrease in accuracy, while in actuality, it could simply reflect a difference in brain extraction that is not decidedly better or worse. Second, significant increases in a pipeline’s accuracy from the addition of an intensity processing procedure are especially noteworthy. This is because the pipeline is showing enhanced accuracy despite bias possibly favoring its control, suggesting that improvements in accuracy from the added intensity processing were precisely in line with the manual corrections or were robust enough to transcend the inherent bias potentially in the gold-standard extractions.

4.2. Future Work

This analysis sets the groundwork for many potential future investigations. Our results suggest that N3 intensity inhomogeneity correction can significantly influence the accuracy and reliability of brain extraction outcomes; however, these effects are highly dependent on the specific procedure of N3 implemented, the MRI sequence, and the brain extraction configuration. While this exploration limited N3 to one and four iterations, it revealed a complex relationship between intensity inhomogeneity correction and brain extraction. This relationship is particularly pronounced in the contrast in outcomes between T1-MPRAGE and T2-FLAIR MRIs and BET configurations. Given this variability, future work should delve into a more granular investigation of intensity inhomogeneity correction beyond the current scope. This could include expanded testing of different N3 iterations using an initial brain mask with N3 (like those required for WS and KDE normalization) or performing intensity normalization after intensity inhomogeneity correction and vice versa. These approaches could enhance the uniformity and contrast of MRI images, potentially improving the accuracy of the subsequent brain extraction. Furthermore, BET|B,f0.1 has long been accepted as the best-performing BET configuration for patient data, such that BET|B,f0.1 is the configuration for the initial extraction in the optiBET script [20]. However, our results suggest that other configurations, namely, BET|R and BET|B,f0.2,g0.3, are better suited for TBI data. This could warrant a more systematic investigation into the optimal BET parameters for TBI data and how incorporating these different configurations into the optiBET script affects brain extraction.

Altogether, the current and proposed investigations contribute to the creation of TBI-optimized brain extraction methods and, ultimately, TBI-optimized MRI preprocessing pipelines, ultimately enhancing the accuracy and efficiency of TBI MRI data analysis. This study has yielded new insights into the interplay between intensity normalization, inhomogeneity correction, and brain extraction. Our findings offer TBI researchers critical guidance for constructing MRI preprocessing pipelines that are better equipped to handle the complex and challenging nature of TBI data. By understanding these interactions, researchers can enhance the accuracy and reliability of brain extraction outcomes, ultimately improving the quality of data available for further TBI-related analyses.