User:Mim.cis/sandbox/Derivation of Landmark LDDMM Conditions

A diffeomorphic mapping system is a system designed to map, manipulate, and transfer information which is stored in many types of spatially distributed medical imagery.

Diffeomorphic mapping is the underlying technology for mapping and analyzing information measured in human anatomical coordinate systems which have been measured via Medical imaging. Diffeomorphic mapping is a broad term that actually refers to a number of different algorithms, processes, and methods. It is attached to many operations and has many applications for analysis and visualization. Diffeomorphic mapping can be used to relate various sources of information which are indexed as a function of spatial postiion as the key index variable. Diffeomorphisms are by their latin root structure preserving transformations, which are in turn differentiable and therefore smooth, allowing for the calculation of metric based quantities such as arc length and surface areas. Spatial location and extents in human anatomical coordinate systems can be recorded via a variety of Medical imaging modalities, generally termed multi-modal medical imagery, providing either scalar and or vector quantities at each spatial location. Examples are scalar T1 or T2 Magnetic resonance imagery, or as 3x3 diffusion tensor matrices Diffusion MRI and Diffusion-weighted imaging, to scalar densities associated to Computed Tomography (CT),or functional imagery such as temporal data of functional magnetic resonance imaging and scalar densities such as Positron emission tomography (PET).

The acronym LDDMM, standing for large deformation diffeomorphic metric mapping is often used to refer to a specific suite of algorithms used for diffeomorphic mapping and manipulating dense imagery based on diffeomorphic metric mapping within the academic discipline of Computational anatomy, to be distinguised from it's precursor based on diffeomorphic mapping. The distinction between the two is that diffeomorphic metric maps satisfy the property that the length associated to their flow away from the identity induces a metric on the group of diffeomorphisms, which in turn induces a metric on the orbit of Shapes and Forms within the field of Computational Anatomy. The study of shapes and forms with the metric of diffeomorphic metric mapping is called Diffeomorphometry.

Computational anatomy is a subdiscipline within the broader field of Neuroinformatics within Bioinformatics and Medical imaging. The first algorithm for dense image mapping via diffeomorphic metric mapping was Beg's LDDMM^[1]^[2] for volumes and Joshi's landmark matching for point sets with correspondence,^[3]^[4] with LDDMM algorithms now available for computing diffeomorphic metric maps between non-corresponding landmarks^[5] and landmark matching intrinsic to spherical manifolds,^[6] curves,^[7] currents and surfaces,^[8]^[9]^[10] tensors,^[11] varifolds,^[12] and time-series.^[13]^[14]^[15] The term LDDMM was first established as part of the National Institutes of Health supported Biomedical Informatics Research Network.^[16]

In a more general sense, diffeomorphic mapping is a term that describes any solution that registers or builds correspondences between dense coordinate systems in Medical imaging by ensuring the solutions are diffeomorphic. There are now many codes organized around diffeomorphic registration^[17] including ANTS,^[18] DARTEL,^[19] DEMONS,^[20] StationaryLDDMM^[21] as examples of actively used computational codes for constructing correspondences between coordinate systems based on dense images.

The distinction between diffeomorphic metric mapping forming the basis for LDDMM and the earliest methods of diffeomorphic mapping is the introduction of a Hamilton principle of least-action in which large deformations are selected of shortest length corresponding to geodesic flows. The lengths of these geodesics form the metric; diffeomorphic mapping in general does not correspond to any metric formulation.

History of Development

Diffeomorphic mapping 3-dimensional information across coordinate systems is central to high-resolution Medical imaging and the area of Neuroinformatics within the newly emerging field of bioinformatics. Diffeomorphic mapping 3-dimensional coordinate systems as measured via high resolution dense imagery has a long history in 3-D beginning with Computed Axial Tomography (CAT scanning) in the early 80's by the University of Pennsylvania group led by Ruzena Bajcsy,^[22] and subsequently the Ulf Grenander school at Brown University with the HAND experiments.^[23]^[24] In the 90's there were several solutions for image registration which were associated to linearizations of small deformation and non-linear elasticity.^[25]^[26]^[27]^[28]^[29]

The central focus of the sub-field of Computational anatomy (CA) within medical imaging is mapping information across anatomical coordinate systems at the 1 millimeter morphome scale. In CA mapping of dense information measured within Magnetic resonance image (MRI) based coordinate systems such as in the brain has been solved via inexact matching of 3D MR images one onto the other. The earliest introduction of the use of diffeomorphic mapping via large deformation flows of diffeomorphisms for transformation of coordinate systems in image analysis and medical imaging was by Christensen et. al.^[17]^[30] and Trouve.^[31] The introduction of flows, which are akin to the equations of motion used in fluid dynamics, exploit the notion that dense coordinates in image analysis follow the Lagrangian and Eulerian equations of motion. This model becomes more appropriate for cross-sectional studies in which brains and or hearts are not necessarily deformations of one to the other. Methods based on linear or non-linear elasticity energetics which grows with distance from the identity mapping of the template, is not appropriate for cross-sectional study. Rather, in models based on Lagrangian and Eulerian flows of diffeomorphisms, the constraint is associated to topological properties, such as open sets being preserved, coordinates not crossing implying uniqueness and existence of the inverse mapping, and connected sets remaining connected. The use of diffeomorphic methods grew quickly to dominate the field of mapping methods post Christensen's original paper, with fast and symmetric methods becoming available.^[19]^[32]

Such methods are powerful in that they introduce notions of regularity of the solutions so that they can be differentiated and local inverses can be calculated. The disadvantages of these methods is that there was no associated global least-action property which could score the flows of minimum energy. This contrasts the geodesic motions which are central to the study of Rigid body kinematics and the many problems solved in Physics via Hamilton's principle of least action. In 1998, Dupuis, Grenander and Miller^[33] established the conditions for guaranteeing the existence of solutions for dense image matching in the space of flows of diffeomorphisms. These conditions require an action penalizing kinetic energy measured via the Sobolev norm on spatial derivatives of the flow of vector fields.

The Large Deformation Diffeomorphic Metric Mapping (LDDMM) code that Faisal Beg derived and implemented for his PhD at Johns Hopkins University^[34] developed the earliest algorithmic code which solved for flows with fixed points satisfying the necessary conditions for the dense image matching problem subject to least-action. Computational anatomy now has many existing codes organized around diffeomorphic registration^[17] including ANTS,^[18] DARTEL,^[19] DEMONS,^[35] LDDMM,^[2] StationaryLDDMM^[21] as examples of actively used computational codes for constructing correspondences between coordinate systems based on dense images.

These large deformation methods have been extended to landmarks without registration via measure matching,^[36] curves,^[37] surfaces,^[38] dense vector^[39] and tensor ^[40] imagery, and varifolds removing orientation.^[41]

The Diffeomorphism Orbit Model in Computational Anatomy

Deformable shape in Computational Anatomy (CA)^[42]^[43]^[44]^[45]is studied via the use of diffeomorphic mapping for establishing correspondences between anatomical coordinates in Medical Imaging. In this setting, three dimensional medical images are modelled as a random deformation of some exemplar, termed the template $I_{temp}$ , with the set of observed images element in the random orbit model of CA for images $I\in {\mathcal {I}}\doteq \{I=I_{temp}\circ \varphi ,\varphi \in Diff_{V}\}$ . The template is mapped onto the target by defining a variational problem in which the template is transformed via the diffeomorphism used as a change of coordinate to minimize a squared-error matching condition between the transformed template and the target.

The diffeomorphisms are generated via smooth flows $\phi _{t},t\in [0,1]$ , with $\varphi \doteq \phi _{1}$ , satisfying the Lagrangian and Eulerian specification of the flow field associated to the ordinary differential equation,

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

,

with $v_{t},t\in [0,1]$ the Eulerian vector fields determining the flow. The vector fields are guaranteed to be 1-time continuously differentiable $v_{t}\in C^{1}$ by modelling them to be in a smooth Hilbert space $v\in V$ supporting 1-continuous derivative.^[46] The inverse $\phi _{t}^{-1},t\in [0,1]$ is defined by the Eulerian vector-field with flow given by

:

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

(Inverse Transport flow)

To ensure smooth flows of diffeomorphisms with inverse, the vector fields ${\mathbb {R} }^{3}$ must be at least 1-time continuously differentiable in space^[35]^[47] which are modelled as elements of the Hilbert space $(V,\|\cdot \|_{V})$ using the Sobolev embedding theorems so that each element $v_{i}\in H_{0}^{3},i=1,2,3,$ has 3-square-integrable derivatives. Thus $(V,\|\cdot \|_{V})$ embeds smoothly in 1-time continuously differentiable functions.^[35]^[47] The diffeomorphism group are flows with vector fields absolutely integrable in Sobolev norm

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)

In CA the space of vector fields $(V,\|\cdot \|_{V})$ are modelled as a reproducing Kernel Hilbert space (RKHS) defined by a 1-1, differential operator $A:V\rightarrow V^{*}$ determining the norm $\|v\|_{V}^{2}\doteq (Av|v),\ v\in V\ .$ The differential operator is selected so that the Green's kernel the inverse of the operator is continuously differentiable in each variable implying that the vector fields support 1-continuous derivative; see^[48] for the necessary conditions on the norm for existence of solutions.

The Variational Problem of Dense Image Matching and Sparse Landmark Matching

Beg's LDDMM algorithm for image matching

The original large deformation diffeomorphic metric mapping (LDDMM) algorithms of Beg, Miller, Trouve, Younes^[49] was derived taking variations with respect to the vector field parameterization of the group, since $v={\dot {\phi }}\circ \phi ^{-1}$ are in a vector spaces. Variations satisfy the necessary optimality conditions. Beg solved the Control Problem for dense image matching mnimizing endpoint matching term $E(\phi _{1})\doteq {\frac {1}{2}}\|I\circ \phi _{1}^{-1}-J\|^{2}$ simultaneously minimizing the action integral of kinetic energy $\int _{0}^{1}\|v_{t}\|_{V}^{2}dt$ , using the vector fields $v_{t},t\in [0,1]$ via the following iterative algorithm optimization problem and iteration.

${\textstyle \min _{v:{\dot {\phi }}=v\circ \phi ,\phi _{0}=id}C(v)\doteq {\frac {1}{2}}\int _{0}^{1}(Av_{t}|v_{t})dt+{\frac {1}{2}}\int _{X}|I\circ \phi _{1}^{-1}-J|^{2}dx}$

(Variational Problem Images)

Update until convergence, $\phi _{t}^{old}\leftarrow \phi _{t}^{new}$ each iteration, with $\phi _{t1}\doteq \phi _{1}\circ \phi _{t}^{-1}$ :

{\begin{cases}&v_{t}^{new}(\cdot )=v_{t}^{old}(\cdot )-\epsilon (v_{t}^{old}-\int _{X}K(\cdot ,y)(I\circ \phi _{t}^{-1old}(y)-J\circ \phi _{t1}^{old}(y))\nabla (I\circ \phi _{t}^{-1old}(y))|D\phi _{t1}^{old}(y)|dy),t\in [0,1]\\&{\dot {\phi }}_{t}^{new}=v_{t}^{new}\circ \phi _{t}^{new},t\in [0,1]\end{cases}}

(Beg-LDDMM-iteration)

This implies that the fixed point at $t=0$ satisfies $\mu _{0}^{*}=Av_{0}^{*}=(I-J\circ \phi _{1}^{*})\nabla I|D\phi _{1}^{*}|$ , which in turn implies it satisfies the Conservation equation given by the Endpoint Matching Condition according to $Av_{t}^{*}=(D\phi _{t}^{*-1})^{T}Av_{0}^{*}\circ \phi _{t}^{*-1}|D\phi _{t}^{*-1}|$ ^[50]^[51]

Joshi's original LDDMM landmark matching

The landmark matching problem has endpoint ${\frac {1}{2}}\sum _{i}\|\phi _{1}(x_{i})-y_{i}\|^{2}$ was solved originally corresponding to the following:

\min _{v:{\dot {\phi }}_{t}=v_{t}\circ \phi _{t}}C(v)\doteq {\frac {1}{2}}\int \|v_{t}\|_{V}^{2}dt+{\frac {1}{2}}\sum _{i}\|\phi _{1}(x_{i})-y_{i}\|^{2}

;

Update until convergence, $\phi _{t}^{old}\leftarrow \phi _{t}^{new}$ each iteration, with $\phi _{t1}\doteq \phi _{1}\circ \phi _{t}^{-1}$ :

{\begin{cases}&v_{t}^{new}(\cdot )=v_{t}^{old}(\cdot )-\epsilon (v_{t}^{old}+\sum _{i}K(\cdot ,\phi _{t}^{old}(x_{i}))(D\phi _{t1})^{oldT}|_{\phi _{t}^{old}(x_{i})}(y_{i}-\phi _{1}^{old}(x_{i})),t\in [0,1]\\&{\dot {\phi }}_{t}^{new}=v_{t}^{new}\circ \phi _{t}^{new},t\in [0,1]\end{cases}}

(Beg-LDDMM-iteration)

with necessary conditions become $Av_{t}=-\sum _{i}(D\phi _{t1})^{T}|_{\phi _{t}(x_{i})}(y_{i}-\phi _{1}(x_{i}))\delta _{\phi _{t}(x_{i})}$ .

Figure depicts LDMM dense image matching. Top row shows transport of the image under the flow

I\circ \phi _{t}^{-1}

; middle row shows sequence of vector fields

v_{t},

t=0,1/5,2/5,3/5,4/5,1; bottom row shows the sequence of grids under

\phi _{t}.

.

Proof of Gradient and Necesssary Maximizer Conditions

The Calculus of variations was used in Beg^[49]^[52] to derive the iterative algorithm as a solution which when it converges satisfies the necessary maximizer conditions given by the necessary conditions for a first order variation. This also gives the update rule for the gradient algorithm requiring the variation of the endpoint with respect to a first order variation of the vector field. We do this by calculating the directional derivative calculating the Gâteaux derivative. These variations were calculated iin Beg's original paper^[49] as well as in.^[53]^[54]

Taking the first order variation in the vector fields $v+\epsilon \delta v$ using the chain rule requires the computation of three variations. The variation of $\phi ^{-1}$ generalizes the matrix perturbation of the inverse via $(\phi +\epsilon \delta \phi \circ \phi )\circ (\phi ^{-1}+\epsilon \delta \phi ^{-1}\circ \phi ^{-1})=id+o(\epsilon )$ giving $\delta \phi ^{-1}\circ \phi ^{-1}=-(D\phi _{1}^{-1})\delta \phi$ .

To express the variation in terms of $\delta v$ , use $\delta \phi _{1}=(D\phi _{1})_{|\phi _{1}^{-1}}\int _{0}^{1}(D\phi _{t})_{|\phi _{1}^{-1}}^{-1}(\delta v_{t})_{\phi _{t}\circ \phi _{1}^{-1}}dt$ solving Lie bracket ${\frac {d}{dt}}\left(\delta \phi _{|\phi }\right)=(Dv)_{|\phi }\delta \phi _{|\phi }+\delta v_{|\phi }$ . i:.

This gives

{\frac {d}{d\epsilon }}{\frac {1}{2}}\int _{X}|I\circ (\phi ^{-1}+\epsilon \delta \phi ^{-1}\circ \phi ^{-1})-J|^{2}dx|_{\epsilon =0}=\int _{X}(I\circ \phi ^{-1}-J)\nabla I|_{\phi ^{-1}}\delta \phi ^{-1}\circ \phi ^{-1}dx

=\int _{X}(I\circ \phi ^{-1}-J)\nabla I|_{\phi ^{-1}}(-D\phi _{1}^{-1})\delta \phi dx

=\int _{X}(I\circ \phi _{1}^{-1}-J)\nabla I|_{\phi _{1}^{-1}}(-D\phi _{1})_{|\phi _{1}^{-1}}^{-1}(D\phi _{1})_{|\phi _{1}^{-1}})\int _{0}^{1}(D\phi _{t})_{|\phi _{1}^{-1}}^{-1}(\delta v_{t})_{|{\phi _{t}\circ \phi _{1}^{-1}}}dtdx

.

Taking the Gâteaux variation substituting $\phi _{t1}\doteq \phi _{1}\circ \phi _{t}^{-1}$ gives the necessary condition for an optimum:

{\frac {d}{d\epsilon }}C(v+\epsilon \delta v)|_{\epsilon =0}=\int _{0}^{1}\int _{X}Av_{t}\cdot \delta v_{t}dxdt-\int _{0}^{1}\int _{X}(I\circ \phi _{1}^{-1}-J)\nabla I|_{\phi _{1}^{-1}}(D\phi _{t})_{|\phi _{1}^{-1}}^{-1}(\delta v_{t})_{|{\phi _{t}\circ \phi _{1}^{-1}}}dxdt

\ \ \ \ \ \ \ \ =\int _{0}^{1}\int _{X}\left(Av_{t}-(I\circ \phi _{t}^{-1}-J\circ \phi _{t1})\nabla I|_{\phi _{t}^{-1}}(D\phi _{t})_{|\phi _{t}^{-1}}^{-1}|D\phi _{t1}|\right)\cdot \delta v_{t}dxdt=0

.

The perturbation in the vector field gives $(\delta \phi _{1})_{|{\phi _{1}}}=(D\phi _{1})\int _{0}^{1}(D\phi _{t})^{-1}\delta v_{t}\circ \phi _{t}dt$ Take the variation in the vector fields $v+\epsilon \delta v$ using the chain rule ${\frac {\partial E(\phi _{1})}{\partial v}}={\frac {\partial E(\phi _{1})}{\partial \phi }}{\frac {\partial \phi _{1}}{\partial v}}$ which gives the first variation

{\frac {d}{d\epsilon }}C(v+\epsilon \delta v)|_{\epsilon =0}=\sum _{i}(\phi _{1}(x_{i})-y_{i})\cdot D\phi _{1}|_{\phi _{1}^{-1}(\phi _{1}(x_{i}))}\int _{0}^{1}(D\phi _{t})_{|\phi _{1}^{-1}(\phi _{1}(x_{i}))}^{-1}\delta v_{t}|_{\phi _{t}\circ \phi _{1}^{-1}(\phi _{1}(x_{i}))}dt

=\int _{0}^{1}\int _{X}\sum _{i}\delta _{\phi _{t}(x_{i})}(x)(\phi _{1}(x_{i})-y_{i})\cdot (D\phi _{1})_{\phi _{t}^{-1}(x)}(D\phi _{t})_{\phi _{t}^{-1}(x)}^{-1}\delta v_{t}(x)dxdt=\int _{0}^{1}\int _{X}\sum _{i}\delta _{\phi _{t}(x_{i})}(y)(D\phi _{t1})_{\phi _{t}(x_{i})}^{T}(\phi _{1}(x_{i})-y_{i})\cdot \delta v_{t}(x)dxdt

Hamiltonian LDDMM for Dense Image Matching

Beg solved the early LDDMM algorithms by solving the variational matching taking variations with respect to the vector fields.^[55] Another solution by Vialard,^[56] reparameterizes the optimization problem in terms of the state $q_{t}\doteq I\circ \phi _{t}^{-1},q_{0}=I$ , for image $I(x),x\in X$ , with the dynamics equation controlling the state by the control given in terms of the advection equation according to ${\dot {q}}_{t}=-\nabla q_{t}\cdot v_{t}$ . The endpoint matching term $E(q_{1})\doteq {\frac {1}{2}}\|q_{1}-J\|^{2}$ gives the variational problem: {NumBlk||

${\begin{matrix}&\ \ \ \ \ \min _{v:{\dot {q}}=v\circ q}C(v)\doteq {\frac {1}{2}}\int _{0}^{1}(Av_{t}|v_{t})dt+{\frac {1}{2}}\int _{{\mathbb {R} }^{3}}|q_{1}(x)-J(x)|^{2}dx\end{matrix}}$ |Advective-State-Image-Matching}}

{\begin{cases}{\text{Hamiltonian Dynamics}}&\ \ \ \ \ \ \ \ \ \ {\dot {q}}_{t}=-\nabla q_{t}\cdot v_{t}\\&\ \ \ \ \ \ \ \ \ \ {\dot {p}}_{t}=-{\text{div}}(p_{t}v_{t}),\ \ \ \ t\in [0,1]\\&\ \ \ \ \ \ \ \ \ \ Av_{t}=\mu _{t}=-p_{t}\nabla q_{t}\\{\text{Endpoint Condition}}&\ \ \ \ \ \ \ \ \ p_{1}=-{\frac {\partial E}{\partial q_{1}}}(q_{1})=-(q_{1}-J)=-(I\circ \phi _{1}^{-1}-J)\\&\ \ \ \ \ \ \ \ \ \ Av_{1}=\mu _{1}=(I\circ \phi _{1}^{-1}-J)\nabla (I\circ \phi _{1}^{-1})\ \ t=1\ .\\{\text{Conserved Dynamics}}&\ \ \ \ \ \ \ \ \ \ p_{t}=-(I\circ \phi _{t}^{-1}-J\circ \phi _{t1})|D\phi _{t1}|\ ,\ \ t\in [0,1]\ .\\\end{cases}}

(Hamiltonian Matching Condition)

Proof of Hamiltonian Dynamics

The Hamiltonian dynamics with advected state and control dynamics $q_{t}=I\circ \phi _{t}^{-1}$ , ${\dot {q}}=-\nabla q\cdot v$ with extended Hamiltonian $H(q,p,v)=(p|-\nabla q\cdot v)-{\frac {1}{2}}(Av|v)$ gives the variational problem^[52]

\min _{p,q,v}C(p,q,v)\doteq (p|{\dot {q}})-\left((p|-\nabla q\cdot v)-{\frac {1}{2}}(Av|v)\right)+E(q_{1})=(p|{\dot {q}})-H(p,q,v)+E(q_{1})\ .

The first variation gives the condition on the optimizing vector field $Av=-p\nabla q$ , with the endpoint condition $p_{1}=-{\frac {\partial E}{\partial q}}(q_{1})$ and dynamics on the Lagrange multipliers determined by the Gatteux derivative conditions $(-{\dot {p}}-\nabla \cdot (pv)|\delta q))=0$ and the state $(\delta p|{\dot {q}}+\nabla q\cdot v)=0$ .