3D Medical Image Registration using SimpleITK

Introduction

Medical Image Registration is the process of aligning two or more images of the same scene so they perfectly overlap. In radiology, this is critical for comparing a patient's current MRI to a previous scan to measure tumor growth, or fusing an MRI (high soft-tissue contrast) with a CT scan (high bone contrast). In this tutorial, we will use SimpleITK, an incredibly robust, open-source C++ library with a Python wrapper, to perform 3D rigid registration.

Architecture Overview

graph LR
    Fixed(Fixed Image - Reference) --> Reg(ImageRegistrationMethod)
    Moving(Moving Image - To be transformed) --> Reg
    
    Reg --> Metric(Metric: Mutual Information)
    Reg --> Opt(Optimizer: Gradient Descent)
    Reg --> Interp(Interpolator: Linear)
    
    Metric -.-> |Evaluate Similarity| Opt
    Opt -.-> |Update Parameters| Transform(3D Euler Transform)
    Transform -.-> |Apply to Moving| Metric
    
    Transform --> Final(Final Transform Parameters)
    Final --> Resampler(Resample Filter)
    Moving --> Resampler
    Fixed --> Resampler
    
    Resampler --> Registered(Registered Moving Image)
        

Prerequisites

• Python 3.8+
• SimpleITK, NumPy, Matplotlib
• Two 3D medical images in NIfTI format (e.g., fixed_image.nii.gz and moving_image.nii.gz)

Step 1: Install Dependencies

pip install SimpleITK numpy matplotlib

Step 2: Load the Images

In registration, we define a Fixed Image (the reference frame that does not move) and a Moving Image (the image we will transform to align with the fixed image).

import SimpleITK as sitk

# Load the fixed and moving images
fixed_image = sitk.ReadImage("fixed_image.nii.gz", sitk.sitkFloat32)
moving_image = sitk.ReadImage("moving_image.nii.gz", sitk.sitkFloat32)

Step 3: Setup the Registration Framework

Registration requires an Optimizer (how we update the transform parameters), a Metric (how we measure similarity between images), and an Initial Transform.

# Initialize the ImageRegistrationMethod
R = sitk.ImageRegistrationMethod()

# Metric: Mutual Information (great for multi-modal registration like MRI to CT)
R.SetMetricAsMattesMutualInformation(numberOfHistogramBins=50)
R.SetMetricSamplingStrategy(R.RANDOM)
R.SetMetricSamplingPercentage(0.01)

# Optimizer: Gradient Descent
R.SetOptimizerAsGradientDescent(learningRate=1.0, numberOfIterations=100, 
                                convergenceMinimumValue=1e-6, convergenceWindowSize=10)
R.SetOptimizerScalesFromPhysicalShift()

# Interpolator: Linear (used to evaluate pixels at non-grid positions)
R.SetInterpolator(sitk.sitkLinear)

# Initial Transform: Center the moving image on the fixed image
initial_transform = sitk.CenteredTransformInitializer(
    fixed_image, moving_image, 
    sitk.Euler3DTransform(), 
    sitk.CenteredTransformInitializerFilter.GEOMETRY
)
R.SetInitialTransform(initial_transform, inPlace=False)

Step 4: Execute the Registration and Visualize

Now we run the optimization process. SimpleITK will iteratively adjust the 3D rotation and translation parameters until the Mutual Information metric is maximized. Finally, we plot the images overlaid on top of each other using a color composite to visually verify the alignment.

# Run the registration
final_transform = R.Execute(fixed_image, moving_image)

print(f"Final metric value: {R.GetMetricValue()}")
print(f"Optimizer stop condition: {R.GetOptimizerStopConditionDescription()}")

# Apply the final transform to the moving image to resample it into the fixed image space
resampled_moving_image = sitk.Resample(
    moving_image, fixed_image, final_transform, sitk.sitkLinear, 0.0, moving_image.GetPixelID()
)

# Save the registered image
sitk.WriteImage(resampled_moving_image, "registered_moving_image.nii.gz")

# --- Visualization ---
import matplotlib.pyplot as plt
import numpy as np

# Extract a middle 2D slice for visualization
z_slice = fixed_image.GetDepth() // 2
fixed_arr = sitk.GetArrayViewFromImage(fixed_image)[z_slice, :, :]
resampled_arr = sitk.GetArrayViewFromImage(resampled_moving_image)[z_slice, :, :]

# Normalize to 0-255 for RGB visualization
def normalize(arr):
    return ((arr - arr.min()) / (arr.max() - arr.min()) * 255).astype(np.uint8)

fixed_norm = normalize(fixed_arr)
resampled_norm = normalize(resampled_arr)

# Create a composite RGB image: Red channel = Fixed, Green channel = Registered Moving
composite = np.zeros((fixed_norm.shape[0], fixed_norm.shape[1], 3), dtype=np.uint8)
composite[:,:,0] = fixed_norm      # Red
composite[:,:,1] = resampled_norm  # Green

plt.figure(figsize=(8,8))
plt.imshow(composite)
plt.title("Registration Overlay (Yellow indicates perfect overlap)")
plt.axis('off')
plt.show()

Conclusion

SimpleITK makes highly complex 3D image registration accessible via a straightforward, modular Python API. Whether you are aligning longitudinal scans to track tumor progression or fusing multi-modal datasets (PET/CT) to locate metabolic activity anatomically, mastering ITK's registration framework is an absolutely foundational skill for any medical imaging AI engineer.