Computational anatomy

Computational anatomy is an interdisciplinary field of biology focused on quantitative investigation and modelling of anatomical shapes variability.^[1][2] It involves the development and application of mathematical, statistical and data-analytical methods for modelling and simulation of biological structures.

The field is broadly defined and includes foundations in anatomy, applied mathematics and pure mathematics, machine learning, computational mechanics, computational science, biological imaging, neuroscience, physics, probability, and statistics; it also has strong connections with fluid mechanics and geometric mechanics. Additionally, it complements newer, interdisciplinary fields like bioinformatics and neuroinformatics in the sense that its interpretation uses metadata derived from the original sensor imaging modalities (of which magnetic resonance imaging is one example). It focuses on the anatomical structures being imaged, rather than the medical imaging devices. It is similar in spirit to the history of computational linguistics, a discipline that focuses on the linguistic structures rather than the sensor acting as the transmission and communication media.

In computational anatomy, the diffeomorphism group is used to study different coordinate systems via coordinate transformations as generated via the Lagrangian and Eulerian velocities of flow in ${\displaystyle {\mathbb {R} }^{3}}$. The flows between coordinates in computational anatomy are constrained to be geodesic flows satisfying the principle of least action for the Kinetic energy of the flow. The kinetic energy is defined through a Sobolev smoothness norm with strictly more than two generalized, square-integrable derivatives for each component of the flow velocity, which guarantees that the flows in ${\displaystyle \mathbb {R} ^{3}}$ are diffeomorphisms.^[3] It also implies that the diffeomorphic shape momentum taken pointwise satisfying the Euler-Lagrange equation for geodesics is determined by its neighbors through spatial derivatives on the velocity field. This separates the discipline from the case of incompressible fluids^[4] for which momentum is a pointwise function of velocity. Computational anatomy intersects the study of Riemannian manifolds and nonlinear global analysis, where groups of diffeomorphisms are the central focus. Emerging high-dimensional theories of shape^[5] are central to many studies in computational anatomy, as are questions emerging from the fledgling field of shape statistics. The metric structures in computational anatomy are related in spirit to morphometrics, with the distinction that Computational anatomy focuses on an infinite-dimensional space of coordinate systems transformed by a diffeomorphism, hence the central use of the terminology diffeomorphometry, the metric space study of coordinate systems via diffeomorphisms.

Genesis

At computational anatomy's heart is the comparison of shape by recognizing in one shape the other. This connects it to D'Arcy Wentworth Thompson's developments On Growth and Form which has led to scientific explanations of morphogenesis, the process by which patterns are formed in biology. Albrecht Durer's Four Books on Human Proportion were arguably the earliest works on computational anatomy.^[6][7][8] The efforts of Noam Chomsky in his pioneering of computational linguistics inspired the original formulation of computational anatomy as a generative model of shape and form from exemplars acted upon via transformations.^[9]

Due to the availability of dense 3D measurements via technologies such as magnetic resonance imaging (MRI), computational anatomy has emerged as a subfield of medical imaging and bioengineering for extracting anatomical coordinate systems at the morphome scale in 3D. The spirit of this discipline shares strong overlap with areas such as computer vision and kinematics of rigid bodies, where objects are studied by analysing the groups responsible for the movement in question. Computational anatomy departs from computer vision with its focus on rigid motions, as the infinite-dimensional diffeomorphism group is central to the analysis of Biological shapes. It is a branch of the image analysis and pattern theory school at Brown University^[10] pioneered by Ulf Grenander. In Grenander's general metric pattern theory, making spaces of patterns into a metric space is one of the fundamental operations since being able to cluster and recognize anatomical configurations often requires a metric of close and far between shapes. The diffeomorphometry metric^[11] of computational anatomy measures how far two diffeomorphic changes of coordinates are from each other, which in turn induces a metric on the shapes and images indexed to them. The models of metric pattern theory,^[12][13] in particular group action on the orbit of shapes and forms is a central tool to the formal definitions in computational anatomy.

History

Computational anatomy is the study of shape and form at the morphome or gross anatomy millimeter, or morphology scale, focusing on the study of sub-manifolds of ${\displaystyle {\mathbb {R} }^{3},}$ points, curves surfaces and subvolumes of human anatomy. An early modern computational neuro-anatomist was David Van Essen^[14] performing some of the early physical unfoldings of the human brain based on printing of a human cortex and cutting. Jean Talairach's publication of Talairach coordinates is an important milestone at the morphome scale demonstrating the fundamental basis of local coordinate systems in studying neuroanatomy and therefore the clear link to charts of differential geometry. Concurrently, virtual mapping in computational anatomy across high resolution dense image coordinates was already happening in Ruzena Bajcy's^[15] and Fred Bookstein's^[16] earliest developments based on computed axial tomography and magnetic resonance imagery. The earliest introduction of the use of flows of diffeomorphisms for transformation of coordinate systems in image analysis and medical imaging was by Christensen, Joshi, Miller, and Rabbitt.^[17][18][19]

The first formalization of computational anatomy as an orbit of exemplar templates under diffeomorphism group action was in the original lecture given by Grenander and Miller with that title in May 1997 at the 50th Anniversary of the Division of Applied Mathematics at Brown University,^[20] and subsequent publication.^[9] This was the basis for the strong departure from much of the previous work on advanced methods for spatial normalization and image registration which were historically built on notions of addition and basis expansion. The structure preserving transformations central to the modern field of Computational Anatomy, homeomorphisms and diffeomorphisms carry smooth submanifolds smoothly. They are generated via Lagrangian and Eulerian flows which satisfy a law of composition of functions forming the group property, but are not additive.

The original model of computational anatomy was as the triple, ${\displaystyle ({\mathcal {G}},{\mathcal {M}},{\mathcal {P}})\ ,}$ the group ${\displaystyle g\in {\mathcal {G}}}$, the orbit of shapes and forms ${\displaystyle m\in {\mathcal {M}}}$, and the probability laws ${\displaystyle P}$ which encode the variations of the objects in the orbit. The template or collection of templates are elements in the orbit ${\displaystyle m_{\mathrm {temp} }\in {\mathcal {M}}}$ of shapes.

The Lagrangian and Hamiltonian formulations of the equations of motion of computational anatomy took off post 1997 with several pivotal meetings including the 1997 Luminy meeting^[21] organized by the Azencott^[22] school at Ecole-Normale Cachan on the "Mathematics of Shape Recognition" and the 1998 Trimestre at Institute Henri Poincaré organized by David Mumford "Questions Mathématiques en Traitement du Signal et de l'Image" which catalyzed the Hopkins-Brown-ENS Cachan groups and subsequent developments and connections of computational anatomy to developments in global analysis.

The developments in computational anatomy included the establishment of the Sobelev smoothness conditions on the diffeomorphometry metric to insure existence of solutions of variational problems in the space of diffeomorphisms,^[23][24] the derivation of the Euler-Lagrange equations characterizing geodesics through the group and associated conservation laws,^[25][26][27] the demonstration of the metric properties of the right invariant metric,^[28] the demonstration that the Euler-Lagrange equations have a well-posed initial value problem with unique solutions for all time,^[29] and with the first results on sectional curvatures for the diffeomorphometry metric in landmarked spaces.^[30] Following the Los Alamos meeting in 2002,^[31] Joshi's^[32] original large deformation singular Landmark solutions in computational anatomy were connected to peaked solitons or peakons^[33] as solutions for the Camassa-Holm equation. Subsequently, connections were made between computational anatomy's Euler-Lagrange equations for momentum densities for the right-invariant metric satisfying Sobolev smoothness to Vladimir Arnold's^[4] characterization of the Euler equation for incompressible flows as describing geodesics in the group of volume preserving diffeomorphisms.^[34][35] The first algorithms, generally termed LDDMM for large deformation diffeomorphic mapping for computing connections between landmarks in volumes^[32][36][37] and spherical manifolds,^[38] curves,^[39] currents and surfaces,^[40][41][42] volumes,^[43] tensors,^[44] varifolds,^[45] and time-series^[46][47][48] have followed.

These contributions of computational anatomy to the global analysis associated to the infinite dimensional manifolds of subgroups of the diffeomorphism group is far from trivial. The original idea of doing differential geometry, curvature and geodesics on infinite dimensional manifolds goes back to Bernhard Riemann's Habilitation (Ueber die Hypothesen, welche der Geometrie zu Grunde liegen^[49][50]); the key modern book laying the foundations of such ideas in global analysis are from Michor.^[51]

The applications within medical imaging of computational anatomy continued to flourish after two organized meetings at the Institute for Pure and Applied Mathematics conferences^[52][53] at University of California, Los Angeles. Computational anatomy has been useful in creating accurate models of the atrophy of the human brain at the morphome scale, as well as Cardiac templates,^[54] as well as in modeling biological systems.^[55] Since the late 1990s, computational anatomy has become an important part of developing emerging technologies for the field of medical imaging. Digital atlases are a fundamental part of modern Medical-school education^[56][57] and in neuroimaging research at the morphome scale.^[58][59] Atlas based methods and virtual textbooks^[60] which accommodate variations as in deformable templates are at the center of many neuro-image analysis platforms including Freesurfer,^[61] FSL,^[62] MRIStudio,^[63] SPM.^[64] Diffeomorphic registration,^[18] introduced in the 1990s, is now an important player with existing codes bases organized around ANTS,^[65] DARTEL,^[66] DEMONS,^[67] LDDMM,^[68] StationaryLDDMM,^[69] FastLDDMM,^[70] are examples of actively used computational codes for constructing correspondences between coordinate systems based on sparse features and dense images. Voxel-based morphometry is an important technology built on many of these principles.

The deformable template orbit model of computational anatomy

The model of human anatomy is a deformable template, an orbit of exemplars under group action. Deformable template models have been central to Grenander's metric pattern theory, accounting for typicality via templates, and accounting for variability via transformation of the template. An orbit under group action as the representation of the deformable template is a classic formulation from differential geometry. The space of shapes are denoted ${\displaystyle m\in {\mathcal {M}}}$, with the group ${\displaystyle ({\mathcal {G}},\circ )}$ with law of composition ${\displaystyle \circ }$; the action of the group on shapes is denoted ${\displaystyle g\cdot m}$, where the action of the group ${\displaystyle g\cdot m\in {\mathcal {M}},m\in {\mathcal {M}}}$ is defined to satisfy

${\displaystyle (g\circ g^{\prime })\cdot m=g\cdot (g^{\prime }\cdot m)\in {\mathcal {M}}.}$

The orbit ${\displaystyle {\mathcal {M}}}$ of the template becomes the space of all shapes, ${\displaystyle {\mathcal {M}}\doteq \{m=g\cdot m_{\mathrm {temp} },g\in {\mathcal {G}}\}}$, being homogenous under the action of the elements of ${\displaystyle {\mathcal {G}}}$.

Figure depicting three medial temporal lobe structures amygdala, entorhinal cortex and hippocampus with fiducial landmarks depicted as well embedded in the MRI background.

The orbit model of computational anatomy is an abstract algebra - to be compared to linear algebra - since the groups act nonlinearly on the shapes. This is a generalization of the classical models of linear algebra, in which the set of finite dimensional ${\displaystyle {\mathbb {R} }^{n}}$ vectors are replaced by the finite-dimensional anatomical submanifolds (points, curves, surfaces and volumes) and images of them, and the ${\displaystyle n\times n}$ matrices of linear algebra are replaced by coordinate transformations based on linear and affine groups and the more general high-dimensional diffeomorphism groups.

Shapes and forms

The central objects are shapes or forms in computational anatomy, one set of examples being the 0,1,2,3-dimensional submanifolds of ${\displaystyle {\mathbb {R} }^{3}}$, a second set of examples being images generated via medical imaging such as via magnetic resonance imaging (MRI) and functional magnetic resonance imaging.

Triangulated mesh surfaces depicting subcortical structures amygdala, hippocampus, thalamus, caudate, putamen, ventricles.The shapes are denoted ${\displaystyle m(u),u\in U\subset {\mathbb {R} }^{1}\rightarrow {\mathbb {R} }^{2}}$ represented as triangulated meshes.

The 0-dimensional manifolds are landmarks or fiducial points; 1-dimensional manifolds are curves such as sulcul and gyral curves in the brain; 2-dimensional manifolds correspond to boundaries of substructures in anatomy such as the subcortical structures of the midbrain or the gyral surface of the neocortex; subvolumes correspond to subregions of the human body, the heart, the thalamus, the kidney.

The landmarks ${\displaystyle X\doteq \{x_{1},\dots ,x_{n}\}\subset {\mathbb {R} }^{3}\in {\mathcal {M}}}$ are a collections of points with no other structure, delineating important fiducials within human shape and form (see associated landmarked image). The sub-manifold shapes such as surfaces ${\displaystyle X\subset {\mathbb {R} }^{3}\in {\mathcal {M}}}$ are collections of points modeled as parametrized by a local chart or immersion ${\displaystyle m:U\subset {\mathbb {R} }^{1,2}\rightarrow {\mathbb {R} }^{3}}$, ${\displaystyle m(u),u\in U}$ (see Figure showing shapes as mesh surfaces). The images such as MR images or DTI images ${\displaystyle I\in {\mathcal {M}}}$, and are dense functions ${\displaystyle I(x),x\in X\subset {\mathbb {R} }^{1,2,3}}$ are scalars, vectors, and matrices (see Figure showing scalar image).

Groups and group actions

Showing an MRI section through a 3D brain representing a scalar image ${\displaystyle I(x),x\in {\mathbb {R} }^{2}}$ based on T1-weighting.

Groups and group actions are familiar to the Engineering community with the universal popularization and standardization of linear algebra as a basic model for analyzing signals and systems in mechanical engineering, electrical engineering and applied mathematics. In linear algebra the matrix groups (matrices with inverses) are the central structure, with group action defined by the usual definition of ${\displaystyle A}$ as an ${\displaystyle n\times n}$ matrix, acting on ${\displaystyle x\in {\mathbb {R} }^{n}}$ as ${\displaystyle n\times 1}$ vectors; the orbit in linear-algebra is the set of ${\displaystyle n}$-vectors given by ${\displaystyle y=A\cdot x\in {\mathbb {R} }^{n}}$, which is a group action of the matrices through the orbit of ${\displaystyle {\mathbb {R} }^{n}}$.

The central group in computational anatomy defined on volumes in ${\displaystyle {\mathbb {R} }^{3}}$ are the diffeomorphisms ${\displaystyle {\mathcal {G}}\doteq Diff}$ which are mappings with 3-components ${\displaystyle \phi (\cdot )=(\phi _{1}(\cdot ),\phi _{2}(\cdot ),\phi _{3}(\cdot ))}$, law of composition of functions ${\displaystyle \phi \circ \phi ^{\prime }(\cdot )\doteq \phi (\phi ^{\prime }(\cdot ))}$, with inverse ${\displaystyle \phi \circ \phi ^{-1}(\cdot )=\phi (\phi ^{-1}(\cdot ))=id}$.

Most popular are scalar images, ${\displaystyle I(x),x\in {\mathbb {R} }^{3}}$, with action on the right via the inverse.

${\displaystyle \phi \cdot I(x)=I\circ \phi ^{-1}(x),x\in {\mathbb {R} }^{3}}$.

For sub-manifolds ${\displaystyle X\subset {\mathbb {R} }^{3}\in {\mathcal {M}}}$, parametrized by a chart or immersion ${\displaystyle m(u),u\in U}$, the diffeomorphic action the flow of the position

${\displaystyle \phi \cdot m(u)\doteq \phi \circ m(u),u\in U}$.

Several group actions in computational anatomy have been defined.

Lagrangian and Eulerian flows for generating diffeomorphisms

For the study of rigid body kinematics, the low-dimensional matrix Lie groups have been the central focus. The matrix groups are low-dimensional mappings, which are diffeomorphisms that provide one-to-one correspondences between coordinate systems, with a smooth inverse. The matrix group of rotations and scales can be generated via a closed form finite-dimensional matrices which are solution of simple ordinary differential equations with solutions given by the matrix exponential.

For the study of deformable shape in computational anatomy, a more general diffeomorphism group has been the group of choice, which is the infinite dimensional analogue. The high-dimensional differeomorphism groups used in Computational Anatomy are generated via smooth flows ${\displaystyle \phi _{t},t\in [0,1]}$ which satisfy the Lagrangian and Eulerian specification of the flow fields as first introduced in.,^[17][19][71] satisfying the ordinary differential equation:

Showing the Lagrangian flow of coordinates ${\displaystyle x\in X}$ with associated vector fields ${\displaystyle v_{t},t\in [0,1]}$ satisfying ordinary differential equation ${\displaystyle {\dot {\phi }}_{t}=v_{t}(\phi _{t}),\phi _{0}=id}$.

${\displaystyle {\frac {d}{dt}}\phi _{t}=v_{t}\circ \phi _{t},\ \phi _{0}=id\$ ;}

(Lagrangian flow)

with ${\displaystyle v\doteq (v_{1},v_{2},v_{3})}$ the vector fields on ${\displaystyle {\mathbb {R} }^{3}}$ termed the Eulerian velocity of the particles at position ${\displaystyle \phi }$ of the flow. The vector fields are functions in a function space, modelled as a smooth Hilbert space of high-dimension, with the Jacobian of the flow ${\displaystyle \ D\phi \doteq ({\frac {\partial \phi _{i}}{\partial x_{j}}})}$ a high-dimensional field in a function space as well, rather than a low-dimensional matrix as in the matrix groups. Flows were first introduced^[72][73] for large deformations in image matching; ${\displaystyle {\dot {\phi }}_{t}(x)}$ is the instantaneous velocity of particle ${\displaystyle x}$ at time ${\displaystyle t}$ .

The inverse ${\displaystyle \phi _{t}^{-1},t\in [0,1]}$ required for the group is defined on the Eulerian vector-field with advective inverse flow

${\displaystyle {\frac {d}{dt}}\phi _{t}^{-1}=-(D\phi _{t}^{-1})v_{t},\ \phi _{0}^{-1}=id\ .}$

(Inverse Transport flow)

The diffeomorphism group of computational anatomy

The group of diffeomorphisms is very big. To ensure smooth flows of diffeomorphisms avoiding shock-like solutions for the inverse, the vector fields must be at least 1-time continuously differentiable in space.^[74][75] For diffeomorphisms on ${\displaystyle {\mathbb {R} }^{3}}$, vector fields are modelled as elements of the Hilbert space ${\displaystyle (V,\|\cdot \|_{V})}$ using the Sobolev embedding theorems so that each element has strictly greater than 2 generalized square-integrable spatial derivatives (thus ${\displaystyle v_{i}\in H_{0}^{3},i=1,2,3,}$ is sufficient), yielding 1-time continuously differentiable functions.^[74][75]

The diffeomorphism group are flows with vector fields absolutely integrable in Sobolev norm:

${\displaystyle Diff_{V}\doteq \{\varphi =\phi _{1}:{\dot {\phi }}_{t}=v_{t}\circ \phi _{t},\phi _{0}=id,\int _{0}^{1}\|v_{t}\|_{V}dt<\infty \}\ ,}$

(Diffeomorphism Group)

where ${\displaystyle \|v\|_{V}^{2}\doteq \int _{X}Av\cdot vdx,\ v\in V\ ,}$ with the linear operator ${\displaystyle A}$ mapping to the dual space ${\displaystyle A:V\mapsto V^{*}}$, with the integral calculated by integration by parts when ${\displaystyle Av\in V^{*}}$ is a generalized function in the dual space.

Diffeomorphometry: The metric space of shapes and forms

The study of metrics on groups of diffeomorphisms and the study of metrics between manifolds and surfaces has been an area of significant investigation.^{[28][76][77][78][79][80]} The diffeomorphometry metric measures how close and far two shapes or images are from each other; the metric length is the shortest length of the flow which carries one coordinate system into the other.

Oftentimes, the familiar Euclidean metric is not directly applicable because the patterns of shapes and images don't form a vector space. In the Riemannian orbit model of computational anatomy, diffeomorphisms acting on the forms ${\displaystyle \phi \cdot m\in {\mathcal {M}},\phi \in Diff_{V},m\in {\mathcal {M}}}$ don't act linearly. There are many ways to define metrics, and for the sets associated to shapes the Hausdorff metric is another. The method we use to induce the Riemannian metric is used to induce the metric on the orbit of shapes by defining it in terms of the metric length between diffeomorphic coordinate system transformations of the flows. Measuring the lengths of the geodesic flow between coordinates systems in the orbit of shapes is called diffeomorphometry.

The right-invariant metric on diffeomorphisms

Define the distance on the group of diffeomorphisms

${\displaystyle d_{Diff_{V}}(\psi ,\varphi )=\inf _{v_{t}}\left(\int _{0}^{1}\int _{X}Av_{t}\cdot v_{t}dx\ dt:\phi _{0}=\psi ,\phi _{1}=\varphi ,{\dot {\phi }}_{t}=v_{t}\circ \phi _{t}\right)^{1/2}\$ ;}

(metric-diffeomorphisms)

this is the right-invariant metric of diffeomorphometry,^[11][28] invariant to reparameterization of space since for all ${\displaystyle \phi \in Diff_{V}}$,

${\displaystyle d_{Diff_{V}}(\psi ,\varphi )=d_{Diff_{V}}(\psi \circ \phi ,\varphi \circ \phi )}$.

The metric on shapes and forms

The distance on shapes and forms,^[81]${\displaystyle d_{\mathcal {M}}:{\mathcal {M}}\times {\mathcal {M}}\rightarrow \mathbb {R} ^{+}}$,

${\displaystyle d_{\mathcal {M}}(m,n)=\inf _{\phi \in \operatorname {Diff} _{V}:\phi \cdot m=n}d_{\operatorname {Diff} _{V}}(id,\phi )\$ ;}

(metric-shapes-forms)

the images^[28] are denoted with the orbit as ${\displaystyle I\in {\mathcal {I}}}$ and metric ${\displaystyle ,d_{\mathcal {I}}}$.

The action integral for Hamilton's principle on diffeomorphic flows

In classical mechanics the evolution of physical systems is described by solutions to the Euler–Lagrange equations associated to the Least-action principle of Hamilton. This is a standard way, for example of obtaining Newton's laws of motion of free particles. More generally, the Euler-Lagrange equations can be derived for systems of generalized coordinates. The Euler-Lagrange equation in computational anatomy describes the geodesic shortest path flows between coordinate systems of the diffeomorphism metric. In computational anatomy the generalized coordinates are the flow of the diffeomorphism and its Lagrangian velocity ${\displaystyle \phi ,{\dot {\phi }}}$, the two related via the Eulerian velocity ${\displaystyle v\doteq {\dot {\phi }}\circ \phi ^{-1}}$. Hamilton's principle for generating the Euler-Lagrange equation requires the action integral on the Lagrangian given by

${\displaystyle J(\phi )\doteq \int _{0}^{1}L(\phi _{t},{\dot {\phi }}_{t})dt\$ ;}

(Hamiltonian-Integrated-Lagrangian)

the Lagrangian is given by the kinetic energy:

${\displaystyle L(\phi _{t},{\dot {\phi }}_{t})\doteq {\frac {1}{2}}\int _{X}A({\dot {\phi }}_{t}\circ \phi _{t}^{-1})\cdot ({\dot {\phi }}_{t}\circ \phi _{t}^{-1})dx={\frac {1}{2}}\int _{X}Av_{t}\cdot v_{t}\ dx\ .}$

(Lagrangian-Kinetic-Energy)

Diffeomorphic or Eulerian shape momentum

In computational anatomy, ${\displaystyle Av}$ was first called the Eulerian or diffeomorphic shape momentum^[82] since when integrated against Eulerian velocity ${\displaystyle v}$ gives energy density, and since there is a conservation of diffeomorphic shape momentum which holds. The operator ${\displaystyle A}$ is the generalized moment of inertia or inertial operator.

The Euler–Lagrange equation on shape momentum for geodesics on the group of diffeomorphisms

Classical calculation of the Euler-Lagrange equation from Hamilton's principle requires the perturbation of the Lagrangian on the vector field in the kinetic energy with respect to first order perturbation of the flow. This requires adjustment by the Lie bracket of vector field, given by operator ${\displaystyle ad_{v}:w\in V\mapsto V}$ which involves the Jacobian given by

${\displaystyle ad_{v}[w]\doteq [v,w]\doteq (Dv)w-(Dw)v\in V}$.

Defining the adjoint ${\displaystyle ad_{v}^{*}:V^{*}\rightarrow V^{*},}$ then the first order variation gives the Eulerian shape momentum ${\displaystyle Av\in V^{*}}$ satisfying the generalized equation:

${\displaystyle {\frac {d}{dt}}Av_{t}+ad_{v_{t}}^{*}(Av_{t})=0\ ,\ t\in [0,1]\$ ;}

(EL-General)

meaning for all smooth ${\displaystyle w\in V,}$

${\displaystyle \int _{X}\left({\frac {d}{dt}}Av_{t}+ad_{v_{t}}^{*}(Av_{t})\right)\cdot wdx=\int _{X}{\frac {d}{dt}}Av_{t}\cdot wdx+\int _{X}Av_{t}\cdot ((Dv_{t})w-(Dw)v_{t})dx=0.}$

Computational anatomy is the study of the motions of submanifolds, points, curves, surfaces and volumes. Momentum associated to points, curves and surfaces are all singular, implying the momentum is concentrated on subsets of ${\displaystyle {\mathbb {R} }^{3}}$ which are dimension ${\displaystyle \leq 2}$ in Lebesgue measure. In such cases, the energy is still well defined ${\displaystyle (Av_{t}\mid v_{t})}$ since although ${\displaystyle Av_{t}}$ is a generalized function, the vector fields are smooth and the Eulerian momentum is understood via its action on smooth functions. The perfect illustration of this is even when it is a superposition of delta-diracs, the velocity of the coordinates in the entire volume move smoothly. The Euler-Lagrange equation (EL-General) on diffeomorphisms for generalized functions ${\displaystyle Av\in V^{*}}$ was derived in.^[83] In Riemannian Metric and Lie-Bracket Interpretation of the Euler-Lagrange Equation on Geodesics derivations are provided in terms of the adjoint operator and the Lie bracket for the group of diffeomorphisms. It has come to be called EPDiff equation for diffeomorphisms connecting to the Euler-Poincare method having been studied in the context of the inertial operator ${\displaystyle A=identity}$ for incompressible, divergence free, fluids.^[35][84]

Riemannian exponential (geodesic positioning) and Riemannian logarithm (geodesic coordinates)

In medical imaging and computational anatomy, positioning and coordinatizing shapes are fundamental operations; the system for positioning anatomical coordinates and shapes built on the metric and the Euler-Lagrange equation a geodesic positioning system as first explicated in Miller Trouve and Younes.^[11] Solving the geodesic from the initial condition ${\displaystyle v_{0}}$ is termed the Riemannian-exponential, a mapping ${\displaystyle Exp_{\rm {id}}(\cdot ):V\to Diff_{V}}$ at identity to the group.

The Riemannian exponential satisfies ${\displaystyle Exp_{id}(v_{0})=\phi _{1}}$ for initial condition ${\displaystyle {\dot {\phi }}_{0}=v_{0}}$, vector field dynamics ${\displaystyle {\dot {\phi }}_{t}=v_{t}\circ \phi _{t},t\in [0,1]}$,

• for classical equation diffeomorphic shape momentum ${\displaystyle \int _{X}Av_{t}\cdot w\,dx}$, ${\displaystyle Av\in V}$, then

${\displaystyle \ \ \ {\frac {d}{dt}}Av_{t}+(Dv_{t})^{T}Av_{t}+(DAv_{t})v_{t}+(\nabla \cdot v)Av_{t}=0\$ ;}

• for generalized equation, then ${\displaystyle Av\in V^{*}}$,${\displaystyle w\in V}$,

${\displaystyle \ \ \ \int _{X}{\frac {d}{dt}}Av_{t}\cdot wdx+\int _{X}Av_{t}\cdot ((Dv_{t})w-(Dw)v_{t})dx=0.}$

Computing the flow ${\displaystyle v_{0}}$ onto coordinates Riemannian logarithm,^[11][81] mapping ${\displaystyle Log_{\rm {id}}(\cdot ):Diff_{V}\to V}$ at identity from ${\displaystyle \varphi }$ to vector field ${\displaystyle v_{0}\in V}$;

${\displaystyle Log_{id}(\varphi )=v_{0}\ {\text{initial condition of EL geodesic }}{\dot {\phi }}_{0}=v_{0},\phi _{0}=id,\phi _{1}=\varphi \ .}$

Extended to the entire group they become

${\displaystyle \phi =Exp_{\varphi }(v_{0}\circ \varphi )\doteq Exp_{id}(v_{0})\circ \varphi }$ ; ${\displaystyle Log_{\varphi }(\phi )\doteq Log_{id}(\phi \circ \varphi ^{-1})\circ \varphi }$ .

These are inverses of each other for unique solutions of Logarithm; the first is called geodesic positioning, the latter geodesic coordinates (see exponential map, Riemannian geometry for the finite dimensional version).The geodesic metric is a local flattening of the Riemannian coordinate system (see figure).

Showing metric local flattening of coordinatized manifolds of shapes and forms. The local metric is given by the norm of the vector field ${\displaystyle \|v_{0}\|_{V}}$ of the geodesic mapping ${\displaystyle Exp_{id}(v_{0})\cdot m}$

Hamiltonian formulation of computational anatomy

In computational anatomy the diffeomorphisms are used to push the coordinate systems, and the vector fields are used as the control within the anatomical orbit or morphological space. The model is that of a dynamical system, the flow of coordinates ${\displaystyle t\mapsto \phi _{t}\in \operatorname {Diff} _{V}}$ and the control the vector field ${\displaystyle t\mapsto v_{t}\in V}$ related via ${\displaystyle {\dot {\phi }}_{t}=v_{t}\cdot \phi _{t},\phi _{0}=id.}$ The Hamiltonian view ^{[81] [86] [87] [88][89]} reparameterizes the momentum distribution ${\displaystyle Av\in V^{*}}$ in terms of the conjugate momentum or canonical momentum, introduced as a Lagrange multiplier ${\displaystyle p:{\dot {\phi }}\mapsto (p\mid {\dot {\phi }})}$ constraining the Lagrangian velocity ${\displaystyle {\dot {\phi }}_{t}=v_{t}\circ \phi _{t}}$.accordingly:

${\displaystyle H(\phi _{t},p_{t},v_{t})=\int _{X}p_{t}\cdot (v_{t}\circ \phi _{t})dx-{\frac {1}{2}}\int _{X}Av_{t}\cdot v_{t}dx.}$

This function is the extended Hamiltonian. The Pontryagin maximum principle^[81] gives the optimizing vector field which determines the geodesic flow satisfying ${\displaystyle {\dot {\phi }}_{t}=v_{t}\circ \phi _{t},\phi _{0}=id,}$ as well as the reduced Hamiltonian

${\displaystyle H(\phi _{t},p_{t})\doteq \max _{v}H(\phi _{t},p_{t},v)\ .}$

The Lagrange multiplier in its action as a linear form has its own inner product of the canonical momentum acting on the velocity of the flow which is dependent on the shape, e.g. for landmarks a sum, for surfaces a surface integral, and. for volumes it is a volume integral with respect to ${\displaystyle dx}$ on ${\displaystyle {\mathbb {R} }^{3}}$. In all cases the Greens kernels carry weights which are the canonical momentum evolving according to an ordinary differential equation which corresponds to EL but is the geodesic reparameterization in canonical momentum. The optimizing vector field is given by

${\displaystyle v_{t}\doteq \arg max_{v}H(\phi _{t},p_{t},v)}$

with dynamics of canonical momentum reparameterizing the vector field along the geodesic

${\displaystyle {\begin{cases}{\dot {\phi }}_{t}={\frac {\partial H(\phi _{t},p_{t})}{\partial p}}\\{\dot {p}}_{t}=-{\frac {\partial H(\phi _{t},p_{t})}{\partial \phi }}\\\end{cases}}}$

(Hamiltonian-Dynamics)

The metric on geodesic flows of landmarks, surfaces, and volumes within the orbit

In computational anatomy the submanifolds are pointsets, curves, surfaces and subvolumes which are the basic primitives. The geodesic flows between the submanifolds determine the distance, and form the basic measuring and transporting tools of diffeomorphometry. At ${\displaystyle t=0}$ the geodesic has vector field ${\displaystyle v_{0}=Kp_{0}}$ determined by the conjugate momentum and the Green's kernel of the inertial operator defining the Eulerian momentum ${\displaystyle K=A^{-1}}$. The metric distance between coordinate systems connected via the geodesic determined by the induced distance between identity and group element:

${\displaystyle d_{Diff_{V}}(id,\varphi )=\|Log_{id}(\varphi )\|_{V}=\|v_{0}\|_{V}={\sqrt {2H(id,p_{0})}}}$

Conservation laws on diffeomorphic shape momentum for computational anatomy

Given the least-action there is a natural definition of momentum associated to generalized coordinates; the quantity acting against velocity gives energy. The field has studied two forms, the momentum associated to the Eulerian vector field termed Eulerian diffeomorphic shape momentum, and the momentum associated to the initial coordinates or canonical coordinates termed canonical diffeomorphic shape momentum. Each has a conservation law. The conservation of momentum goes hand in hand with the EL-General. In computational anatomy, ${\displaystyle Av}$ is the Eulerian Momentum since when integrated against Eulerian velocity ${\displaystyle v}$ gives energy density; operator ${\displaystyle A}$ the generalized moment of inertia or inertial operator which acting on the Eulerian velocity gives momentum which is conserved along the geodesic:

${\displaystyle {\begin{matrix}{\text{Eulerian}}&\ \ \ \ {\frac {d}{dt}}\int _{X}Av_{t}\cdot ((D\phi _{t})w)\circ \phi _{t}^{-1})dx=0\ ,\ t\in [0,1].\\&\\{\text{Canonical}}&\ \ \ \ \ \ \ \ \ \ \ {\frac {d}{dt}}\int _{X}p_{t}\cdot ((D\phi _{t})w)dx=0\ ,\ t\in [0,1]\ {\text{ for all}}\ w\in V\ .\end{matrix}}}$

(Euler-Conservation-Constant-Energy)

Conservation of Eulerian shape momentum was shown in^[94] and follows from EL-General; conservation of canonical momentum was shown in^[81]

Geodesic interpolation of information between coordinate systems via variational problems

Construction of diffeomorphic correspondences between shapes calculates the initial vector field coordinates ${\displaystyle v_{0}\in V}$ and associated weights on the Greens kernels ${\displaystyle p_{0}}$. These initial coordinates are determined by matching of shapes, called Large Deformation Diffeomorphic Metric Mapping (LDDMM). LDDMM has been solved for landmarks with and without correspondence^{[32][95][96][97][98]} and for dense image matchings.^[99][100] curves,^[101] surfaces,^[41][102] dense vector^[103] and tensor^[104] imagery, and varifolds removing orientation.^[105] LDDMM calculates geodesic flows of the EL-General onto target coordinates, adding to the action integral ${\displaystyle {\frac {1}{2}}\int _{0}^{1}\int _{X}Av_{t}\cdot v_{t}dxdt}$ an endpoint matching condition ${\displaystyle E:\phi _{1}\rightarrow R^{+}}$ measuring the correspondence of elements in the orbit under coordinate system transformation. Existence of solutions were examined for image matching.^[24] The solution of the variational problem satisfies the EL-General for ${\displaystyle t\in [0,1)}$ with boundary condition.

Dense image matching in computational anatomy

Dense image matching has a long history now with the earliest efforts^[106][107] exploiting a small deformation framework. Large deformations began in the early 1990s,^[18][19] with the first existence to solutions to the variational problem for flows of diffeomorphisms for dense image matching established in.^[24] Beg solved via one of the earliest LDDMM algorithms based on solving the variational matching with endpoint defined by the dense imagery with respect to the vector fields, taking variations with respect to the vector fields.^[99] Another solution for dense image matching reparameterizes the optimization problem in terms of the state ${\displaystyle q_{t}\doteq I\circ \phi _{t}^{-1},q_{0}=I}$ giving the solution in terms of the infinitesimal action defined by the advection equation.^{[11][27][100]}

Diffusion tensor image matching in computational anatomy

Further information: Group actions in computational anatomy and LDDMM

Image showing a diffusion tensor image with three color levels depicting the orientations of the three eigenvectors of the matrix image ${\displaystyle I(x),x\in {\mathbb {R} }^{2}}$, matrix valued image; each of three colors represents a direction.

Dense LDDMM tensor matching^[104][108] takes the images as 3x1 vectors and 3x3 tensors solving the variational problem matching between coordinate system based on the principle eigenvectors of the diffusion tensor MRI image (DTI) denoted ${\displaystyle M(x),x\in {\mathbb {R} }^{3}}$ consisting of the ${\displaystyle 3\times 3}$-tensor at every voxel. Several of the group actions defined based on the Frobenius matrix norm between square matrices ${\displaystyle \|A\|_{F}^{2}\doteq traceA^{T}A}$ . Shown in the accompanying figure is a DTI image illustrated via its color map depicting the eigenvector orientations of the DTI matrix at each voxel with color determined by the orientation of the directions. Denote the ${\displaystyle 3\times 3}$ tensor image ${\displaystyle M(x),x\in {\mathbb {R} }^{3}}$ with eigen-elements ${\displaystyle \{\lambda _{i}(x),e_{i}(x),i=1,2,3\}}$, ${\displaystyle \lambda _{1}\geq \lambda _{2}\geq \lambda _{3}}$.

Coordinate system transformation based on DTI imaging has exploited two actions one based on the principle eigen-vector or entire matrix.

LDDMM matching based on the principal eigenvector of the diffusion tensor matrix takes the image ${\displaystyle I(x),x\in {\mathbb {R} }^{3}}$ as a unit vector field defined by the first eigenvector. The group action becomes

${\displaystyle \varphi \cdot I={\begin{cases}{\frac {D_{\varphi ^{-1}}\varphi I\circ \varphi ^{-1}\|I\circ \varphi ^{-1}\|}{\|D_{\varphi ^{-1}}\varphi I\circ \varphi ^{-1}\|}}&I\circ \varphi \neq 0;\\0&{\text{otherwise.}}\end{cases}}}$

LDDMM matching based on the entire tensor matrix has group action becomes ${\displaystyle \varphi \cdot M=(\lambda _{1}{\hat {e}}_{1}{\hat {e}}_{1}^{T}+\lambda _{2}{\hat {e}}_{2}{\hat {e}}_{2}^{T}+\lambda _{3}{\hat {e}}_{3}{\hat {e}}_{3}^{T})\circ \varphi ^{-1},}$ transformed eigenvectors

${\displaystyle {\begin{aligned}{\hat {e}}_{1}&={\frac {D\varphi e_{1}}{\|D\varphi e_{1}\|}}\ ,\ \ \ {\hat {e}}_{2}={\frac {D\varphi e_{2}-\langle {\hat {e}}_{1},D\varphi e_{2}\rangle {\hat {e}}_{1}}{\sqrt {\|D\varphi e_{2}\|^{2}-\langle {\hat {e}}_{1},D\varphi e_{2}\rangle ^{2}}}}\ ,\ \ \ {\hat {e}}_{3}={\hat {e}}_{1}\times {\hat {e}}_{2}\end{aligned}}}$.

The variational problem matching onto the principal eigenvector or the matrix is described LDDMM Tensor Image Matching.