Resampling an Image onto Another's Physical Space¶
The purpose of this Notebook is to demonstrate how the physical space described by the meta-data is used when resampling onto a reference image.
In [1]:
%matplotlib inline
import matplotlib.pyplot as plt
import SimpleITK as sitk
# If the environment variable SIMPLE_ITK_MEMORY_CONSTRAINED_ENVIRONMENT is set, this will override the ReadImage
# function so that it also resamples the image to a smaller size (testing environment is memory constrained).
%run setup_for_testing
print(sitk.Version())
from myshow import myshow
# Download data to work on
%run update_path_to_download_script
from downloaddata import fetch_data as fdata
OUTPUT_DIR = "Output"
SimpleITK Version: 2.4.0 (ITK 5.4) Compiled: Aug 15 2024 01:21:37
Load the RGB cryosectioning of the Visible Human Male dataset. The data is about 1GB so this may take several seconds, or a bit longer if this is the first time the data is downloaded from the midas repository.
In [2]:
fixed = sitk.ReadImage(fdata("vm_head_rgb.mha"))
Fetching vm_head_rgb.mha
In [3]:
moving = sitk.ReadImage(fdata("vm_head_mri.mha"))
Fetching vm_head_mri.mha
In [4]:
print(fixed.GetSize())
print(fixed.GetOrigin())
print(fixed.GetSpacing())
print(fixed.GetDirection())
(2048, 1216, 220) (-334.0, 159.0, -20.0) (0.33, 0.33, 1.0) (1.0, 0.0, 0.0, 0.0, -1.0, 0.0, 0.0, 0.0, -1.0)
In [5]:
print(moving.GetSize())
print(moving.GetOrigin())
print(moving.GetSpacing())
print(moving.GetDirection())
(256, 256, 33) (-130.0, -161.6, -74.9) (1.01562, 1.01562, 5.0) (1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0)
In [6]:
import sys
resample = sitk.ResampleImageFilter()
resample.SetReferenceImage(fixed)
resample.SetInterpolator(sitk.sitkBSpline3)
resample.AddCommand(
sitk.sitkProgressEvent,
lambda: print(f"\rProgress: {100*resample.GetProgress():03.1f}%...", end=""),
)
resample.AddCommand(sitk.sitkProgressEvent, lambda: sys.stdout.flush())
out = resample.Execute(moving)
Progress: 0.0%...
Progress: 49.5%...
Progress: 92.6%...
Progress: 1.8%...
Progress: 5.1%...
Progress: 3.6%...
Progress: 5.5%...
Progress: 7.3%...
Progress: 9.1%...
Progress: 10.9%...
Progress: 12.7%...
Progress: 14.5%...
Progress: 16.4%...
Progress: 18.2%...
Progress: 20.0%...
Progress: 21.8%...
Progress: 23.6%...
Progress: 25.5%...
Progress: 27.3%...
Progress: 29.1%...
Progress: 30.9%...
Progress: 32.7%...
Progress: 34.5%...
Progress: 36.4%...
Progress: 38.2%...
Progress: 40.0%...
Progress: 41.8%...
Progress: 43.6%...
Progress: 45.5%...
Progress: 47.3%...
Progress: 49.1%...
Progress: 50.9%...
Progress: 52.7%...
Progress: 54.5%...
Progress: 56.4%...
Progress: 58.2%...
Progress: 60.0%...
Progress: 61.8%...
Progress: 63.6%...
Progress: 65.5%...
Progress: 67.3%...
Progress: 69.1%...
Progress: 70.9%...
Progress: 72.7%...
Progress: 74.5%...
Progress: 76.4%...
Progress: 78.2%...
Progress: 80.0%...
Progress: 81.8%...
Progress: 83.6%...
Progress: 85.5%...
Progress: 87.3%...
Progress: 89.1%...
Progress: 90.9%...
Progress: 92.7%...
Progress: 94.5%...
Progress: 96.4%...
Progress: 98.2%...
Progress: 100.0%...
Progress: 100.0%...
Because we are resampling the moving image using the physical location of the fixed image without any transformation (identity), most of the resulting volume is empty. The image content appears in slice 57 and below.
In [7]:
myshow(out)
In [8]:
# combine the two images using a checkerboard pattern:
# because the moving image is single channel with a high dynamic range we rescale it to [0,255] and repeat
# the channel 3 times
vis = sitk.CheckerBoard(
fixed,
sitk.Compose([sitk.Cast(sitk.RescaleIntensity(out), sitk.sitkUInt8)] * 3),
checkerPattern=[15, 10, 1],
)
In [9]:
myshow(vis)
Write the image to the Output directory: (1) original as a single image volume and (2) as a series of smaller JPEG images which can be constructed into an animated GIF.
In [10]:
import os
sitk.WriteImage(vis, os.path.join(OUTPUT_DIR, "example_resample_vis.mha"))
temp = sitk.Shrink(vis, [3, 3, 2])
sitk.WriteImage(
temp, [os.path.join(OUTPUT_DIR, f"r{i:03d}.jpg") for i in range(temp.GetSize()[2])]
)