A Fast Diffeomorphic Image Registration Algorithm

1. A Fast Diffeomorphic Image Registration Algorithm John Ashburner 2008

2. Overview • Parameterization • Scaling and Squaring • Objective Function • Optimization • Group-wise Registration • Comparison with LDDMM

3. Julian Huxley (1932). PROBLEMS OF RELATIVE GROWTH One essential fact about growth is that it is a process of self-multiplication of living substance – i.e. that the rate of growth of an organism growing equally in all its parts is at any moment proportional to the size of the organism. A constant partition of growth intensity between different regions implies constant differences in their rates of growth. Thus any genes controlling relative size of parts will have to exert their action by influencing the rates of processes, …

4. Allometry • Julian Huxley proposed that simpler relationships among lengths, volumes or areas of anatomical structures could be found by working with the logarithms of the measures. exp(u) = φ(1) = ∫u φ(t)dt,where φ(0) = 1 1 t=0

5. Parameterization Diffeomorphic Anatomical Registration Through Exponentiated Lie Algebra Deformations parameterized by a single flow field, which is considered to be constant in time. Not really a proper Lie Group. Often referred to as a one parameter subgroup.

6. Euler Integration • Parameterising the deformation • φ(0)(x) = x • φ(1)(x) = ∫u(φ(t)(x))dt • u is a flow field to be estimated • Scaling and squaring is used to generate deformations. • c.f. matrix exponentiation 1 t=0

7. Euler integration • The differential equation is dφ(x)/dt = u(φ(t)(x)) • By Euler integration φ(t+h) = φ(t) + hu(φ(t)) • Equivalent to φ(t+h) = (x + hu) oφ(t)

8. Simple integration φ(1/8) = x + u/8 φ(2/8) = φ(1/8)oφ(1/8) φ(3/8) = φ(1/8)oφ(2/8) φ(4/8) = φ(1/8)oφ(3/8) φ(5/8) = φ(1/8)oφ(4/8) φ(6/8) = φ(1/8)oφ(5/8) φ(7/8) = φ(1/8)oφ(6/8) φ(8/8) = φ(1/8)oφ(7/8) 7 compositions Scaling and squaring φ(1/8) = x + u/8 φ(2/8) = φ(1/8)oφ(1/8) φ(4/8) = φ(2/8)oφ(2/8) φ(8/8) = φ(4/8)oφ(4/8) 3 compositions Similar procedure used for the inverse. Starts with φ(-1/8) = x - u/8 For (e.g) 8 time steps

9. Scaling and squaring example

10. Deformations at different times

11. Jacobians • Jacobian fields can also be obtained by scaling and squaring. • If warps are composed by: ϕC=ϕB○ϕA then Jacobian matrices are obtained by the following matrix multiplication: JϕC=(JϕB○ϕA) JϕA

12. Jacobian determinants remain positive (almost)

13. See also… • C. Moler and C. van Loan. “Nineteen Dubious Ways to Compute the Exponential of a Matrix, Twenty-Five Years Later”. SIAM Review 45(1):3-49 (2003). • V. Arsigny, O. Commowick, X. Pennec and N. Ayache. “A Log-Euclidean Polyaffine Framework for Locally Rigid or Affine Registration”. Proc. Of the 3rd International Workshop on Biomedical Image Registration (WBIR'06), 2006, pp. 120-127. LNCS vol 4057. Springer-Verlag, Utrecht, NL. • V. Arsigny, O. Commowick, X. Pennec and N. Ayache. “A Log-Euclidean Framework for Statistics on Diffeomorphisms”. Proc. of the 9th International Conference on Medical Image Computing and Computer Assisted Intervention (MICCAI'06), 2006, pp. 924-931. LNCS 4190. Springer-Verlag, Berlin, Germany. • M. Hernandez, M. N. Bossa, and S. Olmos. “Registration of anatomical images using geodesic paths of diffeomorphisms parameterized with stationary vector fields”. IEEE workshop on Math. Meth. in Biom. Image Anal. (MMBIA’07), 2007.

14. Overview • Parameterization • Objective Function • Optimization • Group-wise Registration • Future directions

15. Original objective function Simultaneously minimize the sum of • Likelihood component • Sum of squares difference • ½ ∑i∑k(tk(xi) – μk(φ(1)(xi)))2 • φ(1) parameterized by u • Prior component • A measure of deformation roughness • ½uTHu

16. Likelihood Term • Images assumed to be partitioned into different tissue classes. • E.g., a 3 class registration simultaneously matches: • Grey matter with grey matter • White matter wit white matter • Background (1 – GM – WM) with background

17. Prior Term • ½uTHu • DARTEL has three different models for H • Membrane energy • Linear elasticity • Bending energy • H is very sparse An example H for 2D registration of 6x6 images (linear elasticity)

18. “Membrane Energy” Penalises first derivatives. Sum of squares of the elements of the Jacobian (matrices) of the flow field. Sparse Matrix Representation Convolution Kernel

19. “Bending Energy” Penalises second derivatives. Sparse Matrix Representation Convolution Kernel

20. “Linear Elasticity” • Decompose the Jacobian of the flow field into • Symmetric component • ½(J+JT) • Encodes non-rigid part. • Anti-symmetric component • ½(J-JT) • Encodes rigid-body part. • Penalise sum of squaresof symmetric part. • Trace of Jacobianencodes volume changes.Also penalised.

21. Regularization models “Membrane energy” Images registered using a small deformation approximation “Bending energy”

22. Overview • Parameterization • Objective Function • Optimization • Gauss-Newton • Multi-grid • Group-wise Registration • Comparison with LDDMM

23. Optimization • Uses Gauss-Newton • Requires a matrix solution to a very large set of equations at each iteration u(k+1) = u(k) - (H+A)-1b • b are the first derivatives of objective function • A is a sparse matrix of second derivatives • Computed efficiently, making use of scaling and squaring

24. Computing Derivatives • First derivatives w.r.t. flow field are considered as a vector field when being computed. • Second derivatives are treated as a field of positive definite symmetric matrices when computed. • Approximations involved, as a result of multiple interpolations used for the computations. • Derivatives resulting from evolution from time 0 to time 1 are simply the sum of the derivatives from time 0 to time 0.5 and those from time 0.5 to time 1. • Derivatives from 0.5 to 1 can be simply computed from those from time 0 to time 0.5. • Similarly derivatives from time 0 to 0.5 are the sum of derivatives from 0 to 0.25 and 0.25 to 0.5. • etc

25. Relaxation • To solve Mx = c Split M into E and F, where • E is easy to invert • F is more difficult • If M is diagonally dominant (membrane energy):x(k+1) = E-1(c – F x(k)) • Otherwise regularize (bending or linear elastic energy):x(k+1) = x(k) + (E+sI)-1(c – M x(k)) • Diagonal dominance is when |mii| > Σi≠j |mij|

26. 2nd derivs of prior term 2nd derivs of likelihood term M = H+A = E+F Easy to invert Difficult to invert

27. Relaxation Strategies • Jacobi’s method if not done in place, such that updates are not used immediately. • Easier to implement in MATLAB • Gauss-Siedel when done in place. • Faster convergence and uses less memory. • “Red-black” alternating update scheme is used with membrane energy. • Sweeps in alternating directions are used for updates when regularization is linear-elastic or bending energy. • Both methods fit high frequencies quickly, but low frequencies slowly.

28. Relaxation 128x128

29. Highest resolution Full Multi-Grid Lowest resolution

30. A • Prolongation of low resolution solution to current resolution. • Add this to existing solution. • Perform a few iterations of relaxation. • Restrict residuals down to lower resolution.

31. B • Prolongation of low resolution solution to current resolution. • Add this to existing solution at current resolution. • Perform a few iterations of relaxation. • Prolongation of solution to higher resolution.

32. C • Restrict high resolution residuals to current resolution. • Perform a few iterations of relaxation. • Restrict residuals down to lower resolution.

33. E • Restrict higher resolution residuals to current resolution. • Obtain exact solution by matrix inversion. • Prolongation of solution to higher resolution.

34. See also… • W. H. Press, S. A. Teukolsky, W. T. Vetterling and B. P. Flannery. Numerical Recipes in C (Second Edition). Cambridge University Press, Cambridge, UK. 1992. • Chapter 15, Section 5 explains Gauss-Newton optimization (Levenberg-Marquardt without the regularisation). • Chapter 19, Section 6 explains the basics of multi-grid methods.

35. Overview • Principles • Objective Function • Optimization • Group-wise Registration • Multinomial distributions • Simultaneous registration of GM & WM • Tissue probability map creation • Comparison with LDDMM

36. “Average Shaped” Template • For CA, work in the tangent space of the manifold, using linear approximations. • Average-shaped templates give less bias, as the tangent-space at this point is a closer approximation. • For spatial normalisation of fMRI, warping to a more average shaped template is less likely to cause signal to disappear. • If a structure is very small in the template, then it will be very small in the spatially normalised individuals. • Smaller deformations are needed to match with an average-shaped template. • Smaller errors.

37. Average shaped templates Average on Riemannian manifold Linear Average (Not on Riemannian manifold)

38. Template Generation Initial Average Iteratively generated from 471 subjects. Began with rigidly aligned tissue probability maps. Regularization lighter for later iterations. After a few iterations Final template

39. Multinomial Model • The extended model is multinomial for matching tissue class images. log p(t|μ,ϕ) = ΣjΣk tjk log(μk(ϕj)) t – individual GM, WM and background μ– template GM, WM and background ϕ– deformation • A general purpose template should not have regions where log(μ) is –Inf.

40. ϕ1 t1 ϕ2 μ t2 t5 ϕ3 t3 ϕ5 t4 ϕ4 Generative Model • p(ϕ1,t1, ϕ2,t2, ϕ3,t3,… μ)= p(t1|ϕ1,μ) p(t2|ϕ2,μ) p(t3|ϕ3,μ) … p(μ) • MAP solution obtainedfor template. • Requires p(μ)

41. Laplacian Smoothness Priors on template 2D 3D

42. Template modelled as softmax of a Gaussian process μk(x) = exp(ak(x))/(Σj exp(aj(x))) MAP solution determined for a, by Gauss-Newton optimisation, using multi-grid.

43. Determining amount of regularisation • Matrices too big for Bayesian variance component estimation. • Used cross-validation. • Smooth an image by different amounts, see how well it predicts other images: Nonlinear registered Rigidly aligned • log p(t|μ) = ΣjΣk tjk log(μjk)

44. ML and MAP templates from 6 subjects Nonlinearly Registered Rigidly registered ML MAP log

45. 471 Subject Average

46. 471 Subject Average

47. 471 Subject Average

49. 471 Subject Average