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]:
from __future__ import print_function
%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: 1.2.4 (ITK 4.13)
Compiled: Nov 12 2019 23:05:07

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.sitkBSpline)
resample.AddCommand(sitk.sitkProgressEvent, lambda: print("\rProgress: {0:03.1f}%...".format(100*resample.GetProgress()),end=''))
resample.AddCommand(sitk.sitkProgressEvent, lambda: sys.stdout.flush())
out = resample.Execute(moving)

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,"r{0:03d}.jpg".format(i)) for i in range(temp.GetSize()[2])])
In [ ]: