simpleitk-image-registration

SimpleITK 图像配准与分析

概述

SimpleITK 是 Insight Toolkit（ITK）用于医学图像处理的简化高级接口。它提供原生 Python 支持，可用于配准（刚性、仿射、B-spline、Demons）、分割（阈值处理、区域生长、分水岭、水平集）、滤波（平滑、形态学、梯度）以及对来自 MRI、CT、超声和荧光显微镜的 3D/4D 图像进行重采样。SimpleITK 图像包含物理空间元数据（间距、原点、方向余弦），这对于正确解读解剖结构和进行多模态对齐至关重要。

适用场景

• 对不同时间点的 MRI 体数据进行配准（纵向研究），或将其配准到标准图谱以进行归一化
• 使用 Otsu 阈值处理和形态学清理，从荧光显微镜图像中分割细胞或细胞核
• 将 DICOM 序列（CT、MRI 扫描仪输出）转换为 NIfTI 格式，以便使用 FSL 或 ANTs 进行下游分析
• 应用预先计算的变换，将图像重采样到统一分辨率或视野
• 根据二值标签掩膜计算区域统计信息（体积、平均强度、表面积）
• 使用互信息度量运行多模态配准（例如，将 PET 与 MRI 对齐）
• 如果需要用于神经影像研究的先进微分同胚配准和多图谱标签融合，请改用 ANTs（通过 antspyx）；SimpleITK 更适合无需原生依赖的原生 Python 可脚本化流水线
• 对非体数据的荧光显微镜数据进行 2D 生物图像分析，并使用 regionprops、形态学操作和分水岭时，请改用 scikit-image（scikit-image-processing）

前置条件

• Python 包：SimpleITK>=2.3、numpy、matplotlib
• 可选：SimpleITK-SimpleElastix，用于使用额外的配准算法（Elastix）
• 数据要求：DICOM 序列（CT/MRI）、NIfTI 文件（.nii 或 .nii.gz），或任何 ITK 支持的格式（MetaImage、NRRD、PNG、TIFF 堆栈）
• 环境：Python 3.8+；无需 GPU；典型 3D 体数据建议使用 8 GB RAM

pip install SimpleITK numpy matplotlib

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

快速开始

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")

核心 API

模块 1：图像 I/O

在完整保留元数据的同时，读取和写入 DICOM 序列、NIfTI 及其他格式。

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")

模块 2：图像滤波

高斯平滑、中值滤波、梯度幅值和保边滤波器。

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")

模块 3：图像配准

使用 ImageRegistrationMethod 进行刚性、仿射和可变形（B-spline、Demons）配准。

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.SetMetricAsMattesMutualInformation(numberOfHistogramBins=50)
affine_method.SetMetricSamplingStrategy(affine_method.RANDOM)
affine_method.SetMetricSamplingPercentage(0.01)
affine_method.SetOptimizerAsGradientDescent(
    learningRate=1.0, numberOfIterations=100,
    convergenceMinimumValue=1e-6, convergenceWindowSize=10
)
affine_method.SetShrinkFactorsPerLevel([4, 2, 1])
affine_method.SetSmoothingSigmasPerLevel([2, 1, 0])
affine_method.SmoothingSigmasAreSpecifiedInPhysicalUnitsOn()
affine_init = sitk.CenteredTransformInitializer(
    fixed, moving, sitk.AffineTransform(3),
    sitk.CenteredTransformInitializerFilter.GEOMETRY
)
affine_method.SetInitialTransform(affine_init, inPlace=False)
affine_method.SetInterpolator(sitk.sitkLinear)
affine_transform = affine_method.Execute(fixed, moving)
print(f"Affine complete: metric = {affine_method.GetMetricValue():.4f}")

# B-spline deformable refinement
bspline_method = sitk.ImageRegistrationMethod()
bspline_method.SetMetricAsMattesMutualInformation(numberOfHistogramBins=50)
bspline_method.SetMetricSamplingStrategy(bspline_method.RANDOM)
bspline_method.SetMetricSamplingPercentage(0.01)
bspline_method.SetOptimizerAsGradientDescent(
    learningRate=0.5, numberOfIterations=50,
    convergenceMinimumValue=1e-6, convergenceWindowSize=10
)
bspline_method.SetShrinkFactorsPerLevel([2, 1])
bspline_method.SetSmoothingSigmasPerLevel([1, 0])
bspline_method.SmoothingSigmasAreSpecifiedInPhysicalUnitsOn()

mesh_size = [8] * fixed.GetDimension()
bspline_init = sitk.BSplineTransformInitializer(image1=fixed, transformDomainMeshSize=mesh_size, order=3)
composite = sitk.CompositeTransform(3)
composite.AddTransform(affine_transform)
composite.AddTransform(bspline_init)
bspline_method.SetInitialTransform(composite, inPlace=True)
bspline_method.SetInterpolator(sitk.sitkLinear)

deformable_transform = bspline_method.Execute(fixed, moving)
resampled = sitk.Resample(moving, fixed, deformable_transform, sitk.sitkLinear, 0.0)
sitk.WriteImage(resampled, "subject_deformable_registered.nii.gz")
print("Saved: subject_deformable_registered.nii.gz")

模块 4：分割

Otsu 阈值处理、区域生长、分水岭和形态学后处理。

import SimpleITK as sitk
import numpy as np

# Read fluorescence microscopy image
image = sitk.ReadImage("cells_gfp.nii.gz", sitk.sitkFloat32)

# Gaussian smoothing before thresholding
smoothed = sitk.SmoothingRecursiveGaussian(image, sigma=1.0)

# Otsu thresholding: automatically finds optimal foreground/background threshold
# Returns binary label image: 1=foreground (cells), 0=background
binary = sitk.OtsuThreshold(smoothed, insideValue=0, outsideValue=1, numberOfHistogramBins=200)

# Count segmented voxels
arr = sitk.GetArrayFromImage(binary)
n_foreground = arr.sum()
total = arr.size
print(f"Foreground: {n_foreground} voxels ({100*n_foreground/total:.1f}%)")

sitk.WriteImage(binary, "cells_binary_otsu.nii.gz")
print("Saved: cells_binary_otsu.nii.gz")

import SimpleITK as sitk

image = sitk.ReadImage("mri_brain.nii.gz", sitk.sitkFloat32)
smoothed = sitk.SmoothingRecursiveGaussian(image, sigma=0.5)

# Region growing from a seed point: ConnectedThreshold
# Seeds must be inside the region of interest; lower/upper bound in image intensity units
seed = (128, 128, 64)   # (x, y, z) in image coordinates
lower_threshold = 50.0
upper_threshold = 200.0
region_grown = sitk.ConnectedThreshold(
    smoothed,
    seedList=[seed],
    lower=lower_threshold,
    upper=upper_threshold,
    replaceValue=1
)
print(f"Region growing: {sitk.GetArrayFromImage(region_grown).sum()} voxels segmented")

# ConfidenceConnected: adaptive thresholding based on neighborhood statistics
confidence_seg = sitk.ConfidenceConnected(
    smoothed,
    seedList=[seed],
    numberOfIterations=2,
    multiplier=2.5,       # number of standard deviations from mean
    initialNeighborhoodRadius=2,
    replaceValue=1
)
print(f"Confidence connected: {sitk.GetArrayFromImage(confidence_seg).sum()} voxels")

sitk.WriteImage(region_grown, "region_grown.nii.gz")
sitk.WriteImage(confidence_seg, "confidence_seg.nii.gz")

import SimpleITK as sitk

# Morphological operations for binary mask cleanup
binary = sitk.ReadImage("cells_binary_otsu.nii.gz")

# Fill small holes
filled = sitk.BinaryFillhole(binary)

# Erosion: remove thin protrusions and separate touching objects
eroded = sitk.BinaryErode(filled, kernelRadius=(2, 2, 1), kernelType=sitk.sitkBall)

# Dilation: restore true boundary after erosion
dilated = sitk.BinaryDilate(eroded, kernelRadius=(2, 2, 1), kernelType=sitk.sitkBall)

# Opening = erosion + dilation: removes small objects
opened = sitk.BinaryMorphologicalOpening(binary, kernelRadius=(1, 1, 1), kernelType=sitk.sitkBall)

# Connected component labeling: assign unique integer to each object
labeled = sitk.ConnectedComponent(opened)
n_objects = sitk.GetArrayFromImage(labeled).max()
print(f"Connected components (objects): {n_objects}")

# Remove objects smaller than 100 voxels
relabeled = sitk.RelabelComponent(labeled, minimumObjectSize=100)
n_kept = sitk.GetArrayFromImage(relabeled).max()
print(f"Objects remaining after size filter: {n_kept}")

sitk.WriteImage(relabeled, "cells_labeled.nii.gz")

模块 5：重采样与变换应用

应用变换、重采样到参考网格，以及在数组和图像表示之间进行转换。

import SimpleITK as sitk
import numpy as np

# Resample moving image to match reference grid
reference = sitk.ReadImage("reference_volume.nii.gz")
moving = sitk.ReadImage("subject_volume.nii.gz")

# Load a pre-computed transform
transform = sitk.ReadTransform("rigid_transform.tfm")

# Resample with linear interpolation (use nearest neighbor for label images)
resampled = sitk.Resample(
    moving,
    reference,
    transform,
    sitk.sitkLinear,    # interpolator; use sitk.sitkNearestNeighbor for labels
    0.0,                # default pixel value for regions outside the moving image
    moving.GetPixelID()
)
print(f"Resampled size: {resampled.GetSize()} (matches reference: {reference.GetSize()})")
sitk.WriteImage(resampled, "resampled_to_reference.nii.gz")

import SimpleITK as sitk
import numpy as np

# Convert between SimpleITK image and NumPy array
image = sitk.ReadImage("ct_volume.nii.gz", sitk.sitkFloat32)

# SimpleITK uses (x, y, z) ordering; NumPy gets (z, y, x) ordering
arr = sitk.GetArrayFromImage(image)
print(f"NumPy array shape (z,y,x): {arr.shape}")
print(f"Value range: {arr.min():.1f} – {arr.max():.1f}")

# Apply NumPy operations
clipped = np.clip(arr, -1000, 3000)   # CT Hounsfield unit clipping
normalized = (clipped - clipped.mean()) / clipped.std()

# Convert back to SimpleITK image (preserving spatial metadata)
out_image = sitk.GetImageFromArray(normalized)
out_image.CopyInformation(image)   # copy spacing, origin, direction from original
print(f"Output size: {out_image.GetSize()}, Spacing: {out_image.GetSpacing()}")

# Resample to isotropic 1mm spacing
new_spacing = [1.0, 1.0, 1.0]
orig_spacing = image.GetSpacing()
orig_size = image.GetSize()
new_size = [
    int(round(orig_size[i] * orig_spacing[i] / new_spacing[i]))
    for i in range(3)
]
resampled_iso = sitk.Resample(
    image,
    new_size,
    sitk.Transform(),         # identity transform
    sitk.sitkLinear,
    image.GetOrigin(),
    new_spacing,
    image.GetDirection(),
    0.0,
    image.GetPixelID()
)
print(f"Isotropic resampled size: {resampled_iso.GetSize()}")
sitk.WriteImage(resampled_iso, "isotropic_1mm.nii.gz")

模块 6：统计与测量

标签形状统计信息（体积、表面积、质心、边界框）和各区域的强度统计信息。

import SimpleITK as sitk
import pandas as pd

# Load intensity image and binary label mask
intensity = sitk.ReadImage("fluorescence_cells.nii.gz", sitk.sitkFloat32)
labels = sitk.ReadImage("cells_labeled.nii.gz", sitk.sitkUInt32)

# Shape statistics: geometric measurements per label
shape_filter = sitk.LabelShapeStatisticsImageFilter()
shape_filter.ComputeOrientedBoundingBoxOn()
shape_filter.Execute(labels)

label_ids = shape_filter.GetLabels()
print(f"Number of labeled objects: {len(label_ids)}")

rows = []
for lbl in label_ids:
    rows.append({
        "label": lbl,
        "volume_voxels": shape_filter.GetNumberOfPixels(lbl),
        "volume_mm3": shape_filter.GetPhysicalSize(lbl),
        "centroid_x": shape_filter.GetCentroid(lbl)[0],
        "centroid_y": shape_filter.GetCentroid(lbl)[1],
        "centroid_z": shape_filter.GetCentroid(lbl)[2],
        "elongation": shape_filter.GetElongation(lbl),
        "roundness": shape_filter.GetRoundness(lbl),
    })

shape_df = pd.DataFrame(rows)
print(shape_df.head())
print(f"\nMean volume: {shape_df['volume_mm3'].mean():.1f} mm³")
shape_df.to_csv("cell_shape_stats.csv", index=False)
print("Saved: cell_shape_stats.csv")

import SimpleITK as sitk
import pandas as pd

intensity = sitk.ReadImage("fluorescence_cells.nii.gz", sitk.sitkFloat32)
labels = sitk.ReadImage("cells_labeled.nii.gz", sitk.sitkUInt32)

# Intensity statistics per label
intensity_filter = sitk.LabelIntensityStatisticsImageFilter()
intensity_filter.Execute(labels, intensity)

label_ids = intensity_filter.GetLabels()
rows = []
for lbl in label_ids:
    rows.append({
        "label": lbl,
        "mean_intensity": intensity_filter.GetMean(lbl),
        "std_intensity": intensity_filter.GetSigma(lbl),
        "min_intensity": intensity_filter.GetMinimum(lbl),
        "max_intensity": intensity_filter.GetMaximum(lbl),
        "median_intensity": intensity_filter.GetMedian(lbl),
        "sum_intensity": intensity_filter.GetSum(lbl),
    })

intensity_df = pd.DataFrame(rows)
print(intensity_df.describe().round(2))
intensity_df.to_csv("cell_intensity_stats.csv", index=False)
print(f"\nSaved: cell_intensity_stats.csv ({len(rows)} cells)")

核心概念

物理空间与像素空间

SimpleITK 图像存储空间元数据（以 mm 为单位的间距、原点、方向余弦），用于定义像素坐标如何映射到物理（世界）坐标。Resample 和配准等操作在物理空间中进行，即使图像的体素大小或方向不同，也能确保解剖结构正确对齐。基于 NumPy 数组创建输出图像时，务必使用 CopyInformation() 保留物理元数据。

import SimpleITK as sitk
import numpy as np

image = sitk.ReadImage("mri.nii.gz")
print(f"Pixel space size:     {image.GetSize()}")       # (x, y, z) in voxels
print(f"Physical spacing mm:  {image.GetSpacing()}")    # mm per voxel
print(f"Physical origin mm:   {image.GetOrigin()}")     # world coords of first voxel
print(f"Direction cosines:    {image.GetDirection()}")  # patient orientation matrix

# Convert pixel index to physical point
pixel_idx = (64, 64, 32)
physical_pt = image.TransformIndexToPhysicalPoint(pixel_idx)
print(f"Pixel {pixel_idx} → physical {physical_pt}")

# NumPy array has REVERSED axis order: (z, y, x)
arr = sitk.GetArrayFromImage(image)
print(f"NumPy shape: {arr.shape}")  # (z, y, x)

变换组合

SimpleITK 变换可以使用 CompositeTransform 组合，以应用一系列变换（例如，先进行仿射预对齐，再进行可变形 B-spline 精细调整）。调用 Resample 时，变换会按添加顺序应用，从而便于构建模块化流水线。

常见工作流

工作流 1：荧光显微镜细胞分割

目标：从 3D 荧光显微镜体数据中分割单个细胞，进行形态学清理，并提取每个细胞的测量值。

import SimpleITK as sitk
import pandas as pd
import numpy as np

# 1. Load fluorescence volume
image = sitk.ReadImage("cells_gfp.nii.gz", sitk.sitkFloat32)
print(f"Image size: {image.GetSize()}, Spacing: {image.GetSpacing()}")

# 2. Denoise with Gaussian smoothing
smoothed = sitk.SmoothingRecursiveGaussian(image, sigma=0.8)

# 3. Otsu threshold to get binary foreground mask
binary = sitk.OtsuThreshold(smoothed, insideValue=0, outsideValue=1, numberOfHistogramBins=200)
n_fg = sitk.GetArrayFromImage(binary).sum()
print(f"Foreground voxels after Otsu: {n_fg}")

# 4. Morphological cleanup
filled = sitk.BinaryFillhole(binary)
opened = sitk.BinaryMorphologicalOpening(filled, kernelRadius=(1, 1, 1))

# 5. Label individual connected components (cells)
labeled = sitk.ConnectedComponent(opened)
relabeled = sitk.RelabelComponent(labeled, minimumObjectSize=200)
n_cells = int(sitk.GetArrayFromImage(relabeled).max())
print(f"Detected cells (>200 voxels): {n_cells}")

# 6. Compute shape statistics
shape_filter = sitk.LabelShapeStatisticsImageFilter()
shape_filter.Execute(relabeled)

intensity_filter = sitk.LabelIntensityStatisticsImageFilter()
intensity_filter.Execute(relabeled, image)

rows = []
for lbl in shape_filter.GetLabels():
    rows.append({
        "cell_id": lbl,
        "volume_mm3": shape_filter.GetPhysicalSize(lbl),
        "roundness": shape_filter.GetRoundness(lbl),
        "elongation": shape_filter.GetElongation(lbl),
        "mean_gfp_intensity": intensity_filter.GetMean(lbl),
        "total_gfp_intensity": intensity_filter.GetSum(lbl),
    })

df = pd.DataFrame(rows)
print(df.describe().round(3))

# 7. Save outputs
sitk.WriteImage(relabeled, "cells_labeled.nii.gz")
df.to_csv("cell_measurements.csv", index=False)
print(f"\nSaved: cells_labeled.nii.gz, cell_measurements.csv ({n_cells} cells)")

工作流 2：MRI 脑图谱配准

目标：使用刚性 + 仿射配准，将受试者的 T1 MRI 配准到标准脑图谱，然后将图谱标签转换到受试者空间。

import SimpleITK as sitk
import numpy as np

# 1. Load atlas (fixed) and subject (moving) T1 MRI
atlas_t1 = sitk.ReadImage("mni152_t1.nii.gz", sitk.sitkFloat32)
subject_t1 = sitk.ReadImage("subject_t1.nii.gz", sitk.sitkFloat32)
atlas_labels = sitk.ReadImage("mni152_labels.nii.gz", sitk.sitkUInt16)
print(f"Atlas size: {atlas_t1.GetSize()}, Subject size: {subject_t1.GetSize()}")

# 2. Normalize intensities for better registration convergence
atlas_norm = sitk.RescaleIntensity(atlas_t1, 0.0, 1.0)
subject_norm = sitk.RescaleIntensity(subject_t1, 0.0, 1.0)

# 3. Initialize with center-of-geometry alignment
initial_transform = sitk.CenteredTransformInitializer(
    atlas_norm, subject_norm,
    sitk.AffineTransform(3),
    sitk.CenteredTransformInitializerFilter.GEOMETRY
)

# 4. Configure and run affine registration
reg = sitk.ImageRegistrationMethod()
reg.SetMetricAsMattesMutualInformation(numberOfHistogramBins=50)
reg.SetMetricSamplingStrategy(reg.RANDOM)
reg.SetMetricSamplingPercentage(0.01)
reg.SetOptimizerAsGradientDescent(
    learningRate=1.0, numberOfIterations=200,
    convergenceMinimumValue=1e-6, convergenceWindowSize=10
)
reg.SetOptimizerScalesFromPhysicalShift()
reg.SetShrinkFactorsPerLevel([4, 2, 1])
reg.SetSmoothingSigmasPerLevel([2, 1, 0])
reg.SmoothingSigmasAreSpecifiedInPhysicalUnitsOn()
reg.SetInitialTransform(initial_transform, inPlace=False)
reg.SetInterpolator(sitk.sitkLinear)

final_transform = reg.Execute(atlas_norm, subject_norm)
print(f"Registration metric: {reg.GetMetricValue():.4f}")
print(f"Optimizer: {reg.GetOptimizerStopConditionDescription()}")

# 5. Apply transform to subject T1 (aligned to atlas space)
subject_registered = sitk.Resample(
    subject_t1, atlas_t1, final_transform,
    sitk.sitkLinear, 0.0, subject_t1.GetPixelID()
)

# 6. Invert transform to bring atlas labels to subject space
inverse_transform = final_transform.GetInverse()
labels_in_subject_space = sitk.Resample(
    atlas_labels, subject_t1, inverse_transform,
    sitk.sitkNearestNeighbor, 0, atlas_labels.GetPixelID()
)

# 7. Save results
sitk.WriteImage(subject_registered, "subject_in_atlas_space.nii.gz")
sitk.WriteImage(labels_in_subject_space, "atlas_labels_in_subject_space.nii.gz")
sitk.WriteTransform(final_transform, "subject_to_atlas.tfm")

n_labels = int(sitk.GetArrayFromImage(labels_in_subject_space).max())
print(f"Transferred {n_labels} atlas regions to subject space")
print("Saved: subject_in_atlas_space.nii.gz, atlas_labels_in_subject_space.nii.gz")

工作流 3：DICOM 序列批量转换为 NIfTI

目标：将扫描仪目录中的多个 DICOM 序列转换为 NIfTI，并进行各向同性重采样。

import SimpleITK as sitk
from pathlib import Path

def convert_dicom_series(dicom_dir: str, output_path: str, target_spacing_mm: float = 1.0) -> dict:
    """Convert a DICOM series directory to NIfTI with optional isotropic resampling."""
    series_ids = sitk.ImageSeriesReader.GetGDCMSeriesIDs(dicom_dir)
    if not series_ids:
        return {"error": f"No DICOM series in {dicom_dir}"}

    reader = sitk.ImageSeriesReader()
    reader.SetFileNames(sitk.ImageSeriesReader.GetGDCMSeriesFileNames(dicom_dir, series_ids[0]))
    volume = reader.Execute()
    orig_size = volume.GetSize()
    orig_spacing = volume.GetSpacing()

    # Resample to isotropic voxels
    new_spacing = [target_spacing_mm] * 3
    new_size = [int(round(orig_size[i] * orig_spacing[i] / target_spacing_mm)) for i in range(3)]
    resampled = sitk.Resample(
        volume, new_size, sitk.Transform(),
        sitk.sitkLinear, volume.GetOrigin(),
        new_spacing, volume.GetDirection(),
        0.0, volume.GetPixelID()
    )

    sitk.WriteImage(resampled, output_path)
    return {
        "input": dicom_dir, "output": output_path,
        "orig_size": orig_size, "orig_spacing": orig_spacing,
        "new_size": new_size, "n_series": len(series_ids)
    }

# Batch convert all subject directories
input_root = Path("DICOM_data/")
output_root = Path("NIfTI_converted/")
output_root.mkdir(exist_ok=True)

results = []
for subject_dir in sorted(input_root.iterdir()):
    if subject_dir.is_dir():
        out_file = output_root / f"{subject_dir.name}_t1.nii.gz"
        info = convert_dicom_series(str(subject_dir), str(out_file), target_spacing_mm=1.0)
        results.append(info)
        print(f"  {subject_dir.name}: {info.get('orig_size')} → {info.get('new_size')}")

print(f"\nConverted {len(results)} subjects to {output_root}/")

关键参数

参数	模块	默认值	范围 / 选项	效果
sigma	SmoothingRecursiveGaussian	—	0.5–5.0 mm	高斯平滑宽度；值越大，模糊越明显，降噪效果越强
numberOfHistogramBins	OtsuThreshold, MutualInfo	128	50–256	用于阈值/指标计算的直方图分辨率
numberOfIterations	ImageRegistrationMethod	100	50–500	优化器的最大步数；复杂或精细配准时可增大此值
learningRate	GradientDescent optimizer	1.0	0.01–2.0	每次优化迭代的步长；过大会导致发散
shrinkFactors	Multi-resolution pyramid	[4,2,1]	[8,4,2,1] 列表	每个分辨率层级的图像降采样；值越大，速度越快但结果越粗略
minimumObjectSize	RelabelComponent	0	1–10000 个体素	移除小于此体素数的连通分量
kernelRadius	BinaryErode/Dilate	(1,1,1)	(1,1,1)–(5,5,5)	各轴的形态学核半径，单位为体素
interpolator	Resample	sitk.sitkLinear	sitkLinear、sitkNearestNeighbor、sitkBSpline	插值方法；整数标签图像应使用 NearestNeighbor
multiplier	ConfidenceConnected	2.5	1.5–4.0	接受体素时允许其偏离邻域均值的标准差数
meshSize	BSplineTransformInitializer	[8]*ndim	[4]*ndim–[20]*ndim	B 样条控制点网格；值越大，可变形自由度越高

最佳实践

1. 配准或滤波前始终转换为 float32：许多 ITK 滤波器需要浮点输入。DICOM 中的整数（通常为 sitk.sitkInt16）可能导致无提示错误或精度损失。

   image = sitk.ReadImage("ct.nii.gz", sitk.sitkFloat32)   # explicit cast at load time
   # Or: image_f32 = sitk.Cast(image, sitk.sitkFloat32)
   

2. 从 NumPy 数组创建图像时使用 CopyInformation()：否则，输出图像会丢失物理间距和原点信息，导致无法在查看器中与源图像叠加。

   arr = sitk.GetArrayFromImage(image)
   modified = some_numpy_operation(arr)
   out = sitk.GetImageFromArray(modified)
   out.CopyInformation(image)  # preserve spacing, origin, direction
   

3. 对标签/掩码重采样使用 sitkNearestNeighbor：对整数标签图像进行线性或 B 样条插值会产生小数标签值，从而破坏掩码。二值掩码和多标签分割始终使用最近邻插值。

4. 使用多分辨率配准以提高速度和稳健性：设置 SetShrinkFactorsPerLevel([4,2,1]) 和 SetSmoothingSigmasPerLevel([2,1,0])，防止优化器在大型 3D 体数据上陷入局部极小值。

5. 使用棋盘格对比进行配准结果的目视验证：在信任配准结果前，使用 sitk.CheckerBoard 检查固定图像与配准后的移动图像是否对齐。

   checker = sitk.CheckerBoard(fixed, resampled_moving, [5, 5, 5])
   sitk.WriteImage(checker, "registration_checker.nii.gz")
   

常见用法

用法：多模态配准（MRI 到 CT）

适用场景：对齐强度不兼容的不同成像模态（例如，PET 到 MRI、CT 到 MRI），需要使用互信息作为指标。

import SimpleITK as sitk

# Load images from different modalities
ct = sitk.ReadImage("patient_ct.nii.gz", sitk.sitkFloat32)
mri = sitk.ReadImage("patient_mri.nii.gz", sitk.sitkFloat32)

# Mutual information metric handles multi-modal intensity relationships
reg = sitk.ImageRegistrationMethod()
reg.SetMetricAsMattesMutualInformation(numberOfHistogramBins=100)
reg.SetMetricSamplingStrategy(reg.RANDOM)
reg.SetMetricSamplingPercentage(0.05)
reg.SetOptimizerAsGradientDescent(
    learningRate=1.0, numberOfIterations=150,
    convergenceMinimumValue=1e-6, convergenceWindowSize=10
)
reg.SetOptimizerScalesFromPhysicalShift()
reg.SetShrinkFactorsPerLevel([4, 2, 1])
reg.SetSmoothingSigmasPerLevel([2, 1, 0])
reg.SmoothingSigmasAreSpecifiedInPhysicalUnitsOn()
initial = sitk.CenteredTransformInitializer(
    ct, mri, sitk.AffineTransform(3),
    sitk.CenteredTransformInitializerFilter.GEOMETRY
)
reg.SetInitialTransform(initial, inPlace=False)
reg.SetInterpolator(sitk.sitkLinear)

transform = reg.Execute(ct, mri)
mri_aligned = sitk.Resample(mri, ct, transform, sitk.sitkLinear, 0.0)
sitk.WriteImage(mri_aligned, "mri_aligned_to_ct.nii.gz")
print(f"Multi-modal registration complete. Metric: {reg.GetMetricValue():.4f}")

用法：批量将变换应用于标签掩码

适用场景：完成一幅图像的配准后，将相同变换应用于分割掩码或图谱标签（使用最近邻插值）。

import SimpleITK as sitk
from pathlib import Path

# Load the transform computed during registration
transform = sitk.ReadTransform("subject_to_atlas.tfm")
reference = sitk.ReadImage("atlas_t1.nii.gz")

# Apply to all label files in a directory
label_dir = Path("subject_labels/")
output_dir = Path("atlas_labels/")
output_dir.mkdir(exist_ok=True)

for label_file in sorted(label_dir.glob("*.nii.gz")):
    label = sitk.ReadImage(str(label_file), sitk.sitkUInt16)
    registered_label = sitk.Resample(
        label, reference, transform,
        sitk.sitkNearestNeighbor,    # critical: integer labels need nearest-neighbor
        0, label.GetPixelID()
    )
    out_path = output_dir / label_file.name
    sitk.WriteImage(registered_label, str(out_path))
    n_labels = int(sitk.GetArrayFromImage(registered_label).max())
    print(f"  {label_file.name}: {n_labels} labels → {out_path.name}")

print(f"Batch transform applied to {len(list(label_dir.glob('*.nii.gz')))} files")

用法：Demons 可变形配准

适用场景：对相同模态且形变较小的图像进行快速可变形配准（例如，存在轻微萎缩的纵向脑 MRI）。

import SimpleITK as sitk

fixed = sitk.ReadImage("baseline_mri.nii.gz", sitk.sitkFloat32)
moving = sitk.ReadImage("followup_mri.nii.gz", sitk.sitkFloat32)

# Normalize intensities
fixed_n = sitk.RescaleIntensity(fixed, 0.0, 1.0)
moving_n = sitk.RescaleIntensity(moving, 0.0, 1.0)

# Demons registration filter
demons = sitk.DemonsRegistrationFilter()
demons.SetNumberOfIterations(50)
demons.SetStandardDeviations(1.0)   # displacement field smoothing

displacement_field = demons.Execute(fixed_n, moving_n)
print(f"Demons complete: RMS change = {demons.GetRMSChange():.6f}")

# Apply displacement field transform
displacement_transform = sitk.DisplacementFieldTransform(displacement_field)
warped = sitk.Resample(moving, fixed, displacement_transform, sitk.sitkLinear, 0.0)

sitk.WriteImage(warped, "followup_registered_demons.nii.gz")
sitk.WriteImage(
    sitk.Cast(sitk.VectorMagnitude(displacement_field), sitk.sitkFloat32),
    "deformation_magnitude.nii.gz"
)
print("Saved: followup_registered_demons.nii.gz, deformation_magnitude.nii.gz")

故障排查

问题	原因	解决方案
TypeError: in method 'ImageRegistrationMethod_Execute'	图像像素类型不是浮点型	配准前转换为 float32：sitk.Cast(image, sitk.sitkFloat32)
配准发散（度量值增加）	学习率过高或图像未预先对齐	将 learningRate 降至 0.1；使用 CenteredTransformInitializer 进行初始对齐；检查图像是否重叠
重采样后标签图像出现小数值	掩膜使用了错误的插值器	对整数标签图像进行重采样时，将插值器设为 sitk.sitkNearestNeighbor
RuntimeError: Exception thrown in SimpleITK ReadImage	格式不受支持或 DICOM 已损坏	使用 dcmdump 验证文件；对于 DICOM 序列，使用 ImageSeriesReader.GetGDCMSeriesFileNames()
N4 偏置校正非常慢	图像过大或迭代次数过多	先降采样：校正前将图像重采样为 2 mm 等体素，然后将变换应用于原始分辨率图像
连通分量标签过多或过少	对象大小阈值不正确或预处理不足	调整 RelabelComponent 中的 minimumObjectSize；增大 Gaussian sigma 以合并相接的对象
NumPy 数组形状与预期的 (z,y,x) 不符	轴顺序混淆	SimpleITK 图像的轴顺序为 (x,y,z)；GetArrayFromImage() 返回轴顺序为 (z,y,x) 的 NumPy 数组；使用 arr.transpose(2,1,0) 进行转换
Otsu 阈值分割出的背景过多	背景较亮（例如，共聚焦反射伪影）	使用 sitk.OtsuMultipleThresholds 并设置 numberOfThresholds=2；或者在阈值分割前应用前景掩膜

相关技能

• cellpose-cell-segmentation — 基于荧光显微镜图像的深度学习细胞分割；当经典阈值法无法处理异质染色或相接细胞时使用
• scikit-image-processing — 使用 regionprops、watershed 和滤波器进行 2D 图像分析；适用于非体积的 2D 显微镜图像分析
• napari-image-viewer — 交互式可视化 3D 体数据和标签掩膜；可与 SimpleITK 配合使用，以目视检查配准和分割结果
• pydicom-medical-imaging — 直接访问 DICOM 标签、匿名化和编辑元数据；可与 SimpleITK 配合使用，后者能够重建体数据，但不负责 PHI 管理
• nnunet-segmentation — 医学图像的自动化深度学习分割；当 SimpleITK 经典分割无法满足复杂解剖结构的需求时使用

参考资料

• SimpleITK 文档 — 完整的 API 参考和 Jupyter 笔记本教程
• SimpleITK GitHub — 源代码、示例和问题跟踪器
• SimpleITK Notebooks — 50 多个涵盖配准、分割和输入/输出的 Jupyter 笔记本
• Lowekamp BC 等人（2013）Front Neuroinform — 介绍 SimpleITK 设计和功能的原创论文
• Yaniv Z 等人（2018）Frontiers in Neuroinformatics — 面向生物医学图像处理领域的 SimpleITK