Simpleitk Image Registration

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

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

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

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

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

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

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

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)

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

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

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