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 • x’=T(x|p)=T(x,y|tx,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 • D=scaling factor • M=cost 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 • P=arbitrary point • C=fixed point of transformation • D=scaling factor • Θ=rotation angle • P and C are complex numbers (x+iy) or reiθ • Store derivates of P in Jacobian matrix for optimizer • Rigid if D=1, otherwise similarity transform

19. Affine Transformation • Collinearity is preserved • x’=A x + T • A is a complex matrix of coefficients • With fixed point • x’=A (x–C) + 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 • Θ(V) = θ • Θ(Vn) = nθ • Versor Increment

26. Rigid Transform in 3D • Use quaternions instead of phasors • P’=V*(P-C)+C • P’=V*P+T, T=C-V*C • P=point, V=Versor, T=Translation, C=fixed 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]= [eiπ/2 ,ejπ/2, ekπ/2]

29. Numerical Representation of a Versor • Any right versor can be represented as • v=xi+yj+zk • x2+y2+z2=1 • Any generic versor can be represented in terms of the right versor parallel to its axis and the rotation angle as • V=evθ

30. Similarity Transform in 3-D • Replace versor of rigid transform with quaternion to represent rotation and scale changes • x’=Q*(x-C)+C • x’=Q*x+T, T=C-Q*C

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 • Doesn’t require floating point calculations

33. Linear Interpolation • Computed as the weighted sum of 2n-1 neighbors • n=dimensionality 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 <1

40. Normalized Correlation • -1 added for minimum-seeking optimizers

41. Normalized Correlation

42. Difference Density • Each pixel’s contribution is calculated using bell-shaped function • f(d) has a maximum of 1 at d=0 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)