Computational Anatomy & Statistical Shape Models

1. Computational Anatomy&Statistical Shape ModelsJohn Ashburnerjohn@fil.ion.ucl.ac.ukFunctional Imaging Lab, 12 Queen Square, London, UK.

2. Why? • The Wellcome Trust is keen that there is a translational component to the work in the FIL. • E.g. develop some potentially useful diagnostic stuff. • For proper generative models of brain shape differences. • More accurate spatial normalisation. • More accurate shape characterisations. • To use these models for proper characterisation of population differences. • These may be multivariate. • Join the mainstream. • How do more established fields of biology compare shapes?

7. NeuroImage Volume 23, Supplement 1, Pages CO2-S299 (2004) Mathematics in Brain Imaging Edited by P.M. Thompson, M.I. Miller, T. Ratnanather, R.A. Poldrack and T.E. Nichols • Mapping cortical change in Alzheimer's disease, brain development, and schizophrenia Paul M. Thompson, Kiralee M. Hayashi, Elizabeth R. Sowell, Nitin Gogtay, Jay N. Giedd, Judith L. Rapoport, Greig I. de Zubicaray, Andrew L. Janke, Stephen E. Rose, James Semple et al. • Computational anatomy: shape, growth, and atrophy comparison via diffeomorphisms Michael I. Miller • Geometric strategies for neuroanatomic analysis from MRI James S. Duncan, Xenophon Papademetris, Jing Yang, Marcel Jackowski, Xiaolan Zeng and Lawrence H. Staib • Variational, geometric, and statistical methods for modeling brain anatomy and function Olivier Faugeras, Geoffray Adde, Guillaume Charpiat, Christophe Chefd'Hotel, Maureen Clerc, Thomas Deneux, Rachid Deriche, Gerardo Hermosillo, Renaud Keriven, Pierre Kornprobst et al. • Computational anatomy and neuropsychiatric disease: probabilistic assessment of variation and statistical inference of group difference, hemispheric asymmetry, and time-dependent change John G. Csernansky, Lei Wang, Sarang C. Joshi, J. Tilak Ratnanather and Michael I. Miller • Sequence-independent segmentation of magnetic resonance images Bruce Fischl, David H. Salat, André J.W. van der Kouwe, Nikos Makris, Florent Ségonne, Brian T. Quinn and Anders M. Dale • Expert knowledge-guided segmentation system for brain MRI Alain Pitiot, Hervé Delingette, Paul M. Thompson and Nicholas Ayache • Surface-based approaches to spatial localization and registration in primate cerebral cortex David C. Van Essen • Cortical surface segmentation and mapping Duygu Tosun, Maryam E. Rettmann, Xiao Han, Xiaodong Tao, Chenyang Xu, Susan M. Resnick, Dzung L. Pham and Jerry L. Prince • Cortical cartography using the discrete conformal approach of circle packings Monica K. Hurdal and Ken Stephenson • A framework to study the cortical folding patterns J.-F. Mangin, D. Rivière, A. Cachia, E. Duchesnay, Y. Cointepas, D. Papadopoulos-Orfanos, P. Scifo, T. Ochiai, F. Brunelle and J. Régis • Geodesic estimation for large deformation anatomical shape averaging and interpolation Brian Avants and James C. Gee • Unbiased diffeomorphic atlas construction for computational anatomy S. Joshi, Brad Davis, Matthieu Jomier and Guido Gerig • Statistics on diffeomorphisms via tangent space representations M. Vaillant, M.I. Miller, L. Younes and A. Trouvé • Soliton dynamics in computational anatomy Darryl D. Holm, J. Tilak Ratnanather, Alain Trouvé and Laurent Younes • Implicit brain imaging Facundo Mémoli, Guillermo Sapiro and Paul Thompson

8. Computational anatomy: shape, growth, and atrophy comparison via diffeomorphismsMichael I. Miller

9. ? ? ? ? Training and Classifying Control Training Data Patient Training Data

10. ? ? ? ? Classifying Controls Patients y=f(aTx+b)

11. Support Vector Classifier

12. Support Vector Classifier (SVC) a is a weighted linear combination of the support vectors Support Vector Support Vector Support Vector

13. Some Equations • Linear classification is by y = f(aTx + b) • where a is a weighting vector, x is the test data, b is an offset, and f(.) is a thresholding operation • a is a linear combination of SVs a = Si wixi • So y = f(Si wixiTx + b)

14. Going Nonlinear • Nonlinear classification is by y = f(Si wi(xi,x)) • where (xi,x) is some function of xiand x. • e.g. RBF classification (xi,x) = exp(-||xi-x||2/(2s2)) • Requires a matrix of distance measures (metrics) between each pair of images.

15. Nonlinear SVC

16. What is a Metric? A B • Positive • Dist(A,B) ≥ 0 • Dist(A,A) = 0 • Symmetric • Dist(A,B) = Dist(B,A) • Satisfy triangle inequality • Dist(A,B)+Dist(B,C) ≥ Dist(A,C) C

17. Concise representations • Information reduction/compression • Most parsimonious representation - best generalisation • Occam’s Razor • Registration compresses data • signal is partitioned into • deformations • residuals

18. The Small deformation setting • Most “nonlinear” registration is done in the small-deformation setting. • Involves adding a smooth displacement field to an identity transform: y = x + u • No one-to-one constraint • Inverse from: x = y - u • Can be a poor approximation to the real inverse • Adding and subtracting displacements doesn’t work properly • Smoothing and averaging displacements doesn’t work properly • Not the most parsimonious model • Unrealistic generative model

19. Small deformation: displacements are linear within the Eulerian framework

20. Small-deformation setting

21. Forward transform Small def. approx. to backward transform Backward transform Small def. approx. to forward transform

22. Illustrating some concepts with rotations • Consider a 2D rotation y=Rx, where • This can be formulated as the solution of a differential equation at time t=1 • x1(t) = x2(t) • x2(t) = -x1(t) • or • x(t) = Ax(t), where

23. Flow field for rigid 2D rotation

24. Exponentials • The solution can be obtained by • The exponential is defined as: • There are many ways of computing R from A, but one of the easiest is by scaling and squaring

25. Averaging rotations • It makes no sense to average the rotation matrices themselves. • The result may not be a rotation • The elements of a rotationmatrix lie on a manifold. • Average by minimising thesum of squared distancestangential to the manifold • Distance derived fromvelocity (distance travelled inunit time) • Shortest distances are geodesics,which require constant velocity r12 r11

26. Groups • 2D rotation matrices form a Liegroup under multiplication (SO2). • Group Requirements • Composition of group members is another group member • The members have inverses • There is an identity member • The composition operations are associative • Lie Group requirements • Continuous and differentiable manifold

27. 3D Rotations • 3D rotations defined by • These do not commute (R1R2≠R2R1) • Similarly exp(A1) exp(A2) ≠exp(A2) exp(A1) • Both differ from exp(A1+A2) • Makes life more difficult • Iterative schemes needed for averaging etc.

28. Lie Algebra • A would be known as the Lie algebra of the rotation matrix. • The amount of non-commutativity is measured by the Lie bracket • Results fromcurvature of themanifold • [A,B] = AB-BA eB e2(AB-BA) eA eA eB

29. How much rotation is in a rotation matrix? • Given two rotation matrices, R and S. The relative difference between the rotations can be found by computing C = log (R-1S) and then computing the RMS of C. (12+ 22 + 32)1/2

30. Cartan decomposition • A matrix can be decomposed into A = (A+AT)/2 + (A-AT)/2 • (A+AT)/2 is symmetric • Encodes zooms and shears • (A-AT)/2 is skew symmetric • Encodes rigid rotations • A can be converted into a column vector (a) • There is a matrix L, such that La gives the elements of (A-AT)/2. • The square of the rotation angle is then given by (La)T(La) = aT(LTL)a • This excludes any zooming and shearing from the measure • Similarly – and more usefully - the amount of zooming and shearing can be computed in a way that is independent of the rotations.

31. Nonlinear Registration

32. Mapping

33. Flow field for nonlinear deformation

34. … and the resulting deformation

35. A diffeomorphism and its inverse

36. Diffeomorphisms have curved trajectories (variable velocity) if followed in the Eulerian reference frame (fixed). If followed within the Lagrangian frame (moves over time), they appear to have constant velocity.

37. Partial Differential Equations Model one image as it deforms to match another. x(t) = u(x(t)) x(1) = eu (x(0))

38. Matrix representations of diffeomorphisms x(1) = eUx(0) x(0) = e-Ux(1) For large k eU ≈ (I+U/k)k

39. Compositions Large deformations generated from compositions of small deformations S1 = S1/8oS1/8oS1/8oS1/8oS1/8oS1/8oS1/8oS1/8 Recursive formulation S1 = S1/2oS1/2,S1/2 = S1/4oS1/4,S1/4 = S1/8oS1/8 Small deformation approximation S1/8 ≈ I + U/8

40. The shape metric • Don’t use the straight distance (i.e. √uTu) • Distance = √uTLTLu • What’s the best form of L? • Membrane Energy • Bending Energy • Linear Elastic Energy

41. LTL for “membrane energy”

42. LTLu for “membrane energy” is generated by convolving with

43. LTL for “bending energy”

44. LTLu for “bending energy” is generated by convolving with

45. Registration with different models

46. Consistent registration Register to a mean shaped image A B A B µ C C Totally impractical for lots of scans Problem: How can the distance between e.g. A and B be computed? Inverse exponentiating is iterative and slow.

47. Baker-Campbell-Hausdorff series • Exp-1(Exp(A)Exp(B))= A+B +[A,B]/2 +[A,[A,B]]/12-[B,[A,B]]/12 -[B,[A,[A,B]]]/48-[A,[B,[A,B]]]/48 + … • Where [A,B] is the Lie bracket applied to flow fields Sometimes unstable. Looks like proper nonlinear methods would be impractical. eB e2(AB-BA) eA eA eB

48. Alternative strategy • Assume the manifold is locally flat around some point (the template image. • The results depend on the point on the manifold that is chosen. • Ideally use a mean shape as the template. • Best approximation.

49. Controls Patients y=f(aTx+b) Visualisation • The results of multivariate analyses are difficult to visualise • For linear classifiers, it can be done by “caricaturing” the difference • E.g. for the separation of two groups, it would be possible to show two exaggerated versions of the mean image.