Pre-processing for “Voxel-Based Morphometry”

1. Pre-processing for “Voxel-Based Morphometry” John Ashburner The Wellcome Trust Centre for Neuroimaging 12 Queen Square, London, UK.

2. Contents • Introduction • Segmentation • DARTEL Registration

3. Voxel-based Morphometry • Pre-process the images of lots of subjects, to generate spatially normalised grey matter maps of each subject. • Smooth spatially. • Perform voxel-wise statistics. • Try to interpret the findings in terms of volumetric differences.

4. Segment into different tissue classes Spatially Normalize – with scaling by Jacobian determinant Smooth Spatially Mass-univariate statistical testing Inference via Random Field Theory

5. Smoothing Each voxel after smoothing effectively becomes the result of applying a weighted region of interest (ROI). Before convolution Convolved with a Gaussian Convolved with a circle

6. Mis-register Mis-classify Folding Thinning Mis-register Thickening Mis-classify Possible Explanations for Findings

7. Contents • Introduction • Segmentation • Mixture of Gaussians • Bias correction • Warping to match tissue probability maps • DARTEL Registration

8. Tissue Segmentation • Circularity: • Registration is helped by tissue classification or bias correction. • Tissue classification helped by registration and bias correction. • Bias correction is helped by registration and tissue classification. • The solution is to put everything in the same generative model. • A MAP solution is found by repeatedly alternating among classification, bias correction and registration steps. • Should produce “better” results than simple serial applications of each component.

9. c1 y1 m g c2 y2 s2 a a0 c3 y3 b0 b Ca cI yI Cb A Generative Model • A model of how the data may have been generated, which comprises: • Mixture of Gaussians (MOG) • Bias Correction • Non-linear Inter-subjectRegistration

10. Mixture of Gaussians (MOG) • Tissue classification is based on a Mixture of Gaussians model (MOG), which represents the intensity probability density by a number of Gaussian distributions. Frequency Image Intensity

11. Belonging Probabilities Belonging probabilities are assigned by normalising to one.

12. Non-Gaussian Intensity Distributions • Multiple Gaussians per tissue class allow non-Gaussian intensity distributions to be modelled. • E.g. accounting for partial volume effects

13. Modelling a Bias Field Corrected image Corrupted image Bias Field

14. Tissue Probability Maps • Tissue probability maps (TPMs) are used instead of the proportion of voxels in each Gaussian as the prior. ICBM Tissue Probabilistic Atlases. These tissue probability maps are kindly provided by the International Consortium for Brain Mapping, John C. Mazziotta and Arthur W. Toga.

15. Deforming the Tissue Probability Maps • Tissue probability images are deformed so that they can be overlaid on top of the image to segment.

16. Optimisation • The “best” parameters are those that maximise the log-probability. • Optimisation involves finding them. • Begin with starting estimates, and repeatedly change them so that the objective function decreases each time.

17. Steepest Descent Start Optimum Alternate between optimising different groups of parameters

18. Spatially normalised BrainWeb phantoms (T1, T2 and PD) Tissue probability maps of GM and WM Cocosco, Kollokian, Kwan & Evans. “BrainWeb: Online Interface to a 3D MRI Simulated Brain Database”. NeuroImage 5(4):S425 (1997)

19. Contents • Introduction • Segmentation • DARTEL Registration • Scaling and squaring • Optimisation • Warping GM and WM images to their average

20. 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.

21. 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

22. 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)

23. 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

24. Scaling and squaring example

25. Deformations at different times

26. 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: JϕC=(JϕB○ϕA) JϕA

27. Jacobian determinants remain positive (almost)

28. 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.

29. Contents • Introduction • Segmentation • DARTEL Registration • Scaling and squaring • Optimisation • Warping GM and WM images to their average

30. Multinomial Likelihood Term • 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.

31. 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)

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

33. 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

34. 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|

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

36. Highest resolution Full Multi-Grid Lowest resolution

37. 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.

38. 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.

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

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

41. 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.

42. Contents • Introduction • Segmentation • DARTEL Registration • Scaling and squaring • Optimisation • Warping GM and WM images to their average

43. 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

44. ϕ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(μ) • = p(t1|ϕ1,μ) p(ϕ1) p(t2|ϕ2,μ) p(ϕ2)… p(μ) • MAP solution obtainedfor template. • Requires p(μ)

45. Laplacian Smoothness Priors on template 2D 3D

46. 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.

47. ML and MAP templates from 6 subjects Nonlinearly aligned Rigidly aligned ML MAP log

48. 471 Subject Average

49. 471 Subject Average

50. 471 Subject Average