Nicolas Guigui

Xavier Pennec

Keywords: Spline regression . . . . . . . . . . . . . . . . . 3.5.2 Hotelling tests on velocities . . . . . . . . . . . . . . . . . . Parallel transport, longitudinal studies, mean trajectory, cardiac motion analysis, Schild's ladder, pole ladder, Riemannian manifolds

Transporting the statistical knowledge regressed in the neighbourhood of a point to a different but related place (transfer learning) is important for many applications. In medical imaging, cardiac motion modelling and structural brain changes are two such examples: for a group-wise statistical analysis, subjectspecific longitudinal deformations need to be transported in a common template anatomy.

In geometric statistics, the natural (parallel) transport method is defined by the integration of a Riemannian connection which specifies how tangent vectors

Parallel transport, a central tool in geometric statistics for computational anatomy. Application to cardiac motion modelling are compared at neighbouring points. In this process, the numerical accuracy of the transport method is critical. Discrete methods based on iterated geodesic parallelograms inspired by Schild's ladder were shown to be very efficient and apparently stable in practice. In this chapter, we show that ladder methods are actually second order schemes, even with numerically approximated geodesics.

We also propose a new original algorithm to implement these methods in the context of the Large Deformation Diffeomorphic Metric Mapping (LDDMM) framework that endows the space of diffeomorphisms with a right-invariant RKHS metric.

When applied to the motion modelling of the cardiac right ventricle under pressure or volume overload, the method however exhibits unexpected effects in the presence of very large volume differences between subjects. We first investigate an intuitive rescaling of the modulus after parallel transport to preserve the ejection fraction. The surprisingly simple scaling/volume relationship that we obtain suggests to decouples the volume change from the deformation directly within the LDMMM metric. The parallel transport of cardiac trajectories with this new metric now reveals statistical insights into the dynamics of each disease. This example shows that parallel transport could become a tool of choice for data-driven metric optimization.

Introduction

At the interface of geometry, statistics, image analysis and medicine, Computational Anatomy aims at analysing and modelling the variability of the biological shape of tissues and organs and their dynamics at the population level. The goal is to estimate representative anatomies across diseases, populations, species or ages, to model the organ development across time (growth or ageing), to discover morphological differences between normal and pathological groups and to estimate and correlate the variability with other functional, genetic or structural information. In the context of cardiology, computational anatomy models the cardiac cycle as a sequence of deformations of an anatomy and allows to characterize quantitatively these deformations to compare the impact of diseases on the cardiac function.

The analysis of the organs' shape and deformations often relies on the identification of features to describe locally the anatomy such as landmarks, curves, surfaces, intensity patches, full images, etc. Modelling their statistical distribution in the population requires to first identify point-to-point anatomical correspondences between these geometric features across subjects. This may be feasible for landmark points, but not for curves or surfaces. Thus, one generally considers relabelled point-sets or reparametrized curves/surfaces/images as equivalent objects. With this geometric formulation, shape spaces are the quotient of the original space of features by their reparametrization group. The difficulty is that taking the quotient generally endows the shape space with a non-linear manifold structure, even when we start from features living in a Euclidean space.

For instance, equivalence classes of k-tuples of points under rigid or similarity transformations result in non-linear Kendall's shape spaces (see e.g. Dryden and Mardia 2016 for a recent account on that subject). The quotient of curves, surfaces and higher dimensional objects by their reparametrizations (diffeomorphisms of their domains) produces in general even more complex infinite dimensional shape spaces [START_REF] Bauer | Overview of the Geometries of Shape Spaces and Diffeomorphism Groups[END_REF]. Thus, shapes belong in general to non-linear manifolds, while statistics were essentially developed for linear and Euclidean spaces. This has motivated the use and development of a consistent statistical framework on Riemannian manifolds and Lie groups during the past 25 years, a field called Geometric Statistics in computational anatomy (Pennec, Sommer, and Fletcher 2020).

Deformation-based morphometry (DBM) is also a popular method to study statistical differences of brain anatomies. It consists in analyzing local deformation features of non-linear image registration to a reference [START_REF] Ashburner | Identifying Global Anatomical Differences: Deformationbased Morphometry[END_REF]. For instance, the Jacobian of the transformation encoding local volumes changes can be used to detect the areas that are statistically growing of shrinking in a population. However, DBM analyses the deformation independently at each point of the image. In order to capture spatial correlations, it is more interesting to model the transformations of all the points together, i.e. to lift statistics from the objects (image voxels, curves, surfaces) to the transformations of their embedding space. This powerful idea also allows us to capture jointly the variability of several structures to model their interactions. In the brain, for instance, one can consider together the shape of the cortex, the ventricles, deep grey nuclei and other internal brain structures. We can also include structures that are not surface, such as the sulcal lines encoding the complex folding patterns of the grey matter at the surface of the cortex.

Diffeomorphometry

Following D'Arcy [START_REF] Thompson | On Growth and Form[END_REF] is often assumed that there exists a template shape or image (called an atlas in medical image analysis) that represents the standard (or mean) anatomy, and that the inter-subject variability is encoded by deformations of that template towards the shapes of each subject or their evolution in time. A very desirable aspect of these transformations is to smoothly preserve the spatial organization of the anatomical tissues by avoiding intersections, folds or tears. Simply encoding deformations with the vector space of displacement fields is not sufficient to preserve the topology: one needs to require diffeomorphic transformations (differentiable one-to-one transformations with differentiable inverse).

Lie groups of diffeomorphisms are examples of space that are both infinite dimensional manifolds and Lie groups, and the statistic analyse of shapes through their diffeomorphic deformations has been coined as diffeomorphometry. This approach was pushed forward by (Grenander and M. Miller 1998) and turned into a mathematically grounded framework by endowing the space of diffeomorphisms with a sufficiently regular right invariant metric (Younes 2019;M. Miller, Trouvé, and Younes 2015), leading to the so called Large Deformation Diffeomorphic Metric Mapping (LDDMM) framework, detailed in Section 2.2.

Since the optimization of time-varying velocity fields was computationally intensive, an alternative parameterization of diffeomorphisms with the flow of stationnary velocity fields (SVF) was introduced by [START_REF] Arsigny | A Log-Euclidean Framework for Statistics on Diffeomorphisms[END_REF]. The flow of SVFs generates one-parameter subgroups, which are simply matrix exponentials in matrix Lie groups. Although a number of theoretical difficulties remain when moving to infinite dimensions, very efficient algorithms can be adapted from the matrix case, like the scaling and squaring procedure [START_REF] Higham | The Scaling and Squaring Method for the Matrix Exponential Revisited[END_REF] to compute the group exponential and its Jacobian, or the Baker-Campbell-Hausdorff (BCH) formula to approximate the composition of two deformations directly in the log-domain. This allows the straightforward extension of the classical "demons" image registration algorithm to encode diffeomorphic transformations parametrized by SVFs. Because the inverse of the flow of an SVF is simply parameterized by the opposite of this SVF, a special feature of the log-demons registration framework is to enforce almost seamlessly the very desirable inverse consistency of the registration. The SVF framework was successfully used in a number of application to brain studies [START_REF] Lorenzi | Schilds Ladder for the Parallel Transport of Deformations in Time Series of Images[END_REF][START_REF] Hadj-Hamou | Longitudinal Analysis of Image Time Series with Diffeomorphic Deformations: A Computational Framework Based on Stationary Velocity Fields[END_REF]) as we will se below.

The differential geometric foundations of this very efficient SVF framework were uncovered in [START_REF] Lorenzi | Geodesics, Parallel Transport & One-parameter Subgroups for Diffeomorphic Image Registration[END_REF]): one-parameters subgroups are actually geodesics (in the sense of auto-parallel curves) of the Cartan-Schouten connection, a canonical bi-invariant symmetric affine connection that exists on every Lie group. Whenever a bi-invariant metric exists on the group, this affine connection coincides with the Levi-Civita connection of that bi-invariant metric, in which case the Fréchet mean is left, right and inverse equivariant. However, only direct products of compact and Abelian groups admit a bi-invariant metric, which is very limiting in practice. Thus, Lie groups generally do not admit any bi-invariant metric, in which case the Fréchet mean based on a left-invariant distance is not consistent with right translations nor with inversions. This is already the case for Euclidean motions SE(3) since it is a semi-direct product. In contrast, exponential barycenters of the canonical Cartan-Schouten connection define a non-metric notion of mean on any Lie groups which is automatically equivariant by left, right translation and inversion [START_REF] Pennec | Exponential Barycenters of the Canonical Cartan Connection and Invariant Means on Lie Groups[END_REF]. One can also define covariance matrices and Mahalanobis distances that are consistent with the group operations. Thus, from the points of view of abstract transformations, independently of their action on objects, the Cartan-Schouten connection offers a unique way to define a consistent statistical framework on Lie groups. This drove the interest for statistics on affine connection spaces, a superclass of Riemannian manifolds where non-metric geodesics are defined locally as auto-parallel curves [START_REF] Pennec | Beyond Riemannian Geometry: The Affine Connection Setting for Transformation Groups[END_REF].

Longitudinal models

Once we have decided for the LDDMM or the SVF parameterization of diffeomorphisms, a first and simple approach to tackle longitudinal data is to regress the template and its deformations over time by minimizing the inter-subject registration square-distance to each observation. The template trajectory may be a geodesic in a group of diffeomorphism for very few data points, such as for studying the remodelling of the heart over time in paediatric images [START_REF] Mansi | A Statistical Model of Right Ventricle in Tetralogy of Fallot for Prediction of Remodelling and Therapy Planning[END_REF][START_REF] Mcleod | Statistical Shape Analysis of Surfaces in Medical Images Applied to the Tetralogy of Fallot Heart[END_REF]). With more samples over time, one may fit a deformation trajectory with kernel regression [START_REF] Gerber | Manifold Modeling for Brain Population Analysis[END_REF] or assume a spline structure (Trouvé and Vialard 2012b;[START_REF] Hinkle | Intrinsic Polynomials for Regression on Riemannian Manifolds[END_REF][START_REF] Singh | Splines for diffeomorphisms[END_REF]. This type of method is good for a cross-sectional design, where we have one data point per subject, but is likely to hide the small longitudinal changes of each subject into the large inter-subject variability when we have more than one time-point per subject.

Longitudinal models are preferable to study time-series of shapes in both dynam-ical systems such as the heart and disease progression modelling such as the aging brain. Due to the non-commutativity of the longitudinal and inter-subject deformations, several paradigms were investigated to compare transformations trajectories at different points [START_REF] Durrleman | Sparse Adaptive Parameterization of Variability in Image Ensembles[END_REF]. One choice is to model the reference trajectory in the template space and to make it patient-specific by deforming it using a template to subject transformation valid for all times [START_REF] Rao | Spatial Transformation of Motion and Deformation Fields Using Nonrigid Registration[END_REF][START_REF] Durrleman | Sparse Adaptive Parameterization of Variability in Image Ensembles[END_REF]. This approach can be geometrized by parallel transporting the initial velocity of the template to subject deformation along the reference longitudinal deformation, as proposed by [START_REF] Schiratti | Learning Spatiotemporal Trajectories from Manifold-Valued Longitudinal Data[END_REF]. This was called exp-parallelization for a geodesic reference trajectory, but the method works for any smooth curve. A recent refinement of this model has been proposed to cluster trajectories into distinct classes with different trajectories [START_REF] Debavelaere | Clustering of Longitudinal Shape Data Sets Using Mixture of Separate or Branching Trajectories[END_REF]). However, this does not allow to easily model the variability on the longitudinal trajectories which is observed for instance in the heart. Moreover, exp-parallel curve may be very different from geodesics when the atlas-to-subject transformation is large: on a sphere, for instance, exp-parallel curves are small circles parallel to the reference geodesic great circle, and reduce to a point at the poles.

In this chapter, we prefer to use the other main approach, where we first regress a patient-specific trajectory, and then normalize it with respect to the patient-to-atlas deformation at a reference time to perform the statistical analysis of all trajectories at the population level [START_REF] Lorenzi | Geodesics, Parallel Transport & One-parameter Subgroups for Diffeomorphic Image Registration[END_REF]. It is sometimes argued that there is no canonical reference time for comparing the evolution of brains. This is not the case for the heart, for which the end-diastole (ED) is a natural time reference of the cardiac cycle. The subject-specific trajectory is usually a geodesic for the analysis of structural brain changes with ageing in Alzheimer's disease [START_REF] Hadj-Hamou | Longitudinal Analysis of Image Time Series with Diffeomorphic Deformations: A Computational Framework Based on Stationary Velocity Fields[END_REF]. For the heart, one can use any motion trajectory parametrized by its logarithm (the initial tangent vector or momentum that geodesically shoots to the deformation at a given point of the trajectory). In this setting, one has a curve of tangent vectors encoding the subject-specific deformation over time that needs to be transported to a reference template.

Transporting scalar values (e.g. the intensity of en image) along a diffeomorphic transformation is simply done by resampling that image. However, when one wants to transport differential (vector) quantities, it is not sufficient to resample each of the vector component individually: there needs to be a normalization that also reorient the vector. For instance, a constant initial vector field encodes a translation over time. However, when the inter-subject transformation is a rotation by 90 degrees, that vector needs to be rotated along. In this chapter, we focus on geometric methods based on parallel transport.

For cardiac motion, an average cardiac deformations across patients is a called a statistical motion atlas [START_REF] Young | Computational Cardiac Atlases: From Patient to Population and Back[END_REF][START_REF] Peyrat | Registration of 4D Cardiac CT Sequences Under Trajectory Constraints With Multichannel Diffeomorphic Demons[END_REF][START_REF] Duchateau | A Spatiotemporal Statistical Atlas of Motion for the Quantification of Abnormal Myocardial Tissue Velocities[END_REF]. In this application, the need to normalize deformations is salient in the example of Figure 1, where the deformation between the patients' end-diastolic (ED) and end-systolic (ES) frames are applied without any normalization to another reference shape, the atlas. These result in meshes of poor quality, with irregularities, especially near the valves and on the septum. These abnormalities are salient when considering ventricles with pressure and volume overload as the amplitude of both the intra-patient and patient-to-atlas deformations may be large.

Parallel transport for inter-subject normalisation

Given a manifold equipped with an affine connection, for example a Riemannian manifold with its Levi-Civita connection, one can define the parallel transport of a tangent vector v along a curve γ.

Definition 1.1 (Parallel vector field). Let M be a smooth manifold and ∇ a connection on M. For any curve γ :

[a, b] → M in M , a vector field X along γ is parallel if ∇ γ(t) X(t) = 0. ( 1
)
From the properties of ODEs, one can prove that given γ and v a tangent vector at x = γ(0), there exists a unique parallel vector field X along γ such that X(0) = v. Intuitively, this ODE constrains the parallel vector field to keep its orientation w.r.t the velocity γ of the curve while moving along it. For any t in the domain of definition of γ, X(t) is called the parallel transport of v at y = γ(t) along γ, and written Π y x v. For normalizing the momentum of a longitudinal deformation along an intersubject geodesic diffeomorphism, parallel transport was first proposed in the LD-

m 1 x 1 x 0 m 2 x 2 v v 2 v 1 (a) Schild's ladder m 1 x 1 x 0 m 2 x 2 v v 2 -v 1 (b) Pole ladder
Figure 2: Representation of the two ladder schemes, Schild's (left) and pole ladder (right). The methods consist in iterating the construction of geodesic parallelograms to approximate the parallel transport of a vector v along a curve which is geodesic by arc.

DMM framework using Jacobi fields [START_REF] Younes | Jacobi Fields in Groups of Diffeomorphisms and Applications[END_REF][START_REF] Qiu | Parallel Transport in Diffeomorphisms Distinguishes the Time-Dependent Pattern of Hippocampal Surface Deformation Due to Healthy Aging and the 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]. The method was used in [START_REF] Cury | Spatio-Temporal Shape Analysis of Cross-Sectional Data for Detection of Early Changes in Neurodegenerative Disease[END_REF] for the analysis of the thalamus in fronto-temporal dementia and it was improved by the Fanning Scheme [START_REF] Louis | A Fanning Scheme for the Parallel Transport Along Geodesics on Riemannian Manifolds[END_REF], which is applied to shape analysis of brain structures [START_REF] Louis | Parallel Transport in Shape Analysis: A Scalable Numerical Scheme[END_REF]. In latter works, the authors claim that the Jacobi Field approach is more precise and less expensive computationally than ladder methods. We have proved in Guigui and Pennec 2021 that this is not the case. Moreover, the proposed method implicitly assumes that the same metric is used for both longitudinal and inter-subject deformation. This is an important concern as the cross-sectional inter-subject and longitudinal intra-subject deformations have a fundamentally different nature.

In order to implement a parallel transport algorithm that remains consistent with the numerical scheme used to compute the geodesics, [START_REF] Lorenzi | Schilds 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] proposed to adapt Schild's ladder to image registration with deformations parametrized by SVF. This method relies on iterating the construction of geodesic parallelograms to approximate the parallel transport of a vector v along a geodesic (see Figure 2a). Interestingly, the Schild's ladder implementation appeared to be more stable in practice than the closed-form expression of the symmetric Cartan-Schouten parallel transport on geodesics. This was attributed to the inconsistency of the numerical schemes used for the computation of the transformation Jacobian in the implementation of this exact parallel transport formula.

Shortly after, it was realized that parallel transport along geodesics could exploit an additional symmetry by using the geodesic along which we want to transport as one of the diagonal of the geodesic parallelogram: this gave rise to the pole ladder scheme (Lorenzi and Pennec 2014). In this case the geodesic along which we are transporting is used as diagonal of the parallelograms (see Figure 2b). In this case the geodesic along which we are transporting is used as diagonal of the parallelograms (see Figure 2b). This greatly reduces the number of geodesics to compute (-v 1 does not even need to be computed on Figure 2b). Pole ladder combined with an efficient numerical scheme based on the Baker-Campbell-Hausdorff formula was found to be more stable on simulated and real experiments than the other parallel transport schemes tested on Lie groups. This result and the higher symmetry led to conjecture that pole ladder could be a higher order scheme than Schild's ladder.

The SVF framework was successfully applied in [START_REF] Lorenzi | Disentangling Normal Aging from Alzheimer's Disease in Structural Magnetic Resonance Images[END_REF] to distinguish pathological evolution from normal aging. An efficient computational pipeline is provided by [START_REF] Hadj-Hamou | Longitudinal Analysis of Image Time Series with Diffeomorphic Deformations: A Computational Framework Based on Stationary Velocity Fields[END_REF]) and used in [START_REF] Sivera | Voxel-Based Assessments of Treatment Effects on Longitudinal Brain Changes in the Multidomain Alzheimer Preventive Trial Cohort[END_REF] to model Alzheimers's disease and to analyze the effect of a potential treatment in a clinical trial (see illustration in Fig. 3).

Chapter organization

In this chapter, we address the longitudinal analysis of cardiac motion across patients and diseases. A key problem of this application is the large cyclic nature of cardiac deformations, along with large to very large inter-subject transformations. In order to caption the large cardiac motion, we chose to use the LDDMM framework because a right invariant metric is compatible with a Lagrangian and a Hamiltonian formulation of the mechanical equations, while the SVF framework is rather related to a Eulerian invariance.

The numerical complexity of the implementation of Jacobi fields and the algorithmic simplicity of ladder methods led us to chose the later for inter-subject normalization. However, despite the practical success of ladder methods for longitudinal brain analyses, their numerical accuracy remained essentially unknown beyond the first order. This led Guigui and Pennec 2021 to conduct a careful numerical analysis of Schild's and pole ladders that we summarize in section 2.1. To the best of our knowledge, ladder methods have not been used in the LDDMM framework, for which a BCH-type approximation was missing. We first give an overview of the LD-DMM framework in section 2.2 before proposing a new second-order implementation of pole ladder on LDDMM deformations in section 2.3.

We turn in section 3 to the application of this methodology to the group-wise analysis of cardiac motion across diseases. For the cardiac motions that we analyse in section 3.2, relative volume differences between subjects appear to have an unexpected effect on the parallel transported trajectories. Since this undesirable effect cannot be attributed to the numerical accuracy of the transport scheme, it indicates that we have to revise the choice of the Riemannian metric used for the inter-subject comparison. We first investigate in section 3.3 an intuitive rescaling of the modulus of the momentum after the parallel transport to preserve on average the ejection fraction along the motion trajectory. Regression results over a population of subjects and patients show a surprisingly simple relationship between the scaling and the inter-subject volume ratio. However, modifying the transport equations in an ad-hoc fashion is not satisfactory from a theoretical point. Thus, we investigate in section 3.4 a more satisfactory strategy that decouples the volume change from the deformation directly within the LDMMM metric.

Parallel transport with ladder methods

Numerical methods have been proposed in geometric data processing for the parallel transport of vectors along curves in manifolds. The oldest algorithm is probably Schild's ladder, a general method for the parallel transport along arbitrary curves, introduced in the theory of gravitation in [START_REF] Misner | Gravitation[END_REF] in the spirit of Schild's geometric constructions. Ehlers, Pirani, and Schild 1972 is often cited as reference for Schild's ladder but no mention of the scheme is made in this work. The method extends the infinitesimal transport through the construction of geodesic parallelograms (Figure 4). It is algorithmically interesting since it only requires the computation of geodesics (initial and boundary value problems) without requiring the knowledge of the second order structure of the space (connection or curvature tensors). Kheyfets, W. A. Miller, and Newton 2000 proved that the scheme realizes a first order approximation of the parallel transport for a symmetric connection. This makes sense since the skew-symmetric part of the connection, the torsion, does not impact the geodesic equation. Schild's ladder is nowadays increasingly used in non-linear data processing and analysis to implement parallel transport in Riemannian manifolds. One can cite for instance Lorenzi, Ayache, and Pennec 2011 and the follow-up publications cited in introduction for the parallel transport of deformations in computational anatomy, or Hauberg, Lauze, and Pedersen 2013 for parallel transporting the covariance matrix in Kalman filtering.

Numerical accuracy of Schild's and pole ladders

Despite the use of ladder methods in applications, no results on their convergence were published before the work of Guigui and Pennec 2021 that we summarize here. We give Taylor approximations of the elementary constructions of Schild's ladder and of the pole ladder with respect to the Riemann curvature of the underlying space. This allows to prove that these methods can be iterated to converge with quadratic speed, even when geodesics are approximated by numerical schemes.

Elementary construction of Schild's ladder Recall that Exp denote the Riemannian exponential map, that maps a point and a tangent vector to the point reached at time t = 1 by the unique geodesic with those initial conditions. Log is its inverse defined locally, such that Exp x (Log x (y)) = y. See Pennec, Sommer, and Fletcher 2020, chapter 1 for more details. The construction to parallel transport v ∈ T x M along the geodesic γ with γ(0) = x and γ(0) = w ∈ T x M (such that (v, w) ∈ U x ) is given by the following steps (Figure 4):

1. Compute the geodesics from x with initial velocities v and w until time s = t = 1 to obtain x v and x w . These are the sides of the parallelogram.

2. Compute the geodesic between x v and x w and the midpoint m of this geodesic:

m = Exp xv 1 2 Log xv (x w ) .
This is the first diagonal of the parallelogram.

3. Compute the geodesic between x and m, let a ∈ T x M be its initial velocity.

Extend it beyond m for the same length as between x and m to obtain z:

a = Log x (m); z = Exp x (2a) = x 2a .
This is the second diagonal of the parallelogram.

4. Compute the geodesic between x w and z. Its initial velocity u w is an approximation of the parallel transport of v along the geodesic from x to x w , i.e. u w = Log xw (x 2a ).

Assuming that there exists a convex neighborhood that contains the entire parallelogram, all the above operations are well defined. In the literature, this construction is then iterated along γ without further precision.

Taylor expansion

We can now reformulate this elementary construction in terms of successive applications of the double exponential and the neighboring logarithm maps, given in Pennec 2019. The computations are detailed in the appendix of Guigui and Pennec 2021, and we only report here the result at third order, meaning that all the terms of the form

∇ • R(•, •)• are summarized in the term O(4).
We first obtain a generalized midpoint rule:

2a = w + v + 1 6 R(v, w)(w -v) + O(4) (2)
We notice that this expression is symmetric in v and w, as expected. Furthermore, the deviation from the Euclidean mean of v, w (the parallelogram law) is explicitly due to the curvature. Accounting for this correction term is a key ingredient to reach a quadratic speed of convergence. Now, by propagating this midpoint rule to compute the error e x made by this construction to parallel transport v, e x = Π xw x vu w , we obtain a third order approximation of the Schild's ladder construction:

Theorem 2.1. Let (M, g) be a finite dimensional Riemannian manifold. Let x ∈ M and v, w ∈ T x M sufficiently small. Then the output u of one step of Schild's ladder parallel transported back to x is given by

u = v + 1 2 R(w, v)v + O(4) (3)
This theorem shows that Schild's ladder is in fact a second-order approximation of parallel transport. Furthermore, this shows that splitting the main geodesic into segments of length 1 n and simply iterating this construction n times will in fact sum n error terms of the form R( w i n , v i )v i , hence by linearity of R, the error won't necessarily vanish as n → ∞. To ensure convergence, it is necessary to also scale v in each parallelogram, as detailed in the next paragraph.

Numerical scheme and convergence With the previous notations, let us define schild(x, w, v) = u w ∈ T xw M . We now divide the geodesic γ into segments of length 1 n for some n ∈ N * large enough so that the previous Taylor expansions can be applied to w n and v n . As mentioned before, v needs to be scaled as w. In fact let α ≥ 1 and consider the sequence defined by (see Figure 4) where

v 0 = v v i+1 = n α • schild x i , w i n , v i n α , (4)
x i = γ( i n ) = Exp x i n w , w i = n Log x i (x i+1 ) = Π x i
x w. We now establish the following result, which ensures convergence of Schild's ladder to the parallel transport of v along γ at order at most two.

Theorem 2.2. Let (x i , v i , w i ) (i≤n) be the sequence defined as above. Then ∃τ >

0, ∃β > 0, ∃N ∈ N, ∀n > N , ∥v n -Π xn x v∥ ≤ τ n α + β n 2 .
Moreover, τ is bounded by a bound on the sectional curvature, and β by a bound on the covariant derivative of the curvature tensor.

The same result can be obtained for pole ladder, and in the case where geodesics are obtained by numerical integration of the geodesic equation.

Pole ladder In this case, the main geodesic is used as diagonal of the constructed parallelograms (see Figure 5). This allows to reduce the number of geodesics that need to be computed. Similarly, a Taylor approximation of the construction gives

u = v + 1 12 (∇ w R)(w, v)(5v -w) + (∇ v R)(w, v)(2v -w) + O(5). (5)
And this can be propagated to show the following bound

∥v n -Π mn m v∥ ≤ β n 2 , ( 6
)
that ensures convergence to the parallel transport with a quadratic speed.

Infinitesimal schemes When geodesics are not available in closed form, we replace the Exp map by a fourth-order numerical scheme (e.g. Runge-Kutta (RK)), and the Log is obtained by gradient descent over the initial condition of Exp. It turns out that only one step of the numerical scheme is sufficient to ensure convergence, and keeps the computational complexity reasonable. As only one step of the integration schemes is performed, we are no longer computing geodesic parallelograms, but infinitesimal ones, and thus refer to this variant as infinitesimal scheme.

As in the previous case, we can show

∥Π xn x v -v n ∥ = O 1 n 2 .
This allows to use this scheme in the SVF and LDDMM frameworks. For the later, we first give an overview of the framework in the next section before presenting our parallel transport algorithm.

A short overview of the LDDMM framework

The LDDMM framework encompasses both algorithms for shape matching (also referred to as registration) and a Riemannian geometric structure on the space of deformations, that projects to the space of shapes. The former allows computing shape descriptors, and to parameterize diffeomorphisms. The latter, provides a distance to compare deformations, and an associated notion of parallelism to transport them. We summarize here the formulation of Durrleman, Allassonnière, and Joshi 2013; [START_REF] Durrleman | Morphometry of Anatomical Shape Complexes with Dense Deformations and Sparse Parameters[END_REF] and their implementation in [START_REF] Bône | Deformetrica 4: An Open-Source Software for Statistical Shape Analysis[END_REF], focusing on the case of landmarks. For a thorough treatment of the topic, we refer to Younes 2019.

In this chapter, we restrict to shapes defined as d -1-dimensional submanifolds of R d , d = 1, 2, with the same given topology, and approximated by triangulated meshes, defined by a set of points in R d , called landmarks, and connectivity between those. The set of shapes is denoted S. The model of computational anatomy stipulates that such a shape is fully described by a diffeomorphism φ ∈ Diff(R d ), where the interaction between shapes and diffeomorphisms is encoded by the action of Diff(R d ) on S. This turns the shape space into a homogeneous space with group Diff(R d ). The isotropy subgroup of a given shape is the set of re-parameterization of this shape, i.e. diffeomorphims that leave the shape unchanged, but move particles along the shape. These diffeomorphisms are infinite dimensional nuisance parameters in a statistical analysis, and one feature of LDDMM is to provide methods that are invariant to such nuisance (M. [START_REF] Miller | Hamiltonian Systems and Optimal Control in Computational Anatomy: 100 Years Since D'Arcy Thompson[END_REF].

A natural and efficient computational construction of diffeomorphisms is obtained by flows associated to ordinary differential equations (ODEs)

∂ϕ t (•) = v t [ϕ t (•)],
with the initial condition ϕ 0 = Id. The time-dependent vector field v t can be interpreted as the instantaneous speed of the points during deformation, and must verify certain regularity conditions to ensure that solutions to the ODE are indeed diffeomorphisms. An efficient way to enforce these conditions is to restrict to vector fields obtained by the convolution of a number N c of momentum vectors carried by distinct control points:

v t (x) = Nc k=1 K(x, c (t) k )µ (t) K , ( 7
)
where K is a kernel function, e.g. the Gaussian kernel:

K(x, y) = exp -∥x-y∥ 2 σ 2
. The effect of the choice of kernel is studied in Micheli and Glaunès 2014, but we restrict to the Gaussian kernel in this study. The (closure of the) set of such vector fields forms a reproducing kernel Hilbert space H K , with the associated norm

∥v∥ 2 K = i,j K(c i , c j )µ T i µ j .
We write Diff K the subgroup of diffeomorphisms obtained this way and use matrix notation with a dN c × dN c block matrix K(c, c ′ ) = K(c i , c ′ j )I d ij and flat vectors of size dN c for landmarks sets, velocities and momentum, so that v(x) = K(c, x)µ and ∥v∥ 2 K = µ ⊤ K(c, c)µ. The total cost, or energy of the deformation can be defined as

1 0 ∥v t ∥ 2 K dt.
It can be shown that the momentum vectors that minimize this energy, considering c

(0) k , c (1)
k , k = 1 . . . N c fixed, together with the equation driving the motion of the control points, follow a Hamiltonian system of ODEs:

ċ(t) = K(c (t) , c (t) )µ (t) μ(t) = -∇ 1 K(c (t) , c (t) )µ (t) T µ (t) (8)
A diffeomorphism ϕ 1 is thus uniquely parameterized by the initial conditions c

(0) k , µ (0)
k , k = 1 . . . N c , and a shape registration criterion between a template q and target q can be defined as

C(c, µ) = ∥q -ϕ c,µ 1 (q)∥ 2 2 + α 2 ∥v c,µ 0 ∥ 2 K . ( 9
)
where α is a regularization parameter that penalises large deformations. Minimizing C therefore amounts to finding the transformation that best deforms q to match q.

The gradient of C can be computed through automatic differentiation, to perform gradient descent. Note that in (9), the L 2 norm between landmark positions is used to evaluate the data-attachment term. When corresponding landmarks are not available, fidelity metrics relying on currents, varifolds or normal cycles have been derived and allow to compute metrics between curves or surfaces with different parametrizations. These are reviewed in [START_REF] Charon | Fidelity Metrics between Curves and Surfaces: Currents, Varifolds, and Normal Cycles[END_REF]. Two examples of registration solutions are shown, the first is on simple parametric shapes, a circle and an ellipse (Figure 6). In this case the full trajectory of each landmark is represented (with a color scheme from blue to red). These differ from a linear interpolation especially where the curvature of the shape increases. The second example is on real RV data (Figure 7). The distribution of momentum vectors on the shape is not intuitive, and it is easier to inspect the velocity field evaluated at the control points. The matching is very efficient on these meshes, thanks to the quality of the correspondences between landmarks and to the smoothness of the shapes.

The optimal value of C defines a distance between ϕ 1 and the identity. In fact this distance derives from a right-invariant Riemannian metric on the group Diff K and the path t → ϕ t is a minimizing geodesic for this metric. Integrating the system (8) computes the exponential map of this metric, and is often called shooting in this section, while minimizing (9) approximates the logarithm, and is called registering. By considering the action of diffeomorphisms on shapes, this distance projects to a distance between the shapes q and q, and we have This in fact turns the space of shapes into a homogeneous space, with invariant metric given by the group.

d K (q, q) 2 = inf 1 2 1 0 qt K(c, c) -1 qt dt | q 0 = q, q 1 = q -2.

Ladder methods with LDDMM

Along with a distance, the Riemannian metric provides a notion of parallel transport thanks to its Levi-Civita connection. A Hamiltonian formulation of parallel transport may be given by an ODE that allows to transport a set of momentum vectors along a path of diffeomorphisms ϕ t (Younes 2019, section 13.3.3). Previous implementations of parallel transport of deformations in the context of LDDMM are however based on Jacobi fields [START_REF] Younes | Jacobi Fields in Groups of Diffeomorphisms and Applications[END_REF][START_REF] Cury | Spatio-Temporal Shape Analysis of Cross-Sectional Data for Detection of Early Changes in Neurodegenerative Disease[END_REF], and in particular on the fanning scheme [START_REF] Louis | A Fanning Scheme for the Parallel Transport Along Geodesics on Riemannian Manifolds[END_REF]. Ladder method in contrast, were proposed in the SVF framework. Adapting them to LDDMM is difficult because of the absence of an equivalent of the BCH formula to implement each rung of the ladder. We propose here a new implementation of pole ladder in the context of LDDMM that leverages the lower dimensional parametrization and the Hamiltonian formulation, and is grounded by the results of the previous section. Indeed, the geodesics are not known in closed-form and we are in the case of infinitesimal ladder schemes, requiring the use of fourth-order integration schemes. Our implementation is largely inspired by that of the fanning scheme and leverages the open-source software Deformetrica [START_REF] Bône | Deformetrica 4: An Open-Source Software for Statistical Shape Analysis[END_REF].

Throughout this section, we consider a reference shape q A , and wish to transport the deformation (c, µ) obtained by registration of a shape q 1 on another shape q 0 to q A , that is, along the geodesic between q 0 and q A , parametrized by (c A , ω) and also obtained by registration. We suppose that registration outputs a set of control points such that K(c, c) is a symmetric positive definite matrix.

Recall from the homogeneous space structure that the tangent space at q 0 is identified with the horizontal subspace of the Lie algebra of Diff K , which is here the reproducing Kernel Hilbert space H K that contains all vector fields defined by the Figure 8: Representation of the pole ladder for diffeomorphisms with an odd number of steps n. The exponential maps are computed with a RK4 scheme, and the log by gradient descent. convolution (7). It defines a one-to-one correspondence between the tangent space at q 0 , and its co-tangent space, meaning that we can transport the momentum µ and flow the points c instead of transporting a full vector field v. Similarly, we don't need to flow the full shape q 0 along the ladder, but only the control points.

We describe the pole ladder in the case where the momentum to transport µ and the momentum along which we transport ω are carried by the same initial control points c = c A . If it is not the case, µ is projected to c A by solving a linear system with a Cholesky decomposition of K(c A , c A ):

µ ′ = K(c A , c A ) -1 K(c, c A )µ.
Define rk : T S * × R + → T S * the map that performs one step of the fourthorder Runge-Kutta numerical scheme on the Hamiltonian system (8), and rk 1 the projection on the control points only, and rk -1

x its approximate inverse by gradient descent, such that ∥ rk 1 (rk -1

x (y), h) -y∥ ≤ h 5 . Then, pole ladder is performed by the following algorithm.

Choose n ∈ N and 1 ≤ α ≤ 2, and divide [0, 1] into n intervals of length h = 1 n . Compute the first midpoint c 1 , ω 1 = rk(c A , ω, h

2 ) along the main geodesic, and the first rung q 0 = rk 1 (µ, 1 n α ). Then iterate for i = 1 . . . n (see Figure 8) 1. Compute the momentum α i = n rk -1 c i (q i-1 );

2. Compute the flow of its inverse q i = rk 1 (c i , -α i , h);

3. Compute the next midpoint c i+1 , ω i+1 = rk(c i , ω i , h) except at the last step where h ← h 2 to compute c; Return μ = n α (-1) n rk -1 c (q n ). Two examples of parallel transport solutions are shown, the first on simple parametric shapes, the evolution from a circle to an ellipse is transported to a smaller circle (Figure 9). We retrieve a smaller ellipse as expected. The second example on real right ventricle (RV) mesh of the heart (Figure 10), shows the ED-to-ES deformation transported to a third reference shape, the atlas. The obtained reconstruction is a personalized estimate of an ES frame for the atlas, specific to the considered patient. The transported frame is much more acceptable than those of Figure 1.

Validation Recall that parallel transport is an isometry, so the norm of the transported deformation must equal that of the initial deformation. We use this property as a first step to validate our implementation on a population of simulated 2d-shapes. These are generated by shooting from a circle with random Gaussian momentum vectors. This defines the source shape. Then an affinity is applied to the source shape to define its evolution (i.e. scaling with different coefficients along each axis). This deformation is then estimated by registration and transported to the template. The reconstruction obtained by deforming the circle with the transported deformation (yellow shape on the figure) resembles an ellipse, as expected by the application of the affinity.

For n = 10 rungs, we obtain a root mean squared error (RMSE) of 8.3 × 10 -3 and less than 1.3 ± 0.9% relative error (in absolute value) for 100 samples. We also apply the affinity to the template and register the template to the result to compute the expected momentum after transport. The deviation of the transported momentum to this expected value is measured with the kernel norm of the difference and averaged over the samples. We obtain an error of 0.34 while the fanning scheme achieves 0.33. The two implementations are therefore very similar on this set of shapes. We however find cases where the fanning scheme does not behave as well as the pole ladder for more complex deformations on the heart data (Figure 13 below).

Application to cardiac motion modelling

Spatio-temporal shape analysis is of growing importance in the study of cardiac diseases. In particular, the assessment of cardiac function requires the measurement and analysis of cardiac motion beyond scalar indicators such as volume, ejection fraction, pressure-volume loops or area strain. We focus on the case of the right ventricle (RV), that is of particular interest as it has been shown to have a large capacity to adapt to overload by remodelling [START_REF] Sanz | Anatomy, Function, and Dysfunction of the Right Ventricle: JACC State-of-the-Art Review[END_REF], raising the issue of disentangling the deformation from the initial anatomy.

The Right Ventricle and its diseases

The right ventricle receives the blood from the right-atrium through the tricuspid valve and ejects it to the lungs through the pulmonary artery. Although neglected in the assessment of normal blood circulation of adult patients for decades, its relevance to determine symptoms and outcomes is now well established in several pathologies. Three such pathologies are studied in our database. Pulmonary hypertension (PH) is a severe condition characterized by increased blood pressure in the pulmonary arteries, appearing between the ages of thirty and forty. It is associated with symptoms and premature deaths [START_REF] Moceri | Three-Dimensional Right-Ventricular Regional Deformation and Survival in Pulmonary Hypertension[END_REF] and the most common source of pressure overload in the right-ventricle, resulting in hypertrophy, flattening of the interventricular septum, progressive dilation and dysfunction [START_REF] Sanz | Anatomy, Function, and Dysfunction of the Right Ventricle: JACC State-of-the-Art Review[END_REF].

Atrial septal defect (ASD) is a congenital heart defect in which blood flows between the left and right atrium. This may lead to lower oxygen levels in the blood that supplies the other organs, resulting in cyanosis. It affects the RV with volume overload, marked by dilation and hypertrophy and leftward septal shift, but it may be tolerated for years, as many studies report no affectation of the function of the RV under volume overload [START_REF] Moceri | 3D Right Ventricular Shape and Strain in Congenital Heart Disease Patients with Right Ventricular Chronic Volume Loading[END_REF][START_REF] Sanz | Anatomy, Function, and Dysfunction of the Right Ventricle: JACC State-of-the-Art Review[END_REF].

Tetralogy of Fallot (ToF) is another congenital heart defect characterized by pulmonary stenosis (i.e. narrowing of the pulmonary valve) and an overriding aorta, allowing blood flows between the left and right ventricles, and resulting in hypertrophy of the RV, reduced oxygen level and cyanosis. This is the unique case where the RV is under both volume and pressure overload. Untreated, ToF patients rarely survive to adulthood. It is treated by surgical repair during the first year of life, greatly improving survival, but post-surgery defects such as pulmonary regurgitation may occur, requiring other operations.

We use a database of 3D meshes extracted from 314 echocardiographic sequences from patients examined at the CHU of Nice (see [START_REF] Moceri | Three-Dimensional Right-Ventricular Regional Deformation and Survival in Pulmonary Hypertension[END_REF][START_REF] Moceri | 3D Right Ventricular Shape and Strain in Congenital Heart Disease Patients with Right Ventricular Chronic Volume Loading[END_REF] for details). The meshes were extracted with a commercial software (4D RV Function 2.0, TomTec Imaging Systems, GmbH, DE) with point-to-point correspondences across time and patients. These are formed of 938 points and 1872 triangles. All the shapes were realigned with a subject-specific rigid-body deformation. An atlas was computed from the end diastolic meshes of the control group, after rigid-body alignment. We warmly thank Pamela Moceri and Nicolas Duchateau for collecting and curating the data, and for the helpful interactions that we had during this project.

Two main measures used by the clinicians to evaluate the cardiac function will be used in the sequel to evaluate our modelling pipeline. Ejection fraction (EF) is the relative volume change between end-diastolic (ED) time and end-systolic (ES) time, i.e. during one contraction. It represents the volume of blood pumped by the heart to the body circulation, and is used in clinical routines to evaluate heart failures. Area strain (AS) is the relative area change of each cell of the mesh between ED time and ES time. As it depends on the quality of the mesh, it is usually filtered by computing the mean at each vertex over the neighboring cells, and visualized as a colormap on the mesh, or averaged by regions of the ventricle. It represents the amount of stretching local tissues undergo during the deformation, and is a common descriptor of cardiac motion [START_REF] Kleijn | Three-Dimensional Speckle Tracking Echocardiography for Automatic Assessment of Global and Regional Left Ventricular Function Based on Area Strain[END_REF][START_REF] Folco | Learning Interactions between Cardiac Shape and Deformation: Application to Pulmonary Hypertension[END_REF][START_REF] Moceri | 3D Right Ventricular Shape and Strain in Congenital Heart Disease Patients with Right Ventricular Chronic Volume Loading[END_REF].

Motion normalization with parallel transport

For each subject, the motion of the right ventricle during the systolic phase of the cardiac cycle is encoded by the LDDMM deformation at a discrete set of sampling times between ED and ES. Each time-point is treated as an independent observation for its reorientation by parallel transport into the atlas reference frame. The framework is summarized in Figure 11. For each observation time t i , the subject's frame at time t i is registered to this subject's ED frame, transported along the path between ED and the atlas, and reconstructed by shooting with the transported momenta.

We reproduce on Figure 12 the reconstructions of the ES frame for a patient of each disease group and of the control group. These are obtained by the pole ladder algorithm (right) as well as the fanning scheme (middle). On these examples the two methods produce very similar results, which validates our implementation. Moreover, the obtained reconstructions are smooth (compare with Figure 1) and seem to adapt the deformation pattern of each patient to the atlas.

In most cases, little difference can be observed qualitatively. We find however a few examples where the fanning scheme did not achieve to transport the deformation, and some points diverge or collapse in the reconstructions. This occurs when the original shapes are significantly larger than the atlas, and with strong spherification of the apex, causing large deformations as for the ASD patient of Figure 13. Figure 13: Examples of reconstruction where the fanning scheme did not achieve to transport the deformation, while pole ladder did, although it is not realistic.

Interaction between shape and deformations: a scale problem Niethammer and Vialard 2013 showed that parallel transport in LDDMM does not conserve global properties such as scale or volume changes. This is visible in our previous experiment (Figure 9), where the anisotropy of the original evolution (from circle to ellipse) is 1.4 whereas the one obtained after parallel transport is 2.3. In the case of cardiac deformations, the magnitude of the temporal deformation is comparable to that of the subject's reference to atlas deformation, and substantial volume changes are observed. Thus the lack of scale-invariance is crucial. Furthermore, there is a significant correlation between the magnitude of the systolic deformation and the ED volume (ρ = 0.42 in our dataset). Additionally, there are significant differences of ED volumes by disease groups (Figure 15a). As the parallel transport is isometric, this deformation may be too large for the atlas. Examples of the obtained ES frames are shown on Figure 14 for patients whose RV volume is greater than that of the atlas, resulting in non-realistic ES frames.

Patient deformation

Fanning Scheme Pole Ladder (ours) ASD n°28 To quantify this phenomenon, we measure the EF and the AS before and after transport and report the RMSE in Table 1. We observe a large relative RMSE EF, representing 31% of the mean value, and equal to the standard deviation of the data. Furthermore, plotting the absolute value of the EF error by disease group, we retrieve a strong correlation between the error and the EF (Figure 15c). This effect is quite undesirable for a normalization procedure!

In conclusion, parallel transport was able to reorient the deformation to a different frame, but was not sufficient to properly normalize it in the presence of large volume differences. Since we control the numerical accuracy of the parallel transport method thanks to 2.1, the error has to be attributed to the choice of the transport method or to the choice of the Riemannian metric used on deformations. In order to achieve proper normalization for a broader range of volume changes and to preserve the relative volume changes, we investigate in Subsection 3.3 a first intuitive ad-hoc strategy where we modulate the amplitude of the parallel transported vector. However, this amounts to state that the transport method and the Riemannian metric are not anymore consistent together. To remove this inconsistency with the general geometric statistics framework, we then investigate in Section 3.4 a modification of the metric on diffeomorphisms that decouples the volume change from the deformation.

An intuitive rescaling of LDDMM parallel transport

In [START_REF] Guigui | Numerical Accuracy of Ladder Schemes for Parallel Transport on Manifolds[END_REF], we hypothesises that the amplitude of the ED to atlas deformation, i.e. the one along which we transport, acts on the transported deformation by scaling to construct a normalization strategy invariant to initial volume changes.

Hypothesis We now propose to introduce a rescaling step after the parallel transport in our framework of Figure 11. More precisely, instead of using R t = Exp A (Π A ED Log ED (S t )) as the subject specific reconstruction of time-frame S t , we introduce a parameter λ > 0 and use

R t (λ) = Exp A λΠ A ED Log ED (S t ) .
where Π A ED is the parallel transport map along the geodesic that joins the ED frame to the atlas. Based on our observation of the results of parallel transport (Figures 12 and14) and of the relationship between the ejection fraction error and the disease group (Figure 15c), we hypothesise that this rescaling should be patient-specific and depend on the volume of the ED frame. Furthermore, rescaling should not be necessary if the volume of the ED frame matches that of the atlas. The relevant quantity is thus the ratio Vol ED / Vol A .

Criterion and Estimation of λ Given a patient's sequence, we now need a criterion to find λ. As discussed in the introduction, a clinically relevant quantity that should be conserved is the ejection fraction, as it is a straightforward scalar indicator of cardiac function. A perfect preservation of the ejection fraction between each time point and the ED frame constraints the volume of each reconstruction:

Vol new (t) = Vol A Vol ED Vol St .
A straightforward loss function for λ is therefore

L(λ) = k i=1 (Vol(R t i (λ)) -Vol new (t i )) 2 .
For each patient, we then minimize this loss function with respect to λ to find the patient-specific scaling parameter. We solve this problem by gradient descent, where the gradient is computed by automatic differentiation. The parameter λ is regularized to be close to 1 to avoid poor solutions. The scaled parallel transport achieves an error that is not related to the disease group.

Results

To validate the method, we first check that the ejection fraction is effectively conserved. The errors on ejection fraction and area strain are reported in Table 2. We indeed observe a three-fold decrease of the RMSE on the ejection fraction. Additionally, the distribution of this error is no longer related to the disease group (Figure 16). Visually, the reconstructions obtained by the scaled parallel transport are more realistic (Figure 17). Interestingly, it is possible to obtain a low error on the EF after the scaled transport, but this does not preserve the area strain (AS). This shows that although these quantities are related, they carry different information and the AS depends on the initial shape, itself related to the pathology.

Relationship between λ and Vol ED As expected, the scaling coefficient is closely related to the ED volume. We use a linear regression to identify a linear relation between log(λ) and log(Vol A / Vol ED ). This is displayed on Figure 18. The relationship seems highly significant as we measure a coefficient of determination R 2 = 0.92. The measured slope is α = 0.64 and corresponds to the relationship

λ = β Vol A Vol ED α . ( 10
)
where β is the exponential of the intercept, but is found close to 1, with β = 1.08 In order to better understand this relationship, we reproduce the experiment on 2d simulated shapes, with the same experimental set-up as in the validation experiment. It is striking that we also find a highly significant relationship with R 2 = 0.99 and almost the same slope α = 0.67 and β = 1.02. The values measured for α are close to 2/3 and surprisingly do not seem to depend on the dimension of the ambient space.

Our experiments thus give strong evidence in favour of relation ( 10), and raise many questions. First is α a constant, or does it depend on the data, the dimension, the kernel, or the deformation model? From our results it seems that it does not depend on the data or the dimension. Performing the same experiments with different kernels could bring further information, this is ongoing work. Similarly, experiments in the SVF framework could be performed to test whether this effect is specific to LDDMM or more general. These experiments could also clarify whether β significantly differs from one.

Furthermore, this effect suggests that all the experiments that have been performed with parallel transport should be revised with a scaling step. Most of these in the medical imaging field concern the evolution of brain images [START_REF] Cury | Spatio-Temporal Shape Analysis of Cross-Sectional Data for Detection of Early Changes in Neurodegenerative Disease[END_REF][START_REF] Lorenzi | Disentangling Normal Aging from Alzheimer's Disease in Structural Magnetic Resonance Images[END_REF]. However, longitudinal deformations of the brain are very different from that of the right ventricle, and we don't expect the scaling step to be necessary. Indeed, brain deformations are characterized by very localized strong Jacobian determinant, but very little overall volume changes. A more relevant type of data to reproduce this effect could be to evaluate tumor evolution during treatment, or the evolution of more plastic organs such as the liver, whose shape is affected by all the liver diseases, whether viral, metabolic or due to alcohol.

Finally, a last hypothesis is that by scaling the parallel transport of the LDDMM metric, we are approximating the parallel transport of a direct product metric with a volume preserving factor and a metric on the volume changes as will be introduced in the next section. This relation would translate in setting the metric on the volume part.

Changing the metric to preserve relative volume changes

Although modifying the parallel transport equations works astonishingly well from the practical point of view, it is not satisfactory from a theoretical point of view since it decorrelates the transport from the metric used on the deformations. We investigate in this section a more satisfactory strategy that decouples the volume change from the deformation directly within the LDMMM metric. This strategy is based on the scale invariant metric proposed by Niethammer and Vialard 2013 and it was implemented thanks to the precious advice of François-Xavier Vialard. More generally, they show that a Riemannian metric on the space of shapes preserves a global non-degenerate function f (such as volume) if and only if it can be decomposed into a product metric on Im(f ) and f -1 (0) (see Niethammer and Vialard 2013, Theorem 41).

Model We apply their result as follows. Restricting to the space of shapes that are embeddings of the circle in R 2 or of the sphere in R 3 , let M = Emb(S d , R d+1 ) be the set of shapes, endowed with a metric g. As we are interested in relative volume changes, we consider the function

f = log Vol Vol A : M → R (so that df = d Vol Vol ),
and let M 0 = f -1 (0), the space of shapes whose volume equals the volume of the atlas Vol A . Then, by Niethammer and Vialard 2013, Theorem 41, (M, g) can be decomposed into a direct product of Riemannian manifolds (M, g) = (R, dt 2 ) × (M 0 , g 0 ) where g 0 is a Riemannian metric on the submanifold M 0 . This metric can be constructed by choosing a projection π : M → M 0 and the restriction of g to M 0 . However, there is no canonical projection and two schemes are proposed in [START_REF] Niethammer | Riemannian Metrics for Statistics on Shapes : Parallel Transport and Scale Invariance[END_REF] by gradient flow: we follow the flow of (d Vol) # where # depends on the choice of metric on M , in this case the LDDMM metric.

by scaling: we centre the shape around 0 and divide all the landmarks by Vol

1/3
q . This choice depends on the centre (barycenter) of the shape which may be unnatural.

In both cases, the framework is slightly modified by applying this projection first to all the time frames, so that their volume matches that of the atlas. Then the previous framework of Figure 11 is applied with volume preserving geodesics and parallel transport. This corresponds to the parallel transport on the factor M 0 of M = R × M 0 . Finally, the (Euclidean) transport on R corresponds to applying the vector log

Vol S t
Vol ED to 0, so that the new volume is Vol new (t) = Vol A Vol ED Vol St as in the previous section, and the ejection fraction is preserved by construction. This volume is obtained either by scaling of by gradient flow, consistently with the projection.

Implementation We now give more details on the computation of geodesics and our implementation in the LDDMM framework of Section 2.2. Recall that tangent spaces are described by the Hilbert space H K of vector fields obtained by convolution. As we restrict to M 0 = f -1 (0), a vertical tangent subspace at any shape q is defined as the set of vector fields of H K that are volume preserving, i.e. ker(df q ) = {v ∈ H K |d Vol q ((v(q)) = 0}, To avoid any confusion, we distinguish the linear form d Vol q from its representation in the canonical basis ∂ Vol q , such that d Vol q (v(q)) = v(q) ⊤ ∂ Vol q . The orthogonal projection on V q (which depends on the LDDMM metric) can be computed by

π q : T q M -→ V q v -→ v -d Vol q (v)⃗ n ( 11
)
where

⃗ n = (d Volq) # ∥(d Volq) # ∥ 2 .
Note that the map # associates a tangent vector to the linear form d Vol such that d Vol(v) = ⟨(d Vol q ) # , v⟩ K , so that (d Vol q ) # = K(•, q)∂ Vol q . Define also the dual π * q : T * q M → V * q to π q by the relation µ(π q (v)) = π * q (µ)(v), i.e. in matrix notations µ ⊤ π q (v) = v ⊤ π * q (µ). From (11), we obtain

π * q : T * q M -→ V * q (c, µ) -→ c, µ - µ ⊤ K(c, q)∂ Vol q ∂ Vol ⊤ q K(q, q)∂ Vol q d Vol q . ( 12
)
Geodesics A geodesic between q 0 , q 1 ∈ Emb 1 (S 2 , R 3 ) is solution of:

inf 1 0 ∥v t ∥ 2 K dt (13)
under the constraints q = v t (q(t)) where v t ∈ V q(t) . The constraints can be relaxed to q = π q(t) (v t (q(t)) where v t ∈ V . This allows solving Problem (13) by solutions of Hamiltonian equations given by H(µ, q) = 1 2 < π * q (µ), K(q)π * q (µ) >. In order to implement these geodesics in the LDDMM framework presented in Section 2.2, we consider that both the control points and the landmark points are part of the system, and write the Hamiltonian

H(µ, q, c) = 1 2 µ T K(c, c)µ - (µ T K(c, q)∂ Vol q ) 2 2∥(d Vol q ) # ∥ 2 K .
The corresponding Hamiltonian system can readily replace (8) for registration and in the geodesics of the pole ladder constructions. We can thus conduct the same experiment as in the previous section with this metric.

Results

We first validate our implementation by measuring the area strain and ejection fraction errors, and complete Table 2 into Table 3. Similarly, we reproduce the boxplot of Figure 15c with the two projection and the previous volume constrained scaling strategy. As in the previous section, there is no relation between the error and the disease group. We show some reconstructions on Figure 19, along the reconstructions obtained by the scaled parallel transport of the previous section. The two metrics, and the two projection methods result in different reconstructions. The projection by gradient flow seems less stable than the two other methods, as the thinning near the valve is exaggerated on the first row, and the apex of the second row is unrealistically spherical. The values of area strain obtained by this normalization procedure were included in the analysis of Di Folco et al. that leveraged dimensionality reduction methods to compare the different normalization procedures and explore the interactions between shape and deformations. Shape descriptors were computed by point-wise differences between the ED mesh and the atlas after a rigid-body alignment, or alignment by a similarity. By comparing the latent spaces obtained with normalized and unnormalized data, they show that both methods reduce the bias introduced by the alignment of the ED shapes [START_REF] Folco | Investigation of the Impact of Normalization on the Study of Interactions Between Myocardial Shape and Deformation[END_REF]). Now we come to the original goal of this chapter: performing group-wise analysis of the systolic deformation. We have proposed several ways to normalize each individual time deformation, we now need to summarize the complete systolic trajectory. We propose to use a geodesic -or higher order-regression. This is detailed in the next section dedicated to the downstream analyses.

Analysis of the normalized deformations

In order to compactly represent the full systolic trajectory in the same space and to proceed with linear statistics, we propose to fit a regression model in the shape space as considered in the registration framework of Section 2.2. It estimates a geodesic path between two shapes (the equivalent of a uniform motion). For a trajectory such as the contraction of the cardiac RV, one may expect to find second-order dynamics, making a (first-order) geodesic regression ill-suited. We thus propose to use the second-order model defined in Trouvé and Vialard 2012a to account for the motion of the RV during systole. It is similar to a geodesic regression [START_REF] Fishbaugh | Geodesic Shape Regression with Multiple Geometries and Sparse Parameters[END_REF] in the shape space, with an additional acceleration term in the model of the trajectory. We therefore compare it to the simpler geodesic regression.

Geodesic and Spline regression

The model of (Trouvé and Vialard 2012a) introduces a second-order terms u t can be interpreted as a random external forces smoothly perturbing the trajectory around a mean geodesic which modifies the continuous-time system of equations (8) as follows: ∀t ∈ [0, 1], ċk

(t) = j K(c (t) k , c (t) j )µ (t) j μk (t) = -j ∇ 1 K(c (t) k , c (t) j )µ (t) T k µ (t) j + u (t) k (14)
Results If we consider a discrete sequence of observation times t 1 = 0, t 2 , . . . t d = 1 and configurations x t 1 , . . . , x t d , one seeks to find the path ϕ t that minimizes the new cost

C S (c, µ, u t ) = 1 α 2 d d i=1 ∥x t i -ϕ t i (x t 0 )∥ 2 2 + 1 0 ∥u (t) ∥ 2 dt + ∥v c,µ 0 ∥ 2 K . (15)
In practice, the ODEs ( 8) and ( 14) are discretized in n time steps and an integration method such as Euler or Runge-Kutta is used. We define all the patients trajectories between t = 0 and t = 1, and use the same discretization for all the patients to ensure that u 0 , . . . , u Nc are estimated at corresponding times. Along with µ (0) , these are estimated by gradient descent as in the case of registration. Setting u 0 = . . . = u Nc = 0 at all times recovers a geodesic trajectory. We use a kernel bandwidth σ = 15 in all the experiments, and 60 control points for all the deformations of the atlas. The initial control points are fixed for the entire dataset so that the initial momenta can be compared consistently. They have been optimized to register the atlas on all the transported ES frames. Visually we only notice very slight differences between the geodesic and spline regressions (Figure 20). We do observe that the fit is slightly less precise in the intermediate frames for the geodesic regression, this is confirmed by measuring the overall data attachment term (left term of (15)) that is reported in Table 4, regardless of the normalization method. The spline regression therefore yields a more faithful representation of the normalized cardiac deformations than the geodesic regression, at the cost of a larger set of parameters, encompassing the external timedependent forces u t whose interpretation and analysis is difficult. Indeed, the size of the acceleration term u t is (d -1) × N c × 3 = 1620.

Hotelling tests on velocities

The sequences are normalized by one of the three proposed methods: scaled parallel transport (SPT), volume-preserving parallel transport (VPPT) with gradient projection and VPPT with scaling projection, and summarized by a regression with one of the two deformation models: geodesic or spline. In this section, we focus on the group mean velocity fields, and show that these differences are mainly of magnitude, with different orientations at a few points along the free wall. The differences observed near the tricuspid valve mainly reflect the difference of magnitude of the deformations, and one should be cautious before drawing further conclusions as the quality of the mesh may vary near this region. However it is interesting to notice that very little differences are observed for the ASD group, which corroborates previous results [START_REF] Moceri | 3D Right Ventricular Shape and Strain in Congenital Heart Disease Patients with Right Ventricular Chronic Volume Loading[END_REF]. Similarly, only small differences are observed on the septum for the PHT group, showing that the shape differences usually observed on the PHT group have been filtered out. This makes the differences observed on the free wall interesting and other markers such as the circumferential strain will be studied to confirm these effects. These experiments also highlight the difference between the PHT and ToF group as significant differences are localized on the inferior part of the free wall and on the inlet (area of the tricuspid valve) for the Tof group while they are distributed across the whole shape for the PHT group. The same visualization is provided for the mean acceleration term by disease group, across the group-wise mean trajectory (Figure 22). Recall that these can be interpreted as external forces. They seem to increase the movement towards contraction at the beginning of motion while slowing it down towards the end. These terms thus provide more insights into the dynamics of each disease.

Similarly to the analysis of the velocities, the mean external force of the ASD group is not significantly different from that of the control group. This shows that the volume differences reported in Figure 15a were successfully normalized by the scaled parallel transport. For the ToF and PHT group, many points have significantly different values of acceleration, and those differences are shown by red arrows. They are distributed across the whole RV for the PHT group while their locations vary in time for the ToF group, and concentrate around the apex towards the end of the deformation. For the ToF group, we do find differences on the infandibular (area of the pulmonary valve) at times ED + 1 and ED + 2 as expected from the consequences of the surgery. This is in favor of our method on the analysis of shape deformations, as these differences were not observed with the traditional analysis of strain maps. As in the case of velocities, the differences mainly reflect amplitude differences, especially near ED and ES times. They thus reflect the late contraction of the RV and the distribution of the differences across the ventricle and in time reflect a longer and asynchronous contraction introduced by the disease. Indeed, the deformation of the control group is more uniformly distributed on the shape, and the introduction of heterogeneity across the shape has been previously reported and associated with the disease, and is known to be a factor of arrhythmia with direct consequences on survival.

Conclusion

In this chapter, we first summarized recent results using parallel transport for geometric statistics in computational anatomy. We proposed in particular a new implementation of the pole ladder scheme. These results guarantee that parallel transport algorithms are well behaved, enabling their use as a normalization step for longitudinal shape data. This further allows to evaluate the underlying metric and the normalization model beyond the numerical stability. Where the differences are significant, control points and differences with the control group are in red, while coloured arrows represent the mean acceleration for the disease group.

As for the velocity fields, there are no differences for the ASD group. For the Tof group, the differences are localized, but these sites vary with time, while for the PHT group, differences are distributed on the entire shape.

On the application side, the results on data of the motion of the cardiac right ventricle exhibited a strong bias due to the volume differences at a reference time point. Two strategies were investigated to correct for this bias: scaling the parallel transport, and using a metric that decomposes volume changes from volume preserving deformations. Both methods were successful at preserving the ejection fraction, and reduced the bias due to volume differences. When scaling the parallel transport, we discovered a significant relationship between the scaling parameter and the initial volume ratio, however its interpretation must be understood further. This will be investigated in future work, together with the possible mathematical relation to the volume preserving metric. The statistical analysis of the normalized deformations was meaningful and coherent with previous knowledge. This framework is readily usable to evaluate the impact of a treatment or a surgery on the deformation. It can also be used to simulate new cardiac sequences for a given shape from a population of cardiac sequences, by performing e.g. a principal component analysis of the momentum vectors of the spline or geodesic regression, and shooting along the principal modes, or by sampling coefficients from a multivariate normal law to mix the principal components.

This study is however not sufficient to decide which method should be preferred. This comes back to the more general question of choosing the metric. As the numerical properties of the parallel transport algorithm are now well grounded, parallel transport can now be used as a proxy to evaluate the pertinence of the metric used. In this implementation of the LDDMM framework, the choice of metric boils down to the choice of kernel, but the framework itself is questioned. As data-driven methods have gained popularity in this field, more papers have proposed strategies to learn the metric from data [START_REF] Vialard | Spatially-Varying Metric Learning for Diffeomorphic Image Registration: A Variational Framework[END_REF][START_REF] Louis | Riemannian Geometry Learning for Disease Progression Modelling[END_REF][START_REF] Niethammer | Metric Learning for Image Registration[END_REF]. These solutions are attractive to remove the arbitrary aspect of the choice of metric. Yet, beyond the computational complexity, these approaches face the problem of defining 'natural' criteria to assess a model. It will be interesting to test the parallel transport normalization method with a metric learned from data.

Figure 1 :

 1 Figure 1: Example of a systolic deformation estimated on a patient and applied to the atlas without normalization (middle) and with normalization by the method of Section 3.3 (right). The ED frame is a red point cloud, while the ES frame is the blue mesh. The ES frames obtained without normalization show irregularities and cannot be used in downstream analyses.

Figure 3 :

 3 Figure 3: Statistics on diffeomoprhisms with the Cartan-Schouten connection to model the normal component of the aging trajectory and the additional component specific to Alzheimer's disease. Longitudinal geodesic deformations regressed in the sequence of images of each subject are parallel transported along the subject-toreference deformation at baseline. In this common space, a linear model in the tangent space of diffeomorphisms estimates the mean trajectory for the two clinical conditions. Derived from original images and illustrations of Marco Lorenzi and Raphaël Sivera.

Figure 4 :

 4 Figure 4: Schild's ladder procedure to parallel transports the vector v along the sampled curve with initial velocity w. Top: First rug of the ladder using an approximate geodesic parallelogram. Bottom: The method is iterated with rungs at each point sampled along the curve. Figures reproduced from Guigui and Pennec 2021

Figure 5 :

 5 Figure 5: Top: Elementary construction of the pole ladder with the previous notations (left), and new notations (right), in a normal coordinate system at m. Bottom: Iterated scheme. Figures reproduced from Guigui and Pennec 2021.

Figure 6 :

 6 Figure6: Example of registration and optimal path between a circle (source) and an ellipse (target) with 60 landmarks, 15 control points (green) and corresponding momemtum vectors (cyan) and kernel width 1. The progression from blue to red shows the time trajectories of the landmarks under the minimizing geodesic deformation.

  Figure 7: Example of the registration between ED (blue wireframe) and ES (red points) for a control right-ventricle. There are 938 landmarks, 34 control points with kernel width 10mm (the overall ventricle at ED is about 55mm high). Note how momentum vectors are noisy due to strong interactions between control points (left).The velocity field at the control points offers a better summary of the deformation (right).

Figure 9 :

 9 Figure 9: Parallel transport of the evolution from a circle to an ellipse along the deformation to a smaller circle.

Figure 10 :

 10 Figure 10: Example ED-to-ES deformation transported along the ED-to-atlas path.The source frames are represented by a blue wireframe, the final frames by red landmarks, the control points by green landmarks and the initial velocities at the control points by orange arrows. The transported frame is much more acceptable than those of Figure 1.

Figure 11 :

 11 Figure 11: Framework using registration and parallel transport to normalize cardiac deformations.

Figure 12 :

 12 Figure 12: Examples of reconstructions of the ES frame (blue meshes) after parallel transport of the ED-to-ES deformation (point cloud to mesh) along the deformation from ED to atlas. Parallel transport is computed by the fanning scheme (middle) and our implementation of pole ladder (right). In most cases, little difference can be observed visually.

Figure 14 :

 14 Figure 14: Examples of reconstructions of the ES frame (blue meshes) after parallel transport of the ED to ES deformation (point cloud to mesh) along the deformation from ED to atlas. The subjects belong to the ASD group, hence show large volume changes that result in unrealistic ES frames.

Figure 15 :

 15 Figure 15: End-diastolic volume by disease group (left), original ejection fraction (middle) and its modification by parallel transport (right).

Figure 16 :

 16 Figure16: Alteration of the ejection fraction per disease group. The scaled parallel transport achieves an error that is not related to the disease group.

OriginalFigure 17 :

 17 Figure 17: Examples of reconstructions of the ES frame (blue meshes) after parallel transport of the ED-to-ES deformation (point cloud to mesh) along ED-to-atlas deformation. The large volume changes between subjects and the atlas result in unrealistic ES frames (middle column). This problem is well addressed by the scaling strategy (right column).

Figure 18 :

 18 Figure 18: Scaling Parameter λ with respect to the ED volume.

Figure 19 :

 19 Figure 19: Examples of reconstructions of the ES frame (blue meshes) after parallel transport of the ED-to-ES deformation (point cloud to mesh) along the deformation from ED to atlas. Comparison of the three methods. The gradient flow seems less stable (first two rows) than the projection by scaling. The scaling method results in more vertical movement of the base than the two other methods. This may be due to the centering step required before scaling.

Figure 20 :

 20 Figure 20: Fit of the normalized sequence of an ASD patient by geodesic (top) and spline (bottom) regression. Fits are the purple transparent meshes, normalized data are the blue meshes with edges. Upon visual inspection, there is barely any difference between the two fits.

Figure 22 :

 22 Figure 22: Results of the Hotelling tests between accelerations.Where the differences are significant, control points and differences with the control group are in red, while coloured arrows represent the mean acceleration for the disease group. As for the velocity fields, there are no differences for the ASD group. For the Tof group, the differences are localized, but these sites vary with time, while for the PHT group, differences are distributed on the entire shape.

Table 1 :

 1 Validation of parallel transport with EF and AS.

		Original Values RMSE PT
	AS EF	-0.24 ± 0.08 0.42 ± 0.13	0.18 0.13

Table 3 :

 3 Validation of the volume preserving parallel transport (VPPT), with the projection performed either by gradient flow or by scaling.

		Original Values PT SPT grad VPPT scaling VPPT
	AS EF	-0.24 ± 0.08 0.18 0.13 0.42 ± 0.13 0.13 0.04	0.13 0.02	0.09 0.02

Table 4 :

 4 Comparison of the mean data attachment term for each method.

	SPT grad VPPT scaling VPPT
	Geodesic 1384 ± 952 1659 ± 906 Spline 549 ± 413 681 ± 390	1483 ± 864 554 ± 324

Acknowledgements This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant G-Statistics agreement No 786854). This work has been supported by the French government, through the 3IA Côte d'Azur Investments in the Future project managed by the National Research Agency (ANR) with the reference number ANR-19-P3IA-0002.

the result obtained with the scaled parallel transport and spline regression. This procedure results in a set of descriptive parameters for each individual that lie in a linear space: the co-tangent space at the atlas. For more stable descriptors of the deformation, we compute the initial velocity vectors at the initial control points with v = K(c, c)µ. We now have all the ingredients to perform statistics on cardiac deformations described by the initial velocities.

We compare the impact of the different diseases by performing pairwise Hotelling tests between each disease group and the control group. One statistic is computed at each control point, that tests if the mean velocity at that point is significantly different. A Bonferroni correction is then applied for multiple testing across the full set of control points, to maintain type I error risk at α = 0.05.

The results for the velocities of the spline regression with the scaled PT are displayed Figure 21 for the ASD, ToF and PHT groups and show significant differences between each disease and the control group. The differences are superimposed with Where the differences are significant, control points and difference with the control group are in red, while coloured arrows represent the mean velocity field for the disease group. All the arrows have been scaled by a factor 2.3 for visualization purposes. The color map of the meshes reflects the norm of the velocity field at that point, if it is significantly different from the control group. As for the area strain, there is little difference for the ASD group, showing that the volume differences have well been filtered out by the normalization. Moreover the differences for the PHT and ToF group reflect deformations with less amplitude than the control group, as observed on the ejection fraction.