Medical Image Preprocessing

A CT or MRI scan arrives as a stack of stored integers with a header, not as a tensor ready for a network. Medical image preprocessing is everything between the two: reading the header to learn where each voxel lies in the patient and what its value means, bringing every scan to a common orientation and voxel spacing, resampling the annotations with them, putting intensities on a comparable scale, and cropping away what carries no information. None of these steps is learned, yet each changes what a model sees, and a mistake in one, such as a flipped axis, labels resampled with the wrong interpolator or statistics computed over air, is silent: the code runs and the numbers come out wrong.

The steps depend on the modality. Hounsfield units give CT a calibrated scale, so a fixed clipping range means roughly the same tissue in every scan; MRI intensities have no absolute scale and must be normalized per image. Both modalities share the geometry: voxels have a physical size, often larger between slices than within them, and a position and orientation in a patient coordinate system that the file format records.

Purpose

Convolutional networks operate on voxel grids and ignore the physical spacing (Isensee et al., 2020). Treating an input as an array, they assume that a step of one voxel means the same distance everywhere, that the first array axis points the same way in every training case, and that equal values mean similar things. Stored scans satisfy none of these by default. Pérez-García et al. (2021) illustrate the problem with an MRI and a brain parcellation of the same size, 181 × 181, that are correctly superimposed only by software that reads the spatial metadata; ignoring the origin, the orientation or the spacing misaligns them, and the authors expect training on such pairs without preprocessing to give poor results. Preprocessing makes the array view consistent with the physical one, so that the model can learn anatomy rather than acquisition settings.

When to Use It

Every learned method on 3D medical images needs some of these steps, and the questions to ask are which and how. Spatial standardization (orientation and spacing) matters whenever cases come from different protocols or scanners, and whenever predictions must be returned on the original grid for measurement. Intensity normalization matters most for MRI, whose values depend on the sequence, the coil and the scanner. Bias-field correction matters, on our reading, for MRI analyses that rely on a tissue having the same intensity across the image. For 2D projection images and whole-slide images the geometry is simpler, and the corresponding problems (detector processing, stain variation) belong to X-ray radiography and whole-slide imaging.

Preprocessing is not a cure for domain shift. Per-image normalization removes a global gain and offset but, by our reasoning, not differences in noise, blur or shading; in a study that Medical Computer Vision describes, a classifier still identified which of two studies each brain MRI came from after intensity whitening.

How It Works

Reading DICOM

Clinical scanners store and exchange images in DICOM. The Image Plane Module describes an image as a two-dimensional plane placed in three-dimensional space (NEMA, 2026, C.7.6.2), so a volume is assembled from the geometry of its slices:

• Image Position (Patient) gives the x, y and z coordinates, in millimetres, of the centre of the first transmitted pixel, the upper left-hand corner of the image.
• Image Orientation (Patient) gives the direction cosines of the first row and then of the first column, each a unit vector, the two orthogonal.
• Pixel Spacing gives the physical distance between pixel centres as a pair: the spacing between adjacent rows first, then the spacing between adjacent columns.
• Slice Thickness is the nominal thickness of the slice, and Spacing Between Slices, when present, the centre-to-centre distance between adjacent slices. A gap between slices shows up as a spacing larger than the thickness.

The standard maps the pixel in column $i$ and row $j$, both counted from zero, to patient coordinates by (NEMA, 2026, Equation C.7.6.2.1-1)

$$\begin{bmatrix} P_x \\ P_y \\ P_z \\ 1 \end{bmatrix} = \begin{bmatrix} X_x \Delta_i & Y_x \Delta_j & 0 & S_x \\ X_y \Delta_i & Y_y \Delta_j & 0 & S_y \\ X_z \Delta_i & Y_z \Delta_j & 0 & S_z \\ 0 & 0 & 0 & 1 \end{bmatrix} \begin{bmatrix} i \\ j \\ 0 \\ 1 \end{bmatrix},$$

where $\mathbf{S}$ is the Image Position, $\mathbf{X}$ and $\mathbf{Y}$ are the row and column direction cosines, $\Delta_i$ is the column spacing (the second value of Pixel Spacing) and $\Delta_j$ the row spacing (the first). The Patient-Based Coordinate System is right-handed: for a patient described as a biped, x increases toward the patient’s left, y toward the posterior and z toward the head. This is the convention usually called LPS. The standard also defines Slice Location, but only relative to an unspecified, implementation-specific reference point, so it is a poor key for ordering slices. It follows from the equation that the position of a slice along the volume’s third axis is the projection of its Image Position onto the normal $\mathbf{X} \times \mathbf{Y}$, which gives a reliable order and the actual spacing between slices (our reasoning).

Stored values are not yet physical values. The Modality LUT Module maps stored values (SV) to output units either by a lookup table or by a linear rescale, output $= m \cdot \text{SV} + b$, with Rescale Slope $m$ and Rescale Intercept $b$, but not both (NEMA, 2026, C.11.1). For CT, the CT Image Module requires the output units to be Hounsfield units when the image is an original image, not a localizer, and not a multi-energy acquisition (NEMA, 2026, C.8.2.1). The output of the rescale may be signed even when the stored pixels are unsigned (C.11.1.1.1), so code that skips the rescale can read air as a positive number rather than about −1000 HU.

CT images are often not rectangular in their native form. A device may pad the area outside the circular reconstruction perimeter with a Pixel Padding Value that, by the standard’s rule, occurs nowhere in the native image, so that display software can leave it out of the dynamic range; the value is given in stored units, before the rescale (NEMA, 2026, C.7.5.1.1.2). It also belongs out of intensity statistics and foreground masks.

For MR images, the MR Image Module fixes one sample per pixel, 16 bits allocated and a MONOCHROME1 or MONOCHROME2 photometric interpretation, and records the acquisition (NEMA, 2026, C.8.3.1): Scanning Sequence (SE spin echo, IR inversion recovery, GR gradient recalled, EP echo planar, RM research mode), Repetition Time and Echo Time in milliseconds, Inversion Time (required for inversion recovery), Echo Train Length, Imaging Frequency in MHz, Imaged Nucleus and the nominal Magnetic Field Strength in tesla. These attributes let a pipeline group series by contrast and field strength instead of trusting free-text series descriptions.

Reading NIfTI

NIfTI-1 was designed for exchanging data between functional MRI analysis packages, whose authors (AFNI, BrainVoyager, FSL and SPM) committed to support it, and research pipelines often convert DICOM series to it. The format was adapted from ANALYZE 7.5 and stores a volume either as one .nii file or as a .hdr and .img pair (NIfTI, 2026). Its header, 348 bytes long, holds the array dimensions (dim), the grid spacings (pixdim) with their units, a data scaling, and two ways to relate voxel indices $(i, j, k)$ to spatial coordinates $(x, y, z)$ of voxel centres (Cox, nifti1.h):

• Method 2, the qform (used when qform_code > 0): a rotation stored as a quaternion, scaled by the pixdim spacings, plus the offsets qoffset_x, qoffset_y, qoffset_z. Because a quaternion encodes only proper rotations, a factor $q_\text{fac} = \pm 1$, stored in pixdim[0], multiplies the third axis so that left-handed voxel grids can be represented. The header comments intend this method for the coordinates the scanner reported, or for those rotated to a fiducial orientation.
• Method 3, the sform (used when sform_code > 0): a general affine matrix stored row by row in srow_x, srow_y and srow_z, intended for coordinates in some standard space, such as one reached by registration.
• Method 1, plain scaling by pixdim with no orientation, exists only for compatibility with ANALYZE files and applies when qform_code is 0.

The codes name the space: 1 for scanner-based anatomical coordinates, 2 for coordinates aligned to another file or to anatomical truth, 3 for Talairach, 4 for MNI 152 and, since 2019, 5 for another template space. If scl_slope is nonzero, stored values are scaled to “true” values by y = scl_slope * x + scl_inter, the counterpart of DICOM’s rescale, and cal_min and cal_max give a display range, the counterpart of a window.

NIfTI’s axes differ from DICOM’s. In Methods 2 and 3, +x points to the subject’s right, +y anterior and +z superior, a right-handed system. The header’s own comments describe the conversion: DICOM’s coordinate system is rotated by 180° about the z-axis relative to NIfTI’s, so converting between them means negating x and y. A DICOM Image Position $(p_x, p_y, p_z)$ becomes the offsets $(-p_x, -p_y, p_z)$, and the first two columns of the rotation are the DICOM row and column direction cosines with their x and y components negated. DICOM does not encode the direction of the third axis simply, the comments add: its sign can be inferred from the positions of successive slices, a method that, in Cox’s words, “occasionally fails”. Comparing the affine after conversion with the DICOM geometry is therefore a useful check.

Coordinate systems and orientation

A preprocessing pipeline juggles voxel indices, patient coordinates in millimetres, and the axis order of arrays in memory. SimpleITK indexes images in x, y, z order (SimpleITK documentation), whereas the NumPy arrays it exchanges are ordered z, y, x, as the example’s code shows. SimpleITK and NiBabel, which Pérez-García et al. (2021) name as the frameworks medical images are commonly read with, follow the two formats. SimpleITK describes an image by its origin (the physical position of the voxel with all-zero index), spacing, size and a direction cosine matrix whose columns give the directions of the axes, with units assumed to be millimetres for historical reasons tied to DICOM; an image with the identity direction is in LPS orientation as DICOM defines it (SimpleITK documentation). NiBabel returns the NIfTI affine, which maps voxel indices to RAS+ coordinates (NiBabel documentation).

Orientation is summarized by three letters, one per array axis, and here the naming is a trap. SimpleITK’s DICOMOrient labels each axis by the direction in which its values increase, so its “LPS” is DICOM’s convention. NiBabel calls the NIfTI system “RAS+” precisely because some authors use “RAS” for the opposite reading, axes starting rather than ending on the right, anterior and superior side, which would make it “LPI” in NiBabel’s terms (NiBabel documentation). A three-letter code is meaningful only with its library’s definition.

Reorienting to a standard code, for example with DICOMOrient(image, "LPS") in SimpleITK or as_closest_canonical in NiBabel, permutes and flips array axes and updates the header so that every voxel keeps its physical position, as the SimpleITK documentation states for DICOMOrient and the example checks. It changes the array a network sees, not the anatomy. A transform that is estimated in voxel indices of one scan, a bounding box, or a prediction saved without the source header loses this connection; medical image registration works in physical coordinates for the same reason.

Resampling and interpolation

Resampling evaluates an image on a new grid of points in physical space. In SimpleITK’s terms it needs the image, a resampling grid (size, origin, spacing and direction, or a reference image), a transform from the grid’s coordinates to the image’s, and an interpolator; points that fall outside the image are set to a default value (SimpleITK documentation). For preprocessing, the transform is usually the identity and the grid has the target spacing. Image warping describes the interpolators themselves.

Images. Linear interpolation is SimpleITK’s default and, in its documentation’s words, “reasonable for ordinary medical images”. nnU-Net uses third-order splines for images (Isensee et al., 2020); SimpleITK’s sitkBSpline is also third order. Thick slices are the hard case: the acquisition averaged the signal over the slice profile, and no interpolator can restore detail that was never sampled; resampling a 4 mm slice stack to 0.8 mm makes a finer array, not a finer image. Going the other way, to a coarser spacing, discards samples, and by the sampling theorem detail finer than the new spacing must be smoothed away first or it aliases (sampling and aliasing); Resample itself only interpolates and does not smooth (our reading of its documentation).

Labels. Interpolating a label map as if it were an image creates values between labels: halfway between background 0 and class 3, linear interpolation produces 1.5, which rounds to a class that is not there. Nearest-neighbour interpolation is, as the SimpleITK documentation puts it, the only one that introduces no new labels. The other standard approach, used by nnU-Net, converts the label map to one binary channel per class, interpolates each channel linearly and takes the class with the largest value at each voxel (Isensee et al., 2020). Images and labels must be resampled to exactly the same grid, or the network learns from misaligned pairs.

Choosing the target spacing. A coarse target spacing loses detail; a fine one makes volumes large and leaves a fixed-size patch less context. For strongly anisotropic data, interpolating many new slices between thick ones creates artifacts. nnU-Net’s rules, summarized below, are one worked answer.

Keeping the extent. A new grid of $n' = n s / s'$ samples should cover the same physical extent as the old one. Since DICOM, NIfTI and SimpleITK all place the origin at the centre of the first voxel, keeping the origin while changing the spacing shifts the grid; the extent is kept by moving the first centre by $(s' - s)/2$ along each axis (our reasoning from the definitions).

Intensity normalization by modality

CT. After the DICOM rescale, values are in Hounsfield units, and a fixed mapping applies to every scan. Clipping to a range and rescaling it is windowing used as preprocessing; Hounsfield Units defines the DICOM window formula and explains why tissue values still shift between scanners and tube voltages. In a CT segmentation pipeline that Pérez-García et al. (2021) describe as a use case, for example, intensities are clipped to [−1000, 1000] and divided by a constant. nnU-Net instead learns the range from the data: it clips to the 0.5th and 99.5th percentiles of foreground intensities in the training set and standardizes with the foreground mean and standard deviation of the whole training set, the same for every image, because CT values are quantitative (Isensee et al., 2020).

MRI. Without a fixed scale, each image is normalized by its own statistics. The simplest scheme is the z-score, $x' = (x - \mu)/\sigma$, with $\mu$ and $\sigma$ computed over the image or over a mask. Reinhold et al. (2019) describe and compare seven schemes for brain MRI, all computed inside a brain mask: the z-score; three that divide by a white-matter intensity found by fuzzy c-means, a Gaussian mixture or a kernel density estimate; piecewise linear histogram matching to a standard histogram learned from landmarks at intensity percentiles, after Nyúl and Udupa; WhiteStripe, a z-score using the white-matter peak of the smoothed histogram and the standard deviation of a narrow band of intensities around it, taken as normal-appearing white matter; and RAVEL, which further removes variation estimated from cerebrospinal fluid intensities across registered images. Several of these depend on brain anatomy; a z-score depends only on the mask.

Bias field. MR images also show smooth intensity nonuniformity, a multiplicative shading across the field of view. Sled et al. (1998) introduced the nonparametric nonuniform intensity normalization method, N3, which, as Tustison and Gee (2009) note, does not require a prior tissue model. Tustison and Gee’s ITK implementation, N4ITK, works on log intensities, where the field becomes additive, alternates between sharpening the intensity histogram and fitting a smooth B-spline field to the residual, and finally divides the input image by the estimated field (Tustison & Gee, 2009). The N4 filter in ITK and SimpleITK, published by Tustison et al. (2010), is described in its documentation as a variation of N3 with an improved B-spline fitting that allows several resolutions and with an update in which the residual bias field is continually updated. The SimpleITK documentation of the filter notes that, because the work is done in log space, negative values and values below 1 can produce poor results, that the authors recommend running the correction on a downsampled image, and that a mask restricts the voxels used to estimate the field.

Cropping and foreground masks

A scan often contains air or padding around the body. Cropping to the foreground reduces memory and computation, and computing intensity statistics inside a foreground mask keeps them from depending on how much air the field of view contains. nnU-Net crops every case to its nonzero region; the authors report that this had no effect on most of their datasets but substantially reduced the size of brain datasets such as Brain Tumor and MSLes. When cropping reduces the average size by a quarter or more, nnU-Net computes the per-image z-score inside the nonzero mask only and ignores the zero voxels around it (Isensee et al., 2020). For CT, DICOM’s padding value marks pixels outside the reconstruction perimeter; for MRI, a threshold or a brain-extraction tool provides the mask. Labels must be cropped with the same bounding box, and the box recorded, so that predictions can be put back in place.

How to Apply It

Example: orientation, resampling and normalization of a synthetic scan

The script below builds a synthetic thick-slice scan with known ground truth and walks through three steps. The phantom is defined in physical LPS millimetres: an ellipsoidal body (intensity 300), an organ (600) containing a small lesion (900), and a second, separate organ (450), with label values 1, 2 and 3 for the three structures. It is “acquired” on a grid of 100 × 100 × 20 voxels with spacing 0.8 × 0.8 × 4.0 mm, each slice averaging eight samples of the phantom across its 4 mm thickness, with Gaussian noise of standard deviation 15 (seed 0) and a direction matrix whose first axis runs toward the patient’s right. Part 1 writes the scan as NIfTI with SimpleITK and reads the header back with NiBabel. Part 2 resamples to isotropic 0.8 mm spacing with three image interpolators and three ways of treating labels, and scores them against the phantom sampled directly on the new grid: root-mean-square error inside the body for images, and the Dice coefficient per structure for labels, together with a count of voxels labelled organ or lesion within 10 mm of the second organ’s centre, where neither exists. Part 3 acquires a second scan of the same phantom with intensities 1.6 times higher and a field of view of 120 mm instead of 80 mm (seed 1), and compares the mean organ and lesion values under three normalizations, the foreground mask coming from Otsu’s threshold.

All parameters were set before the first run. After it, two printed lines were changed, and no parameter or measured value: the affine is printed with negative zeros shown as zeros, and the lesion check prints continuous voxel indices instead of rounded ones, because the lesion centre lies halfway between two voxel centres along the first axis and the rounded indices (55 and 45) obscured the relation between them. Every other printed value is the same as in the first run. The script runs in about 2 seconds with SimpleITK 2.5.7, NiBabel 5.4.2, NumPy 2.5.4 and scikit-image 0.26.0, and gives identical output on repeated runs.

import os
import numpy as np
import SimpleITK as sitk
import nibabel as nib
from skimage.filters import threshold_otsu

print("SimpleITK", sitk.Version.VersionString(), "| nibabel", nib.__version__, "| NumPy", np.__version__)
os.makedirs("out", exist_ok=True)

# Phantom in physical (DICOM/ITK LPS) millimetres: (centre, semi-axes, intensity, label)
SHAPES = [((0, 0, 0), (35, 28, 36), 300, 0),     # body
          ((-8, 0, 2), (18, 14, 20), 600, 1),    # organ
          ((-4, 2, 6), (5, 5, 5), 900, 2),       # lesion inside the organ
          ((22, 4, -10), (7, 7, 7), 450, 3)]     # second organ

def phantom(p, gain=1.0):
    """Intensity and label at physical points p (..., 3)."""
    inten = np.zeros(p.shape[:-1]); lab = np.zeros(p.shape[:-1], np.uint8)
    for c, a, val, l in SHAPES:
        inside = (((p - np.array(c)) / np.array(a)) ** 2).sum(-1) <= 1
        inten[inside] = gain * val; lab[inside] = l
    return inten, lab

def grid_points(size, spacing, origin, direction, dz=0.0):
    """Physical points of every voxel centre, array order (z, y, x, 3)."""
    k, j, i = np.meshgrid(*[np.arange(n) for n in size[::-1]], indexing="ij")
    idx = np.stack([i, j, k + dz / spacing[2]], -1) * np.array(spacing)
    return np.array(origin) + idx @ np.array(direction).reshape(3, 3).T

def acquire(size, spacing, origin, direction, gain, seed):
    """Thick-slice acquisition: average over the 4 mm slice profile, add noise."""
    offsets = np.arange(-1.75, 2.0, 0.5)
    img = np.mean([phantom(grid_points(size, spacing, origin, direction, d), gain)[0]
                   for d in offsets], axis=0)
    img += np.random.default_rng(seed).normal(0, 15, img.shape)
    lab = phantom(grid_points(size, spacing, origin, direction))[1]
    out = []
    for arr in (img.astype(np.float32), lab):
        im = sitk.GetImageFromArray(arr)
        im.SetSpacing(spacing); im.SetOrigin(origin); im.SetDirection(direction)
        out.append(im)
    return out

D = (-1, 0, 0, 0, 1, 0, 0, 0, 1)          # first axis runs toward the patient's right
img, lab = acquire((100, 100, 20), (0.8, 0.8, 4.0), (39.6, -39.6, -38.0), D, 1.0, 0)

# --- 1. Orientation from the header -------------------------------------------
sitk.WriteImage(img, "out/scan.nii.gz")
nii = nib.load("out/scan.nii.gz")
print("\n1. header")
print("SimpleITK orientation:", sitk.DICOMOrientImageFilter.GetOrientationFromDirectionCosines(img.GetDirection()))
print("nibabel axis codes:   ", nib.aff2axcodes(nii.affine), "zooms", nii.header.get_zooms())
print("qform_code", int(nii.header["qform_code"]), "sform_code", int(nii.header["sform_code"]))
print("NIfTI affine (RAS+):\n", np.round(nii.affine, 2) + 0.0)
lps = sitk.DICOMOrient(img, "LPS")
ras = nib.as_closest_canonical(nii)
print("lesion centre at index", img.TransformPhysicalPointToContinuousIndex((-4.0, 2.0, 6.0)), "->",
      lps.TransformPhysicalPointToContinuousIndex((-4.0, 2.0, 6.0)), "after DICOMOrient('LPS')")
print("nibabel canonical codes", nib.aff2axcodes(ras.affine))

# --- 2. Resampling to isotropic spacing ---------------------------------------
new_sp = (0.8, 0.8, 0.8)
new_size = [int(round(n * s / t)) for n, s, t in zip(img.GetSize(), img.GetSpacing(), new_sp)]
# keep the physical extent: shift the first centre by half the old minus half the new voxel
shift = [(t - s) / 2 for s, t in zip(img.GetSpacing(), new_sp)]
new_origin = img.TransformContinuousIndexToPhysicalPoint([sh / s for sh, s in zip(shift, img.GetSpacing())])

def resample(im, interp, default=0.0):
    return sitk.Resample(im, new_size, sitk.Transform(), interp, new_origin, new_sp,
                         im.GetDirection(), default, im.GetPixelID())

truth_img, truth_lab = phantom(grid_points(new_size, new_sp, new_origin, D))
body = phantom(grid_points(new_size, new_sp, new_origin, D))[0] > 0
print("\n2. resampling", img.GetSize(), img.GetSpacing(), "->", tuple(new_size), new_sp)
for name, it in [("nearest", sitk.sitkNearestNeighbor), ("linear", sitk.sitkLinear), ("B-spline", sitk.sitkBSpline)]:
    r = sitk.GetArrayFromImage(resample(img, it))
    print(f"image {name:9s} RMSE in body {np.sqrt(np.mean((r - truth_img)[body] ** 2)):6.1f}")

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

near_organ3 = np.linalg.norm(grid_points(new_size, new_sp, new_origin, D) - np.array([22, 4, -10]), axis=-1) < 10
lab_f = sitk.Cast(lab, sitk.sitkFloat32)
onehot = np.stack([sitk.GetArrayFromImage(resample(sitk.Cast(lab == c, sitk.sitkFloat32), sitk.sitkLinear))
                   for c in range(4)])
results = {"nearest": sitk.GetArrayFromImage(resample(lab, sitk.sitkNearestNeighbor)),
           "linear+round": np.rint(sitk.GetArrayFromImage(resample(lab_f, sitk.sitkLinear))).astype(np.uint8),
           "one-hot linear": np.argmax(onehot, axis=0).astype(np.uint8)}
for name, r in results.items():
    d = [dice(r == c, truth_lab == c) for c in (1, 2, 3)]
    spurious = int(np.sum(np.isin(r, (1, 2)) & near_organ3))
    print(f"labels {name:14s} Dice organ {d[0]:.3f} lesion {d[1]:.3f} organ3 {d[2]:.3f} | spurious voxels {spurious}")

# --- 3. Intensity normalization ------------------------------------------------
imgB, labB = acquire((150, 150, 20), (0.8, 0.8, 4.0), (59.6, -59.6, -38.0), D, 1.6, 1)
print("\n3. normalization (mean of organ / lesion)")
for scan, (im, lb) in {"A": (img, lab), "B": (imgB, labB)}.items():
    x = sitk.GetArrayFromImage(im).astype(np.float64); y = sitk.GetArrayFromImage(lb)
    fg = x > threshold_otsu(x)
    schemes = {"raw": x,
               "z whole volume": (x - x.mean()) / x.std(),
               "z in foreground": (x - x[fg].mean()) / x[fg].std()}
    print(f"scan {scan}: foreground fraction {fg.mean():.3f}")
    for name, v in schemes.items():
        print(f"   {name:16s} organ {v[y == 1].mean():8.2f}   lesion {v[y == 2].mean():8.2f}")

SimpleITK 2.5.7 | nibabel 5.4.2 | NumPy 2.5.4

1. header
SimpleITK orientation: RPS
nibabel axis codes:    ('R', 'P', 'S') zooms (np.float32(0.8), np.float32(0.8), np.float32(4.0))
qform_code 1 sform_code 1
NIfTI affine (RAS+):
 [[  0.8   0.    0.  -39.6]
 [  0.   -0.8   0.   39.6]
 [  0.    0.    4.  -38. ]
 [  0.    0.    0.    1. ]]
lesion centre at index (54.5, 52.0, 11.0) -> (44.5, 52.0, 11.0) after DICOMOrient('LPS')
nibabel canonical codes ('R', 'A', 'S')

2. resampling (100, 100, 20) (0.8, 0.8, 4.0) -> (100, 100, 100) (0.8, 0.8, 0.8)
image nearest   RMSE in body   46.7
image linear    RMSE in body   41.0
image B-spline  RMSE in body   39.4
labels nearest        Dice organ 0.958 lesion 0.827 organ3 0.911 | spurious voxels 0
labels linear+round   Dice organ 0.946 lesion 0.564 organ3 0.805 | spurious voxels 1936
labels one-hot linear Dice organ 0.958 lesion 0.827 organ3 0.911 | spurious voxels 0

3. normalization (mean of organ / lesion)
scan A: foreground fraction 0.286
   raw              organ   591.37   lesion   846.18
   z whole volume   organ     2.98   lesion     4.52
   z in foreground  organ     2.25   lesion     4.54
scan B: foreground fraction 0.127
   raw              organ   946.70   lesion  1351.11
   z whole volume   organ     4.54   lesion     6.63
   z in foreground  organ     2.26   lesion     4.55

Header. SimpleITK’s direction code and NiBabel’s axis codes agree, RPS and (R, P, S), because both name the direction in which each axis increases. The matrices disagree in sign exactly as nifti1.h prescribes: SimpleITK holds the direction $(-1, 1, 1)$ and origin $(39.6, -39.6, -38)$ in LPS, while the NIfTI affine has $+0.8$ and $-0.8$ in its first two rows and the offsets $(-39.6, 39.6, -38)$ in RAS+. SimpleITK wrote both a qform and an sform with code 1, scanner coordinates. Reorienting to LPS moved the lesion centre from index 54.5 to 44.5 along the first axis, that is from $i$ to $99 - i$: the axis was flipped and the origin moved, while the physical point stayed where it was. NiBabel’s canonical form flips the second axis instead, to reach RAS.

Images. Twenty slices became a hundred. The ranking is the expected one, with B-spline interpolation slightly ahead of linear and both ahead of nearest neighbour, but the differences are small next to the error that remains. Most of that error comes from the acquisition, by our reasoning, since noise of standard deviation 15 could not produce it alone: each 4 mm slice averaged across organ boundaries, and interpolation between such slices cannot put the boundaries back. Part of the gain of linear and B-spline interpolation, on our reading, is that they average the noise of neighbouring slices, which nearest neighbour copies.

Labels. Linear interpolation of the label values followed by rounding created 1,936 voxels of organ and lesion around the second organ, where the boundary between 0 and 3 passes through the values 1 and 2. Those false lesion voxels lower the lesion’s Dice coefficient by about a third, because the lesion is small. Nearest neighbour and one-hot linear interpolation create no such voxels, and here they give identical results. That is a property of this example rather than a coincidence: only the slice axis was resampled, and along a single axis the class with the largest linearly interpolated weight is always the class of the nearer slice, which is nearest-neighbour interpolation. When several axes are resampled the two methods can differ.

Normalization. The raw organ values differ by the factor of 1.6 between the scans. A z-score over the whole volume does not remove it: scan B’s wider field of view adds air, the foreground falls from 29% to 13% of the volume, and the organ ends up at 4.54 standard deviations in scan B against 2.98 in scan A. Computed inside the foreground mask, the z-score puts the organ and the lesion at the same values in both scans, to within noise. The statistic depends on what it is computed over, which is why nnU-Net restricts it to the nonzero region when cropping removed a large share of the volume. Here the scans differ only by a gain; real scanners also differ in noise, shading and blur, which no global normalization removes.

Variants

nnU-Net’s rules. nnU-Net derives its preprocessing from a fingerprint of the training set: cropping to the nonzero region, a target spacing at the median of the training spacings with a 10th-percentile exception for the coarse axis of anisotropic datasets, third-order splines for images and one-hot linear interpolation for labels, nearest neighbour along the coarse axis of anisotropic cases, and the CT and MRI normalizations described above (Isensee et al., 2020); its article gives the rules in full and compares them with the code. In the paper’s worked example, the cardiac MRI dataset ACDC, the median image is 9 × 237 × 256 voxels at 10 × 1.56 × 1.56 mm, and the out-of-plane target spacing for the 3D network becomes 5 mm, the 10th percentile. The authors say this reduces interpolation artifacts substantially. Why, they state both ways round in consecutive sentences: more images upsampled for training and downsampled for export, then the reverse; only the first fits a target below the median (our reasoning). The 2D configuration does not resample across slices at all.

Toolkits. SimpleITK and NiBabel read and write the formats and provide resampling and reorientation; registration toolkits built on ITK are listed in Medical Image Registration. TorchIO wraps them for PyTorch training and separates spatial transforms, applied to images and label maps alike, from intensity transforms, applied to images only (its random affine transform, for example, always interpolates label maps with nearest neighbour); its preprocessing transforms include reorientation to a canonical orientation, resampling, cropping or padding, rescaling, z-normalization with an optional mask, and histogram standardization (Pérez-García et al., 2021). As of October 2026 these are open-source Python packages; the versions used here are those of the example.

Where preprocessing ends. Random rotations, scalings, intensity changes and simulated artifacts belong to data augmentation, applied anew in each training step, whereas preprocessing is applied once and identically at training and inference. Registration of one scan to another or to an atlas, which goes beyond resampling on a fixed grid, is the subject of medical image registration.

Best Practices

• Keep the header with the data. Carry origin, spacing and direction (or the affine) through every step, and write predictions with the header of the image they were made on, resampled back to its grid.
• Apply the rescale, then mask padding. For CT, convert to Hounsfield units with the Rescale Slope and Intercept, and exclude the Pixel Padding Value before computing any statistic.
• Reorient to one code, explicitly. Choose a target orientation, state which library’s convention names it, and check a known landmark after conversion.
• Resample labels with nearest neighbour or one-hot interpolation, to exactly the image’s grid, and smooth before downsampling images.
• Compute statistics where the signal is. Per-image MRI normalization inside a foreground or brain mask; CT statistics, if learned, from the training set only, so that the test set does not leak into them.
• Apply identical preprocessing at training and inference, and store its parameters (target spacing, clipping range, normalization statistics) with the model.

Pitfalls

• Ignoring Pixel Spacing and slice spacing. Treating voxels as cubes distorts shapes by the anisotropy factor, five in the example, and makes volumes and distances wrong.
• Sorting slices by file name or Slice Location. Neither is a reliable spatial key; the Image Position projected on the slice normal is.
• Mixing conventions. Applying a DICOM direction matrix to NIfTI coordinates without negating x and y swaps left with right and front with back, a 180° rotation about the head-to-foot axis.
• Taking stored values for physical ones. DICOM’s rescale and NIfTI’s scl_slope and scl_inter turn stored integers into meaningful values.
• Linear label interpolation and whole-volume statistics, which the example shows creating labels that do not exist and making a z-score depend on the field of view.
• Bias-field correction of values with zeros or negatives. N4 works in log space and needs positive input.

Evidence and Impact

Direct comparisons of preprocessing choices are fewer than comparisons of architectures, and they are specific to their tasks. Reinhold et al. (2019) synthesized T2-weighted and FLAIR images from T1-weighted brain MRI of 18 healthy subjects from the Kirby-21 dataset, nine for training and nine for testing, all resampled to 1 mm³, corrected with N4 and affinely registered. They compared three synthesis methods (polynomial regression, random forest regression and a U-Net) by normalized cross-correlation, mean SSIM and mutual information, and concluded that normalization substantially improved synthesis, the network most. In T1-to-T2 synthesis, though, they observed the random forest to be qualitatively more robust to unnormalized input, while the network and the polynomial regression failed. The network on unnormalized images did not converge in the first 100 epochs and was trained for 400 instead, and its output still failed on a test image whose histogram differed from those of most training images. No normalization scheme was consistently better than the others in their statistical tests, and the authors note that their data, from one cohort on one scanner, were a best case.

nnU-Net makes the same argument at the level of a whole pipeline: its authors argue that configuration choices such as resampling and normalization, more than architectural changes, decide segmentation performance, and they encode those choices as rules (nnU-Net gives the evidence). Spatial metadata errors, finally, need no experiment: the misalignment that Pérez-García et al. (2021) show when origin, orientation or spacing is ignored is a geometric fact, and so is the effect of label interpolation in the example above.