SimpleITK Image Registration and Analysis

Overview

SimpleITK is a simplified, high-level interface to the Insight Toolkit (ITK) for medical image processing. It provides Python-native access to registration (rigid, affine, B-spline, Demons), segmentation (thresholding, region growing, watershed, level sets), filtering (smoothing, morphology, gradients), and resampling for 3D/4D images from MRI, CT, ultrasound, and fluorescence microscopy. SimpleITK images carry physical space metadata (spacing, origin, direction cosines) which is critical for correct anatomical interpretation and multi-modal alignment.

When to Use

• Registering MRI volumes across timepoints (longitudinal studies) or to a standard atlas for normalization
• Segmenting cells or nuclei from fluorescence microscopy using Otsu thresholding with morphological cleanup
• Converting DICOM series (CT, MRI scanner output) to NIfTI format for downstream analysis with FSL or ANTs
• Applying pre-computed transforms to resample images to a common resolution or field of view
• Computing region statistics (volume, mean intensity, surface area) from binary label masks
• Running multi-modal registration (e.g., aligning PET to MRI) using mutual information metrics
• Use ANTs (via antspyx) instead when you need state-of-the-art diffeomorphic registration with multi-atlas label fusion for neuroimaging research; SimpleITK is better for Python-native scriptable pipelines without native dependencies
• Use scikit-image (scikit-image-processing) instead for 2D bioimage analysis with regionprops, morphological operations, and watershed on non-volumetric fluorescence microscopy data

Prerequisites

• Python packages: SimpleITK>=2.3, numpy, matplotlib
• Optional: SimpleITK-SimpleElastix for additional registration algorithms (Elastix)
• Data requirements: DICOM series (CT/MRI), NIfTI files (.nii or .nii.gz), or any ITK-supported format (MetaImage, NRRD, PNG, TIFF stacks)
• Environment: Python 3.8+; no GPU required; 8 GB RAM recommended for typical 3D volumes

pip install SimpleITK numpy matplotlib

# For additional Elastix-based registration algorithms:
pip install SimpleITK-SimpleElastix

Quick Start

import SimpleITK as sitk

# Read a NIfTI file, apply Gaussian smoothing, and save
image = sitk.ReadImage("brain_t1.nii.gz")
print(f"Size: {image.GetSize()}, Spacing: {image.GetSpacing()}")

smoothed = sitk.SmoothingRecursiveGaussian(image, sigma=1.0)

# Otsu threshold to create a brain mask
mask = sitk.OtsuThreshold(smoothed, 0, 1, 200)
print(f"Voxels in mask: {sitk.GetArrayFromImage(mask).sum()}")

sitk.WriteImage(mask, "brain_mask.nii.gz")
print("Saved brain_mask.nii.gz")

Core API

Module 1: Image I/O

Reading and writing DICOM series, NIfTI, and other formats with full metadata preservation.

import SimpleITK as sitk

# Read a NIfTI file
image = sitk.ReadImage("subject_t1.nii.gz", sitk.sitkFloat32)
print(f"Size (x,y,z): {image.GetSize()}")
print(f"Spacing (mm): {image.GetSpacing()}")
print(f"Origin:       {image.GetOrigin()}")
print(f"Direction:    {image.GetDirection()}")

# Write as compressed NIfTI
sitk.WriteImage(image, "output.nii.gz")
print("Saved output.nii.gz")

import SimpleITK as sitk
import os

# Read a DICOM series from a directory
dicom_dir = "DICOM/series_001/"
series_ids = sitk.ImageSeriesReader.GetGDCMSeriesIDs(dicom_dir)
print(f"Found {len(series_ids)} DICOM series")

reader = sitk.ImageSeriesReader()
reader.SetFileNames(sitk.ImageSeriesReader.GetGDCMSeriesFileNames(dicom_dir, series_ids[0]))
reader.MetaDataDictionaryArrayUpdateOn()
reader.LoadPrivateTagsOn()
volume = reader.Execute()

print(f"DICOM volume size: {volume.GetSize()}")
print(f"Pixel spacing:     {volume.GetSpacing()}")

# Save the 3D volume as NIfTI
sitk.WriteImage(volume, "ct_volume.nii.gz")
print("DICOM series → ct_volume.nii.gz")

Module 2: Image Filtering

Gaussian smoothing, median filtering, gradient magnitude, and edge-preserving filters.

import SimpleITK as sitk
import numpy as np

image = sitk.ReadImage("fluorescence_cells.nii.gz", sitk.sitkFloat32)

# Gaussian smoothing — reduces noise before segmentation
smoothed = sitk.SmoothingRecursiveGaussian(image, sigma=1.5)

# Median filter — removes salt-and-pepper noise (preserves edges better than Gaussian)
median_filtered = sitk.Median(image, [3, 3, 3])

# Gradient magnitude — highlights edges/boundaries
gradient = sitk.GradientMagnitude(smoothed)

arr = sitk.GetArrayFromImage(gradient)
print(f"Gradient range: {arr.min():.2f} – {arr.max():.2f}")
print(f"Mean gradient:  {arr.mean():.4f}")

sitk.WriteImage(smoothed, "smoothed.nii.gz")
sitk.WriteImage(gradient, "gradient.nii.gz")

import SimpleITK as sitk

image = sitk.ReadImage("ct_volume.nii.gz", sitk.sitkFloat32)

# Normalize intensity to [0, 1] range using RescaleIntensity
rescaled = sitk.RescaleIntensity(image, outputMinimum=0.0, outputMaximum=1.0)

# Histogram equalization — improves contrast for registration
equalized = sitk.AdaptiveHistogramEqualization(rescaled)

# N4 bias field correction for MRI (removes B1 field inhomogeneity)
# Cast to float32 for bias correction
image_f32 = sitk.Cast(image, sitk.sitkFloat32)
mask_otsu = sitk.OtsuThreshold(image_f32, 0, 1, 200)
corrected = sitk.N4BiasFieldCorrection(image_f32, mask_otsu)

print("Applied: rescaling, histogram equalization, N4 bias correction")
sitk.WriteImage(corrected, "bias_corrected.nii.gz")

Module 3: Image Registration

Rigid, affine, and deformable (B-spline, Demons) registration using ImageRegistrationMethod.

import SimpleITK as sitk

# Load fixed (reference/atlas) and moving (subject to align) images
fixed = sitk.ReadImage("atlas_t1.nii.gz", sitk.sitkFloat32)
moving = sitk.ReadImage("subject_t1.nii.gz", sitk.sitkFloat32)

# Set up rigid registration
registration_method = sitk.ImageRegistrationMethod()

# Similarity metric: Mattes mutual information (works for same-modality)
registration_method.SetMetricAsMattesMutualInformation(numberOfHistogramBins=50)
registration_method.SetMetricSamplingStrategy(registration_method.RANDOM)
registration_method.SetMetricSamplingPercentage(0.01)

# Optimizer: gradient descent with line search
registration_method.SetOptimizerAsGradientDescent(
    learningRate=1.0, numberOfIterations=100,
    convergenceMinimumValue=1e-6, convergenceWindowSize=10
)
registration_method.SetOptimizerScalesFromPhysicalShift()

# Multi-resolution pyramid: 3 levels → faster convergence
registration_method.SetShrinkFactorsPerLevel(shrinkFactors=[4, 2, 1])
registration_method.SetSmoothingSigmasPerLevel(smoothingSigmas=[2, 1, 0])
registration_method.SmoothingSigmasAreSpecifiedInPhysicalUnitsOn()

# Initialize with center of geometry
initial_transform = sitk.CenteredTransformInitializer(
    fixed, moving,
    sitk.Euler3DTransform(),
    sitk.CenteredTransformInitializerFilter.GEOMETRY
)
registration_method.SetInitialTransform(initial_transform, inPlace=False)
registration_method.SetInterpolator(sitk.sitkLinear)

# Execute registration
final_transform = registration_method.Execute(fixed, moving)
print(f"Optimizer stop: {registration_method.GetOptimizerStopConditionDescription()}")
print(f"Final metric:   {registration_method.GetMetricValue():.4f}")

# Apply transform and save
resampled = sitk.Resample(
    moving, fixed, final_transform,
    sitk.sitkLinear, 0.0, moving.GetPixelID()
)
sitk.WriteImage(resampled, "subject_registered.nii.gz")
sitk.WriteTransform(final_transform, "rigid_transform.tfm")
print("Saved: subject_registered.nii.gz, rigid_transform.tfm")

import SimpleITK as sitk

# Deformable B-spline registration for non-linear alignment
fixed = sitk.ReadImage("atlas_t1.nii.gz", sitk.sitkFloat32)
moving = sitk.ReadImage("subject_t1.nii.gz", sitk.sitkFloat32)

# Start with affine pre-registration
affine_method = sitk.ImageRegistrationMethod()
affine_method.SetMetricA