Transcript Overview of Medical Image Registration

Overview of Medical
Image Registration
J. Michael Fitzpatrick,
Department of Electrical Engineering and Computer Science
Vanderbilt University, Nashville, TN
CS/EECE 359, Fall, 2007
Acknowledgements
Benoit M. Dawant, PhD, EECS
Robert L. Galloway, PhD, BME
NSF
William C. Chapman, MD, Surgery
NIH
Jeannette L. Herring, PhD, EECS
Jim Stefansic, PhD, Psychology
Diane M. Muratore, MS, BME
Matthew Wang, PhD, IBM
David M. Cash, MS, BME
Jay B. West, PhD, Accuray, Inc.
Steve Hartman, MS, BME
Derek L. G. Hill, PhD Kings College
W. Andrew Bass, BME
Calvin R. Maurer, Jr., PhD, Stanford U.
Computed Tomography (1972)
Siemens CT Scanner
(Somatom AR)
3D Cross-sectional Image
“voxels”
(“volume elements”)
Magnetic Resonance Imaging
GE MR
Scanner
(Signa
1.5T)
Positron Emission Tomography
GE PET Scanner
Physician has 3 or more views.
CT
(bone)
MR
(wet tissue)
PET
(biological
activity)
Combining multiple images
requires image registration
Image Registration: Definition
Determination of corresponding points in two
different views of the same object
Motion relative to the scanners
can be three-dimensional.
Slice orientations vary widely.
transverse
sagittal
coronal
Views may be very different.
But all orientations and all views
can be combined if we have the 3D
point mapping.
Combining Registered Images
= “Image Fusion”
CT
CT + MR
MR
PET
MR + PET
Rigid Registration: Definition
Rigid Registration = Registration using a
“rigid” transformation
Rigid Transformation
Distances between all points remain constant.
Rigid
6 degrees of freedom
Non-rigid
Nonrigid Transformations
can be very complex!
[Thompson, 1996]
Registration Methods
• “Retrospective” methods


Match anatomical features: e.g., surfaces
Maximize similarity of intensity patterns
• “Prospective” methods


Non-invasive: Match skin markers
Invasive: Match bone-implanted markers
Mutual Information:
An Example of
Intensity-Based* Registration
*Also known as Voxel-Based Registration
MR intensity
2D Intensity Histogram (Hill94)
MR
CT intensity
CT
Misregistration Blurs It
0 cm
2 cm
MR
CT
MR
PET
Hill, 1994
5 cm
Mutual Information
(Viola, Collignon, 1996)
•
•
•
•
•
•
A measure of histogram sharpness
Most popular “intensity” method
Assumes a search method is available
Stochastic, multiresolution search common
Requires a good starting pose
May not find global optimum
Example: Mutual Information
Studholme,
Hill, Hawkes,
1996,
“Automated
3D registration
of MR and CT
images of the
head”, MIA,
1996
Goal for
st
1
Day
The Iterative Closest-Point
Algorithm:
An Example of
Surface-Based Registration
Iterative Closest-Point Method
(Besl and McKay, 1992)
•
•
•
•
•
•
Minimizes a positive distance function
Most popular surface method
Assumes surfaces have been delineated
Guaranteed to converge
Requires a good starting pose
May not find global optimum
Start with two surfaces
Reorient one (somehow)
Reorient one (somehow)
Reorient one (somehow)
Pick points on moving surface
Pick points on moving surface
Remove moving surface
Points become proxy for surface
Find closest points on stationary surface
Got here 10/4/2005
Measure the total distance
Remove stationary surface
Points become proxy for surface
Register point sets (rigid)
Register point sets (rigid)
Restore stationary surface
Find (new) closest points
Find (new) closest points
Remove stationary surface
Remove stationary surface
Register Points
Register Points, and so on…
Iterative Closest-Point Algorithm:
• Find closest points
• Measure total distance
• Register points
Example: ICP for Head
Dawant et al.
Example: ICP for Vertebra
Muratore,
Herring,
Dawant,
Galloway
Got to here in CS 395, 8/28/03
ICP requires surface delineation, which is a
problem in Image Segmentation
Example:
Level Set
Segmentation
(Dawant
et al.)
http://www.vuse.vanderbilt.edu/~dawant/levelset_examples/
The Fiducial Marker:
An Example of…
(1) Point-Based Registration
(2) Prospective Methods
(3) Image-guided Therapy
Image-Guided Surgery
Just another image registration problem.
One view is
an image....
...and the other is
the patient.
Acustar
™
Allen, Maciunas,
Fitzpatrick, and
Galloway
1986-1996
(J&J  Z-Kat)
Posts
are implanted
into the skull.
[Maurer, et al., TMI, 1997]
Acustar
™
Allen, Maciunas,
Fitzpatrick, and
Galloway
1986-1996
(J&J  Z-Kat)
Liquid in
marker
shows up
in image
[Maurer, et al., TMI, 1997]
Divot cap is
localizable
in OR
Acustar
™
Allen, Maciunas,
Fitzpatrick, and
Galloway
1986-1996
(J&J  Z-Kat)
Marker center and cap center
occupy the same position relative
[Maurer,
et al.,
TMI, 1997]
to the
post
Acustar
™
Allen, Maciunas,
Fitzpatrick, and
Galloway
1986-1996
(J&J  Z-Kat)
Marker center and cap center
occupy the same position relative
[Maurer,
et al.,
TMI, 1997]
to the
post
Point-based, Rigid Registration
Find
Findall
corresponding
“fiducial” points
Rigid
transformation
View 1
= “Space” 1
Align
corresponding
fiducials
View 2
= “Space” 2
“targets” are
also aligned
What to Optimize?
• Mean-square “Fiducial Registration Error” (FRE2)

Known as the “Orthogonal Procrustes Problem” in
statistics since 1950s.
• Robust estimators (median, M-estimators)

Less sensitive to “outliers”
Color key: Major problems solved, Negligible work
Minimization of FRE2
(Shönemann, Farrell, 1966)
•
•
•
•
•
•
Minimizes a positive distance function
Most popular point method
Assumes points have been localized
Guaranteed to converge
Does not require a good starting pose
Always finds global optimum
Sum of Squares: Step 1
Center the points:
x   xi N
y   yi N
" Centered" points :
~
xi  xi  x ; ~yi  yi  y
Centered
y
x
Step 2 (Shönemann, Farrell, 1966)
Centered
Determine the Rotation:
H  ~
xi ~
y it
SVD : H  UΛV ,
where
t
U tU  V tV  I
Λ  diag( 1 , 2 , 3 )
1  2  3  0
R  V D U t , where

D  diag 1, 1, det(VU t )

Centered
and Rotated
Finding Points = “Localization”
CT
MR-T1
MR-T2
Competition: Acustar v. Leibinger
Competition: Acustar v. Leibinger
Accuracy
Measures of Registration Error
Fiducial
Localization
Error (FLE)
View 2
View 1
Fiducial Registration
Error (FRE)
Registered Views
Target Registration
Error (TRE)
Acustar, 3mm Slice CT-physical
TRE = 0.3 to 0.7 mm
Registration of Head Images:
The State of the Art
Retrospective Median
Maximum (Acustar)
Best CT-MR :
0.6 mm
3.0 mm
(0.5 mm)
Poor CT-MR:
5.4 mm
61 mm
(0.5 mm)
Best PET-MR: 2.5 mm
6.0 mm
(1.7 mm)
Poor PET-MR: 5.3 mm
15 mm
(1.7 mm)
And how do we know?…
Retrospective Image Regstration Evaluation
1995-2007
External site
Vanderbilt
Access: 150+ participants in 20 countries
Evaluation: 57 participants in 17 countries
Goal for
nd
2
Day
Error Theory for Minimization
of
Mean-square FRE
Start with Assumptions
about FLE
Independent, normal, isotropic, zero mean
Space 1
Space 2
“Effective” FLE
FLE
FLE1
FLE 2
Space 1
FLE
2
Space 2
 FLE  FLE
2
1
2
2
FRE Statistics: Sibson 79
Approximate Solution:


To O FLE 2 :
FRE  FLE  / 3 N ,
2
2
2
where the  2 is chi - square
with 3 N  6 DOF.
FRE 2  (1  2 /N ) FLE 2
FRE Statistics: Sibson 1979
Approximate Solution:


To O FLE 22 :
FRE  cFLE
 ,  / 3N ,
2
Configuration doesn’t
matter!
2
2
2
22
where
the

is
where the
is chi
chi -- square
square
with
with 33 N
N
 66 DOF.
DOF.
22
22
FRE

(
1

2
/N
)
FLE
FRE  (1  2 /N ) FLE
TRE statistics, 1998
Approximate Solution:
d1
d3
Principal axes
1
2
 TRE 
N
d2
[Fitzpatrick, West,
Maurer, TMI, ’98]
Configuration does
matter.
  d12 d 22 d 32  
2


1



/
3

FLE

  2

2
2 
f2
f3  
  f1
Marker Placement
TRE
for
FLE
of
1mm
1mm
2mm
[West et al.,
Neurosurgery,
April, 2001]
4mm
3mm
2mm
1mm
FRE = 1mm
Probability density
A distribution would be better
<TRE2>
95% level
TRE2
And what about direction?
TRE statistics, 2001
Approximate Solution:
TRE1
TRE2
TRE3
TRE i  N 0,  i ,
where i  1,2,3
denote independen t
components
along orthogonal
directions .
[West and Fitzpatrick., TMI, Sep 2001]
Some Remaining Problems
Isotropic Scaling
[Actually now solved: Batchelor, West, Fitzpatrick,
Proc. of Med. Im. Undstnd. & Anal. ,
Jul 2002]
Anisotropic Scaling
(Iterative Solution Only)
View 1
Register M points sets
simultaneously
View 2
;
The “Generalized”
Procrustes Problem
View 3
View M
(Iterative Solution Only)
Spatial Weighting
(Iterative Solution Only)
Other Unsolved Problems
• What is the statistical effect on TRE of
dropping or adding a fiducial?
• Does anisotropy in FLE always, sometimes, or
never makes TRE worse?
• How do we configure markers on a given
surface so as to minimize TRE over a given
region?
• Is there a correlation between FRE and TRE?
Other Unsolved Problems (cont.)
• Extenion to perspective transformations.
• Extension to surface matching.
End of Talk
Additional slides follow
Categories within error prediction
•
•
•
•
•
•
•
Number of point sets: Two or more
Scaling: Isotropic or anisotropic
Point-wise weighting: equal or unequal
Anisotropic weighting
Cost function: squared error or other
Point-wise FLE: equal or unequal
Spatial FLE: isotropic or anisotropic...
Key: Approximate,
Negligible progress
Anisotropic Scaling
Problem Statement:
Iterative Algorithm:
R, t = rotation, translation
wi2= point weighting
S = diag( sx , sy , sz )
Search for S that
minimizes
min FRE ( R, t , xi)
2
R, t
Given {xi yi wi} find R, t, S
to minimize mean FRE2
N
 (1 / N ) w RS xi  t  yi
2
i
i
for xi  Sxi .
sz
2
sx
Search
space
sy
Scaling: Anisotropic II
Problem Statement:
Iterative Algorithm:
R, t = rotation, translation
wi2= point weighting
S = diag( sx , sy , sz )
Given {xi yi wi} find R, t, S
to minimize mean FRE2
N
 (1 / N ) w SRxi  t  yi
2
i
i
2
II is (surprisin gly)
much harder tha n I :
Search for S , R, t that
minimizes
FRE 2 ( S, R, t )
Spatial Weighting
Problem Statement:
Iterative Algorithm:
Koschat & Swane [91]
R, t = rotation, translation
wi2= point weighting
S = diag( sx , sy , sz )
A = diag( ax , ay , az )
Chu & Trendafilo v [96]
Partial Solution:
Given {xi yi wi} find R, t, S
to minimize mean FRE2
N
 (1 / N ) w A( RS xi  t  yi )
2
i
i
2
Batchelor and
Fitzpatrick [2000]
Generalized Procrustes Problem
Cost function
Iterative method
(only)
Add Isotropic Scaling
Approximate Solution:
To O(  ) :
2
FRE  FLE  / 3 N .
2
2
2
 2 has 3 N  7 DOF
FRE2 = sum of
squared
fiducial
registration
errors
FRE
2
7 

2
 1 
 FLE
 3N 
FRE: Generalized + Scaling
Approximate Solution:
To O(  2 ) : FRE 2 
FLE 2  2 3 N M  1.
 2 has (M - 1)(3N - 7) DOF.
FRE2 = sum of
squared
fiducial
registration
errors
FRE
2
7 

2
 1 
FLE

 3N 
TRE statistics with scaling
Approximate Solution:
To O  2 
TRE 2
TRE2 =
target
registration
error
West and Fitzpatrick [2001]
Applications of TRE
Statistics
Error Bounds
Surgical
Paths
Radiation
Isodose
Contours
Probe Design
IREDs are fiducials
FLE
Tip = “target”
TRE
Fiducial-Specific FRE
FRE
2
i
 FLE  TRE x i 
2
x1
x2
2
Poor fiducial
alignment tends
to occur where
target registration
is good!!
Four Solution Methods
( All work equally well [Eggert91]! )
t
~
~
t is easy, and R depends only on H   w x yi .
2
i i
1. SVD : R  VU , where H  UDV .
t
t


t 1 / 2
2. Square Root : R  H HH
.
3. Quaternion :
(a) Form 4x4 matrix Q from elements of H .
(b) Find eigenvecto r q of Q with largest eigenvalue .
(c) Elements of R are quadratic in q1 , q2 , q3 , q4 .
t
4. Dual - number quaternion : [Walker91] .
Generalized Procrustes Problem
(We’ve already done it for M=2.)
Problem Statement:
Illustration:
For M sets of N points each,
{x (i m ) | i, m  1 N, 1 M },
find {R ( m ) , S ( m ) , t ( m ) | m  1 M }
to minimize
M
M
N
 (m)  (k )
  w xi  xi
m k m i
2
2
i
 (m)
(m) (m) (m)
(m)
where x i  R S x i  t
Generalized Procrustes Problem
Iterative Algorithm:
 (m)
Start with xi  xi( m ) .
Illustration:
Iterate this :
M
 (m)
xˆ i  (1 / M ) x i .
m
Find R ( m ) , S ( m ) , t ( m ) to minimize
*
2
N
 (m)
 w xi  xˆ i , m  1 ,, M ,
2
i
i
 (m)
where xi  R ( m ) S ( m ) xi( m )  t ( m ) .
until changes are neglible.
*Subject to S(m)
normalization
Approximation Method
(due to Sibson, 1979)
(1) Assume small FLEs :   1 ,  2 .
 x i  z i  1f i
y i  R (true) z i  t (true)   2 g i
2 Note that output same if
 y i  x i   2 g i .
R (true)  I , t (true)  0 :
Approximation Method (cont.)
(3) Expand in , dropping higher order 
R  I  R (1)
t  O ( )
FRE 2  O ( 2 )
TRE 2  O ( 2 )
FRE Statistics
Approximate Solution:
Problem Statement:
Given : {z i } for i  1, N ,
To O(  2 ) :
 1 and  2
FRE  
2
Find : Statistics for
FRE  min 1
R, t N
2
N
 Rx
i
 t  yi
2
2

2
1

2
2

2
/ N,
where the  2 is chi - square
distribute d with 3N  6
degrees of freedom.
i
FRE 2  (1  2 /N ) FLE 2
TRE statistics with scaling
Approximate Solution:
Problem Statement:
Given : {x i } for i  1, N ,
To O  2 
 1 , and  .
TRE 2
Find : Statistics for
West, Fitzpatrick,
and Batchelor [2001]
TRE 2  Rs x  t  y , when
FRE  min 1
R, t , s N
2
N
 Rs x
i
i
 t  yi
2
What do “solved” and “unsolved”
mean?
• “Solved”, working definition:



Reduced to solving algebraic equations
Iterative algorithm that converges to solution
Approximate solution accurate to O ()
• “Unsolved”:

Not solved
Point-wise weighting: Equal or Unequal
(We’ve just looked at this one.)
Problem Statement:
Solution:
R, t = rotation, translation
wi2= point weighting
Given {xi yi wi} find R, t
to minimize mean FRE2
N
 (1 / N ) w Rxi  t  yi
2
i
i
2
See previous
slides again!
1. Performing a Registration
a.k.a. The “Orthogonal Procrustes Problem”
Problem Statement:
xi = point in “from” set; yi = point in “to” set.
t = translation vector.
R = 3x3 rotation matrix (therefore RtR = I ).
Rxi + t
Given {xi yi wi } find R, t
2
to
minimize
mean
FRE
y
xi
i
N
 (1 / N ) w Rxi  t  yi
2
i
2
i
( usually wi=1)
2. Predicting Registration Error
Input--•fiducial positions
r
•target position,
r
•FLE distribution
View 2
View 1
Output---
Output---
statistics for
FRE
statistics for
TRE
Registered Views
Isotropic Scaling
Problem Statement:
Solution:
R, t = rotation, translation
wi2= point weighting
s = isotropic scaling
Given {xi yi wi} find R, t, s
to minimize mean FRE2
N
 (1 / N ) w Rs xi  t  yi
2
i
i
2
(1) Set s  1
Find R, t
(2) s 
w2 R ~
x ~
y

~
~
w
x

x

i
i
2
i
i
i
i
i