Marco Lorenzi

Xavier Pennec

Discrete Ladders for Parallel Transport in Transformation Groups with an Affine Connection Structure

Introduction

The analysis of complex information in medical imaging and in computer vision often requires to represent data in suitable manifolds in which we need to compute trajectories, distances, means and statistical modes. An important example is found in computational anatomy, which aims at developing statistical models of the anatomical variability of organs and tissues. Following D'Arcy Thompson [START_REF] Wentworth | On growth and form by D'Arcy Wentworth Thompson[END_REF], we can assume that the variability of a given observation of the population is encoded by a spatial deformation of a template shape or image (called an atlas in computational anatomy). The analysis of deformations thus enables the understanding of phenotypic variation and morphological traits in populations.

Deformations of images and shapes are usually estimated by image registration, which goal is to estimate the geometric transformation that best superimposes two or more images. Among the numerous methods used for registering medical images, the rich mathematical setting of diffeomorphic non-linear registration is particularly appealing since it provides elegant and grounded methods for atlas building [START_REF] Joshi | Landmark matching via large deformation diffeomorphisms[END_REF], group-wise [START_REF] Bossa | Contributions to 3D diffeomorphic atlas estimation: Application to brain images[END_REF], and longitudinal statistical analysis of deformations [START_REF] Avants | Spatiotemporal normalization for longitudinal analysis of gray matter atrophy in frontotemporal dementia[END_REF][START_REF] Davis | Population shape regression from random design data[END_REF][START_REF] Durrleman | Spatiotemporal atlas estimation for developmental delay detection in longitudinal datasets[END_REF]. In particular, temporal evolutions of anatomies can be modeled by transformation trajectories in the space of diffeomorphisms. However, developing population-based models requires reliable methods for comparing different trajectories.

Among the different techniques proposed so far [START_REF] Rao | Spatial trasformation of motion and deformation fields using nonrigid registration[END_REF][START_REF] Bossa | On changing coordinate systems for longitudinal tensor-based morphometry[END_REF][START_REF] Durrleman | Spatiotemporal atlas estimation for developmental delay detection in longitudinal datasets[END_REF], parallel transport represents a promising method which relies on a solid mathematical background. At the infinitesimal level, a trajectory is a tangent vector to a transformation. Parallel transport consists in transporting the infinitesimal deformation vector across the manifold by preserving its properties with respect to the space geometry. It is one of the fundamental operations of differential geometry which enables to compare tangent vectors, and thus the underlying trajectories, across the whole manifold.

For this reason, parallel transport in transformation groups is currently an important field of research with applications in medical imaging to the development of spatio-temporal atlases for brain images [START_REF] Lorenzi | Mapping the effects of Aβ 1-42 levels on the longitudinal changes in healthy aging: hierarchical modeling based on stationary velocity fields[END_REF], the study of hippocampal shapes [START_REF] Qiu | Parallel transport in diffeomorphisms distinguish the time-dependent pattern of hippocampal surface deformation due to healthy aging and dementia of the Alzheimer's type[END_REF], or cardiac motion [START_REF] Ardekani | Cardiac motion analysis in ischemic and non-ischemic cardiomyopathy using parallel transport[END_REF]. In computer vision one also finds applications to motion tracking and more generally to statistical analysis [START_REF] Wei | Learning deformations with parallel transport[END_REF][START_REF] Subbarao | Rutgers The State University of New Jersey New Brunswick[END_REF][START_REF] Twining | Metrics, connections, and correspondence: the setting for groupwise shape analysis[END_REF][START_REF] Charpiat | Learning Shape Metrics based on Deformations and Transport[END_REF].

Even though the notion of parallel transport comes in somehow intuitively, its practical implementation requires the precise knowledge of the space geometry, in particular of the underlying connection. This is not always easy, especially in infinite dimensions such as in the setting of diffeomorphic image registration. Moreover, parallel transport is a continuous operation involving the computation of (covariant) derivatives which, from the practical point of view, might lead to numerical issues concerning numerical stability and robustness. These issues are related to the unavoidable approximations arising when continuous energy functionals and operators are discretized on grids, especially concerning the evaluation of derivatives through finite difference schemes.

The complexity and limitations deriving from the direct computation of continuous parallel transport methods can be alleviated by considering discrete approximations. In the '70s of the past century [START_REF] Misner | Gravitation[END_REF] proposed a scheme for performing the parallel transport with a very simple geometrical constructions. This scheme was called Schild's Ladder since it was in the spirit of the work of the theoretical physicist Alfred Schild's. The computational interest of Schild's ladder resides in its generality, since it enables the transport of vectors in manifolds by computing geodesics only. This way, the implementation of covariant derivatives is not required anymore, and we can concentrate the implementation effort on the geodesics only.

We recently showed that numerical schemes derived from the Schild's ladder can be effectively applied in the setting of diffeomorphic image registration, by appropriately taking advantage of the underlying geometrical setting [START_REF] Lorenzi | Schild's Ladder for the parallel transport of deformations in time series of images[END_REF][START_REF] Lorenzi | Geodesics, parallel transport & one-parameter subgroups for diffeomorphic image registration[END_REF]. Based on this experience, we believe that discrete transport methods represent promising and powerful techniques for the analysis of transformations due to their simplicity and generality. Bearing in mind the applicative context of the development of transport techniques, this chapter aims to illustrate the principles of discrete schemes for parallel transport in smooth groups equipped with affine connection or Riemannian structures.

The chapter is structured as follows. In Section 2 we provide fundamental notions of finite-dimensional Lie Groups and Riemannian geometry concerning affine connection and covariant derivative. These notions are the basis for the continuous parallel transport methods defined in Section 3. In Section 4 we introduce the Schild's ladder. After detailing its construction and mathematical properties, we derive from it the more efficient Pole ladder construction in Section 4.2. These theoretical concepts are then contextualized and discussed in the applicative setting of diffeomorphic image registration, in which some limitations arise when considering infinite dimensional Lie groups (Section 5). Finally, after illustrating numerical implementations of the aforementioned discrete transport methods, we show in Section 6 their effectiveness when applied to the very practical problem of statistical analysis of the group-wise longitudinal brain changes in Alzheimer's disease.

This chapter summarizes and contextualizes a series of contributions concerning the theoretical foundations of diffeomorphic registration parametrized by stationary velocity fields (SVFs) and discrete transport methods. The application of Schild's ladder in the context of image registration was first presented in [START_REF] Lorenzi | Schild's Ladder for the parallel transport of deformations in time series of images[END_REF]. In [START_REF] Lorenzi | Geodesics, parallel transport & one-parameter subgroups for diffeomorphic image registration[END_REF], we highlighted the geometrical basis of the SVFs registration setting by investigating its theoretical foundations and its connections with Lie group theory and affine geometry. These insights allowed to define a novel discrete transport scheme, the pole ladder, which optimizes the Schild ladder by taking advantage of the geometrical properties of the SVF setting [START_REF] Lorenzi | Efficient parallel transport of deformations in time series of images: from Schild's to pole ladder[END_REF].

Basics of Lie Groups

We recall here the theoretical notions of Lie group theory and affine geometry, that will be extensively used in the following sections.

A Lie group G is a smooth manifold provided with an identity element id, a smooth associative composition rule (g, h) ∈ G × G → gh ∈ G and a smooth inversion rule g → g -1 which are both compatible with the differential manifold structure. As such, we have a tangent space T g G at each point g ∈ G. A vector field X is a smooth function that maps a tangent vector X| g to each point g of the manifold. The set of vector fields (the tangent bundle) is denoted T G. Vector fields can be viewed as the directional (or Lie) derivative of a scalar function φ along the vector field at each point:

∂ X φ| g = d dt φ(ϕ t )| t=0
, where ϕ t is the flow of X and ϕ 0 = g. Composing directional derivatives ∂ X ∂ Y φ leads in general to a second order derivation. However, we can remove the second order terms by subtracting ∂ Y ∂ X φ (this can be checked by writing these expression in a local coordinate system). We obtain the Lie bracket that acts as an internal multiplication in the algebra of vector fields:

[X, Y](φ) = ∂ X ∂ Y φ -∂ Y ∂ X φ.
Given a group element a ∈ G, we call left translation L a the composition with the fixed element a on the left: La : g ∈ G → ag ∈ G. The differential DL a of the left translation maps the tangent space T g G to the tangent space T ag G. We say that a vector field X ∈ T (G) is left invariant if it remains unchanged under the action of the left translation: DL a X|g = X| ag . The sub-algebra of left-invariant vector fields is closed under the Lie bracket and is called the Lie algebra g of the Lie group. Since a left-invariant vector field is uniquely determined by its value at identity through the one-to-one map X| g = DL g X, the Lie algebra can be identified to the tangent space at the identity T id G. One should notice that any smooth vector field can be written as a linear combination of left-invariant vector fields with smooth functional coefficients.

Left-invariant vector fields are complete in the sense that their flow ϕ t is defined for all time. Moreover, this flow is such that ϕ t (g) = gϕ t (id) by left invariance. The map X → ϕ 1 (id) of g into G is called Lie group exponential and denoted by exp. In particular, the group exponential defines the one-parameter subgroup associated to the vector X and has the following properties:

• ϕ t (id) = exp(tX), for each t ∈ R;

• exp((t + s)X) = exp(tX) exp(sX), for each t, s ∈ R.

In finite dimension, it can be shown that the Lie group exponential is a diffeomorphism from a neighborhood of 0 in g to a neighborhood of id in G.

For each tangent vector X ∈ g, the one parameter subgroup exp(tX) is a curve that starts from identity with this tangent vector. One could question if this curve could be seen as a geodesic like in Riemannian manifolds. To answer this question, we first need to define what are geodesics. In a Euclidean space, straight lines are curves which have the same tangent vector at all times. In a manifold, tangent vectors at different times belong to different tangent spaces. When one wants to compare tangent vectors at different points, one needs to define a specific mapping between their tangent spaces: this is the notion of parallel transport. There is generally no way to define globally a linear operator Π h g :

T g G → T h G which is consistent with composition (i.e. Π h g • Π g f = Π h f
). However, specifying the parallel transport for infinitesimal displacements allows integrating along a path, thus resulting in a parallel transport that depend on the path. This specification of the parallel transport for infinitesimal displacements is called the (affine) connection.

Affine Connection Spaces

An affine connection on G is an operator which assigns to each X ∈ T (G) a linear mapping ∇ X : T (G) → T (G) such that, for each vector field X, Y ∈ T (G), and smooth function f, g ∈ C ∞ (G, R)

∇ f X+gY = f ∇ X + g∇ Y (Linearity); (1)
∇ X (f Y) = f ∇ X (Y) + (Xf )Y (Leibniz rule). (2)
The affine connection is therefore a derivation on the tangent space which infinitesimally map them from one tangent plane to another. The connection give rise to two very important geometrical objects: the torsion and curvature tensors. The torsion quantifies the failure to close infinitesimal geodesic parallelograms:

T (X, Y) = ∇ X Y -∇ Y X -[X, Y],
while the curvature measures the local deviation of the space from being flat, and is defined as

R(X, Y)Z = ∇ X ∇ Y Z -∇ Y ∇ X Z -∇ [X,Y] Z.
Once the manifold is provided with a connection, it is possible to generalize the notion of "straight line": a vector field X is parallel along a curve γ(t) if ∇ γ(t) X = 0 for each t. A path γ(t) on G is said to be straight or geodesic if ∇ γ γ = 0, i.e. if its tangent vector remains parallel to itself along the path.

In a local coordinate system, the geodesic equation is a second order differential equation. Thus, given a point p ∈ G and a vector X ∈ T p G, there exist a unique geodesic γ(t, p, X) that passes through p with velocity X at the instant t = 0 [START_REF] Postnikov | Geometry VI: Riemannian Geometry[END_REF]. We define therefore the Affine exponential as the application exp : G × T (G) → G given by exp p (X) = γ(1, p, X).

If, as in the Euclidean case, we want to associate to the straight lines the property of minimizing the distance between points, we need to provide the group G with a Riemannian manifold structure, i.e. with a metric operator g on the tangent space. In this case there is a unique connection, called the Levi Civita connection, which, for each X, Y, Z ∈ T (G):

• Preserves the metric, i.e. the parallel transport along any curve connecting f to g is an isometry:

g(X, Y) g = g(Π f g X, Π f g Y) f .
• Is torsion free:

∇ X Y -∇ Y X = [X, Y],
thus the parallel transport is symmetric with respect to the Lie bracket.

By choosing the Levi Civita connection of a given Riemannian metric, the affine geodesics are the length minimizing paths (i.e. classical Riemannian geodesics). However, given a general affine connection, there may not exist any Riemannian metric for which affine geodesics are length minimizing.

Cartan-Schouten Connections

Given an affine connection ∇ and a vector X on T id G, we can therefore define two curves on G passing through id and having X as tangent vector, one given by the Lie group exponential exp and the other given by the affine exponential exp id . When do they coincide?

The connection ∇ on G is left-invariant if, for each left translation L a (a ∈ G) and any vector fields X and Y, we have ∇ DLaX (DL a Y) = DL a ∇ X (Y). Using two left invariant vector fields X, Ỹ ∈ g generated by the tangent vectors X, Y ∈ T id G, we see that ∇ X Ỹ is itself a left-invariant vector field generated by its value at identity. Since a connection is completely determined by its action on the left-invariant vector fields (we can recover the connection on arbitrary vector fields using Eq. (1,2) from their decomposition on the Lie Algebra), we conclude that each left-invariant connection ∇ is uniquely determined by a product α (symmetric bilinear operator) on

T id G through α(X, Y ) = ∇ X Ỹ id .
Notice that such a product can be uniquely decomposed into a commutative part α ′ = 1 2 (α(X, Y ) + α(Y, X)) and a skew symmetric part α ′′ = 1 2 (α(X, Y ) -α(Y, X)). The symmetric part specifies the geodesics (i.e. the parallel transport of a vector along its own direction) while the skew-symmetric part specifies the torsion which governs the parallel transport of a vector along a transverse direction (the rotation around the direction of the curve if we have a metric connection with torsion).

Following [START_REF] Postnikov | Geometry VI: Riemannian Geometry[END_REF], a left-invariant connection ∇ on a Lie group G is a Cartan-Schouten connection if, for any tangent vector X at the identity, the oneparameter subgroups and the affine geodesics coincide, i.e. exp(tX) = γ(t, id, X) . We can see that a Cartan connection satisfies α(X, X) = 0 or, equivalently, is purely skew-symmetric.

The one-dimensional family of connections generated by α(X, Y ) = λ[X, Y ] obviously satisfy this skew-symmetry condition. Moreover, the connections of this family are also invariant by right translation [START_REF] Pennec | Exponential Barycenters of the Canonical Cartan Connection and Invariant Means on Lie Groups[END_REF], thus invariant by inversion also since they are already left invariant. This make them particularly interesting since they are fully compatible with all the group operations.

In this family, three connections have special curvature or symmetric properties and are called the canonical Cartan-Schouten connections [START_REF] Cartan | On the geometry of the group-manifold of simple and semi-simple groups[END_REF]. The zero curvature connections given by λ = 0, 1 (with torsion T = -[ X, Ỹ] and T = [ X, Ỹ] respectively on left invariant vector fields) are called left and right Cartan connections. The choice of λ = 1/2 leads to average the left and right Cartan connections. It is called the symmetric (or mean) Cartan connection. It is torsion-free, but has curvature R( X, Ỹ) Z = -1 4 [ X, Ỹ], Z . As a summary, the three canonical Cartan connections of a Lie group are (for two left-invariant vector fields):

∇ X Ỹ = 0
Left (Torsion, Flat);

∇ X Ỹ = 1 2 X, Ỹ Symmetric (Torsion-Free, Curved); ∇ X Ỹ = X, Ỹ Right (Torsion, Flat).
Since the three canonical Cartan connections only differ by torsion, they share the same affine geodesics which are the left and right translations of one parameter subgroups. In the following, we call them group geodesics. However, the parallel transport of general vectors along these group geodesics is specific to each connection as we will see below.

Riemannian Setting: Levi Civita Connection

Given a metric < X, Y > on the tangent space at identity of a group, one can propagate this metric to all tangent spaces using left (resp. right) translation to obtain a left-(resp. right-) invariant Riemannian metric on the group. In the left-invariant case we have < DL a X, DL a Y > a =< X, Y > and one can show [START_REF] Kolev | Groupes de Lie et mécanique[END_REF] that the Levi Civita connection is the left-invariant connection generated by the product

α(X, Y ) = 1 2 [X, Y ] - 1 2 (ad * (X, Y ) + ad * (Y, X)),
where the operator ad * is defined by < ad * (Y, X), Z >=< [X, Z], Y > for all X, Y, Z ∈ g. A similar formula can be established for right-invariant metrics using the algebra of right-invariant vector fields .

We clearly see that this left-invariant Levi Civita connection has a symmetric part which make it differ from the Cartan symmetric connection α(X, Y ) =

1 2 [X, Y ].
In fact, the quantity ad * (X, X) specifies the rate at which a left invariant geodesic and a one parameter subgroup starting from the identity with the same tangent vector X deviates from each-other. More generally, the condition ad * (X, X) = 0 for all X ∈ g turns out to be a necessary and sufficient condition to have a bi-invariant metric [START_REF] Postnikov | Geometry VI: Riemannian Geometry[END_REF]. It is important to notice that geodesics of the left-and right-invariant metrics differ in general as there do not exists bi-invariant metrics even for simple groups like the Euclidean motions [START_REF] Pennec | Exponential Barycenters of the Canonical Cartan Connection and Invariant Means on Lie Groups[END_REF]. However, right invariant geodesics can be easily obtained from the left invariant one through inversion: if φ(t) is a left invariant geodesic joining identity to the transformation φ 1 , then φ -1 (t) is a right-invariant geodesic joining identity to φ -1 1 .

Continuous Methods for Parallel Transport

After having introduced the theoretical bases of affine connection spaces, in this section we detail the theoretical relationship between parallel transport of tangent vectors and respectively Cartan-Schouten and Riemannian (Levi Civita) connections .

Cartan-Schouten Parallel Transport

For the left Cartan connection, the unique fields that are covariantly constant are the left-invariant vector fields, and the parallel transport is induced by the differential of the left translation [START_REF] Postnikov | Geometry VI: Riemannian Geometry[END_REF], i.e. Π L : T p G → T q G is defined as

Π L (X) = DL qp -1 X. (3)
One can see that the parallel transport is actually independent of the path, which is due to the fact that the curvature is null: we are in a space with absolute parallelism. Similarly, the right-invariant vector fields are covariantly constant with respect to the right invariant connection only. As above, the parallel transport is given by the differential of the right translation

Π R (X) = DR p -1 q X, (4)
and we have an absolute parallelism as well.

Finally, the parallel transport for the symmetric Cartan connection is given by the infinitesimal alternation of the left and right transports. However, as there is curvature, it depends on the path: it can be shown [START_REF] Helgason | Differential Geometry, Lie groups, and Symmetric Spaces[END_REF] that the parallel transport of X along the geodesic exp(tY ) is:

Π S (X) = DL exp( 1 2 Y ) DR exp( 1 2 Y ) X. (5)

Riemannian Parallel Transport

In the Riemannian setting the parallel transport with the Levi Civita connection can be computed by solving a system of PDEs which locally depend on the associated metric (the interested reader can refer to [START_REF] Carmo | Riemannian Geometry[END_REF] for a more comprehensive description of the parallel transport in Riemannian geometry). Let x i be a local coordinate chart x i with ∂ i = ∂ ∂xi a local basis of the tangent space. The tangent vector to the curve γ is γ = i v i ∂ i . It can be easily shown that a vector Y = i y i ∂ i is parallel transported along γ with respect to the affine connection ∇ iff

∇ γ Y = k   ij Γ k ij v j y k + γ(y k )   ∂ k = 0, (6)
with the Christoffel symbols of the connection being defined from the covariant derivative

∇ ∂i ∂ j = k Γ k ij ∂ k . Let us consider the local expression for the metric tensor g lk =< ∂ ∂x l , ∂ ∂x k >.
Thanks to the compatibility condition, the Christoffel symbols of the Levi Civita connection can be locally expressed via the metric tensor leading to

Γ k ij = 1 2 l ∂ ∂x i g jl + ∂ ∂x j g li + ∂ ∂x l g ij g lk , with (g lk ) = (g lk ) -1 .
Thus, in the Riemannian setting the covariant derivative is uniquely defined by the metric, and the parallel transport thus depends from the path γ and from the local expression of the metric tensor g ij (formula 6).

Discrete Methods for Parallel Transport

In Section 3 we showed that the parallel transport closely depends on the underlying connection, and that thus it assumes very specific formulations depending on the underlying geometry. In this section we introduce discrete methods for the computation of the parallel transport which do not explicitly depend on the connection and only make use of geodesics. Such techniques could be applied more generally when working on arbitrary geodesic spaces.

Schild's Ladder

Schild's ladder is a general method for the parallel transport, introduced in the theory of gravitation in [START_REF] Misner | Gravitation[END_REF] after Schild's similar constructions [START_REF] Schild | Tearing geometry to pieces: More on conformal geometry[END_REF]. The method infinitesimally transports a vector along a given curve through the construction of geodesic parallelograms (Figure 1). The Schild's ladder provides a straightforward method to compute a second order approximation of the parallel transport of a vector along a curve using geodesics only.

Let M a manifold and C a curve parametrized by the parameter τ with ∂C ∂τ | T0 = u, and A ∈ T P0 M , a tangent vector on the curve at the point P 0 = C(0). Let P 1 be a point on the curve relatively close to P 0 , i.e. separated by a sufficiently small parameter value τ . The Schild's ladder computes the parallel transport of A along the curve C as follows:

1. Define a curve on the manifold parametrized by a parameter σ passing through the point P 0 with tangent vector ∂ ∂σ | P0 = A. Chose a point P 2 on the curve separated by P 0 by the value of the parameters σ. The values of the parameters σ and τ should be chosen in order to construct this step of the ladder within a single coordinate neighborhood.

2. Let l be the geodesic connecting P 2 = l(0) and P 1 = l(λ), we choose the "middle point" P 3 = l(λ/2). Now, let us define the geodesic r connecting the starting point P 0 and P 3 parametrized by ρ such that P 3 = r(ρ).

Extending the geodesic at the parameter 2ρ we reach the point P 4 . We can now compute the geodesic curve connecting P 1 and P 4 . The vector A ′ tangent to the curve at the point P 1 is the parallel translation of A along C.

3. If the distance between the points P 0 and P 1 is large, the above construction can be iterated for a sufficient number of small steps.

The algorithmic interest of the Schild's ladder is that it only relies on the computation of geodesics. Although the geodesics on the manifold are not sufficient to recover all the information about the space properties, such as the torsion of the connection, it has been shown that the Schild's ladder implements the parallel transport with respect to the symmetric part of the connection of the space [START_REF] Kheyfets | Schild's Ladder parallel transport for an arbitrary connection[END_REF]. An intuitive view of that point is that the construction of the above diagram is commutative and can be symmetrized with respect to the points P 1 and P 2 . If the original connection is symmetric, then this procedure provides a correct linear approximation of the parallel transport of vectors.

Pole Ladder

We proposed in [START_REF] Lorenzi | Efficient parallel transport of deformations in time series of images: from Schild's to pole ladder[END_REF] a different construction for the parallel transport of vectors based on geodesics parallelograms. If the curve C is geodesic, then it can be itself one of the diagonals and the Schild's ladder can therefore be adapted by requiring the computation of only one new diagonal of the parallelogram. We define in this way a different ladder scheme, that we name "pole ladder" since its geometrical construction recalls the type of ladders with alternating or symmetric steps with respect to a central pole axis which is used on poles.

We now prove that the pole ladder is correctly implementing the parallel transport. In the diagram of Figure 2, the parallel transport of the tangent vector v = Ċ along the geodesic C is specified by the geodesic equation v + Γ k ij v i v j = 0 using the Christoffel symbols Γ k ij (x). In a sufficiently small neighborhood the relationships can be linearized to give

v k (t) = v k (0) -tΓ k ij (x(0))v i (0)v j (0) + O(t 2 ),
and by integrating:

x k (t) = x k (0) + tv k (0) - t 2 2 Γ k ij (x(0))v i (0)v j (0) + O(t 3 ).
By renormalizing the length of the vector v so that C(-1)

= P 0 , C(0) = M and C(1) = Q 0 (and denoting Γ k ij = Γ k ij (M ))
, we obtain the relations:

P 0 k = M k -v k M - 1 2 Γ k ij v i M v j M + O( v 3 ), Q 0 k = M k + v k M - 1 2 Γ k ij v i M v j M + O( v 3 ).
Similarly, we have along the second geodesic:

P 1 k = M k -u k M - 1 2 Γ k ij u i M u j M + O( u 3 ), Q 1 k = M k + u k M - 1 2 Γ k ij u i M u j M + O( u 3 ).
Now, to compute the geodesics joining P 0 to P 1 and Q 0 to Q 1 , we have to use a Taylor expansion of the Christoffel symbols Γ k ij around the point M . In the following, we indicate the coordinate according to which the quantity is derived by the index after a comma:

Γ k ij,a = ∂ a Γ k ij : Γ k ij (P 0 ) = Γ k ij + Γ k ij,a (-v k M - 1 2 Γ k ij v i M v j M ) + 1 2 Γ k ij,ab v a M v b M + O( v 3 ).
However, the Christoffel symbols are multiplied by a term of order O( A 2 ), so that only the first term will be quadratic and all others will be of order 3 with respect to A and v M . Thus, the geodesics joining P 0 to P 1 and Q 0 to Q 1 have equations:

P k 1 = P k 0 + A k - 1 2 Γ k ij A i A j + O(( A + v M ) 3 ), Q k 1 = Q k 0 + B k - 1 2 Γ k ij B i B j + O(( B + v M ) 3 ).
Equating P k 1 in the previous equations gives

u k M + 1 2 Γ k ij u i M u j M = v k M -A k + 1 2 Γ k ij (v i M v j M + A i A j ) + O(( B + v M ) 3 ).
Solving for u as a second order polynomial in v M and A gives

u k = v k M -A k + 1 2 (Γ k ij + Γ k ji )A i v j M + O(( A + v M ) 3 ). Now equating Q k 1 in the previous equations gives B k - 1 2 Γ k ij B i B j = -A k + (Γ k ij + Γ k ji )A i v j M + 1 2 Γ k ij A i A j + O(( A + v M ) 3 ).
Solving for B k as a second order polynomial in v M and A gives:

B k = -A k + (Γ k ij + Γ k ji )A i v j + O(( A + v M ) 3 ). ( 7
)
To verify that this is the correct formula for the parallel transport of A, let us observe that the field A(x) is parallel in the direction of

v j if ∇ V A = 0, i.e. if ∂ v A k + Γ k ij A i v j = 0, which means that A k (x + ǫv) = A k -ǫΓ k ij A i v j + O(ǫ 2 ). If the connection is symmetric, i.e. if Γ k ij = Γ k ji , equation (4.
2) shows that the pole ladder leads to

B k ≃ -A k + 2Γ k ij A i v j .
Thus the pole ladder is realizing the parallel transport for a length ǫ = 2 (remember that our initial geodesic was defined from -1 to 1).

We have thus demonstrated that the vector -B of Figure ( 2) is a second order approximation of the transport of A. In order to optimize the number of time steps we should evaluate the error in equation ( 4.2) at high orders on A and v M . The computation is not straightforward and involves a large number of terms, thus preventing the possibility to synthesize a useful result. However, we believe that the dependency on A is more important that the one on v M , and that we could obtain larger time steps provided that A is sufficiently small.

Diffeomorphic Medical Image Registration

We now describe a practical context in which the previous theoretical insights find useful application. For this purpose, we describe here the link between the theory described in Section 2 and the context of computational anatomy, in particular through the diffeomorphic non-linear registration of time series of images.

Longitudinal and Cross-Sectional Registration Settings

Modeling the temporal evolution of the tissues of the body is an important goal of medical image analysis for understanding the structural changes of organs affected by a pathology, or for studying the physiological growth during the life span. For such purposes we need to analyze and compare the observed anatomical differences between follow-up sequences of anatomical images of different subjects. Non-rigid registration is one of the main instruments for modeling anatomical differences from images. The aim of non-linear registration is to encode the observed structural changes as deformation fields densely represented in the image space, which represent the warping required to match the observed differences. This way, the anatomical changes can be modeled and quantified by analyzing the associated deformations.

We can identify two distinct settings for the application of non-linear registration: longitudinal and cross-sectional. In the former, non-linear registration estimates the deformation field which explains the longitudinal anatomical (intra-subject) changes that usually reflect biological phenomena of interest, like atrophy or growth. In the latter, the deformation field accounts for the anatomical differences between different subjects (inter-subject), in order to match homologous anatomical regions. These two settings are profoundly different: the cross-sectional setting does not involve any physical or mechanical deformations and we might wish to compare different anatomies with different topologies. Moreover, inter-subject deformations are often a scale of magnitude higher than the ones characterizing the usually subtle variations of the longitudinal setting.

In case of group-wise analysis of longitudinal deformations, longitudinal and cross-sectional settings must be integrated in a consistent manner. In fact, the comparison of longitudinal deformations is usually performed after normalizing them in a common reference frame through the inter-subject registration, and the choice of the normalization method might have a deep impact on the following analysis. In order to accurately identify longitudinal deformations in a common reference frame space, a rigorous and reliable normalization procedure need thus to be defined.

Normalization of longitudinal deformations can be done in different ways, depending on the analyzed feature. For instance, the scalar Jacobian determinant of longitudinal deformations represents the associated local volume change, and can be compared by scalar resampling in a common reference frame via intersubject registration. This simple transport of scalar quantities is the basis of the classical deformation/tensor based morphometry techniques [START_REF] Ashburner | Voxel-based morphometry -the methods[END_REF][START_REF] Riddle | Characterizing changes in MR images with color-coded Jacobians[END_REF]. However, transporting the Jacobian determinant is not sufficient to reconstruct a deformation in the Template space.

If we consider vector-values characteristics of deformations instead of scalar quantities, the transport is not uniquely defined anymore. For instance, a simple method of transport consists in reorienting the longitudinal intra-subject displacement vector field by the Jacobian matrix of the subject-to-reference deformation. Another intuitive method was proposed by [START_REF] Rao | Spatial trasformation of motion and deformation fields using nonrigid registration[END_REF] and uses the transformation conjugation (change of coordinate system) in order to compose the longitudinal intra-subject deformation with the subject-to-reference one. As pointed out in [START_REF] Bossa | On changing coordinate systems for longitudinal tensor-based morphometry[END_REF], this practice could potentially introduce variations in the transported deformation and relies on the inverse consistency of the estimated deformations, which can raise numerical problems for large deformations.

In the geometric setting, when we have a Riemannian or affine manifold structure for the space our deformations, one would like to use a normalization which is consistent with the manifold structure. This requirement naturally raise parallel transport as the natural tool for normalizing measurements at different points. In order to elaborate along this idea, we first need to describe the geometric structures on diffeomorphic deformations.

A first formulation of diffeomorphic registration was proposed with the "Large Deformation Diffeomorphic Metric Mapping (LDDMM)" setting [START_REF] Trouvé | Diffeomorphisms groups and pattern matching in image analysis[END_REF][START_REF] Beg | Computing Large Deformation Metric Mappings via Geodesic Flows of Diffeomorphisms[END_REF].

In this framework the images are registered by minimizing the length of the trajectory of transformations in the space of diffeomorphism, once specified an opportune right invariant metric. The solution is the endpoint of the flow of a time-varying velocity field, which is a geodesic parametrized through the Riemannian exponential. The LDDMM deformations are thus Riemannian (metric) geodesics, which are also geodesics of the Levi Civita connection.

Since LDDMM is generally computationally intensive, a different diffeomorphic registration method was later proposed with the stationary velocity field (SVF) setting [3]. In this case the diffeomorphisms are parametrized by stationary velocity fields, in opposition to the time varying velocity fields of the LDDMM framework, through the Lie group exponential. The restriction to stationary velocity fields simplifies the registration problem and provides efficient numerical schemes for the computation of deformations. This time the flow associated to SVFs is a one-parameter subgroup, which is a geodesic with respect to the Cartan-Schouten connections. One-parameter subgroups are generally not metric geodesics, since there do not exist any left and right invariant metric on non-compact and non-commutative groups.

In both the LDDMM and SVF settings, the longitudinal deformation is encoded by the initial tangent velocity field. The transport of longitudinal deformations can be then naturally formulated as the parallel transport of tangent vectors along geodesics according to the underlying connection, i.e. the Levi Civita connection in LDDMM, and the canonical symmetric Cartan-Schouten connection in the SVF setting.

A Glimpse of Lie Group Theory in Infinite Dimension

In Section 2.2, we derived the equivalence of one-parameter subgroups and the affine geodesics of the canonical Cartan connections in a finite dimensional Lie group. In order to use such a framework for diffeomorphisms, we have to generalize the theory to infinite dimensions. However, defining infinite dimensional Lie groups is raising much more difficulties. This is in fact the reason why Lie himself restricted to finite dimensions. The theory was developed since the 70ies and is now an active field of research. We refer the reader to the recent books [START_REF] Khesin | The Geometry of Infinite Dimensional Lie groups[END_REF][START_REF] Younes | Shapes and diffeomorphisms[END_REF] for more details on this theory and to [START_REF] Schmid | Infinite dimensional Lie groups with applications to mathematical physics[END_REF] for a good overview of the problems and applications.

The basic construction scheme is to consider an infinite dimensional manifold endowed with smooth group operations. Such a Lie group is locally diffeomorphic to an infinite-dimensional vector space which can be a Fréchet space (a locally convex space which is complete with respect to a translation invariant distance), a Banach space (where the distance comes from a norm) or a Hilbert space (where the norm is derived from a scalar product). We talk about Fréchet, Banach or Hilbert Lie groups, respectively. Extending differential calculus from R n to Banach and Hilbert spaces is straightforward, but this is not so simple for Fréchet spaces. In particular, the dual of a Fréchet space need not be Fréchet, which means that some extra care must be taken when defining differential forms. Moreover, some important theorems such as the inverse function theorem hold for Banach spaces but not necessarily for Fréchet spaces.

For instance, the set Diff k (M) of C k diffeomorphisms of a compact manifold M is a Banach manifold and the set of Sobolev H s diffeomorphisms Diff s (M) is a Hilbert manifold (if s > dimM/2). However, these are no-classical "Lie groups" since one loses derivatives when differentiating the composition and inversion maps. To obtain the complete smoothness of the composition and inversion maps, one has to go to infinity, but the Banach structure is lost in the process [39, p.12] and we are left with Diff ∞ (M) being only a Fréchet Lie group. Some additional structure can be obtained by considering the sequence of Diff k (M) spaces as a succession of dense inclusions as k goes to infinity: this the Inverse Limit of Banach (ILB)-Lie group setting. Likewise, the succession of dense inclusions of Sobolev H s diffeomorphisms give rise to the Inverse Limit of Hilbert (ILH)-Lie group setting.

As the diffeomorphisms groups considered are Fréchet but not Banach, the usual setting of infinite dimensional Lie groups is the general framework of Fréchet manifolds. This implies that many of the important properties which are true in finite dimension do not hold any more for general infinite dimensional Lie groups [START_REF] Schmid | Infinite-Dimensional Lie Groups and Algebras in Mathematical Physics[END_REF].

First, there is no implicit or inverse function theorem (except Nash-Moser type theorems.) This implies for instance that the log-map (the inverse of the exponential map) may not be smooth even if the differential of the exponential map is the identity.

Second, the exponential map is not in general a diffeomorphism from a neighborhood of zero in the Lie algebra onto a neighborhood of the identity in the group. This means that it cannot be used as a local chart to work on the manifold. For instance in Diff s (M), in every neighborhood of the identity there may exist diffeomorphisms which are not the exponential of an H s vector field. A classical example of the non-surjectivity of the exponential map is the following function in Diff(S 1 ) [START_REF] Milnor | Remarks on infinite-dimensional Lie groups[END_REF]:

f n,ǫ (θ) = θ + π/n + ǫ sin 2 (nθ). ( 8
)
This function can be chosen as close as we want to the identity by opportunely dimensioning ǫ and θ. However, it can be shown that it cannot be reached by any one-parameter subgroup, and therefore the Lie group exponential is not a local diffeomorphisms of Diff(S 1 ). This example is quite instructive and shows that this theoretical problem might actually be a very practical advantage: the norm of the k-th derivative of f n,ǫ is exploding when k is going to infinity, which shows that we would rather want to exclude this type of diffeomorphisms from the group under consideration.

Riemannian Structure and Stationary Velocity Fields

In the Large Deformation Diffeomorphic Metric Mapping (LDDMM) framework [START_REF] Younes | Shapes and diffeomorphisms[END_REF], a different construction is leading to a more restricted subgroup of diffeomorphisms which is more rational from the computational point of view. One first chooses a Hilbert norm on the Lie Algebra which turn it into an admissible Hilbert space. Admissible means that it can be embedded into the space of vector fields which are bounded and vanishing at infinity, as well as all the first order derivatives. Typically, this is a Sobolev norm of a sufficiently high order. Then, one restricts to the subgroup of diffeomorphisms generated by the flow of integrable sequences of such vector fields for a finite time. To provide this group with a Riemannian structure, a right invariant metric is chosen. A first reason for choosing right translation is that it is simply a composition which does not involve a differential operator as for the left translation. A second reason is that the resulting metric on the group is generating an invariant metric on the object space with a right action. One can show that the group provided with this right-invariant metric is a complete metric space [START_REF] Trouvé | Diffeomorphisms groups and pattern matching in image analysis[END_REF][START_REF] Younes | Shapes and diffeomorphisms[END_REF]: the choice of the norm on the Lie algebra is specifying the subgroup of diffeomorphisms which are reachable, i.e. which are at a finite distance.

In the SVF setting, the fact that the flow in an autonomous ODE allows us to generalize efficient algorithms such as the scaling and squaring algorithm: given an approximation exp(δY ) = id + δY for small vector fields δY , the exponential of a SVF Y can be efficiently and simply computed by recursive compositions:

exp(Y ) = exp Y 2 • exp Y 2 = exp Y 2 n 2 n
.

A second algorithm is at the heart of the efficiency of the optimization algorithms with SVFs: the Baker-Campbell-Hausdorff (BCH) formula [START_REF] Bossa | Contributions to 3D diffeomorphic atlas estimation: Application to brain images[END_REF] tells us how to approximate the log of the composition:

BCH(X, δY ) = log(exp(X) • exp(δY )) = X + δY + 1 2 [X, δY ] + 1 12 [X, [X, δY ]] + . . . .
In order to have a well-posed space of deformations, we need to specify on which space is modeled the Lie algebra, as previously. This is the role of the regularization term of the SVF registration algorithms [START_REF] Vercauteren | Symmetric Logdomain diffeomorphic registration: A Demons-based approach[END_REF][START_REF] Hernandez | Registration of anatomical images using paths of diffeomorphisms parameterized with stationary vector field flows[END_REF] or of the spline parametrization of the SVF in [START_REF] Ashburner | A fast diffeomorphic image registration algorithm[END_REF][START_REF] Modat | Log-Euclidean free-form deformation[END_REF]: this restricts the Lie algebra to the subalgebra of sufficiently regular velocity fields. The subgroup of diffeomorphisms considered is then generated by the flow of these stationary velocity fields and their finite composition. So far, the theoretical framework is very similar to the LDDMM setting and we can see that the diffeomorphisms generated by the one-parameter subgroups (the exponential of SVFs) all belong to the group considered in the LDDMM setting, provided that we model the Lie algebra on the same admissible Hilbert space. As in finite dimension, the affine geodesics of the Cartan connections (group geodesics) are metric-free (the Hilbert metric is only used to specify the space on which is modeled the Lie Algebra) and generally differ from the Riemannian geodesics of LDDMM.

It is well known that the subgroup of diffeomorphisms generated by this Lie algebra is significantly larger than what is covered by the group exponential. Indeed, although our affine connection space is geodesically complete (all geodesics can be continued for all time without hitting a boundary), there is no Hopf-Rinow theorem which state that any two points can be joined by a geodesic (metric completeness). Thus, in general, not all the elements of the group G may be reached by the one-parameter subgroups. An example in finite dimension is given by SL(2).

However, this might not necessarily results into a problem in the image registration context since we are not interested in recovering "all" the possible diffeomorphisms, but only those which lead to admissible anatomical transformations. For instance, the diffeomorphism on the circle defined above at Eq. ( 8) cannot be reached by any one-parameter subgroup of S 1 . However, since

lim k→∞ f n,ǫ H k → ∞,
this function is not well behaved from the regularity point of view, which is a critical feature when dealing with image registration.

In practice, we have a spatial discretization of the SVF (and of the deformations) on a grid, and the temporal discretization of the time varying velocity fields by a fixed number of time steps. This intrinsically limits the frequency of the deformation below a kind of "Nyquist" threshold, which prevents these diffeomorphisms to be reached anyway both by the SVF and by the "discrete" LDDMM frameworks. Therefore, it seems more importance to understand the impact of using stationary velocity fields in registration from the practical point of view, than from the theoretical point of view, because we will have necessarily to deal with the unavoidable numerical implementation and relative approximation issues.

Parallel Transport in Diffeomorphic Registration

Continuous and discrete methods for the parallel transport provided in Sections 3,4 can be applied in the diffeomorphic registration setting, once provided the appropriate geometrical context (Section 5). In this section we discuss and illustrate practical implementations of the parallel transport in diffeomorphic registration, with special focus on the application of the ladder schemes exposed in Section 4.

Continuous vs Discrete Transport Methods

As illustrated in abstract form in [START_REF] Arnold | Mathematical methods of classical mechanics[END_REF], parallel transport can be approximated infinitesimally by Jacobi fields. Following this intuition, a computational method for the parallel transport along geodesics of diffeomorphisms provided with a right invariant metric was proposed in the LDDMM context by [START_REF] Younes | Jacobi fields in groups of diffeomorphisms and applications[END_REF]. This framework enables to transport diffeomorphic deformations of point supported and image data, and it was applied to study the hippocampal shape changes in Alzheimer's disease [START_REF] Qiu | Parallel transport in diffeomorphisms distinguish the time-dependent pattern of hippocampal surface deformation due to healthy aging and dementia of the Alzheimer's type[END_REF][START_REF] Qiu | Time sequence diffeomorphic metric mapping and parallel transport track time-dependent shape changes[END_REF]. Although it represents a rigorous implementation of the parallel transport, it comes to the price of the computationally intense scheme. More importantly, it is limited to the transport along geodesics of the considered right invariant metric, and does not allow to specify different metrics for longitudinal and inter-subject registrations. While from the theoretical point of view parallel transporting along a generic curve can be approximated by the parallel transport on piecewise geodesics, the effectiveness of the above methods was shown only on LDDMM geodesics, and no general computational schemes were provided. The parallel transport in the SVF setting was investigated in [START_REF] Lorenzi | Geodesics, parallel transport & one-parameter subgroups for diffeomorphic image registration[END_REF], in which explicit formula for the parallel transport with respect to the standard Cartan-Schouten connections (left, right and symmetric) in the case of finite dimensional Lie groups were derived. Then it was proposed to seamlessly apply these formulas in the infinite dimensional case of the diffeomorphic registration of images. Although further investigations would be needed to better understand the impact of generalizing to infinite dimensions the concepts defined for the Lie Group theory in finite dimension, practical examples of parallel transport of longitudinal diffeomorphisms in synthetic and real images with respect to the Cartan-Schouten connections showed to be an effective and simple way to transport tangent SVFs. In particular the parallel transport of the left, right and symmetric Cartan connection defined in Equations ( 3), (4), and ( 5) was directly applied to longitudinal deformations and compared against the discrete transport implemented with the Schild's ladder. These experiments highlighted the central role of the numerical implementation on stability and accuracy of the methods. For instance, the left and symmetric Cartan transports were much less stable than the right one because they involve the computation of the Jacobian matrix, computed here with standard finite differences. More robust numerical schemes to compute differential operators on discrete image grids are definitely required to compare them on a fair basis.

Discrete Ladders: Application to Image Sequences

Let I i (i = 1 . . . n) be a time series of images with the baseline I 0 as reference. Consider a template image T 0 , the aim of the procedure is to compute the image T i in order to define the transport of the sequence I 0 , . . . , I i in the reference of T 0 . In the sequel, we focus on the transport of a single image I 1 .

We assume that a well posed Riemannian metric is given on the space of images. This could be L 2 , H k or the metric induced on the space of images by the action of diffeomorphisms of a well chosen right-invariant metric (LDDMM).

Schild's ladder can be naturally translated in the image context (Algorithm 1), by requiring the computation of two diagonal geodesics.

The pole ladder is similar to the Schild's one, with the difference of explicitly using as a diagonal the geodesic C which connects I 0 and T 0 (Algorithm 2). This is an interesting property since, given C, it requires the computation of only one additional geodesic, thus the transport of time series of several images is based on the same baseline-to-reference curve C (Figure 3). computation time and resources if we compute them with registration. Moreover, since we look for regular transformations of the space, the registration is usually constrained to be smooth and the perfect match of correspondent intensities in the registered images is not possible. For instance, the definition of I 1 2 using the forward deformation on I 1 or the backward from T 0 would lead to different results. Since we work in computational anatomy with deformations, it seems more natural to perform the parallel transport directly in the group of diffeomorphisms

Effective Ladders Within the SVF Setting

Given a pair of images I i , i ∈ {0, 1}, the SVF framework parametrizes the diffeomorphism ϕ required to match the reference I 0 to the moving image I 1 by a SVF u. The velocity field u is an element of the Lie Algebra g of the Lie group of diffeomorphisms G, i.e. an element of the tangent space at the identity T id G. The diffeomorphism ϕ belongs to the one parameter subgroup ϕ = exp(tu) generated by the flow of u. We can therefore define the paths in the space of the diffeomorphisms from the one parameter subgroup parametrization l(λ) = exp(λ • u).

Figure ( 4) illustrate how we can take advantage of the stationarity properties of the one-parameter subgroup in order to define the following robust scheme:

1. Let I 1 = exp(u) * I 0 . 2. Compute v = argmin v∈G E (T 0 • exp(-v/2), I 0 • exp(v/2))
, where E is a generic registration energy functional to be minimized.

The half space image I 1 2 can be defined in terms of v/2 as exp(-v/2) * T 0 or exp(v/2) * I 0 . While from the theoretical point of view the two images are identical, the choice of one of them, or even their mean, introduces a bias in the construction. The definition of the half step image can be bypassed by relying on the symmetric construction of the parallelogram.

3. The transformation from I 1 to I 1 2 is ρ = exp(v/2) • exp(-u) and the symmetry leads to exp (Π(u)

) = exp(v/2) • ρ -1 = exp(v/2) • exp(u) • exp(-v/2).
The transport of the deformation ϕ = exp(u) can be therefore obtained through the conjugate action operated by the deformation parametrized by v/2.

Since the direct computation of the conjugation by composition is potentially biased by the spatial discretization, we propose a numerical scheme to more robustly evaluate the transport directly in the Lie Algebra.

BCH Formula for the Conjugate Action

The Baker Campbell Hausdorff (BCH) formula was introduced in the SVF diffeomorphic registration in [START_REF] Bossa | Contributions to 3D diffeomorphic atlas estimation: Application to brain images[END_REF] and provides an explicit way to compose diffeomorphisms parametrized by SVFs by operating in the associated Lie Algebra. More specifically, if v, u are SVFs, then exp(v) • exp(u) = exp(w) with r r r exp(u)

✲ ✑ ✑ ✑ ✑ ✑ ✑ ✑ ✸ ✑ ✑ ✑ ✑ ✑ ✑ ✑ ✰ exp(-Π(u)) ✛ r ❆ ❆ ❆ ❆ ❆ ❆ ❆ ❆ ❆ ❆ ❆ ❆ ❆ ❆ ❑ ❆ ❆ ❆ ❆ ❆ ❑ r ✑ ✑ ✑ ✑ ✑ ✑ ✑ ✑ ✑ ✑ ✑ ✑ ✑ I 0 I 1 T ′ 1 T 0 I 1 2 ρ ρ exp(v/2) exp(-v/2)
Figure 4: Ladder with the one parameter subgroups. The transport exp(Π(u)) is the deformation exp

(v/2) • exp(u) • exp(-v/2) w = BCH(v, u) = v + u + 1 2 [v, u] + 1 12 [v, [v, u]] -1 12 [u, [v, u
]] + . . .. In particular, for small u, the computation can be truncated to any order to obtain a valid approximation for the composition of diffeomorphisms. Applying the truncate BCH to the conjugate action leads to

Π BCH (u) ≃ u + [v/2, u] + 1 2 [v/2, [v/2, u]]. (9)
To establish this formula, let consider the following second order truncation of the BCH formula

BCH((v/2, u) ≃ v/2 + u + 1 2 [v/2, u] + 1 12 [v/2, [v/2, u]] - 1 12 [u, [v/2, u]].
The composition

Π BCH (u) = BCH (v/2, BCH(u, -v/2)) is Π(u) v = v/2 + BCH(u, -v/2) A + 1 2 [v/2, BCH(u, -v/2)] B + 1 12 [v/2, [v/2, BCH(u, -v/2)] C - 1 12 [BCH(u, -v/2), [v/2, BCH(u, -v/2)]] D .
3. Let u = w.

4. Iterate the steps 2 and 3 n times.

The BCH formula allows to perform the transport directly in the Lie algebra and avoids multiple exponentiation and interpolations, thus reducing the bias introduced by the numerical approximations. Moreover, this method preserves the original "ladder" formulation, operated along the inter-subject geodesic exp(tv).

In fact it iterates the construction of the ladder along the path exp(tv) over small steps of size exp( v n ). The stability of the proposed method critically depends from the initial scaling step n, which determines the step-size of the numerical scheme. Ideally the step-size should depend on the curvature, and be therefore small enough in order to minimize the error in case of highly curved space. For this purpose, given the image domain Ω, we define a global scaling factor n in order to guarantee that the given SVF stays sufficiently close to 0, i.e. in order to satisfy the global condition max x∈Ω v(x) /n < δ, with δ = 0.5 * voxel size. This condition ensures reasonably small SVFs, and thus enables the iterative construction of the parallelogram in small neighborhoods.

Pole Ladder for Estimating Longitudinal Changes in Alzheimer's Disease

We provide here an application of the pole ladder for the estimation of a group-wise model of the longitudinal changes in a group of patients affected by Alzheimer's disease (AD). AD is a neurodegenerative pathology of the brain, characterized by the co-occurrence of different phenomena, starting from the deposition of amyloid plaques and neurofibrillary tangles, to the development of functional loss and finally to cell deaths [START_REF] Jack | Hypothetical model of dynamic biomarkers of the Alzheimer's pathological cascade[END_REF]. In particular brain atrophy detectable from magnetic resonance imaging (MRI) is currently considered as a potential outcome measure for the monitoring of the disease progression. Structural atrophy was shown to strongly correlate with cognitive performance and neuropsychological scores, and characterizes the progression from pre-clinical to pathological stages [START_REF] Jack | Hypothetical model of dynamic biomarkers of the Alzheimer's pathological cascade[END_REF]. For this reason, the development of reliable atlases of the pathological longitudinal evolution of the brain is of paramount importance for improving the understanding of the pathology. A preliminary approach to the group-wise analysis of longitudinal morphological changes in AD consists in performing the longitudinal analysis after the subject-to-template normalization [START_REF] Chetelat | Using voxel-based morphometry to map the structural changes associated with rapid conversion to MCI[END_REF][START_REF] Thompson | Dynamics of gray matter loss in Alzheimer's disease[END_REF]. A key issue here is the different nature of the changes occurring at the intra-subject level, which reflects the biological phenomena of interest, and the changes across different subjects, which are usually large and not related to any biological process. In fact, the inter-subject variability is a scale of magnitude higher than the more subtle longitudinal subject-specific variations. To provide a more sensitive quantification of the longitudinal dynamics, the intra-subject changes should be modeled independently from the subject-to-template normalization, and only transported in the common reference for statistical analysis afterward. Thus, novel techniques such as the parallel transport of longitudinal deformations might lead to better accuracy and precision for the modeling and quantification of longitudinal pathological brain changes.

Data analysis and Results

Images corresponding to the baseline I 0 and the one-year follow-up I 1 scans were selected for 135 subjects affected by Alzheimer's disease. For each subject i, the pairs of scans were rigidly aligned. The baseline was linearly registered to a reference template and the parameters of the transformation were applied to I i 1 . Finally, for each subject, the longitudinal changes were measured by non-linear registration using the LCC-Demons algorithm [START_REF] Lorenzi | LCC-Demons: a robust and accurate symmetric diffeomorphic registration algorithm[END_REF].

The resulting deformation fields ϕ i = exp(v i ) were transported with the pole ladder (BCH scheme) in the template reference along the subject-to-template deformation. The group-wise longitudinal progression was modeled as the mean of the transported SVFs v i . The areas of significant longitudinal changes were investigated by one-sample t-test on the group of log-Jacobian scalar maps corresponding to the transported deformations, in order to detect the areas of measured expansion/contraction significantly different from zero.

For the sake of comparison, the one sample t-statistic was tested on the subject specific longitudinal log-Jacobian scalar maps warped into the template space along the subject-to-template deformation. This is the classical transport used in tensor's based morphometry studies [START_REF] Ashburner | Symmetric diffeomorphic modeling of longitudinal structural MRI[END_REF].

Figure 5 shows a detail from the mean SVF from the transported one-year longitudinal trajectories. The field flows outward from the ventricles indicates a pronounced enlargement. Moreover, we notice an expansion in the temporal horns of the ventricles as well as a consistent contracting flow in the temporal areas. The same effect can be statistically quantified by evaluating the areas where the log-Jacobian maps are statistically different from zero. The areas of significant expansion are located around the ventricles and spread in the CSF areas, while a significant contraction is appreciable in the temporal lobes, hippocampi, parahippocampal gyrus and in the posterior cingulate. The statistical result is in agreement with the one provided by the simple scalar interpolation of the longitudinal subject specific log-Jacobian maps. In fact we do not experience any substantial loss of localization power by transporting SVFs instead of scalar log-Jacobian maps. However by parallel transporting we preserve also the multidimensional information of the SVFs that, as experienced in [START_REF] Lorenzi | Mapping the effects of Aβ 1-42 levels on the longitudinal changes in healthy aging: hierarchical modeling based on stationary velocity fields[END_REF], potentially leads to more powerful voxel-by-voxel comparisons than the ones obtained with univariate tests on scalars.

Finally, figure 6 shows that the apparent volume changes associated to the average trajectory computed with the pole ladder describe biologically plausible dynamics of longitudinal atrophy. We notice that the estimated one-year longitudinal trajectory is associated to local volume changes ranging from +12% for the expansion of the ventricles, to 5% for the volume loss of the hippocampi.

We recall that the apparent expansion of the CSF areas is detectable thanks to the diffeomorphic registration constraint. In fact, since the deformation is spatially smooth, the apparent volume loss (e.g. brain atrophy) detectable in the image is modeled as voxels shrinkage, which is associated to the enlargement of the surrounding areas (e.g. ventricles).

Conclusions

This chapter illustrates the principles of parallel transporting in transformation groups, with particular focus on discrete transport methods. The use of discrete transport is motivated from both theoretical and practical point of view. In fact, discrete methods such as the pole ladder are based only on the computation of geodesics, and thus they do not require the explicit knowledge of the connection Figure 6: Apparent relative volume changes encoded by the average longitudinal trajectory computed with the pole ladder (Figure 5). The trajectory describes a pattern of apparent volume gain in the CSF areas, and of apparent volume loss in temporal areas and around the ventricles.

of the space. This is a rather interesting characteristic that enables to employ the ladder without requiring the design of any additional tool outside geodesics.

From the practical point of view, discrete methods can alleviate the numerical problems arising from the discretization of continuous functional on finite grids, and thus provide feasible and numerically stable alternatives to continuous transport approaches. The application shown in Section 6.4 is a promising example of the potential of such approaches when applied to challenging problems such as the estimation of longitudinal atlases in diffeomorphic registration.

As shown in Section 4 the construction of the ladder holds in sufficiently small neighborhoods. From the practical point of view this is related to the choice of an appropriate step size for the iterative scheme proposed in Section 6, and future studies are required in order to investigate the impact of the step size from the numerical point of view.

Finally, future studies aimed to directly compare discrete versus continuous approaches might shed more light on the theoretical and numerical properties of different methods of transport.

Figure 1 :

 1 Figure 1: (1) The transport of the vector A along the curve C is performed by the Schild's ladder by (2) the construction of geodesic parallelograms in a sufficiently small neighborhood. (3) The construction is iterated for a sufficient number of neighborhoods.

Figure 2 :

 2 Figure 2: The pole ladder parallel transports the vector A along the geodesic C. Contrarily to the Schild's ladder, it only requires to compute one diagonal geodesic.

Figure 3 : 2 . 3 .Algorithm 2 2 . 3 . 1 2 2 . 4 .

 323223124 Figure 3: Geometrical schemes in the Schild's ladder and in the pole ladder. By using the curve C as diagonal, the pole ladder requires the computation of half times of the geodesics (blue) required by the Schild's ladder (red).

Figure 5 :

 5 Figure 5: One year structural changes for 135 Alzheimer's patients. A) Mean of the longitudinal SVFs transported in the template space with the pole ladder. We notice the lateral expansion of the ventricles and the contraction in the temporal areas. B) T-statistic for the correspondent log-Jacobian values significantly different from 0 (p < 0.001 FDR corrected). C) T-statistic for longitudinal log-Jacobian scalar maps resampled from the subject to the template space. Blue color: significant expansion, Red color: significant contraction. The figure is reproduced from[START_REF] Lorenzi | Efficient parallel transport of deformations in time series of images: from Schild's to pole ladder[END_REF]

The second order truncation of the four terms is:

From the additive and anticommutative properties of the Lie bracket, adding the four terms leads to [START_REF] Bossa | Contributions to 3D diffeomorphic atlas estimation: Application to brain images[END_REF].

Iterative Computation of the ladder

Once defined the formula for the computation of the ladder, we need a consistent scheme for the iterative construction along trajectories. We recall that the transport by geodesic parallelograms holds only if both sides of the parallelogram are sufficiently small, which in our case means that both longitudinal and intersubject vectors must be small. This is not the case in practice, since the intersubject deformation is usually very large. By definition, the ladder requires to scale down vectors to a sufficiently small neighborhood, in order to correctly approximate the transport by parallelograms.

From the theoretical point of view, the degree of approximation of the ladder is approximately proportional to the curvature of the space of deformations. This can be seen by the higher order terms that we dropped off in the proof of Section 4.2, which are all derivatives of the Christoffel symbols. While on a linear space the ladder is the exact parallel transport, when working on curved spaces the error resulting from the non-infinitesimal geodesic parallelogram is proportional to the distance between the points.

From the numerical point of view, we notice that Formula (9) requires the computation of the Lie brackets of the velocity fields. Lie brackets involve the differentiation of the vector which is usually computed on images by finite differences, and which are know to be very sensitive to noise and to be unstable in case of large deformations.

For all these reasons we propose the following iterative scheme based on the properties of SVFs. To provide a sufficiently small vector for the computation of the conjugate we observe that exp

The conjugation can then be recursively computed in the following way:

1. Scaling step. Find n such that v/n is small.

Ladder

Step. Compute