Volumetric Intersubject Registration for fMRI Analysis

1. Volumetric Intersubject RegistrationJohn AshburnerWellcome Department of Imaging Neuroscience,12 Queen Square, London, UK.

2. Intersubject registration for fMRI • Inter-subject averaging • Increase sensitivity with more subjects • Fixed-effects analysis • Extrapolate findings to the population as a whole • Mixed-effects analysis • Standard coordinate system • e.g., Talairach & Tournoux space

3. Typical overview of fMRI analysis Statistical Parametric Map Design matrix fMRI time-series Motion Correction Smoothing General Linear Model Parameter Estimates Spatial Normalisation Anatomical Reference

4. Overview • Part I: General Inter-subject registration • Spatial transformations • Affine • Global nonlinear • Local nonlinear • Objective functions for registration • Likelihood Models • Mean squared difference • Information Theoretic measures • Prior Models • Part II: The Segmentation Method in SPM5

5. Image Registration • Registration - i.e. Optimise the parameters that describe a spatial transformation between the source and reference (template) images • Transformation - i.e. Re-sample according to the determined transformation parameters

6. A Mapping from one image to another Need x, y and z coordinates in one image that correspond to those of another

7. Affine Transforms • Rigid-body transformations are a subset • Parallel lines remain parallel • Operations can be represented by: x’ = m11x + m12y + m13z + m14 y’ = m21x + m22y + m23z + m24 z’ = m31x + m32y + m33z + m34 • Or as matrices: Y=Mx

8. 2D Affine Transforms • Translations by tx and ty • x’ = x + tx • y’ = y + ty • Rotation around the origin by  radians • x’ = cos() x + sin() y • y’ = -sin() x + cos() y • Zooms by sx and sy • x’ = sx x • y’ = sy y • Shear • x’ = x + h y • y’ = y

9. 2D Affine Transforms • Translations by tx and ty • x’ = 1 x + 0 y + tx • y’ = 0 x + 1 y + ty • Rotation around the origin by  radians • x’ = cos() x + sin() y + 0 • y’ = -sin() x + cos() y + 0 • Zooms by sx and sy: • x’ = sx x+ 0 y + 0 • y’ = 0 x+ sy y+ 0 • Shear • x’ = 1 x + h y+ 0 • y’ = 0 x + 1 y+ 0

10. Polynomial Basis Functions As used by Roger Woods’ AIR Software

11. Cosine Transform Basis Functions As used by SPM software

12. SPM Spatial Normalisation • Begin with affine registration • Refine with some non-linear registration Affine registration Non-linear registration

13. Accuracy of Automated Volumetric Inter-subject Registration Hellier et al. Inter subject registration of functional and anatomical data using SPM.MICCAI'02 LNCS 2489 (2002) Hellier et al. Retrospective evaluation of inter-subject brain registration.MIUA (2001)

14. Local Basis Functions • More detailed deformations use lots of basis functions with local support. • Local support means that the basis functions are mostly all zero • Faster computations

15. Simple addition of displacements Notice that there is no longer a one-to-one mapping

16. Generating large one-to-one deformations Y1 Y2 = Y1  Y1 Y3 = Y1  Y2 Y4 = Y1  Y3 Y5 = Y1  Y4 Y6 = Y1  Y5 Y7 = Y1  Y6 Y8 = Y1  Y7 The principle behind the one-to-one mappings of viscous fluid registration

17. Faster to repeatedly square the deformation Y1 Y2 = Y1  Y1 Y4 = Y2  Y2 Y8 = Y4  Y4 Note that this is analogous to computing a matrix exponential (c.f. Lie Groups and exponential mappings)

18. Y16 = Y8  Y8

19. One-to-One Mappings • One-to-one mappings break down beyond a certain scale • The concept of a single “best” mapping may be meaningless at higher resolution Pictures taken from http://www.messybeast.com/freak-face.htm

20. Optimisation • Optimisation involves finding some “best” parameters according to an “objective function”, which is either minimised or maximised • The “objective function” is often related to a probability based on some model Most probable solution (global optimum) Objective function Local optimum Local optimum Value of parameter

21. Bayes Rule • Most registration procedures maximise a joint probability of the deformation (warp) and the images (data). • P(Warp,Data) = P(Warp | Data) x P(Data) = P(Data | Warp) x P(Warp) • In practice, this can be by minimising • -log P(Warp,Data) = -log P(Data | Warp) -log P(Warp) Prior Likelihood

23. Mean Squared Difference Objective Function • Assumes one image is a warped version of the other with Gaussian noise added… • P(fi|t) = (2ps2)-1/2 exp(-(fi-gi(t))2/(2s2)) so • -log P(fi |t) = (fi-gi(t))2 /(2s2) + 1/2 log(2ps2) • Assumes that voxels are independent... • P(f1,f2,…,fN,...) = P(f1) P(f2) … P(fN) so • -log P(f1,f2,…,fN) = ((f1-g1(t))2 + (f2-g2(t))2 +…+ (fN-gN(t))2)/(2s2) + 1/2 N log(2ps2)

24. Information Theoretic Approaches • Used when there is no simple relationship between intensities in one image and those of another

25. Joint Probability Density • Intensities in one image predict those of another. • Joint probability often represented by a histogram.

26. Mutual Information • MI= ab P(a,b) log2 [P(a,b)/( P(a) P(b) )] • Related to entropy: MI = -H(a,b) + H(a) + H(b) • H(a) = -a P(a) log P(a) da • H(a,b) = -a P(a,b) logP(a,b) da

27. More Joint Probabilities 4x256 Joint Histograms 4x256 Joint Histograms

28. Joint Probabilities generated from Tissue Probability Maps Rather than using an image of discrete values, multiple images showing which voxels are in which class can be used. 4x256 Joint Histogram These can be constructed from an average of many subjects

29. Priors enforce “smooth” deformations • Membrane Energy • Bending Energy • Linear Elastic Energy

30. Priors enforce “smooth” deformations • The form of prior determines how the deformations behave in regions with no matching information

31. Overview • Part I: General Inter-subject registration • Part II: The Segmentation Method in SPM5 • Modelling intensities by a Mixture of Gaussians • Bias correction • Tissue Probability Maps to assist the segmentation • Warping the tissue probability maps to match the image

32. Traditional View of Pre-processing • Brain image processing is often thought of as a pipeline procedure. • One tool applied before another etc... • For example… Original Image Skull Strip Non-uniformity Correct Extract Brain Surfaces Classify Brain Tissues

33. Segmentation in SPM5 • Uses a generative model, which involves: • Mixture of Gaussians (MOG) • Bias Correction Component • Warping (Non-linear Registration) Component c1 y1 m g c2 y2 s2 a a0 c3 y3 b0 b Ca cI yI Cb Ashburner & Friston.Unified Segmentation. NeuroImage 26:839-851 (2005).

34. Gaussian Probability Density • If intensities are assumed to be Gaussian of mean mk and variance s2k, then the probability of a value yi is:

35. Non-Gaussian Probability Distribution • A non-Gaussian probability density function can be modelled by a Mixture of Gaussians (MOG): Mixing proportion - positive and sums to one

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

37. Mixing Proportions • The mixing proportion gk represents the prior probability of a voxel being drawn from class k - irrespective of its intensity. • So:

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

39. Probability of Whole Dataset • If the voxels are assumed to be independent, then the probability of the whole image is the product of the probabilities of each voxel: • A maximum-likelihood solution can be found by minimising the negative log-probability:

40. Modelling a Bias Field • A bias field is included, such that the required scaling at voxel i, parameterised by b, is ri(b). • Replace the means by mk/ri(b) • Replace the variances by (sk/ri(b))2

41. c1 y1 m g c2 y2 s2 a a0 c3 y3 b0 b Ca cI yI Cb Modelling a Bias Field • After rearranging... y yr(b) r(b)

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

43. c1 y1 m g c2 y2 s2 a a0 c3 y3 b0 b Ca cI yI Cb “Mixing Proportions” • Tissue probability maps for each class are included. • The probability of obtaining class k at voxel i, given weights g is then:

44. c1 y1 m g c2 y2 s2 a a0 c3 y3 b0 b Ca cI yI Cb Deforming the Tissue Probability Maps • Tissue probability images are deformed according to parameters a. • The probability of obtaining class k at voxel i is then:

45. The Extended Model • By combining the modified P(ci=k|q) and P(yi|ci=k,q), the overall objective function (E) becomes: The Objective Function

46. Optimisation • The “best” parameters are those that minimise this objective function. • Optimisation involves finding them. • Begin with starting estimates, and repeatedly change them so that the objective function decreases each time.

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

48. Schematic of optimisation Repeat until convergence… Hold g, m, s2 and a constant, and minimise E w.r.t. b - Levenberg-Marquardt strategy, using dE/db and d2E/db2 Hold g, m, s2 and b constant, and minimise E w.r.t. a - Levenberg-Marquardt strategy, using dE/da and d2E/da2 Hold a and b constant, and minimise E w.r.t. g, m and s2 -Use an Expectation Maximisation (EM) strategy. end

49. Levenberg-Marquardt Optimisation • LM optimisation is used for nonlinear registration (a) and bias correction (b). • Requires first and second derivatives of the objective function (E). • Parameters a and b are updated by • Increase l to improve stability (at expense of decreasing speed of convergence).

50. Expectation Maximisation is used to update m, s2 and g • For iteration (n), alternate between: • E-step: Estimate belonging probabilities by: • M-step: Set q(n+1) to values that reduce: