Neuroimaging Analysis

Neuroimaging analysis turns brain scans, mostly MRI, into measurements: the volume of the hippocampus or of white matter lesions, the thickness of the cortex at each point of its surface, the response of each voxel to a task, the integrity of white matter tracts. Most of these measurements are made not by a person looking at the images but by long automated pipelines in which nearly every step is a vision task: removing everything that is not brain, correcting intensity shading, aligning each brain to a common template, labelling tissues and structures in 3D, reconstructing the cortical surface. The field built some of the most widely used software in medical image computing, and its large shared datasets, ADNI and UK Biobank among them, are where many methods for medical image segmentation and classification were developed and where their failures were found.

Problem

The questions are mostly about groups and about change. Does cortical thickness differ between patients and controls? How fast does the brain shrink in one condition compared with normal ageing? Which regions respond to a stimulus? Which differences between individuals in brain structure relate to differences in cognition? Each question needs the same measurement made consistently on hundreds or thousands of scans from many scanners and sites, and often on the same person over years. Manual delineation cannot scale to that, and it is not consistent between raters either; automated pipelines can, provided their errors do not depend on the scanner, the site or the group being studied.

The software this article describes is built for research. FreeSurfer’s licence, for example, states that the software is for research purposes only, has not been reviewed or approved by the FDA or any other agency, and that clinical applications are neither recommended nor advised (FreeSurfer licence, version 1.0). Clinical neuroradiology and its regulated software are a different setting, with the questions of clinical validation that apply to any medical AI.

Visual Data

Structural MRI. The workhorse is a 3D T1-weighted volume at about 1 mm isotropic resolution, in which white matter is brighter than grey matter and cerebrospinal fluid (CSF) is dark; MRI explains why the contrast depends on the sequence. UK Biobank’s brain protocol, run on a 3 T Siemens Skyra in about 35 minutes, gives the flavour of a modern research protocol: a T1-weighted MPRAGE at 1.0 × 1.0 × 1.0 mm (4 min 54 s), a T2-weighted FLAIR at 1.0 × 1.0 × 1.05 mm, susceptibility-weighted imaging, diffusion MRI at 2 mm and resting and task functional MRI at 2.4 mm (Littlejohns et al., 2020, Table 2). FLAIR suppresses the CSF signal and shows white matter lesions as bright; susceptibility weighting shows iron and small bleeds.

Functional MRI is a 4D series: a whole-brain volume every repetition time (735 ms in UKB’s protocol) for several minutes, whose signal follows blood-oxygenation-level-dependent (BOLD) contrast. A voxel’s time course, not its brightness, carries the information. Diffusion MRI acquires many volumes, each sensitized to water diffusion along one direction and at one strength, the b-value; UKB acquires 50 directions on each of two shells, at b = 1000 and 2000 s/mm² (Littlejohns et al., 2020). Diffusion and functional volumes are usually read out with echo-planar imaging, which MRI shows to be distorted along the phase-encoding direction by field inhomogeneity.

Formats and organization. Scanners store DICOM; research pipelines convert to NIfTI, whose header carries the voxel-to-world transform (Medical Image Preprocessing explains both). The Brain Imaging Data Structure (BIDS) standardizes how a study’s files are named and arranged: NIfTI images, JSON sidecar files with acquisition metadata and tab-separated tables of participants and events, in one directory per subject, with names such as sub-01_T1w.nii.gz (Gorgolewski et al., 2016). Its authors developed it because no widely adopted standard existed for organizing MRI data, which made sharing and automated processing difficult. A head MRI also contains the face, so shared structural scans are usually defaced: UKB’s default T1 and T2 files are defaced, and in 99.5% of cases the defacing mask does not overlap the brain mask (Littlejohns et al., 2020).

Vision Tasks Involved

• Segmentation: brain versus non-brain (skull stripping), tissue classes, subcortical structures, lesions and tumors (medical image segmentation, tumor segmentation).
• Registration: of a person’s scans to each other (functional to structural, follow-up to baseline) and of every person to a template or atlas (medical image registration).
• Intensity correction: bias-field correction and normalization (medical image preprocessing).
• Surface reconstruction: meshes of the grey–white and grey–CSF boundaries, from which cortical thickness and surface-based registration follow (polygon meshes are the representation).
• Classification and regression from images: diagnosis, prognosis or age prediction (medical image classification).
• Detection and quality control: lesions, artefacts, failed processing steps.

The Structural Pipeline

FreeSurfer’s citation page summarizes its own stream, and the summary reads as a list of the standard steps: motion correction and averaging of several T1 volumes, removal of non-brain tissue, an automated Talairach transformation, segmentation of subcortical white matter and deep grey matter structures, intensity normalization, tessellation of the grey–white boundary, topology correction, and surface deformation to place the grey–white and grey–CSF borders, followed by surface inflation, registration to a spherical atlas and cortical parcellation (FreeSurferMethodsCitation). Other suites order the steps differently or merge them, but the same operations recur.

Skull stripping

Brain extraction removes the scalp, skull, eyes and neck so that later steps see only brain and the CSF around it. Methods differ in how they find the boundary:

• Deformable surfaces. FSL’s BET deletes non-brain tissue from an image of the whole head by fitting a surface mesh that starts as a sphere at the head’s centre of gravity, with half the estimated head radius; its main parameter is a fractional intensity threshold between 0 and 1 (default 0.5), where smaller values give larger brain outlines (FSL documentation, BET). The documentation’s options for robust centre estimation, eye clean-up and bias reduction show the method’s usual failure points: neck in the field of view, residual eye and optic nerve voxels, and shading.
• Hybrids with atlases. FreeSurfer’s classical stream uses a hybrid of a watershed and surface deformation (FreeSurferMethodsCitation); fMRIPrep uses a Nipype implementation of ANTs’ antsBrainExtraction.sh, which its documentation describes as atlas-based (fMRIPrep documentation).
• Learned models trained on synthetic data. SynthStrip, part of FreeSurfer since version 7.3.0, is a deep network trained on images synthesized from segmentation maps, so that it does not depend on a particular contrast or resolution; its authors’ dataset of full-head images with ground-truth brain masks spans 622 MRI, CT and PET scans, from infants to glioblastoma patients (SynthStrip page). FreeSurfer 8.0 made it part of the default stream (release notes).

The output is a binary mask whose errors propagate: brain cut away is lost to every later measurement, and dura or bone left inside can be classified as grey matter.

Bias-field correction

A smooth multiplicative shading from coil sensitivity and transmit-field inhomogeneity makes the same tissue brighter in one part of the head than another (MRI). Medical Image Preprocessing describes N3 and N4, which estimate the field from the image alone. Neuroimaging suites also build the correction into segmentation: FSL’s FAST estimates the bias field while it segments and outputs both the field and a corrected image, and SPM’s segmentation saves an “intensity nonuniformity corrected” image (FSL documentation, FAST; SPM segmentation tutorial). The two problems are coupled, since the field is easier to estimate when tissue classes are known and classes are easier to find when the field is removed. The example below shows what an uncorrected field does to tissue volumes.

Registration to a template

To compare brains, each is mapped to a common space. Templates are average brains: the MNI152 family (fMRIPrep uses MNI152NLin2009cAsym as its default spatial-standardization reference, and FAST can initialize from MNI152 tissue probability maps) and Talairach space, which FreeSurfer’s classical stream computes (fMRIPrep documentation; FSL documentation, FAST; FreeSurferMethodsCitation). The mapping is a rigid or affine alignment followed by a deformable one; medical image registration covers the transforms, similarity measures, diffeomorphic methods such as ANTs SyN, and the largest evaluation of brain registration algorithms, by Klein et al. (2009). Templates also carry labels: an atlas registered to a new brain gives a segmentation, and several atlases with label fusion give a better one, the multi-atlas approach whose history organ segmentation traces from the papers of Rohlfing and colleagues in 2003 and 2004, and which its BTCV organizers described as very successful for brain and cranial structures.

The cortex is a folded sheet, and folds differ between people. FreeSurfer therefore registers cortical surfaces through a spherical atlas based on each person’s folding pattern, rather than relying only on volume registration (FreeSurferMethodsCitation).

Tissue segmentation

Classifying brain voxels as grey matter, white matter or CSF is one of the field’s longest-studied tasks and the base of volumetry. Intensity alone almost separates the three classes on a T1-weighted image, which is why the classical methods are mixtures of distributions with spatial priors:

• FAST (FSL) uses a hidden Markov random field model fitted by expectation–maximization, starting by default from a k-means segmentation, and outputs partial volume maps (the fraction of each tissue in each voxel) as well as hard labels; three classes is the default, two when grey–white contrast is poor, four when strong lesions form a class of their own (FSL documentation, FAST).
• SPM segments with tissue probability maps of six classes (grey matter, white matter, CSF, non-brain soft tissue, skull and other), and the same model yields the bias-corrected image and a deformation field to template space (SPM segmentation tutorial).

The mixture view is the one Gaussian mixture models and expectation–maximization describe; the spatial priors and the Markov field handle noise and the overlap of class intensities at boundaries. Voxels at a tissue boundary contain both tissues, so a hard label over- or under-counts each class, and FreeSurfer’s SynthSeg page notes that its reported volumes are computed from soft segmentations and so differ from counts of labelled voxels (SynthSeg page).

Subcortical segmentation, cortical surfaces and parcellation

Beyond tissue classes, pipelines label structures. FreeSurfer’s volumetric segmentation (aseg) labels subcortical structures such as the hippocampus, amygdala, caudate, putamen and ventricles; its atlas was built from 40 subjects scanned with the same MP-RAGE sequence, in four groups of ten: young, middle-aged, healthy older adults and older adults with Alzheimer’s disease (FreeSurferMethodsCitation). FSL’s FIRST is a model-based segmentation and registration method for subcortical structures such as the putamen and caudate nucleus, which needs a good-quality T1-weighted image (FSL documentation).

The cortex is reconstructed as two surfaces, the grey–white boundary and the pial surface, placed where intensity changes most between classes. Cortical thickness is then the closest distance from the grey–white boundary to the grey–CSF boundary at each vertex of the mesh, which is not limited to the voxel grid; FreeSurfer’s citation page states that its thickness measurements have been validated against histological analysis and manual measurements. The cortex is then parcellated into regions by gyral and sulcal structure, in FreeSurfer with methods that its citation page attributes to Desikan et al. (2006) and to Fischl and colleagues’ 2004 paper on automatically parcellating the cerebral cortex (FreeSurferMethodsCitation). Parcellations turn a brain into a table of regional volumes, thicknesses and areas, the “imaging-derived phenotypes” that UK Biobank releases in thousands (Littlejohns et al., 2020).

Learned and contrast-agnostic methods

Convolutional networks now do many of these steps directly, much faster than the classical streams. FreeSurfer’s own move is the clearest case: version 8.0 (February 2025) brought SynthSeg, SynthStrip and SynthMorph into recon-all, deep learning tools the release notes describe as making the stream more robust and much faster (SynthMorph is not needed for recon-all to run, but its nonlinear registration can map to MNI152), and added learned segmentations of regions that often cause surface errors, such as the major venous sinuses (FreeSurfer release notes). SynthSeg is described by its page as segmenting brain scans of any contrast and resolution without retraining, robust to bias field, skull, intensity normalization and the absence of registration because of aggressive augmentation during training, with output at 1 mm regardless of input resolution, and with optional cortical parcellation, quality control and intracranial volume estimation (SynthSeg page). The common idea, training on images synthesized from label maps with randomized contrast, addresses the brittleness to acquisition that Medical Computer Vision documents.

Lesions and Tumors

White matter lesions. In multiple sclerosis and in small vessel disease, lesions appear bright on FLAIR. They also break tissue segmentation, since a T1 pipeline sees them as grey matter or CSF; FAST’s documentation suggests a fourth class when lesions are strong, FSL includes BIANCA, a supervised k-nearest-neighbour classifier that outputs each voxel’s probability of being a white matter hyperintensity, and a lesion-filling tool, and FreeSurfer offers WMH-SynthSeg for scans with lesions or from low-field scanners (FSL documentation; SynthSeg page). Lesion counts and loads are measured by detection and segmentation together, so lesion detection metrics, as well as overlap, matter.

The MICCAI 2016 MS lesion segmentation challenge shows how hard the task is (Commowick et al., 2018). Its 53 patients came from three French sites on four scanners (Siemens, Philips and GE; three 3 T, one 1.5 T) under a common protocol; each was delineated by seven experts, and a consensus built with LOP STAPLE served as reference. (The abstract counts four centres, one per scanner; the text and the scanner table give three hospitals.) The 15 training cases came from three scanners, so one scanner appeared only in the 38 test cases, and the organizers ran all thirteen methods themselves on their platform, so that teams could not tune parameters on test images. Every method scored below every expert: the best mean Dice was 0.59 against 0.67 for the lowest expert, and lesion detection F1 ranged from 0.13 to 0.49 against at least 0.66 for the experts. On the unseen scanner, the average detection F1 over teams fell to 0.22, against 0.32–0.39 on the other centres. Fusing all teams’ outputs with LOP STAPLE came closer to the experts on segmentation scores but less so on detection. One test case had an empty consensus, though five of the seven experts had marked between one and eight lesions in it, an illustration of the disagreement that tumor segmentation discusses for annotators in general.

Brain tumors. Glioma segmentation in multi-contrast MRI, with its nested whole tumor, tumor core and enhancing tumor regions, the BraTS annotation protocol and its ranking, is the main case study of tumor segmentation. Tumors also disturb the standard pipeline: a template registration assumes that a brain is a deformed template, which a mass that displaces tissue violates, and skull strippers and tissue classifiers trained on healthy brains can fail around them. SynthStrip’s dataset includes glioblastoma patients (SynthStrip page).

Functional and Diffusion MRI in Outline

Functional MRI. Preprocessing aligns the volumes of a run to correct head motion, corrects slice timing and the distortion of echo-planar images (field unwarping), registers the run to the person’s structural scan and maps it to a template; fMRIPrep defines its “minimal preprocessing” as motion correction, field unwarping, normalization, bias-field correction and brain extraction, built from FSL, ANTs, FreeSurfer and AFNI tools (fMRIPrep documentation). Task analysis then fits a linear model to each voxel’s time course,

$$y(t) = a\,x(t) + b + e(t),$$

where $x(t)$ is the stimulus time course convolved with a haemodynamic response function that blurs and delays it, $b$ the baseline and $e(t)$ the error; the parameter estimate divided by its standard error gives a $t$ statistic per voxel, and contrasts between parameters compare conditions (FSL documentation, FEAT user guide and GLM overview). FEAT smooths each volume spatially by default with a 5 mm kernel, which reduces noise as long as the activation is larger than the kernel. Resting-state analysis instead relates time courses between regions (functional connectivity) or decomposes them, for example by independent component analysis into spatial and temporal components, as FSL’s MELODIC does. The vision content is mostly in the preprocessing: registration across modalities, distortion correction and brain masks.

Diffusion MRI. The diffusion tensor model represents diffusion in each voxel by a symmetric 3 × 3 matrix, which needs diffusion encoding along at least six non-collinear directions, after correction for eddy-current distortions and head motion (Soares et al., 2013). With eigenvalues $\lambda_1, \lambda_2, \lambda_3$ of the tensor, the mean diffusivity and the fractional anisotropy are

$$\mathrm{MD} = \frac{\lambda_1 + \lambda_2 + \lambda_3}{3}, \qquad \mathrm{FA} = \sqrt{\frac{3}{2}}\,\frac{\sqrt{(\lambda_1 - \bar\lambda)^2 + (\lambda_2 - \bar\lambda)^2 + (\lambda_3 - \bar\lambda)^2}}{\sqrt{\lambda_1^2 + \lambda_2^2 + \lambda_3^2}},$$

with $\bar\lambda = \mathrm{MD}$; FA runs from 0 for isotropic diffusion to 1 for diffusion along a single direction, and the principal eigenvector gives the direction of fastest diffusion, which colour-coded FA maps show (Soares et al., 2013). Tractography follows these directions from voxel to voxel to trace white matter pathways, stopping at low FA or sharp turns. A single tensor cannot represent two fibre populations crossing in one voxel, which leads to wrong directions and truncated tracts; models with more directions or several compartments, such as high angular resolution methods and spherical deconvolution, address this some of them at the cost of longer acquisitions because they need more encoding directions (Soares et al., 2013). FSL’s diffusion tools cover the chain from distortion correction (TOPUP, EDDY) to tensor fitting (DTIFIT), probabilistic tractography (BedpostX, ProbtrackX) and group analysis of FA (TBSS) (FSL documentation).

Software

The main suites, as of November 2026 from their official pages:

Suite	Maintainer	Scope (from its documentation)	Latest release	Licence
FreeSurfer	Laboratory for Computational Neuroimaging, Martinos Center	structural stream (skull stripping, bias correction, segmentation, surfaces, parcellation), longitudinal stream, FS-FAST for fMRI, TRACULA for diffusion	8.2.0 (18 March 2026)	FreeSurfer licence, research use
FSL	Analysis group, Oxford Centre for Integrative Neuroimaging	structural (BET, FAST, FIRST, SIENA), registration (FLIRT, FNIRT), task and resting fMRI (FEAT, MELODIC), diffusion (TOPUP, EDDY, DTIFIT, TBSS), statistics	6.0.7.23 (28 July 2026)	free for non-commercial use, with some components open source
SPM	Functional Imaging Laboratory, UCL	spatial normalization, segmentation, statistical parametric maps; fMRI, PET, SPECT, EEG and MEG	26.10.01 on GitHub (13 October 2026)	GPL-2.0 (repository)
ANTs	ANTsX	registration, template construction, N4, Atropos segmentation, cortical thickness, brain extraction, multi-atlas label fusion	v2.6.5 (January 2026)	Apache-2.0
fMRIPrep	NiPreps	fMRI preprocessing assembled from FSL, ANTs, FreeSurfer and AFNI, on BIDS input	25.2.5 (10 March 2026)	Apache-2.0

SPM’s software page still named 25.01 as the latest version when it was read, while its GitHub repository had released 26.10.01; the table gives the newer release. FreeSurfer’s wiki calls it the structural MRI analysis software of choice for the Human Connectome Project, and UK Biobank’s processing pipeline is built mainly around FSL with FreeSurfer and other packages (Littlejohns et al., 2020). SPM began as Karl Friston’s software for PET analysis, made available in 1991, and SPM’94 was its first major revision (SPM website). Two release notes bear on reproducibility: FreeSurfer 8.2.0 notes that different combinations of build and execution operating systems can give different results, and FreeSurfer 7.4.0 that one topology-correction step could differ slightly with the number of threads (FreeSurfer release notes).

Datasets and Access

Large neuroimaging datasets are shared under data use agreements; none is a download without terms. The summaries below come from the projects’ pages and papers.

• ADNI (Alzheimer’s Disease Neuroimaging Initiative) is a longitudinal, multi-centre, observational study whose goal is to validate biomarkers for Alzheimer’s disease clinical trials, with clinical, imaging, genetic and biomarker data (structural, functional and molecular brain imaging among them). Data are shared through the LONI Image and Data Archive; access requires accepting the ADNI Data Use Agreement, an online application with institutional affiliation and proposed use, review by the Data Sharing and Publications Committee, an annual update, and adherence to its publication policy (ADNI data pages).
• UK Biobank set out in 2014 to re-invite 100,000 of its 500,000 participants, aged 40–69 at recruitment, for brain, cardiac and abdominal MRI, DXA and carotid ultrasound, and to re-image at least 10,000 of them. The resource is open to bona fide researchers at academic and commercial institutions who register and apply for a specified, health-related project, and researchers can request the derived phenotypes or the scans themselves; results and derived variables are expected to be returned to the resource (Littlejohns et al., 2020). UK Biobank’s own website refused automated access when this article was prepared, so its current figures are not given here.
• OASIS (Open Access Series of Imaging Studies) aims to make neuroimaging data freely available to the scientific community. OASIS-1 is a cross-sectional set of 416 subjects aged 18 to 96, each with three or four T1-weighted scans from one session, 100 of those over 60 with very mild to moderate Alzheimer’s disease; OASIS-2 is longitudinal (150 subjects aged 60 to 96, 373 sessions); OASIS-3 is a retrospective compilation of 1,378 participants and 2,842 MR sessions with several contrasts and many FreeSurfer segmentations. Its data use terms forbid attempts to identify participants, any facial reconstruction or recognition from the images, and commercial use without written authorization (OASIS website).
• HCP Young Adult (Human Connectome Project) has 1,200 subjects aged 22–35 (the project ran 2010–2016; its final 1200 Subjects release has data from 1,206 participants collected in 2012–2015), with lifespan and disease studies following the same style of protocol. Open-access data (all imaging and most behavioural data) require an account and acceptance of the open-access data use terms; restricted data such as family structure and age by year require a separate application. The HCP notes that its released data are not considered de-identified, because combinations of restricted data might identify individuals (HCP pages).
• Challenge datasets pair images with expert references for one task: MRBrainS for tissue classes, the MICCAI 2016 MS lesion data, and BraTS for gliomas (tumor segmentation); their designs appear below.

Many methods papers use OASIS, ADNI or both; Medical Image Classification describes how the subject-level structure of these datasets leads to leakage when slices of one brain fall on both sides of a split.

Example: Skull Stripping, a Bias Field and Tissue Volumes

The script builds a synthetic T1-like head of 80 × 96 × 80 voxels at 2 mm, with scalp, a dark skull, a CSF layer, two CSF-filled ventricles, a grey matter ribbon whose inner boundary is folded by a sinusoidal pattern, and white matter, each with a known label. It renders the head with a little blur and Rician noise (the magnitude of complex Gaussian noise, as MRI explains), once without and once with a smooth multiplicative bias field that ranges from 0.79 to 1.25 inside the skull, a range chosen for illustration rather than taken from a scanner. A simple skull stripper (a threshold between skull and CSF, one erosion, the largest connected component, one dilation and hole filling, the mathematical morphology toolkit) gives the intracranial mask, and three-class multi-level Otsu thresholding inside the mask gives CSF, grey and white matter. The biased image is also corrected with SimpleITK’s N4 filter, estimated on the mask at half resolution. Volumes are compared with the true label volumes, overall and in each half of the head along $x$, the axis of the field’s largest linear term. The script loops over noise seeds 0–4 and prints the ranges. It runs in a few seconds with SimpleITK 2.5.7, NumPy 2.5.4, SciPy 1.18.1 and scikit-image 0.26.0.

All parameters were fixed before the first run. After it, which printed the whole-brain table and seed ranges, the per-half grey matter measurement was added, because the whole-brain volume errors had stayed within 2.5% in every condition while the per-class Dice fell; the labels of the halves were then renamed from “left” and “right” to “x<0” and “x>0”, since the phantom has no anatomical orientation. Nothing else changed.

import numpy as np
import SimpleITK as sitk
from scipy import ndimage as ndi
from skimage.filters import threshold_multiotsu

SHAPE, SP = (80, 96, 80), 2.0                 # voxels, isotropic spacing in mm
VOX_ML = SP**3 / 1000.0                        # millilitres per voxel
NAMES = {3: "CSF", 4: "GM", 5: "WM"}

def phantom():
    """Label map: 0 air, 1 scalp, 2 skull, 3 CSF, 4 grey matter, 5 white matter."""
    z, y, x = [(np.arange(n) - n / 2 + 0.5) * SP for n in SHAPE]
    Z, Y, X = np.meshgrid(z, y, x, indexing="ij")
    rho = np.sqrt((X / 60) ** 2 + (Y / 76) ** 2 + (Z / 62) ** 2)
    r = np.sqrt(X**2 + Y**2 + Z**2) + 1e-9
    theta, phi = np.arccos(Z / r), np.arctan2(Y, X)
    folds = np.sin(6 * phi) * np.sin(5 * theta)          # crude "gyri" on the GM/WM boundary
    lab = np.zeros(SHAPE, np.uint8)
    lab[rho < 1.19] = 1
    lab[rho < 1.11] = 2
    lab[rho < 1.03] = 3
    lab[rho < 0.97] = 4
    lab[rho < 0.80 + 0.06 * folds] = 5
    for sx in (-1, 1):                                    # two ventricles
        v = ((X - sx * 9) / 6) ** 2 + ((Y + 4) / 22) ** 2 + ((Z - 6) / 9) ** 2
        lab[v < 1] = 3
    return lab, (X, Y, Z)

def bias_field(X, Y, Z):
    u, v, w = X / 80, Y / 96, Z / 80
    return np.exp(0.18 * u + 0.12 * v - 0.15 * w + 0.12 * (u**2 - v**2) + 0.10 * u * w)

def render(lab, bias, seed):
    mean = np.array([0, 90, 5, 35, 75, 110], float)       # T1-like: WM > GM > CSF, dark skull
    img = ndi.gaussian_filter(mean[lab], 0.6) * bias
    rng = np.random.default_rng(seed)
    noise = rng.normal(0, 3, (2,) + SHAPE)
    return np.abs(img + noise[0] + 1j * noise[1])         # magnitude image (Rician noise)

def skull_strip(img):
    m = ndi.binary_erosion(img > 20)                      # 6-connected, 1 voxel
    cc, n = ndi.label(m)
    sizes = ndi.sum(m, cc, range(1, n + 1))
    m = cc == (1 + np.argmax(sizes))                      # keep the largest component
    return ndi.binary_fill_holes(ndi.binary_dilation(m))

def n4(img, mask):
    im = sitk.GetImageFromArray(img.astype(np.float32)); im.SetSpacing((SP,) * 3)
    mk = sitk.GetImageFromArray(mask.astype(np.uint8)); mk.CopyInformation(im)
    f = sitk.N4BiasFieldCorrectionImageFilter()
    f.SetMaximumNumberOfIterations([50, 50, 50, 50])
    f.Execute(sitk.Shrink(im, [2] * 3), sitk.Shrink(mk, [2] * 3))
    log_b = sitk.GetArrayFromImage(f.GetLogBiasFieldAsImage(im))
    return img / np.exp(log_b), log_b

def dice(a, b):
    return 2 * np.sum(a & b) / (np.sum(a) + np.sum(b))

def segment(img, mask):
    t1, t2 = threshold_multiotsu(img[mask], classes=3)
    seg = np.zeros(SHAPE, np.uint8)
    seg[mask & (img < t1)] = 3
    seg[mask & (img >= t1) & (img < t2)] = 4
    seg[mask & (img >= t2)] = 5
    return seg

lab, (X, Y, Z) = phantom()
B = bias_field(X, Y, Z)
icv = lab >= 3
true_vol = {k: np.sum(lab == k) * VOX_ML for k in NAMES}
print("bias field over the intracranial volume: %.2f to %.2f" % (B[icv].min(), B[icv].max()))
print("true volumes (mL): " + "  ".join(f"{NAMES[k]} {true_vol[k]:.0f}" for k in NAMES))

errs = {c: [] for c in ("no bias", "bias", "bias + N4")}
gm_half = {c: {} for c in errs}
for seed in range(5):
    imgs = {"no bias": render(lab, 1.0, seed), "bias": render(lab, B, seed)}
    m_b = skull_strip(imgs["bias"])
    imgs["bias + N4"], log_b = n4(imgs["bias"], m_b)
    if seed == 0:
        r = np.corrcoef(log_b[m_b], np.log(B[m_b]))[0, 1]
        print(f"N4 log-field vs true log-field inside the mask: r = {r:.3f}")
        print(f"{'condition':<10} {'mask Dice':>9} " + " ".join(f"{NAMES[k] + ' mL (err)':>16}" for k in NAMES)
              + " " + " ".join(f"{'Dice ' + NAMES[k]:>8}" for k in NAMES))
    for c, img in imgs.items():
        mask = m_b if c == "bias + N4" else skull_strip(img)
        seg = segment(img, mask)
        vol = {k: np.sum(seg == k) * VOX_ML for k in NAMES}
        e = {k: 100 * (vol[k] / true_vol[k] - 1) for k in NAMES}
        errs[c].append([e[k] for k in NAMES])
        for side, half in (("x<0", X < 0), ("x>0", X > 0)):
            g = 100 * (np.sum((seg == 4) & half) / np.sum((lab == 4) & half) - 1)
            gm_half[c].setdefault(side, []).append(g)
        if seed == 0:
            print(f"{c:<10} {dice(mask, icv):9.3f} "
                  + " ".join(f"{vol[k]:8.0f} ({e[k]:+5.1f}%)" for k in NAMES) + " "
                  + " ".join(f"{dice(seg == k, lab == k):8.3f}" for k in NAMES))
print("volume error range over noise seeds 0-4 (%):")
for c, v in errs.items():
    v = np.array(v)
    print(f"  {c:<10} " + "  ".join(f"{NAMES[k]} {v[:, i].min():+5.1f} to {v[:, i].max():+5.1f}"
                                     for i, k in enumerate(NAMES)))
print("GM volume error in each half of the head, range over seeds 0-4 (%):")
for c, d in gm_half.items():
    print(f"  {c:<10} " + "  ".join(f"{side} {min(v):+5.1f} to {max(v):+5.1f}" for side, v in d.items()))

Output:

bias field over the intracranial volume: 0.79 to 1.25
true volumes (mL): CSF 223  GM 472  WM 600
N4 log-field vs true log-field inside the mask: r = 0.631
condition  mask Dice     CSF mL (err)      GM mL (err)      WM mL (err) Dice CSF  Dice GM  Dice WM
no bias        0.999      221 ( -1.0%)      474 ( +0.5%)      599 ( -0.1%)    0.987    0.994    0.998
bias           0.997      225 ( +0.7%)      465 ( -1.4%)      602 ( +0.4%)    0.946    0.943    0.970
bias + N4      0.997      218 ( -2.4%)      464 ( -1.6%)      609 ( +1.6%)    0.977    0.974    0.982
volume error range over noise seeds 0-4 (%):
  no bias    CSF  -1.2 to  -1.0  GM  +0.5 to  +0.5  WM  -0.1 to  -0.0
  bias       CSF  +0.7 to  +1.3  GM  -1.8 to  -0.8  WM  -0.3 to  +0.5
  bias + N4  CSF  -2.4 to  -2.3  GM  -1.6 to  -0.9  WM  +1.1 to  +1.6
GM volume error in each half of the head, range over seeds 0-4 (%):
  no bias    x<0  +0.5 to  +0.5  x>0  +0.4 to  +0.6
  bias       x<0  +1.5 to  +2.7  x>0  -5.0 to  -4.1
  bias + N4  x<0  -1.1 to  -0.5  x>0  -2.3 to  -1.4

Skull stripping is barely affected: the dark skull leaves a wide intensity gap between CSF and bone, which a field of this size does not close, so the mask Dice stays near 1 in every condition. Tissue segmentation is affected, and in a way the whole-brain volumes hide. With the field, a single pair of thresholds is too low for the bright side of the head and too high for the dark side, so white matter on the dark side is labelled grey matter and grey matter on the bright side is labelled white matter. The Dice of every class falls, yet the whole-brain grey matter volume changes by only a percentage point or two, because the two errors partly cancel. Split the head in two and they separate: grey matter is overestimated on the dark half and underestimated by 4–5% on the bright half. A group comparison of regional volumes, or a longitudinal one in which the head sits differently in the coil, would inherit that error.

N4 recovers most of the regional error and the Dice, but not all of it, and it does not improve the whole-brain volumes, which now carry a small shift of their own (CSF lower, white matter higher). Its estimated field correlates only moderately with the true one; on our reading, the B-spline fit and histogram sharpening absorb part of the tissue contrast near the folded boundary into the field. The lesson for pipelines is the one the MRBrainS organizers drew in their design: bias correction is a step with its own errors, so they corrected all challenge data with one method (SPM8) before release, to keep differences in bias correction out of the comparison of segmentation methods (Mendrik et al., 2015). The phantom is far simpler than a brain (one contrast, no partial-volume model in the classifier, a smooth field of chosen shape), so the numbers show the mechanism, not the size of the effect in real scans.

Evaluation and Results with Their Designs

Neuroimaging methods are judged in two ways: against manual references on a few scans, and by the reliability and validity of the measurements they produce on many.

Overlap with a manual reference. MRBrainS (Mendrik et al., 2015) segmented grey matter, white matter and CSF on 3 T scans of 20 older adults (aged 65–80, mean 71) from a cohort without invalidating stroke or other brain disease, chosen to have varying degrees of atrophy and white matter lesions, acquired on one Philips scanner with T1, T1 inversion recovery and FLAIR sequences. Trained research assistants segmented every axial slice manually, and three experts checked and corrected the result; there was one reference per scan, so inter-observer variability could not be measured. Five scans were released for training and fifteen kept for testing. Methods were evaluated by Dice, the 95th-percentile Hausdorff distance and absolute volume difference over grey matter, white matter, CSF, brain and intracranial volume, and ranked on the grey matter, white matter and CSF scores. Eleven workshop methods and three freeware packages, run with default settings except where the authors set parameters, were compared; the top two methods tied on score and were separated by the standard deviation of their ranks, and the best method reached a mean grey matter Dice of 84.7% and white matter Dice of 88.4%. Among the freeware packages, FreeSurfer (v5.3.0) was the most robust to white matter lesions, which the authors attribute to its atlas prior, while its grey matter Dice (77.4%) was lower than the other packages’, possibly because it was evaluated on the thick-slice T1 rather than the high-resolution scan it had processed. SPM12 with T1 and FLAIR performed well among the packages on the thick-slice scans. The authors caution that the results hold for one scanner, one protocol and one group of elderly subjects.

Lesions. The MS lesion challenge above (Commowick et al., 2018) is the counterpart for lesions, with seven raters per case, a held-out scanner and organizer-run methods.

Registration. Klein et al. (2009) compared 14 nonlinear registration algorithms on 80 manually labelled brains by the overlap of propagated labels; medical image registration gives the design and results.

Reproducibility of findings. Measurement accuracy is not the same as a reproducible scientific result, and two large studies measured the gap. Marek et al. (2022) used three of the largest datasets (ABCD, n = 11,874; HCP, n = 1,200; UK Biobank, n = 35,735), observational data analysed by resampling, to estimate brain-wide association effects, correlations between cortical thickness or resting-state connectivity and behavioural phenotypes (41 measures of demographics, cognition and mental health in ABCD). In the rigorously denoised ABCD sample (n = 3,928) the median univariate effect was $|r| = 0.01$, the top 1% of associations exceeded $|r| = 0.06$, and the largest that replicated out of sample was $|r| = 0.16$. At n = 25, about the median size of neuroimaging studies, the 99% confidence interval of a univariate correlation was $r \pm 0.52$, so small studies produce inflated effects and replication failures; functional measures, cognitive tests and multivariate methods gave more robust effects than structural measures, mental health questionnaires and univariate methods. In a study by Botvinik-Nezer et al. (2020), 70 independent teams analysed the same fMRI dataset to test the same nine hypotheses; no two teams chose identical workflows, and the results of their hypothesis tests varied considerably, even among teams whose statistical maps were highly correlated at intermediate stages (abstract).

Diagnosis from images. Classifiers of Alzheimer’s disease from MRI are a frequent deep learning application, and a clear case of inflated results: Wen et al. (2020) and Yagis et al. (2021), whose designs Medical Image Classification gives, found near-perfect accuracy from slice-level splits and much lower accuracy when subjects were kept apart. Site effects are as strong: in the study by Glocker et al. (2019) that Medical Computer Vision describes, a classifier identified which of two studies a brain MRI came from with about 97% accuracy after skull stripping, bias correction, registration and intensity whitening.

Failure Modes and Risks

• Silent pipeline failures. A skull strip that leaves dura, a surface that crosses into the skull, a registration that converged to the wrong solution: each produces plausible numbers. Quality control by visual review or automated checks is part of the method, and FreeSurfer documents manual editing for both its basic and longitudinal streams (FreeSurferWiki). SynthSeg includes an automated quality control output (SynthSeg page).
• Errors that depend on position, scanner or group. As the example shows, a spatially varying error can leave totals unchanged and bias regions; a lesion or a tumor breaks assumptions of healthy anatomy; and acquisition differences between sites survive preprocessing (Glocker et al., 2019, in Medical Computer Vision).
• Version and platform dependence. Results change between major versions (FreeSurfer 8 gives different recon-all output from version 7) and can change with the operating system (FreeSurfer release notes). FreeSurfer’s release notes advise analysing all subjects with a single version, and recording the version and platform makes a study repeatable.
• Analytic flexibility and small samples. Many defensible pipelines give different conclusions (Botvinik-Nezer et al., 2020), and small samples give inflated effects (Marek et al., 2022).
• Leakage. Slices or visits of one subject in both training and test sets inflate classifier accuracy (Medical Image Classification).
• Re-identification. Structural scans contain the face; defacing, data use terms that forbid facial reconstruction (OASIS) and the HCP’s statement that its data are not considered de-identified show that privacy is a property of the whole dataset, not only of the images.

Deployment Considerations

• Organize data as BIDS and process with versioned containers. BIDS input lets pipelines such as fMRIPrep run with minimal manual parameters, and containers fix the software version; FreeSurfer 8.2.0 adds a file that restricts the operating systems on which recon-all runs, which, with a file of allowed versions, the notes say goes a long way towards assuring repeatability (fMRIPrep documentation; FreeSurfer release notes).
• Use the longitudinal streams for change over time. FreeSurfer’s longitudinal stream builds an unbiased within-subject template and initializes several processing steps of each time point from it, which its citation page says significantly increases reliability and statistical power (FreeSurferMethodsCitation).
• Report measurements with their method. A volume from soft segmentation, from counted voxels, or from a surface is a different number, and intracranial volume estimates depend on the method too (SynthSeg page; FreeSurfer release notes).
• Test on the scanners, populations and pathologies of the intended use, and, for anything aimed at patient care, with the evidence that clinical validation describes; FreeSurfer’s licence states that clinical applications are neither recommended nor advised (FreeSurfer licence).
• Respect the access terms. ADNI, UK Biobank, OASIS and HCP each require registration and acceptance of terms, and some restrict commercial use, redistribution or the reporting of restricted variables.