Riemannian metric and Lie bracket in computational anatomy

Computational anatomy (CA) is the study of shape and form in medical imaging. The study of deformable shapes in CA relies on the high-dimensional diffeomorphism group $\operatorname{Diff}_V$ , which generates orbits of the form $\mathcal{M} \doteq \{ \varphi \cdot m \mid \varphi \in \operatorname{Diff}_V \}$ . In CA the orbit is usually treated as a smooth Riemannian manifold, because at every point of the shape manifold $m \in \mathcal{M}$ an inner product on the tangent space induces a norm $\| \cdot \|_m$ that varies smoothly over the manifold. This is achieved by treating the diffeomorphism group $\varphi \in \operatorname{Diff}_V$ as a Riemannian manifold, with norm $\| \cdot \|_\varphi$ on the tangent space at $\varphi \in \operatorname{Diff}_V$ . The action of the diffeomorphism group then induces a norm and a metric on the orbit $m \in \mathcal{M}$ .

The diffeomorphism group generated by Lagrangian and Eulerian flows

Diffeomorphisms in computational anatomy are constructed as flows $\varphi_t, t \in [0,1]$ that satisfy the Lagrangian and Eulerian velocity-field specifications and are generated by the ordinary differential equation

\frac{d}{dt} \varphi_t = v_t \circ \varphi_t , \ \varphi_0 = \operatorname{id};

where $v_t = \dot \varphi_t \circ \varphi_t^{-1}, t \in [0,1]$ is the Eulerian vector field $v \doteq (v_1,v_2,v_3)$ in $\mathbb{R}^3$ . The inverse of the flow satisfies

\frac{d}{dt} \varphi_t^{-1} = -(D \varphi_t^{-1}) v_t, \ \varphi_0^{-1} = \operatorname{id},

and for a flow in $\mathbb{R}^3$ the $3 \times 3$ Jacobian matrix is $D\varphi \doteq \left(\frac{\partial \varphi_i}{\partial x_j}\right)$ .

To ensure a smooth flow of diffeomorphisms with smooth inverses, the vector fields on $\mathbb{R}^3$ must be at least once continuously differentiable in space.^[1]^[2] By the Sobolev embedding theorems they are modelled as elements of a Hilbert space $(V, \| \cdot \|_V)$ , so that each element $v_i \in H_0^3, i=1,2,3,$ with three square-integrable derivatives means that $(V, \| \cdot \|_V)$ embeds smoothly in the space of once continuously differentiable functions. The diffeomorphism group is the set of flows generated by vector fields that are absolutely integrable in the Sobolev norm:

\operatorname{Diff}_V \doteq \{\varphi=\varphi_1: \dot \varphi_t = v_t \circ \varphi_t , \varphi_0 = \operatorname{id}, \int_0^1 \|v_t \|_V \,dt < \infty \} \ .

Riemannian orbit model

Shapes in computational anatomy are studied through diffeomorphic mappings that establish correspondences between anatomical coordinate systems. In this framework, three-dimensional medical images are modelled as diffeomorphic transforms of a template $I_{temp}$ , so that an observed image is an element of the CA random orbit model. For images this gives $I \in \mathcal{I} \doteq \{ I = I_{temp} \circ \varphi, \varphi \in \operatorname{Diff}_V \}$ , and for charts representing submanifolds it gives $\mathcal{M} \doteq \{ \varphi \cdot m_{temp} : \varphi \in \operatorname{Diff}_V \}$ .

Riemannian metric

Orbits of shape and form in computational anatomy are generated by group action, $\mathcal{M} \doteq \{ \varphi \cdot m : \varphi \in \operatorname{Diff}_V \}$ . They are made into Riemannian orbits by introducing a metric associated with each point and its tangent space. A metric is defined on the group, and it induces a metric on the orbit. In computational anatomy the metric on the tangent space at each element $\varphi \in \operatorname{Diff}_V$ is taken as

\| \dot \varphi \|_\varphi \doteq \| \dot \varphi \circ \varphi^{-1} \|_V=\| v \|_V

where vector fields are modelled as elements of a Hilbert space with Hilbert-space norm $(V, \| \cdot \|_V)$ . The space $V$ is modelled as a reproducing kernel Hilbert space (RKHS) defined by a one-to-one differential operator $A: V \rightarrow V^*$ . For $\sigma(v) \doteq Av \in V^*$ , a distribution or generalised function, the linear form $(\sigma\mid w) \doteq \int_{\mathbb{R}^3} \sum_i w_i(x) \sigma_i(dx)$ determines the norm and inner product. For $v \in V$ they are defined by

\langle v , w \rangle_V \doteq \int_X A v \cdot w \, dx, \ \| v\|_V^2 \doteq \int_X A v \cdot v \, dx, \ v,w \in V \ .

The integrals are computed by integration by parts, with $Av$ a generalised function in the dual space $V^*$ . The differential operator is chosen so that the Green's kernel associated with its inverse is smooth enough to support once continuously differentiable vector fields.

Right-invariant metric on the diffeomorphism group

The metric on the diffeomorphism group is defined as the distance between pairs of group elements:

d_{\operatorname{Diff}_V}(\psi, \varphi) = \inf_{v_t} \left(\frac{1}{2} \int_0^1 \int_X Av_t \cdot v_t \,dx \ dt : \varphi_0 = \psi, \varphi_1 = \varphi, \dot \varphi_t = v_t \circ \varphi_t \right)^{1/2} \ .

This distance provides a right-invariant metric for diffeomorphometry^[3]^[4], and it is invariant to reparameterisation of space for all $\varphi \in \operatorname{Diff}_V$ :

d_{\operatorname{Diff}_V}(\psi, \varphi) = d_{\operatorname{Diff}_V}(\psi \circ \varphi, \varphi \circ \varphi).

Lie bracket on the diffeomorphism group

The Lie bracket describes the adjustment of the velocity term produced by a perturbation of motion in a curved space. Using Hamilton's principle of least action, the optimised flow is derived as a critical point of the action integral, the integral of kinetic energy. The Lie bracket of vector fields in computational anatomy was first introduced by Miller, Trouvé and Younes^[5]. The derivation computes the perturbation $\delta v$ of the vector field $v^\varepsilon = v + \varepsilon \delta v$ , expressed as the time derivative of a group perturbation and corrected by the Lie bracket of vector fields involving the Jacobian, as distinct from the matrix-group case:

ad_v:V \mapsto V

is defined by

ad_v(w)\doteq (Dv)w - (Dw)v , \ v,w \in V .

Proof: To prove the Lie bracket of vector fields, take a first-order perturbation of the flow at a point $\varphi \in \operatorname{Diff}_V$ .

Taking a first-order perturbation gives $\varphi_t^\varepsilon \doteq (\operatorname{id} + \varepsilon w) \circ \varphi = \varphi+\varepsilon w \circ \varphi$ , with fixed boundary $w_0= w_1=0$ and $\frac{d}{dt} \varphi_t^\varepsilon = v_t^\varepsilon \circ \varphi_t^\varepsilon, \varphi_0^\varepsilon = \operatorname{id}, \varphi_1^\varepsilon = \varphi_1$ , which yields the two equations

$\frac{d}{dt} \varphi_t^\varepsilon = \frac{d}{dt} \varphi_t + \varepsilon \frac{d}{dt}(w_t \circ \varphi_t) =v_t \circ \varphi_t + (D w_t)\circ \varphi_t v_t \circ \varphi_t +o(\varepsilon) \ .$
$\dot \varphi_t^\varepsilon = (v_t + \varepsilon \delta v_t) \circ (\varphi_t + \varepsilon w_t \circ \varphi_t) \simeq v_t \circ \varphi_t +\varepsilon (Dv_t)\circ \varphi_t w_t \circ \varphi_t + \delta v_t \circ \varphi_t +o(\varepsilon) \ .$

Setting the two equations equal gives the vector-field perturbation adjusted by the Lie bracket. The Lie bracket gives the first-order variation of the vector field relative to the flow:

\delta v_t = \frac{d}{d t} w_t - ad_{v_t}(w_t) =\frac{d}{d t} w_t-((Dv_t) w_t - (Dw_t)v_t) \ .

Generalised Euler-Lagrange equation for the diffeomorphic flow metric

The Euler-Lagrange equation can be used to compute the geodesic flow in the group, which is the basis of the metric. The action integral of the kinetic-energy Lagrangian in Hamilton's principle becomes

J(\varphi) \doteq \frac{1}{2}\int_0^1 \| \dot \varphi_t \|_{\varphi_t}^2 \,dt = \frac{1}{2}\int_0^1 \| \dot \varphi_t \circ \varphi_t^{-1} \|_V^2 \,dt =\frac{1}{2}\int_0^1 \int_X A (\dot \varphi_t \circ \varphi_t^{-1})\cdot (\dot \varphi_t \circ \varphi_t^{-1}) \, dx \,dt \ .

The action integral expressed in terms of the vector field corresponds to the integral of kinetic energy:

J(v) \doteq \frac{1}{2} \int_0^1 \| v_t \|_V^2 dt = \frac{1}{2} \int_0^1 \int_X Av_t \cdot v_t \,dx \ dt \ .

The geodesic connections of shortest path in the orbit are defined by Hamilton's principle of least action, which requires the first variation of solutions in the computational anatomy orbit and is based on critical points of the path metric length or energy.

The original derivation of the Euler equation^[5] uses the generalised-function equation when $Av \in V^*$ is a distribution. Using the adjoint operator of the Lie bracket to take the first variation of the action integral, for all smooth $w \in V$ one obtains

\frac{d}{d \varepsilon} J(\varphi^\varepsilon)|_{\varepsilon=0} = \int_0^1 \int_X Av_t \cdot \delta v_t\,dx \, dt = \int_0^1 \int_X Av_t \cdot \left( \frac{d}{d t} w_t -((Dv_t)w-(Dw)v_t) \right)\,dx \,dt.

Using the bracket $ad_v: w \in V \mapsto V$ and $ad_v^*: V^* \rightarrow V^*$ gives

\frac{d}{dt} Av_t + ad_{v_t}^*(Av_t)=0 \ , \ t \in [0,1] \ ,

which means that for all smooth $w \in V$ ,

\int_X \left( \frac{d}{dt} Av_t + ad_{v_t}^*(Av_t) \right) \cdot w\, dx = \int_X \frac{d}{dt} Av_t \cdot w \,dx + \int_X Av_t \cdot \left( (Dv_t)w-(Dw)v_t \right) \,dx =0 .

This equation is the Euler equation when the diffeomorphic shape momentum is a generalised function^[6].

It is called EPDiff, the Euler-Poincaré equation for diffeomorphisms, and it has been studied in fluid mechanics in connection with incompressible fluids and the $L^2$ metric^[7]^[8].

Riemannian exponential for positioning

In the random orbit model of computational anatomy, the entire flow is reduced to an initial condition that forms the coordinate encoding the diffeomorphism and provides a way to position information in the orbit. This was first called the geodesic positioning system by Miller, Trouvé and Younes^[9]. Starting from the initial condition $v_0$ , geodesic positioning relative to the computational anatomy Riemannian metric solves the flow of the Euler-Lagrange equation. Solving the geodesic from the initial condition $v_0$ is called the Riemannian exponential, a map $\operatorname{Exp}_{\operatorname{id}}(\cdot): V \to \operatorname{Diff}_V$ from the identity to the group.

The Riemannian exponential satisfies $\operatorname{Exp}_{\operatorname{id}}(v_0)= \varphi_1$ , with initial condition $\dot \varphi_0 = v_0$ and vector-field dynamics $\dot \varphi_t = v_t \circ \varphi_t, t \in [0,1]$ .

For the classical equation in which the diffeomorphic shape momentum is a smooth vector, $Av_t = \mu_t \,dx$ , with $\int_X \mu_t \cdot w \,dx \ ,w \in V$ , the Euler equation exists in classical density form^[10]:

\frac{d}{dt} \mu_t + (Dv_t)^T \mu_t +(D\mu_t)v_t + ( \nabla \cdot v) \mu_t =0 \ , \ Av_t= \mu_t \, dx \ ;

For the generalised equation $Av \in V^*$ ,

\frac{d}{dt} Av_t + ad_{v_t}^*(Av_t)=0 \ , \ t \in [0,1] \ .

It extends to the whole group as

\varphi= \operatorname{Exp}_\varphi(v_0\circ \varphi) \doteq \operatorname{Exp}_{\operatorname{id}}(v_0) \circ \varphi

Variational problem for matching or registering coordinate-system information

Matching information across coordinate systems is central to computational anatomy. A matching term $E: \varphi \in \operatorname{Diff}_V \rightarrow R^+$ is added to the action integral, and it represents the target endpoint:

C(\varphi) \doteq \int_0^1 \int_X Av_t \cdot v_t \,dx \, dt + E(\varphi_1) \ .

The endpoint term adds a boundary condition to the Euler-Lagrange equation and gives an Euler equation with a boundary term. Taking the variation gives

the necessary geodesic conditions:

\begin{cases} & \dfrac{d}{dt} Av_t + (Dv_t)^T Av_t +(DAv_t)v_t + ( \nabla \cdot v) Av_t =0 \ ; \\ & Av_1 + \frac{\partial E(\varphi)}{\partial \varphi_1} = 0 \end{cases}

Proof: The proof by variational calculus uses the perturbation above and classical variational arguments^[10].

\begin{align} & \int_0^1 \int_X Av_t \cdot \left(\frac{d}{d t} \delta \varphi_t - (Dv_t \delta \varphi_t-D\delta \varphi_t v_t) \right)\,dx \, dt + \int_X \left( \frac{\partial E(\varphi)}{\partial \varphi_1}\cdot \delta \varphi_1 \,dx\right) \\ = {} & -\int_0^1 \int_X \left(\frac{d Av_t}{d t}+ad_{v_t}^*(Av_t)\right) \cdot \delta \varphi_t \,dx \,dt + \int_X \left(Av_1 + \frac{\partial E(\varphi)}{\partial \varphi_1} \right) \cdot \delta \varphi_1 \,dx . \end{align}

Euler-Lagrange geodesic endpoint conditions for image matching

The earliest large deformation diffeomorphic metric mapping (LDDMM) algorithms solved the matching problem associated with images and the registration of landmarks. Landmarks lie in a vector space. Image-matching geodesic equations satisfy the classical dynamical equations with endpoint conditions. The necessary conditions for image-matching geodesics appear as the classical Euler-Lagrange equation with boundary conditions:

\min_{\varphi: \dot \varphi = v_t \circ \varphi_t } C(\varphi) \doteq \frac{1}{2} \int_0^1 \int_X Av_t \cdot v_t \,dx\,dt +\frac{1}{2} \int_X |I \circ \varphi_1^{-1}(x) - J(x) |^2 \, dx

the necessary geodesic conditions:

\begin{cases} & \dfrac{d}{dt} Av_t + (Dv_t)^T Av_t +(DAv_t)v_t + ( \nabla \cdot v) Av_t =0 \ ; \\ & Av_1 =(I \circ \varphi_1^{-1} -J) \nabla (I\circ \varphi_1^{-1}) \end{cases}

Euler-Lagrange geodesic endpoint conditions for landmark matching

The registered landmark-matching problem satisfies the generalised-function dynamical equation with endpoint conditions:

\min_{\varphi: \dot \varphi = v_t \circ \varphi_t } C(\varphi) \doteq \frac{1}{2} \int_0^1 \int_X Av_t \cdot v_t\, dx\,dt +\frac{1}{2} \sum_i (\varphi_1(x_i)-y_i )\cdot (\varphi_1(x_i)-y_i) .

the necessary geodesic conditions:

\begin{cases} & \dfrac{d}{dt} Av_t + ad_{v_t}^*(Av_t)=0 \ , \ t \in [0,1] \ , \\ & Av_1 = \sum_{i=1}^n \delta_{\varphi_1 (x_i)} (y_i-\varphi_1(x_i)) \end{cases}

Proof: The variation $\frac{\partial}{\partial \varphi} E(\varphi)$ requires the variation of the inverse $\varphi^{-1}$ , obtained by generalising the matrix inverse perturbation through $(\varphi + \varepsilon \delta \varphi \circ \varphi)\circ (\varphi^{-1} + \varepsilon \delta \varphi^{-1} \circ \varphi^{-1}) = \operatorname{id} + o(\varepsilon)$ , which gives $\delta \varphi^{-1} \circ \varphi^{-1} =-(D \varphi_1^{-1}) \delta \varphi$ and therefore

\begin{align} & \frac{d}{d \varepsilon} \frac{1}{2} \left. \int_X | I \circ (\varphi^{-1} + \varepsilon \delta \varphi^{-1} \circ \varphi^{-1})-J|^2 \, dx\right|_{\varepsilon =0} \\ = {} & \int_X (I \circ \varphi^{-1} -J) \nabla I|_{\varphi^{-1}} (-D \varphi_1^{-1}) \delta \varphi \, dx \\ = {} & -\int_X(I \circ \varphi_1^{-1} -J) \nabla (I\circ \varphi_1^{-1}) \delta \varphi \, dx. \end{align}

References

^P. Dupuis, U. Grenander, M.I. Miller, Existence of Solutions on Flows of Diffeomorphisms, Quarterly of Applied Math, 1997.
^A. Trouvé. Action de groupe de dimension infinie et reconnaissance de formes. C R Acad Sci Paris Sér I Math, 321(8):1031–1034, 1995.
^Michael I. Miller, Laurent Younes B. Group Actions, Homeomorphisms, and Matching: A General Framework. International Journal of Computer Vision. 2001, 41, 61-84. DOI: 10.1023/a:1011161132514
^MI Miller, A Trouvé, L Younes. Hamiltonian Systems and Optimal Control in Computational Anatomy: 100 Years Since D'Arcy Thompson. Annual review of biomedical engineering. 2015, 17, 447-509. DOI: 10.1146/annurev-bioeng-071114-040601
^MI Miller, A Trouvé, L Younes. Geodesic Shooting for Computational Anatomy. Journal of mathematical imaging and vision. 2006, 24, 209-228. DOI: 10.1007/s10851-005-3624-0
^M.I. Miller, A. Trouvé, L. Younes, Geodesic Shooting in Computational Anatomy, IJCV, 2006.
^R Camassa, DD Holm. An integrable shallow water equation with peaked solitons. Physical review letters. 1993, 71, 1661-1664. DOI: 10.1103/physrevlett.71.1661
^Darryl D. Holm, Jerrold Eldon Marsden, Tudor Stefan Ratiu. The Euler–Poincaré Equations and Semidirect Products with Applications to Continuum Theories. Advances in Mathematics. 1998, 137, 1-81. DOI: 10.1006/aima.1998.1721
^MI Miller, L Younes, A Trouvé. Diffeomorphometry and geodesic positioning systems for human anatomy. Technology. 2014, 2, 36. DOI: 10.1142/s2339547814500010
^