ITK Deformable Registration

1
ITK Deformable Registration
Finite Elements Methods
2
Open Source Disclaimer
Many of the slides in this talk were taken from
the ITK course presented at MICCAI 2003 by
Dr. James Gee (U. Penn)Brian Avants
(U. Penn)Tessa Sundaram (U. Penn)Dr.
Lydia Ng (AllenInstitue)
Of course, any errors are mine.
3
Deformable Registration
Finite Element Methodsfor Deformable
Registration
4
Deformable Transformation
y
y
Transform
x
x
Fixed Image
Moving Image
5
Deformable Transformation
y
y
Transform
x
x
Fixed Image
Moving Image
6
Deformable Transformation
y
x
7
Deformable Transformation
y
x
8
FEM Grid
y
FEM Grid
Resampling Grid
x
9
FEM Grid
y
FEM Grid
x
10
FEM Grid
y
FEM Grid
Computed Deformation
x
11
FEM Grid
y
FEM Grid
Displacements
Forces
x
12
FEM Iterative Linear System
Finite Element Methods
F
Forces
Vector N
U
Vector N
Displacements
K
Matrix NxN
Regularization
13
FEM Iterative Linear System
Finite Element Methods
F
F
U
K
U

?
K
14
FEM Iterative Linear System
N Number of Nodes
N x N
N
N

?
F
U
K
15
FEM Iterative Linear System
Iteratively Solving a Linear System
F
U
K

?
Linearization of a Physical Model
Image based forces
Node Displacements
16
FEM Energy Formulation
Find registration transformation that maximizes
17
FEM Energy Formulation
Benefits

• Intuitive easier to express constraints
• Powerful numerical solutions available
• Optimality of solutions easier to debug

Limitations

• Difficult / overhead to implement

18
Deformable Registration
To solve the deformation, consider only
displacements of the form
19
Deformable Registration
Linear Elements
20
Deformable Registration
a1
Element
f1
Shape Function
21
Deformable Registration
Element
f2
a2
Shape Function
22
Deformable Registration
a3
f3
Element
Shape Function
23
Deformable Registration
u
a3
f3
Element
a1
f1
f2
a2
Shape Functions
24
Deformable Registration
Higher Order Elements
25
Deformable Registration
a1
Element
f1
Shape Function
26
Deformable Registration
a4
f4
Element
Shape Function
27
Deformable Registration
Element
f2
a2
Shape Function
28
Deformable Registration
f5
a5
Element
Shape Function
29
Deformable Registration
a3
f3
Element
Shape Function
30
Deformable Registration
f6
a6
Element
Shape Function
31
Deformable Registration
a4
f4
u
a3
f3
Element
f6
a1
a6
f1
a5
f5
f2
a2
Shape Functions
32
Deformable Registration
Substitute uh into E, then minimizing with
respect to ai
33
BSplines Grid Image Grid
Calculation are made in an Element by Element
basis
34
BSplines Grid Image Grid
Elements are connected at Nodes at which the
displacement is solved
35
BSplines Grid Image Grid
Efficiency is gained by elemental computation
36
BSplines Grid Image Grid
Domain subdivision (Mesh) can be tailored to the
underlying geometry of the image.
37
FEM Solver
Start Iteration Loop

• Begin Loop by making physical assumptions and
  then taking the derivative of the similarity
  metric.
• End loop when the solution stabilizes.

Physical Assumptions
Solve
New Solution
Image Metric Derivative
End Iteration Loop
38
FEM Solver
Start Iteration Loop
K
Physical Assumptions
Solve
New Solution
U
Image Metric Derivative
F
End Iteration Loop
39
FEM Solver
Start Iteration Loop
K
U
Unew
F
If ( Unew Uold) lt e then Stop
40
KUF in Code
itkFEMRegistrationFilterIterativeSolve()
itkFEMSolverAssembleK()
FEMSolver Solve()
FEMSolver AddSolution()
FEMSolverAssembleF() calls FEMImageMetricLoad
Fe()
41
FEM-Based Registration Options

• Element Type
• Triangles
• Quadrilaterals
• Hexahedra
• Tetrahedra

42
FEM-Based Registration Options

• Continuum / Physical Model
• Linear elasticity
• Membrane
• Other specialized

43
FEM-Based Registration Options

• Mesh geometry
• Uniform grid vs. adaptive
• Anatomy-specific mesh

44
FEM-Based Registration Options

• Metric
• Mean square
• Normalized cross-correlation
• Mutual information
• Pattern intensity

45
ITK FEM Library
Introduction to the ITK Finite Element Library
46
ITK FEM Library

• Library for solving general FEM problems
• Object oriented
• C classes are used to
• specify the geometry and behavior of the elements
• apply external forces and boundary conditions
• solve problem and post-process the results

47
ITK FEM Library

• Applications
• Mechanical modeling
• Image registration

48
FEM Basics

• Mesh
• Nodes
• Points in space where solutions are obtained
• Elements
• e.g., 2-D triangular elements
• Loads
• e.g., gravity (body) load
• Boundary conditions
• e.g., nodes fixed in space

49
ITK FEM Elements

• Core of the library is the Element class
• Code is in two functionally independent parts
• Geometry and Physics
• Arbitrarily combined to create new elements
• Problem domain is specified by a mesh

Geometry
Physics
50
Loads

• Classes that apply external forces (loads) to
  elements
• Various types
• Easily extensible

51
Solvers

• Provide functionality to obtain and process the
  solution
• Different solution methods ? different solver
  classes
• Static problems
• Time dependent - dynamic problems

52
Solvers

• Use linear system wrappers to link FEM classes
  to an external numeric library
• Any numeric library can be used to solve the
  systems of linear equations in FEM problems
• VNL and ITPACK currently supported

53
Setting Up a FEM Problem

• Four-step process
• Select element classes
• Discretize problem domain
• Specify boundary conditions
• Specify/Apply external loads
• Two options
• Directly ? create proper objects in code
• Indirectly ? read object definitions from a file

54
Deformable Registration
FEM-Base RegistrationParameters
55
Parameter File Part 1
------------------------------------------------
--------- Parameters for the single- or
multi-resolution techniques --------------------
------------------------------------- 1 Number
of levels in the multi-resolution pyramid (1
single-res) 1 Highest level to use in the
pyramid 1 1 Scaling at lowest level for each
image dimension 8 Number of pixels per
element 1.e5 Elasticity (E) 1.e4 Density
(RhoC) 1. Image energy scaling 4
NumberOfIntegrationPoints 1
WidthOfMetricRegion 25 MaximumIterations
------------------------------- Parameters for
the registration -------------------------------
0 1.0 Similarity metric (0mean sq, 1ncc,
2pattern int, 3MI) 1.0 Alpha 1
DescentDirection 2 DoLineSearch (0never,
1always, 2if needed) 1.e1 TimeStep 1.e-15
Energy Reduction Factor
56
Parameter File Part 2
----------------------------------
Information about the image inputs
---------------------------------- 2
ImageDimension 256 Nx
(image x dimension) 256 Ny
(image y dimension) 128 Nz
(image z dimension - not used if
2D) brain_slice1.mhd ReferenceFileName brain
_slice1warp.mhd TargetFileName
--------------------------------------------------
----------------- The actions below depend on
the values of the flags preceding them. For
example, to write out the displacement fields,
you have to set the value of WriteDisplacementFi
eld to 1. --------------------------------------
----------------------------- 0
UseLandmarks? -
LandmarkFileName brain_result
ResultsFileName (prefix only) 1
WriteDisplacementField? brain_disp
DisplacementsFileName (prefix only) 1
ReadMeshFile? brain_mesh.fem
MeshFileName END
57
Configuring Parameters 1
this-gtDoMultiRes(true) this-gtm_NumLevels
nlev this-gtm_MaxLevel mlev for (jj0
jj lt ImageDimension jj) m_ImageScalingjj
dim for (jj0 jj lt this-gtm_NumLevels
jj) this-gtm_MeshPixelsPerElementAtEachResolu
tion(jj) p this-gtSetElasticity(e, jj)
this-gtSetRho(p, jj) this-gtSetGamma(g, jj)
this-gtSetNumberOfIntegrationPoints(ip, jj)
this-gtSetWidthOfMetricRegion(w, jj)
this-gtSetMaximumIterations(mit, jj)
58
Configuring Parameters 2
this-gtSetDescentDirectionMinimize() or this-
gtSetDescentDirectionMaximize()
this-gtDoLineSearch( n ) // n 0, 1, 2
this-gtSetTimeStep( t ) this-gtSetEnergyReduction
Factor( fbuf )
59
Configuring Parameters 3
this-gtm_ImageSize0 xdim this-gtm_ImageSize1
ydim if (dim 3) this-gtm_ImageSize2
zdim this-gtSetReferenceFile( imgfile1
) this-gtSetTargetFile( imgfile2
) this-gtUseLandmarks( true ) this-gtSetLandmarkF
ile( lmfile ) this-gtSetResultsFile( resfile
) this-gtSetWriteDisplacements( true
) this-gtSetDisplacementsFile( dispfile
) this-gtm_ReadMeshFile true this-gtm_MeshFileN
ame meshfile
60
Deformable Registration
FEM-Based Registration Writing the Code
../ Insight / Examples / Registration /
DeformableRegistration1.cxx
61
Header Declarations
include "itkImageFileReader.h" include
"itkImageFileWriter.h include
"itkFEM.h" include itkFEMRegistrationFilter.h"
62
Type Definitions
typedef itkImagelt unsigned char, 2 gt
fileImageType typedef itkImagelt float, 2 gt
ImageType typedef itkfemElement2DC0LinearQu
adrilateralMembrane
ElementType typedef
itkfemElement2DC0LinearTriangularMembrane

ElementType2 typedef itkfemFEMRegistrationFi
lterlt
ImageType,
ImageType gt
RegistrationType
63
Registering Objects
ElementTypeLoadImplementationFunctionPointer
fp1 itkfemImageMetricLoadImplementa
tionlt ImageLoadType
gtImplementImageMetricLoad DispatcherTypeReg
isterVisitor( (ImageLoadType)0 , fp1
) ElementType2LoadImplementationFunctionPoin
ter fp2 itkfemImageMetricLoadImpleme
ntationlt ImageLoadType
gtImplementImageMetricLoad DispatcherType2Re
gisterVisitor( (ImageLoadType)0 , fp2 )
64
Input / Output
RegistrationTypePointer registration
RegistrationTypeNew()
registration-gtSetConfigFileName( paramname
) registration-gtReadConfigFile()
65
Material and Element Setup
// Create the material properties itkfemMateri
alLinearElasticityPointer m m
itkfemMaterialLinearElasticityNew() m-gtGN
0 m-gtE registration-gtGetElasticity() m-gtA
1.0 // Cross-sectional area m-gth
1.0 // Thickness m-gtI 1.0
// Moment of inertia m-gtnu 0.
// Poisson's ratio m-gtRhoC 1.0
// Density // Create the element
type ElementTypePointer e1
ElementTypeNew() e1-gtm_mat dynamic_castlt
itkfemMaterialLinearElasticity gt( m
) registration-gtSetElement( e1
) registration-gtSetMaterial( m )
66
Running the Registration
registration-gtRunRegistration() registration-gtWr
iteWarpedImage() if ( registration-gtGetWriteDisp
lacements() ) registration-gtWriteDisplacementF
ield( 0 ) // x registration-gtWriteDisplacement
Field( 1 ) // y registration-gtWriteDisplaceme
ntFieldMultiComponent()
67
FEM - Deformable Registration
Example 1
68
Fixed Image
69
Moving Image
70
Registered Image
71
Registered Image
72
FEM - Deformable Registration
Example 2
73
Fixed Image
74
Moving Image
75
Registered Image
76
Registered Image
77
FEM - Deformable Registration
Example 3
78
Fixed Image
79
Moving Image
80
Registered Image
81
Registered Image
82
FEM - Deformable Registration
Example 4Elasticity value was doubled
83
Fixed Image
84
Moving Image
85
Registered Image
86
Registered Image
87
Enjoy ITK !