ITK Registration Methods

1
ITK Registration Methods
Kitware Inc.
2
Overview

• Image Resampling
• Registration Framework
• Multi-Modality
• Multi-Resolution
• Deformable registration

3
Image Resampling
4
Why Resampling ?
Resamplingis the Essence of Intensity-BasedImag
e Registration
5
What is an Image ?
An Image is a sampling of a continuous field
using a discrete grid
6
Image Origin Spacing
Spacing (Sx)
Spacing (Sy)
Origin (Ox,Oy)
7
Image Sampling Grid
Spacing (Sx)
Spacing (Sy)
Origin (Ox,Oy)
8
Image Pixel
Spacing (Sx)
Pixel Value
Pixel Region
Spacing (Sy)
Origin (Ox,Oy)
9
Image Indices
Spacing (Sx)
Pixel Index
0,7
4,7
0,6
0,5
0,4
0,3
0,2
Spacing (Sy)
0,1
0,0
1,0
2,0
3,0
4,0
5,0
Origin (Ox,Oy)
10
Index to Physical Coordinates
Spacing (Sx)
Pixel Index
0,7
4,7
0,6
0,5
P0 Index0 x Spacing0 Origin0
P1 Index1 x Spacing1 Origin1
0,4
Index0 floor( ( P0 - Origin0 ) /
Spacing0 0.5 )
0,3
Index1 floor( ( P1 - Origin1 ) /
Spacing1 0.5 )
0,2
Spacing (Sy)
0,1
0,0
1,0
2,0
3,0
4,0
5,0
Origin (Ox,Oy)
11
Image Region
Spacing (Sx)
Pixel Value
Image Region
Pixel Region
Spacing (Sy)
Origin (Ox,Oy)
12
Image Region
Spacing (Sx)
Region Size
Image Region
3,5
Starting Index
2,3
Spacing (Sy)
Origin (Ox,Oy)
13
Basic Resampling
ResamplingTrivial Cases
14
Sub-Sampling by Half
Spacing (Sx)
Image Region
Spacing (Sy)
Origin (Ox,Oy)
15
Sub-Sampling by Half
New Spacing Sx
New Spacing Sy
New Origin (Ox,Oy)
Origin (Ox,Oy)
16
Super-Sampling by Double
Spacing ( Sx/2 )
Spacing (Sx)
Spacing ( Sy/2 )
Image Region
Spacing (Sy)
Origin (Ox,Oy)
17
Super-Sampling by Double
New Spacing Sx
New Origin (Ox,Oy)
Origin (Ox,Oy)
18
Resampling in ITK
itkResampleImageFilter
19
Resampling in ITK
Origin
Spacing
Region Size
Resample Filter
Region Start
Transform
Interpolator
20
Resample Image Filter
include "itkImage.h" include "itkResampleImageFi
lter.h" include "itkIdentityTransform.h include
"itkLinearInterpolateImageFunction.h" typedef
itkImagelt char, 2 gt ImageType ImageTypePoin
ter inputImage GetImageSomeHow() typedef
itkResampleImageFilterlt ImageType gt
FilterType FilterTypePointer resampler
FilterTypeNew() ImageTypeSizeType
size size0 200 size1 300 ImageTypeIn
dexType start start0 0 start1 0
21
Resample Image Filter
ImageTypePointType origin origin0 10.0
// millimeters origin1 25.5 //
millimeters ImageTypeSpacingType spacing
spacing0 2.0 // millimeters spacing1
1.5 // millimeters resampler-gtSetOutputSpacing(
spacing ) resampler-gtSetOutputOrigin( origin
) resampler-gtSetSize( size ) resampler-gtSetOutp
utStartIndex( start ) resampler-gtSetDefaultPixel
Value( 100 ) resampler-gtSetInput( inputImage )
22
Resample Image Filter
typedef itkLinearInterpolateImageFunctionlt
ImageType,
double gt
InterpolatorType InterpolatorTypePointer
interpolator InterpolatorTypeNew() typedef
itkTranslationTransformlt double, 2 gt
TransformType TransformTypePointer transform
TransformTypeNew() transform-gtSetIdentity()
resampler-gtSetInterpolator( interpolator
) resampler-gtSetTransform( transform
) resampler-gtUpdate() const ImageType
outputImage resampler-gtGetOutput()
23
Basic Registration
24
Coordinate System Conversions
25
Things I will not do
26
Fixed Image Moving Image
27
Selecting Moving Fixed Images
In principle the denomination of Fixed Image
Moving Image is arbitrary
In practice the moving image is the one that will
be resampled into the fixed image coordinate
system
28
Quiz 1
Images from the same patient
Moving Image ? Fixed Image ?
Images provided as part of the project
Retrospective Image Registration Evaluation,
NIH, Project No. 8R01EB002124-03, Principal
Investigator, J. Michael Fitzpatrick, Vanderbilt
University, Nashville, TN.
29
Quiz 2
Images from the same patient
What scale factor ?

1. 2.0
2. 1.0
3. 0.5

Images provided as part of the project
Retrospective Image Registration Evaluation,
NIH, Project No. 8R01EB002124-03, Principal
Investigator, J. Michael Fitzpatrick, Vanderbilt
University, Nashville, TN.
30
Things I will not do
31
Registration in ITK
MultiResolutionRegistrationFramework
ImageRegistrationFramework
PDEBasedRegistration
FEMBasedRegistration
32
Registration Framework
33
Components
Registration Method
FixedImage
Metric
Optimizer
Interpolator
MovingImage
Transform
34
Image Metrics

• Mean Squares
• Normalized Correlation
• Mean Reciprocal Square Difference
• Mutual Information- Viola-Wells- Mattes-
  Histogram based- Histogram normalized

35
Transforms

• Translation
• Scaling
• Rotation
• Rigid3D
• Rigid2D
• Affine
• BSplines
• Splines TPS, EBS, VS

36
Optimizers

• Gradient Descent
• Regular Step Gradient Descent
• Conjugate Gradient
• Levenberg-Marquardt
• One plus One Evolutionary Algorithm

37
Interpolators

• Nearest Neighbor
• Linear
• BSpline

38
Exercise 20
39
Image Registration
40
Image Registration
include itkImageRegistrationMethod.h include
itkTranslationTransform.h include
itkMeanSquaresImageToImageMetric.h include
itkLinearInterpolateImageFunction.h include
itkRegularStepGradientDescentOptimizer.h includ
e itkImage.h include itkImageFileReader.h in
clude itkImageFileWriter.h include
itkResampleImageFilter.h
41
Image Registration
const unsigned int Dimension 2 typedef
unsigned char PixelType typedef itkImagelt
PixelType , Dimension gt
FixedImageType typedef itkImagelt PixelType ,
Dimension gt
MovingImageType typedef itkTranslationTransfo
rmlt double, Dimension gt TransformType typed
ef itkRegularStepGradientDescentOptimizer
OptimizerType typedef
itkLinearInterpolateImageFunctionlt
MovingImageType ,
double gt InterpolatorType typedef
itkMeanSquaresImageToImageMetriclt
FixedImageType , MovingImageType gt
MetricType typedef itkImageRegistratio
nMethodlt FixedImageType
, MovingImageType gt RegistrationType
42
Image Registration
TransformTypePointer transform
TransformTypeNew() OptimizerTypePointer
optimizer OptimizerTypeNew() Interpolato
rTypePointer interpolator
InterpolatorTypeNew() MetricTypePointer
metric MetricTypeNew() Regist
rationTypePointer registrator
RegistrationTypeNew() registrator-gtSetTransfor
m( transform ) registrator-gtSetOptimizer(
optimizer ) registrator-gtSetInterpolator(
interpolator ) registrator-gtSetMetric(
metric ) registrator-gtSetFixedImage(
fixedImageReader-gtGetOutput() ) registrator-gtSetM
ovingImage( movingImageReader-gtGetOutput() )
43
Image Registration
registrator-gtSetFixedImageRegion(
fixedImageReader-gtGetOutput()-gtGetBufferedR
egion() ) typedef RegistrationTypeParameters
Type ParametersType transform-gtSetIdentity()
registrator-gtSetInitialTransformParameters(

transform-gtGetParameters()
) optimizer-gtSetMaximumStepLength( 4.00 )
optimizer-gtSetMinimumStepLength( 0.01
) optimizer-gtSetNumberOfIterations( 100
) optimizer-gtMaximizeOff()
44
Image Registration
try registrator-gtStartRegistrati
on () catch( itkExceptionObject excp
) stdcerr ltlt Error in
registration ltlt stdendl stdcerr
ltlt excp ltlt stdendl transform-gtSetParam
eters(
registrator-gtGetLastTransformParameters() )
45
Image Registration
typedef itkResampleImageFilterlt
FixedImageType , MovingImageType gt
ResamplerType ResamplerType Pointer
resampler ResamplerTypeNew() resampler-gtSet
Transform ( transform ) resampler-gtSetInput(
movingImageReader-gtGetOutput() ) FixedImageType
Pointer fixedImage fixedImageReader-gtGetOut
put() resampler-gtSetOrigin( fixedImage-gtGetOrigin
() ) resampler-gtSetSpacing( fixedImage-gtGetSpacin
g() ) resampler-gtSetSize(
fixedImage-gtGetLargestPossibleRegion()-gtGetSize
() ) resampler-gtUpdate()
46
Tracking theRegistration Process
47
Observing Registration
include itkCommand.h class CommandIteration
public itkCommand public typedef
CommandIteration Self typedef itkCommand
SuperClass typedef itkSmartPointerlt
Self gt Pointer itkNewMacro( Self
) protected CommandIteration()
public typedef itkRegularStepGradientD
escentOptimizer OptimizerType typedef
const OptimizerType
OptimizerPointer
48
Observing Registration

• void Execute( itkObject caller, const
  itkEventObject event ) this-gt
  Execute( (const itkObject ) caller, event
  )
• void Execute( const itkObject caller, const
  itkEventObject event )
  OptimizerPointer optimizer
  dynamic_castlt
  OptimizerPointer gt( caller )
• if( typeid( event ) typeid(
  itkIterationEvent ) )
• stdcout ltlt optimizer-gtGetCurrentIterat
  ion() ltlt
• stdcout ltlt optimizer-gtGetValue() ltlt
  
• stdcout ltlt optimizer-gtGetCurrentPositi
  on() ltlt stdendl

49
Image Registration

• CommandIterationPointer observer
  CommandIterationNew()
• optimizer-gtAddObserver( itkIterationEvent(),
  observer )
• try registrator-gtStartRegistrati
  on ()
• catch( itkExceptionObject excp )
• stdcerr ltlt Error in registration
  ltlt stdendl stdcerr ltlt excp ltlt
  stdendl

50
Image Similarity Metrics
51
Image Metrics
How similar is image A to image B ?
52
Image Metrics
Does Image B matches Image A better than
Image C ?
53
Image Metrics
gtlt
Match( A , B )
Match( A , C )
Image A
Image C
Image B
54
Image Metrics
Match( A , B )Simplest MetricMean Squared
Differences
55
Mean Squared Differences
For each pixel in A
Image A
Image B
Difference( index ) A( index ) B( index )
Sum Difference( index ) 2
Match( A , B ) Sum / numberOfPixels
56
For each pixel in the Fixed Image
57
Image Metrics
FixedImage
Metric
Value
MovingImage
Interpolator
Transform
Parameters
58
Mean Squared Differences
include "itkImage.h" include "itkMeanSquaresImag
eToImageMetric.h" include "itkLinearInterpolateIm
ageFunction.h" include "itkTranslationTransform.h
" typedef itkImagelt char, 2 gt
ImageType ImageTypeConstPointer fixedImage
GetFixedImage() ImageTypeConstPointer
movingImage GetMovingImage() typedef
itkLinearInterpolateImageFunctionlt
ImageType,
double gt
InterpolatorType InterpolatorTy
pePointer interpolator InterpolatorTypeNew()
typedef itkTranslationTransformlt double, 2
gt TransformType TransformTypePointer
transform TransformTypeNew()
59
Mean Squared Differences
typedef itkMeanSquaresImageToImageMetriclt
ImageType, ImageType gt
MetricType MetricTypePointer metric
MetricTypeNew() metric-gtSetInterpolator(
interpolator ) metric-gtSetTransform( transform
) metric-gtSetFixedImage( fixedImage
) metric-gtSetMovingImage( movingImage
) MetricTypeTransformParametersType
translation( Dimension ) translation0
12 translation1 27 double value
metric-gtGetValue( translation )
60
Mean Squared Differences
MetricTypeTransformParametersType translation(
Dimension ) double value2121 for( int dx
0 dx lt 20 dx) for( int dy 0 dy lt
20 dy) translation0
dx translation1 dy valuedxdy
metric-gtGetValue( translation )
61
Evaluating many matches
y
y
Transform
x
x
Fixed Image
Moving Image
62
Plotting the Metric
Mean Squared Differences
Transform Parametric Space
63
Plotting the Metric
Mean Squared Differences
Transform Parametric Space
64
Evaluating many matches
y
y
Transform
(-15,-25) mm
x
x
Fixed Image
Moving Image
65
Plotting the Metric
Mean Squared Differences
Transform Parametric Space
66
Plotting the Metric
Mean Squared Differences
Transform Parametric Space
67
The Best Transform Parameters
Evaluation of thefull parameter spaceis
equivalent to performingoptimization by
exhaustive search
68
The Best Transform Parameters
Very SafebutVery Slow
69
The Best Transform Parameters
Better Optimization Methodsfor
exampleGradient Descent
70
Image Registration Framework
FixedImage
Metric
MovingImage
Interpolator
Optimizer
Transform
Parameters
71
Gradient Descent Optimizer
f( x , y )
S Step
L LearningRate
S L G( x , y )
72
Gradient Descent Optimizer
f( x , y )
S L G( x , y )
73
Gradient Descent Optimizer
f( x , y )
L too large
S L G( x , y )
74
Gradient Descent Optimizer
f( x , y )
L too small
S L G( x , y )
75
Gradient Descent Optimizer
Whats wrong with this algorithm ?
76
Gradient Descent Optimizer
millimeters
S Units ?
f(x,y) Units ?
intensity
G(x,y) Units ?
intensity / millimeters
S L G( x , y )
L Units ?
millimeters2 / intensity
77
Gradient Descent Optimizer
f( x )
S L G( x )
1
1
78
Gradient Descent Optimizer
f( x )
S L G( x )
S Large in high gradients
S Small in low gradients
1
1
79
Gradient Descent Optimizer
f( x )
S L G( x )
Works great with this function
Works badly with this function
80
Gradient Descent Variant
Driving Safe !
81
Regular Step Gradient Descent
f( x )

S D G( x )
If G changes direction
Di Di-1 / 2.0
D1
D1
D2
D1
82
Regular Step Gradient Descent
f( x )

S D G( x )
User Selects D1
User Selects Dstop
User Selects Gstop
D1
D1
D2
D1
83
Optimizers are like a car
Watch while you are driving !
84
Watch over your optimizer
ExampleOptimizer registering an image with
itself starting at (-15mm, -25mm)
85
Plotting the Optimizers Path
Mean Squared Differences
Step Length 1.0 mm
86
Plotting the Optimizers Path
Mean Squared Differences
Step Length 2.0 mm
87
Plotting the Optimizers Path
Mean Squared Differences
Step Length 5.0 mm
88
Plotting the Optimizers Path
Mean Squared Differences
Step Length 10.0 mm
89
Plotting the Optimizers Path
Mean Squared Differences
Step Length 20.0 mm
90
Plotting the Optimizers Path
Mean Squared Differences
Step Length 40.0 mm
91
Watch over your optimizer
ExampleOptimizer registering an image shifted
by (-15mm, -25mm)The optimizer starts at
(0mm,0mm)
92
Plotting the Optimizers Path
Mean Squared Differences
Step Length 1.0 mm
93
Plotting the Optimizers Path
Mean Squared Differences
Step Length 2.0 mm
94
Plotting the Optimizers Path
Mean Squared Differences
Step Length 5.0 mm
95
Plotting the Optimizers Path
Mean Squared Differences
Step Length 10.0 mm
96
Plotting the Optimizers Path
Mean Squared Differences
Step Length 20.0 mm
97
Plotting the Optimizers Path
Mean Squared Differences
Step Length 40.0 mm
98
Multi - Modality
99
Multi-Modality Registration
100
Multiple Image Modalities
Number of pairs
101
Multiple Image Modalities
More possible pairs
102
Mutual Information
103
Mutual Information
104
Exercise 21
105
Image Registration
106
Image Registration
include itkImageRegistrationMethod.h include
itkTranslationTransform.h include
itkMattesMutualInformationImageToImageMetric.h
include itkLinearInterpolateImageFunction.h inc
lude itkRegularStepGradientDescentOptimizer.h i
nclude itkImage.h include itkImageFileReader.h
include itkImageFileWriter.h include
itkResampleImageFilter.h
107
Image Registration
const unsigned int Dimension 2 typedef
unsigned char PixelType typedef itkImagelt
PixelType , Dimension gt
FixedImageType typedef itkImagelt PixelType ,
Dimension gt
MovingImageType typedef itkTranslationTransfo
rmlt double, Dimension gt TransformType typed
ef itkRegularStepGradientDescentOptimizer
OptimizerType typedef
itkLinearInterpolateImageFunctionlt
MovingImageType ,
double gt InterpolatorType typedef
itkMattesMutualInformationImageToImageMetriclt
FixedImageType ,
MovingImageType gt MetricType typedef
itkImageRegistrationMethodlt
FixedImageType , MovingImageType gt
RegistrationType
108
Image Registration
TransformTypePointer transform
TransformTypeNew() OptimizerTypePointer
optimizer OptimizerTypeNew() Interpolato
rTypePointer interpolator
InterpolatorTypeNew() MetricTypePointer
metric MetricTypeNew() Regist
rationTypePointer registrator
RegistrationTypeNew() registrator-gtSetTransfor
m( transform ) registrator-gtSetOptimizer(
optimizer ) registrator-gtSetInterpolator(
interpolator ) registrator-gtSetMetric(
metric ) registrator-gtSetFixedImage(
fixedImageReader-gtGetOutput() ) registrator-gtSetM
ovingImage( movingImageReader-gtGetOutput() )
109
Image Registration
registrator-gtSetFixedImageRegion(
fixedImageReader-gtGetOutput()-gtGetBufferedR
egion() ) typedef RegistrationTypeParameters
Type ParametersType transform-gtSetIdentity()
registrator-gtSetInitialTransformParameters(

transform-gtGetParameters()
) optimizer-gtSetMaximumStepLength( 4.00 )
optimizer-gtSetMinimumStepLength( 0.05
) optimizer-gtSetNumberOfIterations( 200
) optimizer-gtMaximizeOff()
110
Image Registration
metric-gtSetNumberOfHistogramBins( 20 )
// Metric specificmetric-gtSetNumb
erOfSpatialSamples( 10000 ) //
Metric specific try
registrator-gtStartRegistration () catch(
itkExceptionObject excp )
stdcerr ltlt Error in registration ltlt
stdendl stdcerr ltlt excp ltlt
stdendl transform-gtSetParameters(
registrator-gtGetLast
TransformParameters() )
111
Image Registration
typedef itkResampleImageFilterlt
FixedImageType , MovingImageType gt
ResamplerType ResamplerType Pointer
resampler ResamplerTypeNew() resampler-gtSet
Transform ( transform ) resampler-gtSetInput(
movingImageReader-gtGetOutput() ) FixedImageType
Pointer fixedImage fixedImageReader-gtGetOut
put() resampler-gtSetOrigin( fixedImage-gtGetOrigin
() ) resampler-gtSetSpacing( fixedImage-gtGetSpacin
g() ) resampler-gtSetSize(
fixedImage-gtGetLargestPossibleRegion()-gtGetSize
() ) resampler-gtUpdate()
112
Rotation TranslationParameter Scaling
113
Exercise 22
114
Image Registration
115
Image Registration
include itkImageRegistrationMethod.h include
itkCenteredRigid2DTransform.h include
itkMeanSquaresImageToImageMetric.h include
itkLinearInterpolateImageFunction.h include
itkRegularStepGradientDescentOptimizer.h includ
e itkCenteredTransformInitializer.h include
itkImage.h include itkImageFileReader.h incl
ude itkImageFileWriter.h include
itkResampleImageFilter.h
116
Image Registration
const unsigned int Dimension 2 typedef
unsigned char PixelType typedef itkImagelt
PixelType , Dimension gt
FixedImageType typedef itkImagelt PixelType ,
Dimension gt MovingImageType typedef
itkCenteredRigid2DTransformlt double gt
TransformType typedef itk
CenteredTransformInitializerlt

TransformType ,
FixedImageType ,

MovingImageType

gt InitializerType
117
Image Registration
TransformTypePointer transform
TransformTypeNew() InitializerTypePointer
initializer InitializerTypeNew() Opt
imizerTypePointer optimizer
OptimizerTypeNew() InterpolatorTypePointer
interpolator InterpolatorTypeNew() MetricTy
pePointer metric
MetricTypeNew() RegistrationTypePointer
registrator RegistrationTypeNew() regist
rator-gtSetTransform( transform )
registrator-gtSetOptimizer( optimizer )
registrator-gtSetInterpolator( interpolator )
registrator-gtSetMetric( metric
) registrator-gtSetFixedImage(
fixedImageReader-gtGetOutput() ) registrator-gtSetM
ovingImage( movingImageReader-gtGetOutput() )
118
Image Registration
registrator-gtSetFixedImageRegion(
fixedImageReader-gtGetOutput()-gtGetBuffe
redRegion() ) initializer-gtSetTransform (
transform ) initializer-gtSetFixedImage(
fixedImageReader-gtGetOutput()
) initializer-gtSetMovingImage(
movingImageReader-gtGetOutput() ) initializer-gtMo
mentsOn() initializer-gtInitializeTransform() r
egistrator-gtSetInitialTransformParameters(

transform-gtGetParameters() )
119
Image Registration
typedef OptimizerTypeScaleType
OptimizerScalesType OptimizerScalesType
optimizerScales(
optimizer-gtSetMaximumStepLength(
) ) const double translationScale 1.0 /
1000.0 optimizerScales 0 1.0
optimizerScales 1 translationScale
optimizerScales 2 translationScale
optimizerScales 3 translationScale
optimizerScales 4 translationScale
optimizer-gtSetScales( optimizerScales )
120
Image Registration
try registrator-gtStartRegistrati
on () catch( itkExceptionObject excp
) stdcerr ltlt Error in
registration ltlt stdendl stdcerr
ltlt excp ltlt stdendl transform-gtSetParam
eters(
registrator-gtGetLastTransformParameters() )
121
Image Registration
typedef itkResampleImageFilterlt
FixedImageType , MovingImageType gt
ResamplerType ResamplerType Pointer
resampler ResamplerTypeNew() resampler-gtSet
Transform ( transform ) resampler-gtSetInput(
movingImageReader-gtGetOutput() ) FixedImageType
Pointer fixedImage fixedImageReader-gtGetOut
put() resampler-gtSetOrigin( fixedImage-gtGetOrigi
n() ) resampler-gtSetSpacing( fixedImage-gtGetSpaci
ng() ) resampler-gtSetSize(
fixedImage-gtGetLargestPossibleRegion()-gtGetSiz
e() ) resampler-gtUpdate()
122
Multi - Resolution
123
Multi-ResolutionRegistration Framework

• Improve speed
• Accuracy
• Robustness

124
Multi-ResolutionRegistration Framework
125
Multi-ResolutionRegistration Framework

• Flexible framework allows change of
• Parameters
• Components
• between resolution levels

126
Enjoy ITK !