First, we initialise the Python environment to make use of the development version of SimpleITK.
%matplotlib inline
import sys
# first you need to run:
# cd simpleitk && tar -zxvf SimpleITK-0.8.0-STRAW-branch-py2.7-linux-x86_64.egg_FILES.tar.gz
sys.path.insert(0,'simpleitk/SimpleITK-0.8.0-STRAW-branch-py2.7-linux-x86_64.egg_FILES')
sys.path.insert(0,'../src')
import numpy as np
import math as mt
import SimpleITK as sitk
import matplotlib.pyplot as plt
from imagedisplay import myshow
from IPython.display import display, clear_output
The next cell is a simple registration experiment for comparing different (dis)similarity measures. Here, we use the same 2D brain image for both the fixed and the moving image. We translate the moving image along the x-axis and record and plot the (dis)similarity values.
TASK: a) Try out different (dis)similarity measures, b) try out different interpolation methods, c) change the range of translations
#load images
fixedInputFilename = "./data/BrainProtonDensitySliceBorder20.png"
fixedInput = sitk.ReadImage(fixedInputFilename)
#set fixed and moving to the same image
fixed = sitk.VectorIndexSelectionCast(fixedInput,0,sitk.sitkFloat32)
moving = fixed
# the following defines a function that takes two images, an array with 1D translations,
# and a flag defining the (dis)similarity measure
def metric_calc_translation (fixed_img, moving_img, tx_array, metric='SSD'):
# Define the registration object and properties
# (just to calculate the metric value - we are not going to perform optimization)
R = sitk.ImageRegistrationMethod()
R.SetOptimizerAsRegularStepGradientDescent( maxStep=1e-10,
minStep=1e-11,
numberOfIterations=1 )
# Interpolation method
# MAKE A SELECTION HERE
# Linear interpolation
R.SetInterpolator(sitk.sitkLinear)
# Nearest neighbor interpolation
#R.SetInterpolator(sitk.sitkNearestNeighbor)
if metric == 'SSD':
R.SetMetricAsMeanSquares()
title_str = 'Sum of Squared Differences'
elif metric == 'CC':
R.SetMetricAsNormalizedCorrelation(subtractMean=True)
title_str = 'Correlation Coefficient'
elif metric == 'NMI':
R.SetMetricAsMattesMutualInformation(numberOfHistogramBins=50)
title_str = 'Normalized Mutual Information'
else:
raise RuntimeError('undefined similarity metric')
cost = []
for tx in tx_array:
# Create a transformation (translation in x-direction)
t_input=sitk.Transform(fixed.GetDimension(), sitk.sitkTranslation)
t_input.SetParameters( [tx,0] ) # t_input[0]=offset_x t_input[1]=offset_y
# Run the cost function calculation
R.SetTransform(t_input)
R.Execute(fixed,moving)
# Collect value
cost.append(R.GetMetricValue())
# Plot the metric values
fig = plt.figure()
ax = fig.add_subplot(111)
ax.set_title(title_str,fontsize=16)
ax.set_xlabel('Translation in x direction (in mm)',fontsize=16)
ax.set_ylabel('Cost function',fontsize=16)
ax.plot(tx_arr,cost,'r')
# CHANGE CODE HERE: Try different a different translation range
tx_arr = np.linspace(-15.0, 15.0, num=100, endpoint=True)
# Run the experiment
metric_calc_translation (fixed,moving,tx_arr,'SSD')
metric_calc_translation (fixed,moving,tx_arr,'CC')
metric_calc_translation (fixed,moving,tx_arr,'NMI')
Now we explore the effect of Gaussian smoothing on the (dis)similarity values.
TASK: a) try out different filter kernel sizes, b) try out differernt (dis)similarity measures
tx_arr = np.linspace(-15.0, 15.0, num=100, endpoint=True)
# Apply Gaussian smoothing to the fixed and the moving image
moving = sitk.DiscreteGaussian(moving, 2.0)
fixed = sitk.DiscreteGaussian(fixed, 2.0)
# Re-run the little translation experiments
metric_calc_translation (fixed,moving,tx_arr,metric='SSD')
metric_calc_translation (fixed,moving,tx_arr,metric='CC')
metric_calc_translation (fixed,moving,tx_arr,metric='NMI')
In the next cell, we repeat the earlier experiment, but this time the moving image is also rotated.
TASK: Explore how different rotations impact the (dis)similarity values.
# CHANGE CODE HERE: Try different values for the rotation
angle = 3.14/4
# Define display function
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)
# Define the warping function
def warpmyimage(reference,img_source,transformation):
resampler = sitk.ResampleImageFilter()
resampler.SetReferenceImage(reference)
resampler.SetInterpolator(sitk.sitkLinear)
resampler.SetDefaultPixelValue(0.0)
resampler.SetTransform(transformation)
moving = resampler.Execute(img_source)
return moving
# Read the image again
fixedInputFilename = "./data/BrainProtonDensitySliceBorder20.png"
fixedInput = sitk.ReadImage(fixedInputFilename)
fixed = sitk.VectorIndexSelectionCast(fixedInput,0,sitk.sitkFloat32)
# Smooth the input image
fixed = sitk.DiscreteGaussian(fixed, 1.0)
# Move the image center to origin - the rotation will be performed around the point (0,0)
spacing = fixed.GetSpacing();dimension = fixed.GetSize();
fixed.SetOrigin([ -1*spacing[0]*dimension[0]/2,
-1*spacing[1]*dimension[1]/2 ])
# Transformation object for rotation:
# http://www.itk.org/Doxygen/html/classitk_1_1Euler2DTransform.html
# Get the parameters that uniquely define the transform.
# These are typically used by optimizers. There are 3 parameters.
# The first one represents the angle or rotation in radians and the last two represents the translation.
# The center of rotation is fixed.
t_input=sitk.Transform(fixed.GetDimension(), sitk.sitkEuler)
p = list(t_input.GetParameters()); p[0] = angle; t_input.SetParameters(p);
# Resample the moving image based on this transformation
moving = warpmyimage(fixed,fixed,t_input)
cmyshow(fixed,moving, title="Fixed Red, Moving Green", dpi = 50)
# Run the (dis)similarity evaluation for a range of translations
tx_arr = np.linspace(-15.0, 15.0, num=100, endpoint=True)
metric_calc_translation (fixed,moving,tx_arr,metric='SSD')
metric_calc_translation (fixed,moving,tx_arr,metric='CC')
metric_calc_translation (fixed,moving,tx_arr,metric='NMI')
The parameter space has 3 dimensions (1 rotation & 2 translations). To obtain the optimal transformation, it is possible to perform a greedy search which is computationally very expensive. Rather than doing this, in this cell we will try to optimize the cost function with gradient descent algorithm.
TASK: Try out different (dis)similarity measures
# Setting up the registration
R = sitk.ImageRegistrationMethod()
# (Dis)similarity measure
# MAKE A SELECTION HERE
# SSD
R.SetMetricAsMeanSquares()
# CC
#R.SetMetricAsNormalizedCorrelation(subtractMean=True)
# NMI
#R.SetMetricAsMattesMutualInformation(numberOfHistogramBins=50)
R.SetOptimizerAsRegularStepGradientDescent(maxStep=2.0, minStep=.01, numberOfIterations=100 )
R.SetOptimizerScales([100.0, 1.0, 1.0])
R.SetTransform(sitk.Transform(fixed.GetDimension(), sitk.sitkEuler))
R.SetInterpolator(sitk.sitkLinear)
def command_iteration(method) :
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) )
t_output=R.Execute(fixed,moving)
warped_image = warpmyimage(fixed,moving,t_output)
cmyshow(fixed,moving,title="Before registration, Fixed Red, Moving Green", dpi = 50)
cmyshow(fixed,warped_image, title="After registration, Fixed Red, Moving Green", dpi = 50)
0 = 2733.91151295 : (-0.011247367522158776, 0.22291662197915948, 1.6386809994200457) 1 = 2673.6432611 : (-0.02784494462369627, -0.30787160833537575, 2.620245158374898) 2 = 2652.23799383 : (-0.04466714487640808, -1.3696719221827545, 2.8269233625722343) 3 = 2623.32745972 : (-0.05989246324173654, -2.645497419757683, 2.5941872319830974) 4 = 2588.17210585 : (-0.07684988084701647, -3.683785835714133, 2.8096333897396724) 5 = 2556.08101082 : (-0.09504456163161594, -4.323046097155317, 3.339616624953215) 6 = 2526.46169992 : (-0.1148772751382055, -4.389822245673656, 3.090266077240859) 7 = 2496.9389636 : (-0.13419886557392374, -4.458262196789887, 3.602201037401648) 8 = 2464.65479018 : (-0.15331236265558748, -4.546484858186561, 3.0199991350443973) 9 = 2436.46200082 : (-0.17126446269688633, -4.5364819198385264, 3.901543792050715) 10 = 2403.93766927 : (-0.1870942024000613, -4.3634535138929555, 2.6914799176898616) 11 = 2384.52978892 : (-0.20181275596115084, -4.33542982752636, 4.045310218612328) 12 = 2350.02867525 : (-0.20850800321595414, -4.251558968163576, 3.3072707068514484) 13 = 2330.62629003 : (-0.21844299808556833, -4.264698900337555, 3.4203463965016256) 14 = 2310.3349304 : (-0.22843239351339398, -4.227537563176255, 3.393165063095336) 15 = 2289.87683952 : (-0.2383842478865777, -4.1370491939369, 3.3555113729462813) 16 = 2269.51592567 : (-0.24836585564259692, -4.086533688106959, 3.3219964550661585) 17 = 2249.09147532 : (-0.2583084484249927, -4.008663737959035, 3.2486155413454734) 18 = 2229.37078179 : (-0.26829529453501005, -3.9610643689473632, 3.2295538567924615) 19 = 2209.19439241 : (-0.2782677196923177, -3.8929418225677197, 3.2001143681127076) 20 = 2189.46230801 : (-0.2881538173101924, -3.769155143959688, 3.11451298320735) 21 = 2170.00917435 : (-0.2981290112416358, -3.698845357018357, 3.1111079632060967) 22 = 2149.98125238 : (-0.3079723146048562, -3.566137451803283, 2.9949942220314667) 23 = 2130.40952175 : (-0.31792789833086005, -3.4851645765817576, 2.946963812667945) 24 = 2109.76094867 : (-0.3278631549188619, -3.3722876851184913, 2.934096307888813) 25 = 2088.31996755 : (-0.33766184703574376, -3.2352616244465624, 2.7889061336954253) 26 = 2066.89553519 : (-0.34752773391546, -3.0750357759048246, 2.757754276198386) 27 = 2043.60953186 : (-0.3573952116535142, -2.934048707003279, 2.677433187481891) 28 = 2019.00325961 : (-0.3673216977086295, -2.826318124037057, 2.6222709950874084) 29 = 1992.55887746 : (-0.3771563060968801, -2.667347135492692, 2.5354781192350266) 30 = 1965.02069225 : (-0.38703985266319013, -2.5327261362621836, 2.4645397002929834) 31 = 1935.40093104 : (-0.3968423700968658, -2.382753779809364, 2.335641244426209) 32 = 1904.82596132 : (-0.40670121674401577, -2.2191001027548256, 2.3003014387345857) 33 = 1872.62276502 : (-0.41649978079928673, -2.059576681010913, 2.1801610938764395) 34 = 1839.64742593 : (-0.4264058329026601, -1.9272607421342227, 2.145609544392996) 35 = 1804.8557676 : (-0.4362064344212007, -1.7655286764830327, 2.030176283640344) 36 = 1769.92426886 : (-0.44612910282607476, -1.663705665037161, 1.9591932461649593) 37 = 1733.99448463 : (-0.4560220743436341, -1.5288183468232672, 1.9035465309640076) 38 = 1697.55116878 : (-0.4658310752044147, -1.3910085607737581, 1.7662749607061927) 39 = 1660.85821921 : (-0.4757211073495198, -1.2860153777051473, 1.6621160763179337) 40 = 1622.18841768 : (-0.48563709293688767, -1.1781403266659447, 1.590734904609866) 41 = 1582.44454803 : (-0.4955022293969491, -1.0508120984559948, 1.487883805711614) 42 = 1541.18365146 : (-0.5053789790705908, -0.9309321742679558, 1.3872509948652127) 43 = 1498.35098873 : (-0.5152353941588856, -0.7961659400785669, 1.2855215586396926) 44 = 1453.44383635 : (-0.5251310100787123, -0.6995874428056386, 1.1785615215359295) 45 = 1406.06850552 : (-0.5349678644541209, -0.5733175704644434, 1.0504253466711078) 46 = 1356.91194027 : (-0.5448367763786498, -0.4418495371102757, 0.9568201627072082) 47 = 1305.10814411 : (-0.5547469635726093, -0.37906244552805074, 0.8387537452234752) 48 = 1251.2948536 : (-0.5646491869036635, -0.2848788607485886, 0.7358499873846819) 49 = 1195.96976419 : (-0.5745756757352789, -0.2088950170224761, 0.6416446252721015) 50 = 1139.73219132 : (-0.5845331576031029, -0.1645653418322926, 0.5608954555110848) 51 = 1082.94302726 : (-0.5945107575208185, -0.1351196392235611, 0.5008292468918719) 52 = 1025.78447594 : (-0.6044860042233053, -0.1010876130964182, 0.4392958946108472) 53 = 968.681601664 : (-0.6144636217850268, -0.07312977445198171, 0.3785518646275332) 54 = 911.035839854 : (-0.6244527386913546, -0.054523363924487894, 0.33578229829326606) 55 = 852.631707895 : (-0.6344406681429724, -0.03643929975354372, 0.29011371699131716) 56 = 793.630232965 : (-0.6444266438230178, -0.010442538610934058, 0.24399359572637302) 57 = 733.836968468 : (-0.654394570611073, 0.021911957803475674, 0.17079836335776322) 58 = 673.468070426 : (-0.6643888672306525, 0.006505562895288773, 0.2008482185080155) 59 = 612.225153486 : (-0.6743470622766472, 0.04828648776499008, 0.11962095352698242) 60 = 551.330497439 : (-0.6843448721437366, 0.03798640568664381, 0.13783871750785598) 61 = 489.859732602 : (-0.6943038865209173, 0.08692592473552022, 0.06177788773099778) 62 = 428.781644057 : (-0.7042768236831439, 0.035787521574492935, 0.11459935405807775) 63 = 367.713667374 : (-0.7141334030579917, 0.13200782294202537, -0.024037163057761435) 64 = 308.386105125 : (-0.7237961228183785, -0.03588947217526142, 0.17122923603487963) 65 = 252.777899844 : (-0.7325297167077376, 0.26099717781065224, -0.2149076500788816) 66 = 207.807519403 : (-0.7394051396309214, -0.20333126911964927, 0.3433796322810868) 67 = 179.932051662 : (-0.7419669340435101, 0.06601042155871806, 0.008973969143019789) 68 = 146.386820523 : (-0.7469662222225472, 0.05759946290641282, 0.009631847511894505) 69 = 120.285156615 : (-0.7519653396381675, 0.06025472069317848, 0.0006207291291754069) 70 = 95.6918780597 : (-0.7569629614763383, 0.04556519399561809, 0.005308738490941686) 71 = 73.0362684528 : (-0.7619607057634356, 0.04929689106362953, -0.00923756154949932) 72 = 52.7370444187 : (-0.7669412961190282, 0.016210034425622702, 0.01978754853808919) 73 = 35.4709789323 : (-0.7718168973343232, 0.07824074592526409, -0.0720673300747605) 74 = 22.6879328517 : (-0.7753517387228692, -0.1450907293203528, 0.20210782515057268) 75 = 23.2328936626 : (-0.7759892127609961, 0.006274845293913256, 0.013627769630537911) 76 = 13.0345728629 : (-0.7783906821098298, 0.04426345162531626, -0.04456428700198432) 77 = 9.67233493697 : (-0.7799866533677959, -0.0739827863058087, 0.10724672924056293) 78 = 9.99725948149 : (-0.7803086684159862, 0.00021246315255557535, 0.011941283198896036) 79 = 7.18628005511 : (-0.7814379676844662, 0.029694379143778173, -0.032809674909217164) 80 = 6.42269731493 : (-0.7820898291550449, -0.035345406931062584, 0.0517219313637498) 81 = 6.39856713244 : (-0.7822780169783238, 0.0014980186149175165, 0.0048746750328473795) 82 = 5.70534975091 : (-0.7828643071716243, 0.012628266471620924, -0.013699653354601198) 83 = 5.44648579424 : (-0.7832831007840065, -0.014046169376744072, 0.024258792009373916) 84 = 5.37704910385 : (-0.7834010104278911, 0.002475384823260483, 0.0004980280147887269) 85 = 5.21587470203 : (-0.7837109740435225, 0.0034576002169779264, -0.003352078836648201) 86 = 5.12850074197 : (-0.7840098708907106, -0.0011163765082649953, 0.0045376603432755866) 87 = 5.06607748441 : (-0.7842479600370101, 0.008128888671513528, -0.01346853736470751) 88 = 5.04832890833 : (-0.7843005478365795, 0.0005128554329254867, -0.0008795780720095327) 89 = 5.01301803195 : (-0.7844557158449231, 0.00232435599528142, -0.0011761966477929257) 90 = 4.99409854278 : (-0.7846052707956724, -0.0018507826840805227, 0.0005679389753134019) 91 = 4.98071711758 : (-0.7847189126375667, 0.007328820903941795, -0.004975599666719616) ------- Optimizer stop condition: Iteration: 93 Metric value: 4.98124618928
In the next two cells, we are plotting the progress of the (dis)similarity measure during optimisation with two different methods. The first one is using Gradient Descent, the second one is using the Nealder-Mead/Downhill Simplex optimiser (http://en.wikipedia.org/wiki/Nelder%E2%80%93Mead_method).
TASK: Try to explain the two differernt curves
def command_plot_value_start(method,title):
global x, ax, f
x = []
f, ax = plt.subplots()
ax.set_title(title,fontsize=16)
ax.set_xlabel('Optimization Iteration Number',fontsize=16)
ax.set_ylabel('Metric Value',fontsize=16)
def command_plot_value_end(method):
global x, ax
del x
del ax
plt.close()
def command_plot_value(method):
global x, ax, f
v = method.GetMetricValue()
x.append(v)
if( len(ax.lines) > 0 ):
ax.lines.pop(0)
ax.plot(x,'r')
display(f)
titleSSD='SSD Registration -GradientDescent Optimizer'
R.RemoveAllCommands()
R.AddCommand( sitk.sitkIterationEvent, lambda: clear_output() )
R.AddCommand( sitk.sitkStartEvent, lambda: command_plot_value_start(R,titleSSD) )
R.AddCommand( sitk.sitkEndEvent, lambda: command_plot_value_end(R) )
R.AddCommand( sitk.sitkIterationEvent, lambda: command_plot_value(R) )
R.SetTransform(sitk.Transform(fixed.GetDimension(), sitk.sitkEuler))
R.Execute(fixed,moving);
titleNelderMead='SSD Registration - Nelder-Mead Simplex'
R.RemoveAllCommands()
R.AddCommand( sitk.sitkIterationEvent, lambda: clear_output() )
R.AddCommand( sitk.sitkStartEvent, lambda: command_plot_value_start(R,titleNelderMead) )
R.AddCommand( sitk.sitkIterationEvent, lambda: command_plot_value(R) )
R.SetTransform(sitk.Transform(fixed.GetDimension(), sitk.sitkEuler))
R.SetOptimizerAsAmoeba(simplexDelta=3,parametersConvergenceTolerance=0.01,functionConvergenceTolerance=0.01);
R.Execute(fixed,moving);