simpleitk-image-registration

SimpleITK による画像レジストレーションと解析

概要

SimpleITK は、医療画像処理向けの Insight Toolkit（ITK）を簡略化した高レベルインターフェースです。レジストレーション（剛体、アフィン、B スプライン、Demons）、セグメンテーション（しきい値処理、領域拡張、ウォーターシェッド、レベルセット）、フィルタリング（平滑化、形態学的処理、勾配）、および MRI、CT、超音波、蛍光顕微鏡による 3D／4D 画像のリサンプリングを、Python ネイティブに利用できます。SimpleITK の画像には物理空間のメタデータ（間隔、原点、方向余弦）が含まれており、解剖学的に正しく解釈し、マルチモーダル画像を位置合わせするうえで重要です。

使用する場面

• 縦断研究で時点間の MRI ボリュームを位置合わせする、または標準アトラスに位置合わせして正規化する場合
• 蛍光顕微鏡画像から、Otsu しきい値処理と形態学的な後処理を用いて細胞や細胞核をセグメント化する場合
• DICOM シリーズ（CT、MRI スキャナー出力）を NIfTI 形式に変換し、FSL や ANTs で後続の解析を行う場合
• 事前計算済みの変換を適用して、画像を共通の解像度または視野にリサンプリングする場合
• 2 値ラベルマスクから領域の統計量（体積、平均輝度、表面積）を計算する場合
• 相互情報量メトリクスを使ってマルチモーダルレジストレーション（例：PET と MRI の位置合わせ）を行う場合
• 神経画像研究で最先端の微分同相レジストレーションとマルチアトラスラベル融合が必要な場合は、代わりに ANTs（antspyx 経由）を使用してください。ネイティブ依存関係なしで、Python ネイティブのスクリプト可能なパイプラインを構築する場合は SimpleITK が適しています
• regionprops、形態学的処理、ウォーターシェッドを使った非ボリューム蛍光顕微鏡データの 2D バイオ画像解析には、代わりに scikit-image（scikit-image-processing）を使用してください

前提条件

• Python パッケージ：SimpleITK>=2.3、numpy、matplotlib
• オプション：追加のレジストレーションアルゴリズム（Elastix）を使用する場合は SimpleITK-SimpleElastix
• データ要件：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：画像入出力

メタデータを完全に保持したまま、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 スプライン、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 スプラインによる調整）を適用できます。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 フィルターは float 入力を想定しています。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}")

レシピ：ラベルマスクへの変換の一括適用

使用する場面：1 つの画像を位置合わせした後、同じ変換をセグメンテーションマスクやアトラスラベルに適用する場合です（最近傍補間を使用します）。

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'	画像のピクセル型が float ではありません	登録の前に float32 にキャストします：sitk.Cast(image, sitk.sitkFloat32)
登録が発散する（メトリックが増加する）	学習率が高すぎる、または画像が事前に位置合わせされていません	learningRate を 0.1 に下げます。初期位置合わせには CenteredTransformInitializer を使用し、画像が重なっていることを確認します
リサンプリング後のラベル画像に小数値が含まれる	マスクに誤った補間器が使用されています	整数ラベル画像をリサンプリングする際は、補間器に sitk.sitkNearestNeighbor を設定します
RuntimeError: Exception thrown in SimpleITK ReadImage	サポートされていない形式、または破損した DICOM	dcmdump でファイルを確認します。DICOM シリーズには ImageSeriesReader.GetGDCMSeriesFileNames() を使用します
N4 バイアス補正が非常に遅い	画像サイズが大きい、または反復回数が多すぎます	まずダウンサンプリングします。補正前に 2mm 等方解像度へリサンプリングし、その後、元の解像度の画像に変換を適用します
連結成分ラベルが多すぎる、または少なすぎる	オブジェクトサイズのしきい値が不適切、または前処理が不十分です	RelabelComponent の minimumObjectSize を調整します。接触しているオブジェクトを結合するには、ガウシアンのシグマを大きくします
NumPy 配列の形状が想定される (z,y,x) と一致しない	軸の順序が混同されています	SimpleITK 画像は (x,y,z) です。GetArrayFromImage() は (z,y,x) の NumPy 配列を返します。変換には arr.transpose(2,1,0) を使用します
大津のしきい値処理で背景が多くセグメント化される	背景が明るい（例：共焦点画像の反射アーティファクト）	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 管理には対応しないため、SimpleITK と併用します
• nnunet-segmentation — 医用画像の自動ディープラーニングセグメンテーション。複雑な解剖学的構造に対して SimpleITK の従来型セグメンテーションでは不十分な場合に使用します

参考資料

• SimpleITK のドキュメント — API リファレンス一式と Jupyter ノートブックのチュートリアル
• SimpleITK の GitHub — ソースコード、サンプル、課題トラッカー
• SimpleITK ノートブック — 位置合わせ、セグメンテーション、I/O を網羅した 50 以上の Jupyter ノートブック
• Lowekamp BC et al. (2013) Front Neuroinform — SimpleITK の設計と機能を説明する原著論文
• Yaniv Z et al. (2018) Frontiers in Neuroinformatics — 生物医学画像処理コミュニティ向けの SimpleITK