Large deformation diffeomorphic metric mapping

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, and subsequently the Ulf Grenander school at Brown University with the HAND experiments. In the 90's there were several solutions for image registration which were associated to linearizations of small deformation and non-linear elasticity. 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, Rabbitt and Miller and Trouvé. 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. 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 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 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 LDDMM, curves, surfaces, dense vector and tensor imagery, and varifolds removing orientation. ==The diffeomorphism orbit model in computational anatomy==

Deformable shape in computational anatomy (CA)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_\text{temp} \circ \varphi, \varphi \in \operatorname{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 \varphi_t, t \in [0,1], with \varphi \doteq \varphi_1, satisfying the Lagrangian and Eulerian specification of the flow field associated to the ordinary differential equation, : \frac{d}{dt} \varphi_t=v_t \circ \varphi_t, \ \varphi_0={\rm 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. 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-times square-integrable weak-derivatives. Thus (V, \|\cdot \|_V ) embeds smoothly in 1-time continuously differentiable functions. The diffeomorphism group are flows with vector fields absolutely integrable in Sobolev norm {{NumBlk|::| \operatorname{Diff}_V \doteq \{\varphi=\varphi_1: \dot \varphi_t=v_t \circ \varphi_t, \varphi_0={\rm id}, \int_0^1 \|v_t \|_V \, dt |}} ==The variational problem of dense image matching and sparse landmark matching==

LDDMM algorithm for dense image matching 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 \int_{R^3} Av\cdot v \, dx, \ v \in V \ , where the integral is calculated by integration by parts when Av is a generalized function in the dual space 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 for the necessary conditions on the norm for existence of solutions. The original large deformation diffeomorphic metric mapping (LDDMM) algorithms of Beg, Miller, Trouve, Younes 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. Beg solved the dense image matching minimizing the action integral of kinetic energy of diffeomorphic flow while minimizing endpoint matching term according to{{NumBlk|| \min_{v: \dot \phi=v \circ \phi, \phi_0=id} C(v) \doteq \frac{1}{2} \int_0^1 \int_{R^3} Av_t \cdot v_t \, dx \, dt +\frac{1}{2} \int_{R^3}|I\circ \phi_1^{-1}-J|^2 \, dx • Beg's Iterative Algorithm for Dense Image Matching Update until convergence, \phi_t^{old} \leftarrow \phi_t^{new} each iteration, with \phi_{t1} \doteq \phi_1 \circ \phi_t^{-1}:{{NumBlk|| \begin{cases} &v_t^{new}(\cdot)=v_t^{old}(\cdot) - \epsilon(v_t^{old} - \int_{R^3} K(\cdot, y) (I \circ \phi_t^{-1 old}(y) -J \circ \phi_{t1}^{old}(y) ) \nabla(I \circ \phi_t^{-1old}(y)) \\ & \dot \phi_t^{new}=v_t^{new} \circ \phi_t^{new}, t \in [0,1] \end{cases} |}} 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 according to :Av_t^*=(D \phi_t^{*-1})^T Av_0^* \circ \phi_t^{*-1}|D \phi_t^{*-1}| {{Hide| First Order Variation of the Flow and Vector Field for Dense Image and Landmark Matching| The first order variation in the vector fields v+\epsilon \delta v requires 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 the solution to the Lie bracket \frac{d }{d t} \left( \delta \phi_{|\phi} \right)=(Dv)_{|\phi} \delta \phi_{|\phi} + \delta v_{|\phi} giving : \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 • Image Matching: Taking the directional derivative of the image endpoint condition E(\phi)=\int_X|I \circ \phi^{-1} - J|^2 dx 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}}} dt dx . Substituting \phi_{t1} \doteq \phi_1 \circ \phi_t^{-1} gives the necessary condition for an optimum: : \begin{align} \frac{d}{d \epsilon}C(v +\epsilon\delta v)|_{\epsilon=0}&=\int_0^1 \int_X Av_t \cdot \delta v_t\ dx\ dt -\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}}}\ dx\,dt \\ & =\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\ dx\,dt \\ &=0 \end{align} . • Landmark Matching: Take the variation in the vector fields v+\epsilon \delta v of \frac{1}{2} \sum_{i}|\phi_1(x_i)-y_i)|^2 use the chain rule for the perturbation \delta \phi \circ \phi to gives the first variation : \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) dx dt =\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) dx dt }} ==LDDMM Diffusion Tensor Image Matching==

LDDMM matching based on the principal eigenvector of the diffusion tensor matrix takes the image I(x), x \in {\mathbb R}^3 as a unit vector field defined by the first eigenvector. has group action \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 :\begin{align} \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{align}. Dense matching problem onto principle eigenvector of DTI The variational problem matching onto vector image I^\prime (x), x \in {\mathbb R}^3with endpoint : E(\phi_1) \doteq \alpha \int_{{\mathbb R}^3} \|\phi_1 \cdot I -I^\prime \|^2 \, dx+\beta \int_{{\mathbb R}^3} (\|\phi_1 \cdot I\|-\|I^\prime\|)^2 \, dx). becomes : \min_{v: \dot \phi \circ \phi^{-1}} \frac{1}{2} \int_0^1 \int_{R^3} Av_t \cdot v_t dx dt + \alpha \int_{{\mathbb R}^3} \|\phi_1 \cdot I -I^\prime \|^2 \, dx+\beta \int_{{\mathbb R}^3} (\|\phi_1 \cdot I\|-\|I^\prime\|)^2 \, dx \ . Dense matching problem onto DTI MATRIX The variational problem matching onto: M^\prime (x), x \in {\mathbb R}^3 with endpoint : E(\phi_1) \doteq \int_{{\mathbb R}^3} \|\phi_1 \cdot M(x) - M^\prime(x) \|_F^2 dx with \|\cdot \|_F Frobenius norm, giving variational problem {{NumBlk|::| \min_{v: v=\dot \phi \circ \phi^{-1}} \frac{1}{2} \int_0^1 \int_{R^3} Av_t \cdot v_t dx dt +\alpha \int_{{\mathbb R}^3} \|\phi_1 \cdot M(x) - M^\prime(x) \|_F^2 dx |}} ==LDDMM ODF==

High angular resolution diffusion imaging (HARDI) addresses the well-known limitation of DTI, that is, DTI can only reveal one dominant fiber orientation at each location. HARDI measures diffusion along n uniformly distributed directions on the sphere and can characterize more complex fiber geometries by reconstructing an orientation distribution function (ODF) that characterizes the angular profile of the diffusion probability density function of water molecules. The ODF is a function defined on a unit sphere, {\mathbb{S}}^2 . Denote the square-root ODF ({\sqrt{\text{ODF}}} ) as \psi({\bf s}) , where \psi({\bf s}) is non-negative to ensure uniqueness and \int_{{\bf s}\in {\mathbb{S}}^2} \psi^2({\bf s}) d{\bf s}=1. The metric defines the distance between two {\sqrt{\text{ODF}}} functions \psi_1, \psi_2 \in \Psi as : \begin{align} \rho(\psi_1,\psi_2)= \|\log_{\psi_1}(\psi_2)\|_{\psi_1} =\cos^{-1} \langle \psi_1, \psi_2 \rangle= \cos^{-1}\left(\int_{{\bf s}\in {\mathbb{S}}^2} \psi_1({\bf s}) \psi_2({\bf s})d{\bf s}\right), \end{align} where \langle \cdot, \cdot \rangle is the normal dot product between points in the sphere under the \mathrm{L}^2 metric. The template and target are denoted \psi_{\mathrm{temp}}({\bf s},x), \psi_{\mathrm{targ}}({\bf s}, x), {\bf s} \in {{\mathbb{S}}^2} x \in X indexed across the unit sphere and the image domain, with the target indexed similarly. Define the variational problem assuming that two ODF volumes can be generated from one to another via flows of diffeomorphisms \phi_t , which are solutions of ordinary differential equations \dot \phi_t=v_t (\phi_t), t \in [0,1],\phi_0={id}. The group action of the diffeomorphism on the template is given according to \phi_1 \cdot \psi (x)\doteq (D\phi_1) \psi \circ \phi_1^{-1}(x), x \in X , where (D\phi_1) is the Jacobian of the affined transformed ODF and is defined as : \begin{align} (D \phi_1) \psi\circ \phi_1^{-1}(x)=\sqrt{\frac{\det{\bigl(D_{\phi_1^{-1}}\phi_1 \bigr)^{-1}} }{\left\|{\bigl(D_{\phi_1^{-1}}\phi_1 \bigr)^{-1} } {\bf s} \right\|^3} } \quad \psi \left( \frac{(D_{\phi_1^{-1}}\phi_1 \bigr)^{-1} {\bf s}}{\|(D_{\phi_1^{-1}}\phi_1 \bigr)^{-1} {\bf s}\|}, \phi_1^{-1}(x) \right) . \end{align} The LDDMM variational problem is defined as : \begin{align} \min_{v: \dot \phi_t= v_t \circ \phi_t,\phi_0={ id}}\int_0^1 \int_{R^3} Av_t \cdot v_t dx \ dt +\lambda \int_{R^3}\|\log_{(D\phi_1) \psi_{\mathrm{temp}} \circ \phi_1^{-1}(x)}(\psi_{\mathrm{targ}}(x))\|^2_{(D\phi_1) \psi_{\mathrm{temp}}\circ \phi_1^{-1}(x)}dx \end{align} . ==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. Another solution by Vialard, 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=R^3 , 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 \int_{R^3} Av_t \cdot v_t dx dt +\frac{1}{2} \int_{{\mathbb R}^3}|q_1(x)-J(x)|^2dx \end{matrix} |}}{{NumBlk|| \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} {{hidden 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 : \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 . }} ==Software for diffeomorphic mapping==

Software suites containing a variety of diffeomorphic mapping algorithms include the following: • Deformetrica • ANTS Voxel-based morphometry(VBM) • DEMONS • LDDMM ==See also==