Medical Image Registration: Concepts and Implementation

1
Medical Image Registration Concepts and
Implementation

• Feb 28, 2006
• Jen Mercer

2
Registration

• Spatial transform that maps points from one image
  to corresponding points in another image

3
Registration Criteria

• Landmark-based
• Features selected by the user
• Segmentation-based
• Rigidly or deformably align binary structures
• Intensity-based
• Minimize intensity difference over entire image

4
Spatial Transformation

• Rigid
• Rotations and translations
• Affine
• Also, skew and scaling
• Deformable
• Free-form mapping

5
Registration Framework
6
Transforms

• xT(xp)T(x,ytx,ty,?)
• Goal Find parameter values that optimize image
  similarity metric

7
Optimizer

• Often require derivative of image similarity
  metric (S)

8
Jacobian and Image Gradient
9
Identity Transform

• Maps every point to itself
• Only used for testing
• Fixed set (C) set of points that remain
  unchanged by transform

10
Translation Transform

• Fixed set is an empty set

11
Scaling Transform

• Isotropic vs. anisotropic
• Fixed set is the origin of the coordinates

12
Scaling and Translation
13
Rotation Transform

• Fixed set is the origin

14
Rotations in Polar Coordinates
15
Optimization

• Search for value of ? that minimizes cost
  function S
• Gradient descent algorithm
• Update of parameter
• G is the variation from the gradient of the cost
  function
• ? is step length of algorithm

16
Combined Scaling and Rotation

• Dscaling factor
• Mcost function
• Apply transform to a point as

17
Add Translation

• Find fixed point of transformation
• Translation (d) is result of scaling and rotation

18
Scaling, Rotation,Translation

• Parbitrary point
• Cfixed point of transformation
• Dscaling factor
• Trotation angle
• P and C are complex numbers (xiy) or rei?
• Store derivates of P in Jacobian matrix for
  optimizer
• Rigid if D1, otherwise similarity transform

19
Affine Transformation

• Collinearity is preserved
• xA x T
• A is a complex matrix of coefficients
• With fixed point
• xA (xC) C
• A is optimized similar to the scaling factor

20
Quaternions

• Quotient of two vectors
• Q A / B
• Operator that produces second vector
• A Q ? B
• Represents orientation of one vector with respect
  to another, as well as ratio of their magnitudes
• Versor-rotates vector
• Tensor-changes vector magnitude

21
Scalars and Versors

• Quaternion represented by 4 numbers
• Versor
• Direction parallel to axis of rotation
• Rotation angle
• Norm function of rotation angle
• Tensor
• Magnitude

22
Unit Sphere Versor Representation
23
Versor Composition

• Versor definition (vector quotient)
• VCB B / C
• VBA A / B
• VCA A / C
• Versor composition
• VCA VBA ? VCB
• Not communative

24
Versor Addition
25
Optimization of Versors

• Versor exponentiation
• V2 V ? V
• V V1/2 ? V1/2
• T(V) ?
• T(Vn) n?
• Versor Increment

26
Rigid Transform in 3D

• Use quaternions instead of phasors
• PV(P-C)C
• PVPT, TC-VC
• Ppoint, VVersor, TTranslation, Cfixed point
• Transform represented by 6 parameters
• Three numbers representing versor
• Three components of fixed coordinate system

27
Numerical Representation of a Versor

• Right versor

28
Numerical Representation of a Versor

• -i k ? j
• -j i ? k
• -k j ? i
• Set of elementary quaternions
• i,j,k eip/2 , ejp/2, ekp/2

29
Numerical Representation of a Versor

• Any right versor can be represented as
• vxiyjzk
• x2y2z21
• Any generic versor can be represented in terms of
  the right versor parallel to its axis and the
  rotation angle as
• Vev?

30
Similarity Transform in 3-D

• Replace versor of rigid transform with quaternion
  to represent rotation and scale changes
• xQ(x-C)C
• xQxT, TC-QC

31
Image Interpolators

• 2 functions
• Compute interpolated intensity at requested
  position
• Detect whether or not requested position lies
  within moving-image domain

32
Nearest Neighbor

• Uses intensity of nearest grid position
• Computationally cheap
• Doesnt require floating point calculations

33
Linear Interpolation

• Computed as the weighted sum of 2n-1 neighbors
• ndimensionality of image
• Weighting is based on distance between requested
  position and neighbors

34
B-spline Interpolation

• Intensity calculated by multiplying B-spline
  coefficients with shifted B-spline kernels
• Higher spline orders require more pixels to
  computer interpolated value
• Third-order B-spline kernels typically used
  because good tradeoff between smoothness and
  computational burden

35
Metrics

• Scalar function of the set of transform
  parameters for a given fixed image, moving image,
  and transformation type
• Typically samples points within fixed image to
  compute the measure

36
Mean Squares

• Mean squared difference over all the pixels in an
  image
• Intensities are interpolated for the moving image
• For gradient-based optimization, derivative of
  metric is also required

37
Mean Squares

• Optimal value of zero
• Interpolator will affect computation time and
  smoothness of metric plot
• Assumes intensity representing the same
  homologous point is in both images
• Images must be from same modality

38
Mean Squares
Smoothness affected by interpolator
39
Normalized Correlation

• Computes pixel-wise cross-correlation between the
  intensity of the two images, normalized by the
  square root of the autocorrelation of each image
• For two identical images, metric 1
• Misalignment, metric lt1

40
Normalized Correlation

• -1 added for minimum-seeking optimizers

41
Normalized Correlation
42
Difference Density

• Each pixels contribution is calculated using
  bell-shaped function
• f(d) has a maximum of 1 at d0 and minimum of
  zero at d/-infinity
• d is difference in intensity b/w F and M

43
Difference Density

• ? controls the rate of drop off
• Corresponds to the difference in intensity where
  f(d) has dropped by 50

44
Difference Density

• Optimal value is N
• Poor matches small measure values
• Approximates the probability density function of
  the difference image and maximizes its value at
  zero

45
Difference Density

• Width of peak controlled by ?

46
Multi-modal Volume Registration by Maximization
of Mutual InformationWells W, Viola P, Atsumi
H, Nakajima S, Kikinis R
47
Registering Images from Same Modality

• Typical measure of error is sum of squared
  differences between voxels values
• This measure is directly proportional to the
  likelihood that the images are correctly
  registered
• Same measure is NOT effective for images of
  different modalities

48
Relationship Between Images of Different
Modalities

• Example We should be able to construct a
  function F() that predicts CT voxel value from
  corresponding MRI value
• Registration could be evaluated by computing
  F(MR) and comparing it to the CT image
• Via sum of squared differences (or correlation)
• In practice, this is a difficult and
  under-determined problem

49
Mutual Information

• Theory Since MR and CT both describe the
  underlying anatomy, there will be mutual
  information between the two images
• Find the best registration by maximizing the
  information that one image provides about the
  other
• Requires no a priori model of the relationship
• Assumes that max. info. is provided when the
  images are correctly registered

50
Notation

• Reference (fixed) volume u(x)
• Test (moving) volume v(x)
• x coordinates of the volume
• T transformation from coordinate frame of
  reference volume to test volume
• v(T(x)) test volume voxel associated with
  reference volume voxel u(x)

51
Mutual Information

• Defined in terms of entropies
• If there are any dependencies, H(A,B)ltH(A)H(B)

52
Maximizing Mutual Information

• h(v(T(x))) encourages transformations that
  project u into complex parts of v
• Last term of MI eqn contributes when u and v are
  functionally related
• Together, last two terms of MI eqn identify
  transforms that find complexity and explain it
  well

53
Parzen Windowing

• Used to estimate probability density P(z)
• Entropy estimated based on P(z)

54
Finding Maximum of I(T)

• To find maximum of mutual information, calculate
  its derivative
• Derivative of reference volume is 0, b/c not a
  function of T
• Entropies depend on covariance of Parzen window
  functions

55
Stochastic Maximization of Mutual Information

• Similar to gradient descent
• Steps are taken that are proportional to dI/dT
• Repeat
• A ? sample of size NA drawn from x
• B ? sample of size NB drawn from x
• T ? T?(dI/dT)
• ? is the learning rate
• Repeated a fixed of times, or until convergence

56
Stochastic Approximation

• Uses noisy derivative estimates instead of the
  true derivative for optimizing a function
• Authors have found that technique always
  converges to a transformation estimate that is
  close to locally optimal
• NANB50 has been successful
• The noise introduced by the sampling can
  effectively penetrate small local minima

57
MRI-CT Example

• Coarse to fine registration
• Images were smoothed by convolving with binomial
  kernel
• Rigid transform represented by displacement
  vectors and quaternions
• Images were sampled and tri-linear interpolation
  was used
• 5 levels of resolution
• 10000, 5000(4) iterations

58
Initial Condition of MR-CT Registration
59
Final Configuration for MR-CT Registration
60
Initial Condition of MR-PET Registration
61
Final Configuration for MR-PET Registration
62
Application

• Register 2 MRIs of brain (SPGR and T2-weighted)
  to visualize anatomy and tumor
• Create at 3-D model for surgical planning and
  visualization

63
3-D Model
Tumor(green), Vessels(red), Ventricles(blue),
Edema (orange)
64
Correlation

• Conventional correlation aligns two signals by
  minimizing a summed quadratic difference between
  their intensities
• If intensity of one signal is negated, then
  intensities no longer agree, and alignment by
  correlation will fail
• Mutual information is not affected by negation of
  either signal

65
Occlusion

• Correlation is significantly affected by
  occlusion because intensity is substantially
  different
• Occlusion will reduce mutual information at
  alignment
• But mutual information measure degrades
  gracefully when subject to partially occluded
  imagery

66
Comparison to Other Methods

• Many researchers use surface-based methods to
  register MRI and PET imagery
• Need for reliable segmentation is a drawback
• Others use joint entropy to characterize
  registration
• not robust difficulty describing partial
  overlap
• Mutual Information is better because it has a
  larger capture range
• Additional influence from term that rewards for
  complexity in portion of test volume into which
  reference volume is transformed

67
Comparison to Other Methods

• Woods et al. register MR and PET by minimizing
  range of PET values associated with a particular
  MR intensity value
• Effective when test volume distribution is
  Gaussian
• Mutual Information can handle data that is
  multi-modal
• Woods measure is sensitive to noise and outliers

68
Conclusions

• Intensity based techniques work directly with
  volumetric data (vs. segmentation)
• Mutual information does not rely on assumptions
  about nature of imaging modalities
• Have also used this technique to register 3D
  volumetric images to video images of patients