Toward a Comprehensive Framework for the Spatiotemporal Statistical Analysis of Longitudinal Shape Data

Abstract

This paper proposes an original approach for the statistical analysis of longitudinal shape data. The proposed method allows the characterization of typical growth patterns and subject-specific shape changes in repeated time-series observations of several subjects. This can be seen as the extension of usual longitudinal statistics of scalar measurements to high-dimensional shape or image data. The method is based on the estimation of continuous subject-specific growth trajectories and the comparison of such temporal shape changes across subjects. Differences between growth trajectories are decomposed into morphological deformations, which account for shape changes independent of the time, and time warps, which account for different rates of shape changes over time. Given a longitudinal shape data set, we estimate a mean growth scenario representative of the population, and the variations of this scenario both in terms of shape changes and in terms of change in growth speed. Then, intrinsic statistics are derived in the space of spatiotemporal deformations, which characterize the typical variations in shape and in growth speed within the studied population. They can be used to detect systematic developmental delays across subjects. In the context of neuroscience, we apply this method to analyze the differences in the growth of the hippocampus in children diagnosed with autism, developmental delays and in controls. Result suggest that group differences may be better characterized by a different speed of maturation rather than shape differences at a given age. In the context of anthropology, we assess the differences in the typical growth of the endocranium between chimpanzees and bonobos. We take advantage of this study to show the robustness of the method with respect to change of parameters and perturbation of the age estimates.

Appendices

A Differentiation of the Temporal Shape Regression Criterion

1.1 A.1 Matrix Notations

For the sake of simplicity, we introduce matrix notations: $\mathbf{x}_{0}=\{x_{p}\}_{p=1,\ldots ,N}$ denotes the $3N$ vector which is the concatenation of the coordinates of $N$ vertices of the baseline shape $M_0$. Using the notations of Sect. 3.2.2, we denote the moving points $\mathbf{x}(t)$ (resp. the parameterizing vectors $\varvec{\alpha }(t)$) the $3N$ vector: $(x_p(t))_{p=1,\ldots ,N}$ (resp. $(\alpha _p(t))_{p=1,\ldots ,N}$). We denote also $\mathbf{K}^\chi (\mathbf{x}(t),\mathbf{x}(t))$ the $3N$-by-$3N$ block matrix whose block $p,q$ is given by the 3-by-3 matrix $(K^\chi (x_p(t),x_q(t)))$. This matrix is symmetric, positive definite by definition of the kernel $K^\chi$.

Thanks to these notations, the norm of the speed vector $v_t$ is written: $\left\Vert{v_t}\right\Vert_{V^\chi }^2 = \varvec{\alpha }(t)^t \mathbf{K}^\chi (\mathbf{x}(t),\mathbf{x}(t)) \varvec{\alpha }(t)$. For $A$, a function from $\mathbb R ^3$ to $\mathbb R$, we denote by $d_xA$ its Jacobian matrix at point $x$, so that for any vector $V$: $(d_x A)V = (\nabla _x A)^t V$. By extension, $\nabla _{\mathbf{x}}A$ denotes the $3N$ vector $(\nabla _{x_1}A,\ldots ,\nabla _{x_N}A)$.

With these notations, the regression criterion (17) becomes:

$$\begin{aligned} E\left((\varvec{\alpha }(t))_{t\in [0,T]}\right) = \sum _{t_i} A_i(\mathbf{x}_{t_i}) + \int _0^T L^\chi (\mathbf{x}(t),\varvec{\alpha }(t)) dt \end{aligned}$$

(39)

subject that:

$$\begin{aligned} \frac{d\mathbf{x}(t)}{dt} = f(\mathbf{x}(t),\varvec{\alpha }(t)) \quad \text{ with} \quad \mathbf{x}(0) = \mathbf{x}_0 \end{aligned}$$

(40)

where we denote:

$$\begin{aligned} \begin{aligned}&f(\mathbf{x}(t),\varvec{\alpha }(t)) = \mathbf{K}^\chi (\mathbf{x}(t),\mathbf{x}(t))\varvec{\alpha }(t)\\&L^\chi (\mathbf{x}(t),\varvec{\alpha }(t)) = \gamma _\chi \varvec{\alpha }(t)^t f(\mathbf{x}(t),\varvec{\alpha }(t)) \end{aligned} \end{aligned}$$

(41)

For the sake of simplicity, we will denote in the sequel, $f(t)$ and $L^\chi (t)$ instead of $f(\mathbf{x}(t),\varvec{\alpha }(t))$ and $L^\chi (\mathbf{x}(t),\varvec{\alpha }(t))$.

1.2 A.2 Gradient in a matrix form

Let $\delta E$ be a variation of the criterion $E$ with respect to a variation $\delta \varvec{\alpha }(t)$ of the momenta $\varvec{\alpha }(t)$:

$$\begin{aligned} \delta E&= \sum _{t_i} \left(d_{\mathbf{x}(t_i)} A_i\right)\delta \mathbf{x}(t_i) + \int _0^T (\partial _{\mathbf{x}} L^\chi (t))\delta \mathbf{x}(t)\nonumber \\&\quad + (\partial _{\varvec{\alpha }} L^\chi (t))\delta \varvec{\alpha }(t) dt \end{aligned}$$

(42)

where $\delta \mathbf{x}(t)$ denotes the variations of the positions $\mathbf{x}(t)$ with respect to the variations of the momenta $\varvec{\alpha }(t)$. The differentiation of the flow Eq. (40) shows that these variations $\delta \mathbf{x}(t)$ satisfy a linear ODE with source term:

$$\begin{aligned} \frac{d}{dt}\delta \mathbf{x}(t)&= (\partial _{\mathbf{x}} f(t))\delta \mathbf{x}(t)\nonumber \\&\quad + (\partial _{\varvec{\alpha }} f(t))\delta \varvec{\alpha }(t) \quad \text{ with} \quad \delta \mathbf{x}(0) = 0 \end{aligned}$$

(43)

We introduce the flow $R_{ut}$ for $u,t\in [0,T]$ which is solution of the homogeneous equation:

$$\begin{aligned} \frac{d R_{ut}}{dt} = R_{ut}(\partial _{\mathbf{x}} f(t)) \quad \text{ with} R_{tt} = \text{ Id}\end{aligned}$$

(44)

The method of the variations of the parameters leads to the following solution of the ODE:

$$\begin{aligned} \delta \mathbf{x}(t) = \displaystyle \int _0^t R_{ut}\partial _{\varvec{\alpha }} f(u)\delta \varvec{\alpha }(u) du \end{aligned}$$

(45)

In particular, this shows that we can write the variations $\delta \mathbf{x}(t_i)$ as:

$$\begin{aligned} \delta \mathbf{x}(t_i) = \int _0^T R_{tt_i}\partial _{\varvec{\alpha }} f(t)\delta \varvec{\alpha }(t)\mathbf{1}_{\{t\le t_i\}} dt \end{aligned}$$

(46)

where $\mathbf{1}_{\{t\le t_i\}}=1$ if $t\le t_i$ and 0 otherwise.

Now, we can plug these last two equations into (42). Using Fubini’s theorem, which implies that $\int _0^T \int _0^t F(u,t)dudt = \int _0^T \int _u^T F(u,t)dt du = \int _0^T \int _t^T F(t,u)dudt$ for any $L^2$ function $F(u,t)$, this leads to:

$$\begin{aligned} \delta E&= \int _0^T \Biggl (\partial _{\varvec{\alpha }} L^\chi (t) \nonumber \\&+ \underbrace{\Bigl (\sum _i d_{\mathbf{x}(t_i)}A_i R_{tt_i}\mathbf{1}_{\{t\le t_i\}} + \int _t^T \partial _{\mathbf{x}} L^\chi (u) R_{tu} du\Bigr )}_{\eta (t)^t}\partial _{\varvec{\alpha }} f(t) \Biggr )\nonumber \\&\delta \varvec{\alpha }(t) dt\nonumber \\ \end{aligned}$$

(47)

This gives the gradient of $E$ with respect to the $L^2$ metric as:

$$\begin{aligned} \nabla _{\varvec{\alpha }} E(t) = \partial _{\varvec{\alpha }} L^\chi (t)^t + \partial _{\varvec{\alpha }} f(t)^t\varvec{\eta }(t) \end{aligned}$$

(48)

where we denote the auxiliary variable $\varvec{\eta }(t)$:

$$\begin{aligned} \varvec{\eta }(t)&= \sum _i (R_{tt_i})^t \nabla _{\mathbf{x}(t_i)}A_i \mathbf{1}_{\{t\le t_i\}}\nonumber \\&\quad + \int _t^T (R_{tu})^t\partial _{\mathbf{x}} L^\chi (u)^t du \end{aligned}$$

(49)

The auxiliary variable $\varvec{\eta }(u)$ depends on the flows $R_{ut}$ and therefore satisfies an ODE. To make this ODE explicit, we write the inverse flow $R_{ut}$ in integral form. Noticing that $R_{tu}R_{ut} = \text{ Id}$, we have $\frac{d R_{ut}}{du} = -\partial _{\mathbf{x}}f(u) R_{ut}$, which gives in integral form (noticing that $R_{ut}$ and $f$ commute):

$$\begin{aligned} R_{ut} = \text{ Id}+ \int _u^t R_{st}\partial _{\mathbf{x}} f(s) ds. \end{aligned}$$

(50)

Now, we can plug this equation into the definition of $\varvec{\eta }(t)$ in (49). Writing $R_{tt_i}= \text{ Id}+ \int _t^T R_{ut_i}\partial _{\mathbf{x}}f(u)\mathbf{1}_{\{u\le t_i\}} du$ and noticing that for any $L^2$ function $F(u,s)$, the Fubini’s theorem implies that $\int _t^T \int _t^u F(u,s)dsdu = \int _t^T \int _u^T F(s,u) dsdu$, this leads to:

$$\begin{aligned}&\varvec{\eta }(t) = \sum _i \nabla _{\mathbf{x}(t_i)} A_i \mathbf{1}_{\{t\le t_i\}} + \int _t^T \partial _{\mathbf{x}} L^\chi (u)^t +\nonumber \\&\partial _{\mathbf{x}}f(u)^t \underbrace{ \left( \sum _i (R_{ut_i})^t\nabla _{\mathbf{x}(t_i)}A_i\mathbf{1}_{\{u\le t_i\}}\mathbf{1}_{\{t\le t_i\}} + \int _u^T (R_{us})^t\partial _{\mathbf{x}}L^\chi (s)^t ds\right)}_{(\star )} du\nonumber \\ \end{aligned}$$

(51)

Now, we notice that $t\le u$ within the integral, which implies that $\mathbf{1}_{\{t\le t_i\}}\mathbf{1}_{\{u\le t_i\}} = \mathbf{1}_{\{u\le t_i\}}$. Hence, $(\star )$ is exactly equal to $\varvec{\eta }_u$. Therefore, $\varvec{\eta }_t$ is the solution of the integral equation (integrated upstream in time):

$$\begin{aligned} \varvec{\eta }(t)&= \sum _i \nabla _{\mathbf{x}_{t_i}}A_i\mathbf{1}_{\{t\le t_i\}} + \int _t^T \partial _{\mathbf{x}} L^\chi (u)^t\nonumber \\&\quad + \partial _{\mathbf{x}}f(u)^t\varvec{\eta }(u) du \end{aligned}$$

(52)

1.3 A.3 Gradient in Coordinates

Due to the definition of the functions $f$ and $L^\chi$ in (41), we have:

$$\begin{aligned} \partial _{\mathbf{x}}f = (\partial _1 + \partial _2)(\mathbf{K}^\chi (\mathbf{x},\mathbf{x})\varvec{\alpha }) \qquad \qquad \quad \partial _{\varvec{\alpha }}f = \mathbf{K}^\chi (\mathbf{x},\mathbf{x})\nonumber \\ \partial _{\mathbf{x}}L^\chi = \gamma _\chi \varvec{\alpha }^t\left((\partial _1 + \partial _2)\mathbf{K}^\chi (\mathbf{x},\mathbf{x})\varvec{\alpha }\right)\quad \partial _{\varvec{\alpha }}L^\chi = 2\gamma _\chi \varvec{\alpha }^t\mathbf{K}^\chi (\mathbf{x},\mathbf{x}) \nonumber \\ \end{aligned}$$

(53)

Therefore, the gradient of the regression criterion with respect to the $L^2$ metric given in (48) is now equal to:

$$\begin{aligned} \nabla _{\varvec{\alpha }} E(t) = \mathbf{K}^\chi (\mathbf{x}(t),\mathbf{x}(t))\left(2\gamma _\chi \varvec{\alpha }(t) + \varvec{\eta }(t)\right). \end{aligned}$$

The matrix $\mathbf{K}^\chi (\mathbf{x}(t),\mathbf{x}(t))$ is precisely the Sobolev metric induced by the kernel on the set of $L^2$ functions (see Sect. 3.2.2), so that the gradient with respect to this metric is given by:

$$\begin{aligned} \nabla _{\alpha _p} E(t) = 2\gamma _\chi \alpha _p(t) + \eta _p(t) \end{aligned}$$

(54)

The auxiliary variable $\varvec{\eta }(t)$ satisfies the ODE (52), now written as:

$$\begin{aligned} \varvec{\eta }(t)&= \sum _i \nabla _{\mathbf{x}_{t_i}}A_i\mathbf{1}_{\{t\le t_i\}}+ \int _t^T \Bigl ((\partial _1 + \partial _2)\mathbf{K}^\chi (\mathbf{x}(u),\nonumber \\&\mathbf{x}(u))\varvec{\alpha }(u)\Bigr )^t\left(\gamma _\chi \varvec{\alpha }(u) + \varvec{\eta }(u)\right) du \end{aligned}$$

(55)

The $3N$ vector $\nabla _{\mathbf{x}_{t_i}}A_i$ is equal to $(\nabla _{x_1(t_i)} A_i,\ldots ,\nabla _{x_N(t_i)} A_i)$. For generic $3N$ vectors $\mathbf{x},\,\mathbf{y}$ and $\varvec{\alpha }$, the $k$ th coordinate of the $3N$-vector $\mathbf{K}^\chi (\mathbf{x},\mathbf{y})\varvec{\alpha }$ is given as: $(\mathbf{K}^\chi (\mathbf{x},\mathbf{y})\varvec{\alpha })_k = \sum _{p=1}^N K^\chi (x_k,y_p)\alpha _p$. The kernel $K^\chi$ is scalar, namely of the form $K^\chi (x,y)=k^\chi (x,y)\text{ Id}$ for a scalar function $k^\chi$. We have therefore for every $i,j=1,\ldots ,N$:

$$\begin{aligned} \begin{array}{l} \partial _{x_i}\left(\mathbf{K}^\chi (\mathbf{x},\mathbf{y})\varvec{\alpha }\right)_j = \sum \limits _{p=1}^N \alpha _p \left(\nabla _1 k^\chi (x_i,y_p)\right)^t \delta (i-j)\\ \partial _{y_i}\left(\mathbf{K}^\chi (\mathbf{x},\mathbf{y})\varvec{\alpha }\right)_j = \alpha _i\left(\nabla _2 k^\chi (x_j,y_i)\right)^t \end{array} \end{aligned}$$

(56)

Therefore, for a generic $3N$-vector $\varvec{\beta }$, we have:

$$\begin{aligned} \left((\partial _1 + \partial _2)\left(\mathbf{K}^\chi (\mathbf{x},\mathbf{y})\varvec{\alpha }\right)^t\varvec{\beta }\right)_k&= \sum \limits _{p=1}^N \alpha _p^t\beta _k\nabla _1 k^\chi (x_k,y_p)\nonumber \\&+ \alpha _k^t\beta _p\nabla _2k^\chi (x_p,y_k) \end{aligned}$$

(57)

Now, we can apply this equation with $\mathbf{y}= \mathbf{x}$ and $\varvec{\beta }= \gamma _\chi \varvec{\alpha }+ \varvec{\eta }$ and combine it with (55). Noticing that for a symmetric kernel, we have $\nabla _1 k(x,y) = \nabla _2 k(y,x)$, we get eventually the set of ODEs satisfied by the functions $\eta _p(t)$ as given in (19).

B Differentiation of the Spatiotemporal Matching Criterion

1.1 B.1 Matrix Notations

Let $\mathbf{t}_0 = \{t_j\}_{j=1,\ldots ,N_\mathrm{target}}$ be the vector of time-points associated to the target shapes. The 1D diffeomorphism $\psi _u$ changes $\mathbf{t}_0$ into $\mathbf{t}(u) = \{t_j(u)\}_{j=1,\ldots ,N_\mathrm{target}}$ for $u\in [0,1]$. This vector satisfies the ODE: $\tfrac{d \mathbf{t}}{du}(u) = \mathbf{K}^\psi (\varvec{\beta }(u),\varvec{\beta }(u))\mathbf{t}(u)$ with $\mathbf{t}(0)=\mathbf{t}_0$, where $\varvec{\beta }(u)$ is the concatenation of the vectors $\beta _j(u)$ defined in (23), $\mathbf{K}^\psi (\varvec{\beta }(u),\varvec{\beta }(u))$ is the block matrix whose block $(i,j)$ is given by: $K^\psi (\beta _i(u),\beta _j(u))$.

Similarly, we denote $\mathbf{x}_0(t) = \{x_p(t)\}_{p=1,\ldots ,N}$ be the concatenation of the positions of all the points of the source evolution $S(t)$ for any time-point $t$ and $\mathbf{x}_0(\mathbf{t}(1))$ the concatenation of the $\mathbf{x}_0(t_j(1))$ for $j=1,\ldots ,N_\mathrm{target}$. The diffeomorphism $\phi _u$ maps this vector to $\mathbf{x}(u)$, which satisfies the ODE: $\frac{d\mathbf{x}}{du} = \mathbf{K}^\phi (\mathbf{x}(u),\mathbf{x}(u))\varvec{\alpha }(u)$ with initial condition: $\mathbf{x}(0) = \mathbf{x}_0(\mathbf{t}(1))$ (which depends on the final time-points $\mathbf{t}(1)$), where $\varvec{\alpha }(u)$ is the concatenation of the vectors $\alpha _{p,j}(u)$ defined in (22) for $p=1,\ldots ,N$ and $j=1,\ldots ,N_\mathrm{target}$.

Therefore, we can write the matching criterion (4) as:

$$\begin{aligned} E(\varvec{\alpha }(u),\varvec{\beta }(u))&= A(\mathbf{x}(1)) + \int _0^1 L^\phi \left(\mathbf{x}(u),\varvec{\alpha }(u)\right) du\nonumber \\&+ \int _0^1 L^\psi \left(\mathbf{t}(u),\varvec{\beta }(u)\right) du \end{aligned}$$

(58)

subject to:

$$\begin{aligned} \left\{ \displaystyle \begin{array}{lll} \dfrac{d\mathbf{x}(u)}{du} = f(\mathbf{x}(u),\varvec{\alpha }(u))&\text{ with}&\mathbf{x}(0) = \mathbf{x}_0(\mathbf{t}(1))\\ \dfrac{d\mathbf{t}(u)}{du} = g(\mathbf{t}(u),\varvec{\beta }(u))&\text{ with}&\mathbf{t}(0) = \mathbf{t}_0 \end{array}\right. \end{aligned}$$

(59)

where we denote:

$$\begin{aligned} \begin{aligned}&f(\mathbf{x}(u),\varvec{\alpha }(u)) = \mathbf{K}^\phi (\mathbf{x}(u),\mathbf{x}(u))\varvec{\alpha }(u)\\&g(\mathbf{t}(u),\varvec{\beta }(u)) = \mathbf{K}^\psi (\mathbf{t}(u),\mathbf{t}(u))\varvec{\beta }(u)\\&L^\phi (\mathbf{x}(u),\varvec{\alpha }(u)) = \gamma _\phi \varvec{\alpha }(u)^t f(\mathbf{x}(u),\varvec{\alpha }(u))\\&L^\psi (\mathbf{t}(u),\varvec{\beta }(u)) = \gamma _\psi \varvec{\beta }(u)^t g(\mathbf{t}(u),\varvec{\beta }(u)) \end{aligned} \end{aligned}$$

(60)

For the sake of simplicity, we will write in the sequel $f(u),\,g(u),\,L^\phi (u)$ and $L^\psi (u)$ instead of $f(\mathbf{x}(u),\varvec{\alpha }(u))$, $g(\mathbf{t}(u),\varvec{\beta }(u))$, $L^\phi (\mathbf{x}(u),\varvec{\alpha }(u))$ and $L^\psi (\mathbf{t}(u),\varvec{\beta }(u))$ respectively.

1.2 B.2 Gradient in a Matrix Form

Now, let $\delta E$ be a variation of the criterion $E$ induced by a variation of the momenta $\delta \varvec{\alpha }(u)$ and $\delta \varvec{\beta }(u)$:

$$\begin{aligned} \delta E&= (d_{\mathbf{x}(1)}A)\delta \mathbf{x}(1) + \int _0^1 (\partial _\mathbf{x}L^\phi (u))\delta \mathbf{x}(u)\nonumber \\&+ ( \partial _{\varvec{\alpha }} L^\phi (u) )\delta \varvec{\alpha }(u)+ (\partial _\mathbf{t}L^\psi (u))\delta \mathbf{t}(u)\nonumber \\&+ (\partial _{\varvec{\beta }} L^\psi (u))\delta \varvec{\beta }(u) du \end{aligned}$$

(61)

where we denote $\delta \mathbf{x}(u)$ and $\delta \mathbf{t}(u)$ the variations of the path $\mathbf{x}(u)$ and $\mathbf{t}(u)$ induced by the variations of the momenta $\delta \varvec{\alpha }(u)$ and $\delta \varvec{\beta }(u)$. These vectors satisfy the linear ODEs with source term derived from (59):

$$\begin{aligned} \dfrac{d}{du}\delta \mathbf{x}(u)&= (\partial _\mathbf{x}f(u))\delta \mathbf{x}(u)\nonumber \\&+ (\partial _{\varvec{\alpha }} f(u))\delta \varvec{\alpha }(u) \quad \text{ with} \quad \delta \mathbf{x}(0) = \delta \mathbf{x}_0(\mathbf{t}(1))\nonumber \\ \dfrac{d}{du}\delta \mathbf{t}(u)&= (\partial _\mathbf{t}g(u))\delta \mathbf{t}(u)\nonumber \\&+ (\partial _{\varvec{\beta }} g(u))\delta \varvec{\beta }(u) \quad \text{ with} \quad \delta \mathbf{t}(0) = 0 \end{aligned}$$

(62)

We introduce the flows $P_{su}$ and $R_{su}$ for all $s,u\in [0,1]$, which are solution of the homogeneous equations:

$$\begin{aligned} \begin{array}{lll} \dfrac{d}{du}P_{su} = P_{su}(\partial _\mathbf{x}f(u))&\text{ with}&P_{uu} = \text{ Id}\\ \dfrac{d}{du}R_{su} = R_{su}(\partial _\mathbf{t}g(u))&\text{ with}&R_{uu} = \text{ Id}\end{array} \end{aligned}$$

(63)

The method of variations of the parameters leads to the following solution of the ODEs:

$$\begin{aligned} \begin{array}{l} \delta \mathbf{x}(u) = P_{0u}\delta \mathbf{x}(0) + \displaystyle \int _0^u P_{su}\partial _{\varvec{\alpha }} f(s)\delta \varvec{\alpha }(s) ds\\ \delta \mathbf{t}(u) = \displaystyle \int _0^u R_{su}\partial _{\varvec{\beta }} g(s)\delta \varvec{\beta }(s)ds \end{array} \end{aligned}$$

(64)

where the variations of the initial condition $\delta \mathbf{x}(0) = \delta \mathbf{x}_0(\mathbf{t}(1))$ equals:

$$\begin{aligned} \delta \mathbf{x}(0) \!= \!(d_{\mathbf{t}(1)}\mathbf{x}_0)\delta \mathbf{t}(1) \!=\! (d_{\mathbf{t}(1)}\mathbf{x}_0)\int _0^1 R_{u1}\partial _{\varvec{\beta }} g(u)\delta \varvec{\beta }(u)du,\nonumber \\ \end{aligned}$$

(65)

according to (64).

Plugging (64) into (61) leads to the variation of the criterion (noticing that for any $L^2$ function $F(s,u)$ we have that $\int _0^1\int _0^u F(s,u)dsdu = \int _0^1\int _s^1 F(s,u)duds = \int _0^1\int _u^1 F(u,s)dsdu$):

$$\begin{aligned} \delta E&= \Bigl ( \underbrace{ (d_{\mathbf{x}(1)}A)P_{01} + \int _0^1 \partial _\mathbf{x}L^\phi (u)P_{0u}du}_{\varvec{\eta }(0)^t} \Bigr ) \delta \mathbf{x}(0)\nonumber \\&+ \int _0^1 \Bigl ( \partial _{\varvec{\alpha }} L^\phi (u) \nonumber \\&+ \Bigl ( \underbrace{(d_{\mathbf{x}(1)}A)P_{u1} + \int _u^1 \partial _\mathbf{x}L^\phi (s)P_{us}ds}_{\varvec{\eta }(u)^t} \Bigr )\partial _{\varvec{\alpha }} f(u) \Bigr )\nonumber \\&\delta \varvec{\alpha }(u) du+ \int _0^1 \Bigl ( \partial _{\varvec{\beta }} L^\psi (u) + \int _u^1 \partial _\mathbf{t}L^\psi (s)R_{us}ds\nonumber \\&\partial _{\varvec{\beta }}g(u) \Bigr ) \delta \varvec{\beta }(u) du \end{aligned}$$

(66)

Now, we denote,

$$\begin{aligned} \varvec{\eta }(u)^t = (d_{\mathbf{x}(1)}A) P_{u1} + \int _u^1 \partial _\mathbf{x}L^\phi (s) P_{us}ds \end{aligned}$$

(67)

which appears twice in (66) as $\varvec{\eta }(0)$ and $\varvec{\eta }(u)$. Given the expression of $\delta \mathbf{x}(0)$ in (65), we have:

$$\begin{aligned} \delta E&= \int _0^1 \left( \partial _{\varvec{\alpha }} L^\phi (u) + \eta (u)^t\partial _{\varvec{\alpha }} f(u) \right)\nonumber \\&\delta \varvec{\alpha }(u)du + \int _0^1 \Biggl ( \partial _{\varvec{\beta }} L^\psi (u)\nonumber \\&\quad + \Bigl ( \underbrace{\varvec{\eta }(0)^t (d_{\mathbf{t}(1)}\mathbf{x}_0) R_{u1} + \int _u^1 \partial _\mathbf{t}L^\psi (s)R_{us}ds}_{\varvec{\xi }(u)^t} \Bigr )\partial _{\varvec{\beta }} g(u) \Biggr )\nonumber \\&\delta \varvec{\beta }(u)du \end{aligned}$$

(68)

Denoting

$$\begin{aligned} \varvec{\xi }(u)^t = \varvec{\eta }(0)^t (d_{\mathbf{t}(1)}\mathbf{x}_0) R_{u1} + \int _u^1 \partial _\mathbf{t}L^\psi (s)R_{us}ds, \end{aligned}$$

(69)

we end up with the gradient of the criterion with respect to the $L^2$ metric written as:

$$\begin{aligned} \left\{ \begin{array}{l} \nabla _{\varvec{\alpha }} E(u) = \partial _{\varvec{\alpha }} L^\phi (u)^t + \partial _{\varvec{\alpha }} f(u)^t\varvec{\eta }(u)\\ \nabla _{\varvec{\beta }} E(u) = \partial _{\varvec{\beta }} L^\psi (u)^t + \partial _{\varvec{\beta }} g(u)^t\varvec{\xi }(u) \end{array}\right. \end{aligned}$$

(70)

The auxiliary variables $\varvec{\eta }(u)$ and $\varvec{\xi }(u)$ depend on the flows $R_{us}$ and $P_{us}$. Therefore they satisfy a ODE, which we need to make explicit now. The inverse flows are written in integral form as:

$$\begin{aligned} P_{us} \!= \!\text{ Id}\!+\! \int _u^s P_{rs}\partial _\mathbf{x}f(r) dr \quad R_{us} \!=\! \text{ Id}+ \int _u^s R_{rs}\partial _\mathbf{t}g(r)dr,\nonumber \\ \end{aligned}$$

(71)

so that the auxiliary variable $\varvec{\eta }(u)$ satisfies:

$$\begin{aligned} \varvec{\eta }(u)&= \nabla _{\mathbf{x}(1)}A + \int _u^1 \partial _\mathbf{x}L^\phi (s)^t ds + \int _u^1 \partial _\mathbf{x}f(s)^t(P_{s1})^t\nonumber \\&\times (\nabla _{\mathbf{x}(1)}A) ds\!+\! \int _u^1\int _u^s (\partial _\mathbf{x}f(r))^t(P_{rs})^t(\partial _\mathbf{x}L^\phi (s))^t dr ds\nonumber \\ \end{aligned}$$

(72)

where we denote $\nabla _{\mathbf{x}}A = \left(d_{\mathbf{x}}A\right)^t$ for any scalar function $A$.

Since we have for any $L^2$ functions $F(r,s)$, $\int _u^1\int _u^sF(r,s) drds = \int _u^1\int _s^1 F(s,r) dr ds$ by permuting the two integrals, we have:

$$\begin{aligned} \varvec{\eta }(u)&= \nabla _{\mathbf{x}(1)}A + \int _u^1 \partial _\mathbf{x}L^\phi (s)^t +\partial _\mathbf{x}f(s)^t\nonumber \\&\Bigl ( \underbrace{ (P_{s1})^t\nabla _{\mathbf{x}(1)}A + \int _s^1 (P_{sr})^t(\partial _{\mathbf{x}} L^\phi (r))^t dr }_{\varvec{\eta }(s)} \Bigr ) ds \end{aligned}$$

(73)

The term in the brackets is exactly $\varvec{\eta }(s)$, so that the integral equation satisfied by $\varvec{\eta }(u)$ is eventually given by:

$$\begin{aligned} \varvec{\eta }(u) = \nabla _{\mathbf{x}(1)}A + \int _u^1 (\partial _\mathbf{x}L^\phi (s))^t + (\partial _\mathbf{x}f(s))^t\varvec{\eta }(s)ds. \end{aligned}$$

(74)

Similar computations using the integral form of the flow $R_{us}$ leads to the integral equation satisfied by $\varvec{\xi }(u)$:

$$\begin{aligned} \varvec{\xi }(u) = (d_{\mathbf{t}(1)}\mathbf{x}_0)^t\varvec{\eta }(0) + \int _u^1 (\partial _\mathbf{t}L^\psi (s))^t + (\partial _\mathbf{t}g(s))^t\varvec{\xi }(s)ds.\nonumber \\ \end{aligned}$$

(75)

1.3 B.3 Gradient in Coordinates

Given the definition of the functions $f,\,g,\,L^\phi$ and $L^\psi$, we have:

$$\begin{aligned} \begin{aligned}&\partial _{\mathbf{x}}f = (\partial _1 + \partial _2) (\mathbf{K}^\phi (\mathbf{x},\mathbf{x})\varvec{\alpha }) \,\,\,\quad \quad \quad \quad \partial _{\varvec{\alpha }} f = \mathbf{K}^\phi (\mathbf{x},\mathbf{x})\\&\partial _{\mathbf{t}} g = (\partial _1 + \partial _2) (\mathbf{K}^\psi (\mathbf{t},\mathbf{t})\varvec{\beta }) \,\quad \quad \quad \quad \quad \partial _{\varvec{\beta }} g = \mathbf{K}^\psi (\mathbf{t},\mathbf{t})\\&\partial _\mathbf{x}L^\phi = \gamma _\phi \varvec{\alpha }^t \Bigl ( (\partial _1 + \partial _2)\mathbf{K}^\phi (\mathbf{x},\mathbf{x})\alpha \Bigr ) \quad \partial _{\varvec{\alpha }} L^\phi = 2\gamma _\phi \varvec{\alpha }^t\mathbf{K}^\phi (\mathbf{x},\mathbf{x})\\&\partial _\mathbf{t}L^\psi = \gamma _\psi \varvec{\beta }^t\Bigl ( (\partial _1 + \partial _2)\mathbf{K}^\psi (\mathbf{t},\mathbf{t})\varvec{\beta }\Bigr ) \quad \partial _{\varvec{\beta }} L^\psi = 2\gamma _\psi \varvec{\beta }^t\mathbf{K}^\psi (\mathbf{t},\mathbf{t}) \end{aligned} \end{aligned}$$

(76)

so that the gradient with respect to the Sobolev metric (the matrices $\mathbf{K}^\phi (\mathbf{x},\mathbf{x})$ and $\mathbf{K}^\psi (\mathbf{t},\mathbf{t})$ factorize in (70)) is given as:

$$\begin{aligned} \left\{ \begin{array}{l} \nabla _{\varvec{\alpha }} E(u) = 2\gamma _\phi \varvec{\alpha }(u) + \varvec{\eta }(u)\\ \nabla _{\varvec{\beta }} E(u) = 2\gamma _\psi \varvec{\beta }(u) + \varvec{\xi }(u) \end{array}\right. \end{aligned}$$

(77)

where

$$\begin{aligned}&\varvec{\eta }(u) = \nabla _{\mathbf{x}(1)}A + \nonumber \\&\displaystyle \int _u^1 \Bigl ( (\partial _1 + \partial _2)\mathbf{K}^\phi (\mathbf{x}(s),\mathbf{x}(s))\alpha (s) \Bigr )^t(\gamma _\phi \varvec{\alpha }(s) + \varvec{\eta }(s))ds\nonumber \\ \end{aligned}$$

(78)

and

$$\begin{aligned} \varvec{\xi }(u)&= (d_{\mathbf{t}(1)}\mathbf{x}_0)^t\varvec{\eta }(0) + \nonumber \\&\displaystyle \int _u^1 \Bigl ( (\partial _1 + \partial _2)\mathbf{K}^\psi (\mathbf{t}(s),\mathbf{t}(s))\varvec{\beta }(s) \Bigr )^t(\gamma _\psi \varvec{\beta }(s) + \varvec{\xi }(s))\nonumber \\ \end{aligned}$$

(79)

The $3NN_\mathrm{target}$ vector $\nabla _{\mathbf{x}(1)}A$ is the concatenation of the vectors $\nabla _{x_p(t_j(1))}A$ for $p=1,\ldots ,N$ and $j=1,\ldots ,N_\mathrm{target}$. Similarly, the $3NN_\mathrm{target}$ vector is the concatenation of the vectors $\{x_p(t_j(1))\}_{p=1,\ldots ,N,j=1,\ldots ,N_\mathrm{target}}$. $d_{\mathbf{t}(1)}\mathbf{x}_0$ is the $3NN_\mathrm{target}$-by-$N_\mathrm{target}$ matrix: $(d_{t_1(1)}\mathbf{x}_0,\ldots ,d_{t_{N_\mathrm{target}}}\mathbf{x}_0)$. In the vector $d_{t_j(1)}\mathbf{x}_0$ almost every coordinate vanishes except the ones corresponding at the jth block of size $3N$: $(d_{t_j(1)}x_1(t_j(1)),\ldots ,d_{t_j(1)}x_N(t_j(1)))$ (since $d_{t_j(1)}x_p (t_i(1))\!=\!0$ when $i\ne j$). Therefore, we have: $(d_{t_j(1)}\mathbf{x}_0)^t\eta = \sum _{p=1}^N \left(\frac{d x_p}{dt}(t_j(1))\right)^t\eta _{p,j}$, which is the jth coordinate of the $N_\mathrm{target}$ vector $(d_{\mathbf{t}(1)}\mathbf{x}_0)^t\eta$.

Eventually, using the generic expression (57) for scalar kernels $K^\phi (x,y)=k^\phi (x,y)\text{ Id}$ and $K^\psi (x,y) = k^\psi (x,y)\text{ Id}$, the evolution of $\varvec{\eta }(u)$ and $\varvec{\xi }(u)$ in (78) are written in coordinates as in (28) and (29).

C Algorithms