We start by setting up the Python environment and pre-define some useful functions for warping and displaying images.
%matplotlib inline
import sys
sys.path.insert(0,'simpleitk/SimpleITK-0.8.0-STRAW-branch-py2.7-linux-x86_64.egg_FILES')
sys.path.insert(0,'../src')
import SimpleITK as sitk
import numpy as np
import math as mt
import matplotlib.pyplot as plt
from imagedisplay import myshow
from IPython.display import display, clear_output
# Define the warping function
def warpmyimage(reference,img_source,transformation):
stats = sitk.StatisticsImageFilter()
stats.Execute( img_source )
minval = stats.GetMinimum()
resampler = sitk.ResampleImageFilter()
resampler.SetReferenceImage(reference)
resampler.SetInterpolator(sitk.sitkLinear)
resampler.SetDefaultPixelValue(minval)
resampler.SetTransform(transformation)
moving = resampler.Execute(img_source)
return moving
def cmyshow(img1,img2,*args, **kw):
simg1 = sitk.Cast(sitk.RescaleIntensity(img1), sitk.sitkUInt8)
simg2 = sitk.Cast(sitk.RescaleIntensity(img2), sitk.sitkUInt8)
cimg = sitk.Compose(simg1, simg2, simg1/2.+simg2/2.)
myshow(cimg,*args, **kw)
Here, we are running a full multi modal CT-MRI registration. The original images are already aligned, so we apply a random transformation to the moving image initially. The registration method is then trying to recover this transformation.
TASK: a) Modify the code below such that the registration 'works', b) Try out differernt initial transformation parameters
# Multi Modal Image Registration
fixedInputFilename = "./data/ct.nii.gz"
movingInputFilename = "./data/mr.nii.gz"
fixed = sitk.ReadImage(fixedInputFilename)
moving = sitk.ReadImage(movingInputFilename)
# Define a transformation and warp the moving image (this time transformation is in 3D)
# Resample the moving image based on this transformation
# CHANGE CODE HERE: Try different initial transformation parameters
angleZ = -3.14*40/180; tx = -35.00; ty = 41.00;
t_input=sitk.Transform(fixed.GetDimension(), sitk.sitkEuler)
t_input.SetParameters([0,0,angleZ,tx,ty,0])
moving = warpmyimage(fixed,moving,t_input)
global t_output
t_output = sitk.Transform(fixed.GetDimension(), sitk.sitkEuler)
# Multi-scale registration
levels = [8,4]
for level in levels:
fixed_s = sitk.DiscreteGaussian(fixed, 2.0 * level )
moving_s = sitk.DiscreteGaussian(moving, 2.0 * level )
fixed_s = sitk.Shrink(fixed_s, [1,1] * level )
moving_s = sitk.Shrink(moving_s, [1,1] * level )
# Create the registration object and set its properties
R = sitk.ImageRegistrationMethod()
R.SetOptimizerAsRegularStepGradientDescent(maxStep=4.0, minStep=.05, numberOfIterations=400 )
R.SetOptimizerScales([1e10,1e10,100,1,1,1e10])
R.SetInterpolator(sitk.sitkLinear)
R.SetTransform(t_output)
R.SetNumberOfThreads(8)
# Define the (dis)similarity measure, sum of squared differences
R.SetMetricAsMeanSquares()
# To observe the cost function over iterations
def command_iteration(method) :
if ((method.GetOptimizerIteration() % 100) == 0):
print("{0} = {1} : {2}".format(method.GetOptimizerIteration(),
method.GetMetricValue(),
method.GetOptimizerPosition()));
sys.stdout.flush();
def command_end(R):
print("-------")
print("Optimizer stop condition: {0}".format(R.GetOptimizerStopConditionDescription()))
print(" Iteration: {0}".format(R.GetOptimizerIteration()))
print(" Metric value: {0}".format(R.GetMetricValue()))
sys.stdout.flush()
R.RemoveAllCommands()
R.AddCommand( sitk.sitkIterationEvent, lambda: command_iteration(R) )
R.AddCommand( sitk.sitkEndEvent, lambda: command_end(R) )
# Run the registration and warp moving image
t_output = R.Execute(fixed_s,moving_s)
warped_image = warpmyimage(fixed,moving,t_output)
cmyshow(fixed,moving, title="Fixed Red, Moving Green - Before Registration", dpi = 50)
cmyshow(fixed,warped_image, title="Fixed Red, Moving Green - After Registration", dpi = 50)
0 = 994950.373688 : (1.5524825498136799e-18, -1.0743783672051473e-18, 0.0344688981185441, -1.9432108572259512, -0.5855614308451771, 5.3206890344523745e-21) 100 = 728892.181617 : (6.074729805325386e-17, -6.204852149253858e-17, 0.4203838018875175, 2.742206722789201, -8.10792654652087, 3.30865430113541e-19) 200 = 730760.538218 : (9.120673405223583e-17, -8.958944412607805e-17, 0.478735648675651, 15.6966629305593, -12.946335662107089, 4.901385976296533e-19) 300 = 731940.929968 : (1.1420128189407496e-16, -1.0719310294088356e-16, 0.5265456126533368, 26.54655406557398, -16.607703625342882, 5.947013028762683e-19) ------- Optimizer stop condition: Iteration: 400 Metric value: 733250.191871 0 = 751127.259849 : (1.3900337589558087e-16, -1.2200220426175787e-16, 0.5422046347155325, 39.79414426408663, -19.46251874258399, 6.858996938313943e-19) 100 = 752595.694069 : (1.524231146366928e-16, -1.260152078788725e-16, 0.6005596336967433, 44.105999548661025, -21.484659199013276, 7.148927849127011e-19) 200 = 753538.489712 : (1.6831201208302585e-16, -1.3047401816640664e-16, 0.6239653725759613, 49.53917556153984, -22.884842898724724, 7.442845823338585e-19) ------- Optimizer stop condition: Iteration: 266 Metric value: 754152.278213