Parallel transport in shape analysis: a scalable numerical scheme

Publication GSI2017
contenu protégé  Document accessible sous conditions - vous devez vous connecter ou vous enregistrer pour accéder à ou acquérir ce document.
- Accès libre pour les ayants-droit


The analysis of manifold-valued data requires efficient tools from Riemannian geometry to cope with the computational complexity at stake. This complexity arises from the always-increasing dimension of the data, and the absence of closed-form expressions to basic operations such as the Riemannian logarithm. In this paper, we adapt a generic numerical scheme recently introduced for computing parallel transport along geodesics in a Riemannian manifold to finite-dimensional manifolds of diffeomorphisms. We provide a qualitative and quantitative analysis of its behavior on high-dimensional manifolds, and investigate an application with the prediction of brain structures progression.

Parallel transport in shape analysis: a scalable numerical scheme


application/pdf Parallel transport in shape analysis: a scalable numerical scheme Maxime Louis, Alexandre Bône, Benjamin Charlier, Stanley Durrleman
Détails de l'article
contenu protégé  Document accessible sous conditions - vous devez vous connecter ou vous enregistrer pour accéder à ou acquérir ce document.
- Accès libre pour les ayants-droit

Parallel transport in shape analysis: a scalable numerical scheme


Voir la vidéo


493.4 Ko


Creative Commons Aucune (Tous droits réservés)


Sponsors Platine


Sponsors Bronze


Sponsors scientifique





<resource  xmlns:xsi=""
        <identifier identifierType="DOI">10.23723/17410/22572</identifier><creators><creator><creatorName>Stanley Durrleman</creatorName></creator><creator><creatorName>Maxime Louis</creatorName></creator><creator><creatorName>Alexandre Bône</creatorName></creator><creator><creatorName>Benjamin Charlier</creatorName></creator></creators><titles>
            <title>Parallel transport in shape analysis: a scalable numerical scheme</title></titles>
        <resourceType resourceTypeGeneral="Text">Text</resourceType><dates>
	    <date dateType="Created">Thu 8 Mar 2018</date>
	    <date dateType="Updated">Thu 8 Mar 2018</date>
            <date dateType="Submitted">Fri 10 Aug 2018</date>
	    <alternateIdentifier alternateIdentifierType="bitstream">1c5619df3cd4fabd50142c34e69e5f5bbe29df31</alternateIdentifier>
            <description descriptionType="Abstract">The analysis of manifold-valued data requires efficient tools from Riemannian geometry to cope with the computational complexity at stake. This complexity arises from the always-increasing dimension of the data, and the absence of closed-form expressions to basic operations such as the Riemannian logarithm. In this paper, we adapt a generic numerical scheme recently introduced for computing parallel transport along geodesics in a Riemannian manifold to finite-dimensional manifolds of diffeomorphisms. We provide a qualitative and quantitative analysis of its behavior on high-dimensional manifolds, and investigate an application with the prediction of brain structures progression.

Parallel transport in shape analysis: a scalable numerical scheme Maxime Louis†12 , Alexandre Bône†12 , Benjamin Charlier23 , Stanley Durrleman12 , and the Alzheimer’s Disease Neuroimaging Initiative 1 Sorbonne Universités, UPMC Université Paris 06, Inserm, CNRS, Institut du Cerveau et de la Moelle (ICM) – Hôpital Pitié-Salpêtrière, 75013 Paris, France, 2 Inria Paris, Aramis project-team, 75013 Paris, France, 3 Université de Montpellier, France Abstract. The analysis of manifold-valued data requires efficient tools from Riemannian geometry to cope with the computational complexity at stake. This complexity arises from the always-increasing dimension of the data, and the absence of closed-form expressions to basic operations such as the Riemannian logarithm. In this paper, we adapt a generic numerical scheme recently introduced for computing parallel transport along geodesics in a Riemannian manifold to finite-dimensional manifolds of diffeomorphisms. We provide a qualitative and quantitative analysis of its behavior on high-dimensional manifolds, and investigate an appli- cation with the prediction of brain structures progression. 1 Introduction Riemannian geometry is increasingly meeting applications in statistical learning. Indeed, working in flat space amounts to neglecting the underlying geometry of the laws which have produced the considered data. In other words, such a sim- plifying assumption ignores the intrinsic constraints on the observations. When prior knowledge is available, top-down methods can express invariance properties as group actions or smooth constraints and model the data as points in quotient spaces, as for Kendall shape space. In other situations, manifold learning can be used to find a low-dimensional hypersurface best describing a set of observations. Once the geometry has been modeled, classical statistical approaches for constrained inference or prediction must be adapted to deal with structured data, as it is done in [4,5,11,13]. Being an isometry, the parallel transport arises as a natural tool to compare features defined at different tangent spaces. In a system of coordinates, the parallel transport is defined as the solution to an ordinary differential equation. The integration of this equation requires to compute the Christoffel symbols, which are in general hard to compute –e.g. in the case of the Levi-Civita connection– and whose number is cubic in the dimension. The Schild’s ladder [5], later improved into the Pole ladder [7] when transporting along geodesics, is a more geometrical approach which only requires † Equal contributions. the computation of Riemannian exponentials and logarithms. When the geodesic equation is autonomous, the scaling and squaring procedure [6] allows to com- pute exponentials very efficiently. In Lie groups, the Baker-Campbell Haussdorff formula allows fast computations of logarithms with a controlled precision. In such settings, the Schild’s ladder is computationnally tractable. However, no theoretical study has studied the numerical approximations or has provided a convergence result. In addition, in the more general case of Riemannian mani- folds, the needed logarithm operators are often computationally intractable. The Large Deformation Diffeomorphic Metric Mapping (LDDMM) frame- work [1] focuses on groups of diffeomorphisms, for shape analysis. Geodesic trajectories can be computed by integrating the Hamiltonian equations, which makes the exponential operator computationally tractable, when the logarithm remains costly and hard to control in its accuracy. In [12] is suggested a numeri- cal scheme which approximates the parallel transport along geodesics using only the Riemannian exponential and the metric. The convergence is proved in [8]. In this paper, we translate this so-called fanning sheme to finite-dimensional manifolds of diffeomorphisms built within the LDDMM framework [2]. We pro- vide a qualitative and quantitative analysis of its behavior, and investigate a high-dimensional application with the prediction of brain structures progression. Section 2 gives the theoretical background and the detailed steps of the algo- rithm, in the LDDMM context. Section 3 describes the considered application and discusses the obtained results. Section 4 concludes. 2 Theoretical background and practical description 2.1 Notations and assumptions Let M be a finite-dimensional Riemannian manifold with metric g and tangent space norm k · kg. Let γ : t → [0, 1] be a geodesic whose coordinates are known at all time. Given t0, t ∈ [0, 1], the parallel transport of a vector w ∈ Tγ(s)M from γ(t0) to γ(t) along γ will be noted Pt0,t(w) ∈ Tγ(t)M. We recall that this mapping is uniquely defined by the integration from u = t0 to t of the differential equation ∇γ̇(u)Pt0,u(w) = 0 with Pt0,t0 (w) = w where ∇ is the Levi- Civita covariant derivative. We denote Exp the exponential map, and for h small enough we define Jw γ(t)(h), the Jacobi Field emerging from γ(t) in the direction w ∈ Tγ(t)M by: Jw γ(t)(h) = ∂ ∂ε ε=0 Expγ(t) h [γ̇(t) + εw]  ∈ Tγ(t+h)M. (1) 2.2 The key identity The following proposition relates the parallel transport to a Jacobi field [12]: Proposition. For all t > 0 small enough and w ∈ Tγ(0)M, we have: P0,t(w) = Jw γ(0)(t) t + O t2  . (2) Proof. Let X(t) be the time-varying vector field corresponding to the parallel transport of w, i.e. such that Ẋi + Γi klXl γ̇k = 0 with X(0) = w. At t = 0, in normal coordinates the differential equation simplifies into Ẋi (0) = 0. Besides, near t = 0 in the same local chart, the Taylor expansion of X(t) writes Xi (t) = wi + O t2  . Noticing that the ith normal coordinate of Expγ(0) (t [γ̇(t) + εw]) is t(vi 0 + εwi ), the ith coordinate of Jw γ(0)(t) = ∂ ∂ε |ε=0Expγ(0) (t [γ̇(0) + εw]) is therefore twi , and we thus obtain the desired result. Subdividing [0, 1] into N intervals and iteratively computing the Jacobi fields 1 N Jw γ(k/N)( 1 N ) should therefore approach the parallel transport P0,1(w). With an error in O 1 N2  at each step, a global error in O 1 N  can be expected. We propose below an implementation of this scheme in the context of a manifold of diffeomorphisms parametrized by control points and momenta. Its convergence with a rate of O 1 N  is proved in [8]. 2.3 The chosen manifold of diffeomorphisms The LDDMM-derived construction proposed in [2] provides an effective way to build a finite-dimensional manifold of diffeomorphims acting on the d-dimensional ambient space Rd . Time-varying vector fields vt(.) are generated by the convolu- tion of a Gaussian kernel k(x, y) = exp  −kx − yk2 /2σ2  over ncp time-varying control points c(t) = [ci(t)]i, weighted by ncp associated momenta α(t) = [αi(t)]i, i.e. vt(.) = Pncp i=1 k [. , ci(t)] αi(t). The set of such vector fields forms a Repro- ducible Kernel Hilbert Space (RKHS). Those vector fields are then integrated along ∂tφt(.) = vt[φ(.)] from φ0 = Id into a flow of diffeomorphisms. In [10], the authors showed that the kernel- induced distance between φ0 and φ1 –which can be seen as the deformation kinetic energy– is minimal i.e. the obtained flow is geodesic when the control points and momenta satisfy the Hamiltonian equations : ċ(t) = Kc(t)α(t), α̇(t) = − 1 2 gradc(t)  α(t)T Kc(t) α(t) , (3) where Kc(t) is the kernel matrix. A diffeomorphism is therefore fully parametrized by its initial control points c and momenta α. Those Hamiltonian equations can be integrated with a Runge-Kutta scheme without computing the Christoffel symbols, thus avoiding the associated curse of dimensionality. The obtained diffeomorphisms then act on shapes embedded in Rd , such as images or meshes. For any set of control points c = (ci)i∈{1,..,n}, we define the finite-dimensional subspace Vc = span  k(., ci)ξ | ξ ∈ Rd , i ∈ {1, .., n} of the vector fields’ RKHS. We fix an initial set c = (ci)i∈{1,..,n} of distinct control points and define the set Gc = {φ1 | ∂tφt = vt ◦ φt, v0 ∈ Vc , φ0 = Id}. Equipped with Kc(t) as –inverse– metric, Gc is a Riemannian manifold such that Tφ1 Gc = Vc(1), where for all t in [0, 1], c(t) is obtained from c(0) = c through the Hamiltonian equations (3) [9]. 2.4 Summary of the algorithm We are now ready to describe the algorithm on the Riemannian manifold Gc. Algorithm. Divide [0, 1] into N intervals of length h = 1 N where N ∈ N. We note ωk the momenta of the transported diffeomorphism, ck the control points and αk the momenta of the geodesic γ at time k N . Iteratively : (i) Compute the main geodesic control points ck+1 and momenta αk+1, using a Runge-Kutta 2 method. (ii) Compute the control points c±h k+1 of the perturbed geodesics γ±h with initial momenta and control points (αk ±hωk, ck), using a Runge-Kutta 2 method. (iii) Approximate the Jacobi field Jk+1 by central finite difference : Jk+1 = c+h k+1 − c−h k+1 2h . (4) (iv) Compute the transported momenta ω̃k+1 according to equation (2) : Kck+1 ω̃k+1 = Jk+1 h . (5) (v) Correct this value with ωk+1 = βk+1ω̃k+1 +δk+1αk+1, where βk+1 and δk+1 are normalization factors ensuring the conservation of kωkVc = ωT k Kck ωk and of hαk, ωkick = αT k Kck ωk. As step of the scheme is illustrated in Figure 1. The Jacobi field is com- puted with only four calls to the Hamiltonian equations. This operation scales quadratically with the dimension of the manifold, which makes this algorithm practical in high dimension, unlike Christoffel-symbol-based solutions. Step (iv) –solving a linear system of size ncp– is the most expensive one, but remained within reasonable computational time in the investigated examples. Fig. 1: Step of the parallel transport of the vector w (blue arrow) along the geodesic γ. Jw γ is computed by central finite difference with the perturbed geodesics γh and γ−h, integrated with a second-order Runge-Kutta scheme (dot- ted black arrows). A fan of geodesics is formed. In [8], the authors prove the convergence of this scheme, and show that the error increases linearly with the size of the step used. The convergence is guaranteed as long as the step (ii) is performed with a method of order at least two. A first order method in step (iii) is also theoretically sufficient to guarantee convergence. Those variations will be studied in Section 3.3. 3 Application to the prediction of brain structures 3.1 Introducing the exp-parallelization concept Fig. 2: Time-reparametrized exp-parallelization of a reference geodesic model. The black dots are the observations, on which are fitted a geodesic regression (solid black curve, parametrized by the blue arrow) and a matching (leftmost red arrow). The red arrow is then parallel-transported along the geodesic, and exponentiated to define the exp-parallel curve (black dashes). Fig. 3: Illustration of the exp-parallelization concept. Top row: the reference geodesic at successive times. Bottow row: the exp-parallel curve. Blue arrows: the geodesic momenta and velocity field. Red arrows: the momenta describing the initial registration with the target shape and its transport along the geodesic. Exploiting the fanning scheme described in Section 2.4, we can parallel- transport any set of momenta along any given reference geodesic. Figure 2 illustrates the procedure. The target shape is first registered to the reference geodesic : the diffeomorphism that best transforms the chosen reference shape into the target one is estimated with a gradient descent algorithm on the initial control points and momenta [2]. Such a procedure can be applied generically to images or meshes. Once this geodesic is obtained, its initial set of momenta is parallel-transported along the reference geodesic. Taking the Riemannian ex- ponential of the transported vector at each point of the geodesic defines a new trajectory, which we will call exp-parallel to the reference one. As pointed out in [5], the parallel transport is quite intuitive in the context of shape analysis, for it is an isometry which transposes the evolution of a shape into the geometry of another shape, as illustrated by Figure 3. 3.2 Data and experimental protocol Repeated Magnetic Resonance Imaging (MRI) measurements from 71 subjects are extracted from the ADNI database and preprocessed through standard pipeli- nes into affinely co-registered surface meshes of hippocampi, caudates and putam- ina. The geometries of those brain sub-cortical structures are altered along the Alzheimer’s disease course, which all considered subjects finally convert to. Two subjects are successively chosen as references, for they have fully de- veloped the disease within the clinical measurement protocol. As illustrated on Figure 2, a geodesic regression [3] is first performed on each reference subject to model the observed shape progression. The obtained trajectory on the chosen manifold of diffeomorphisms is then exp-parallelized into a shifted curve, which is hoped to model the progression of the target subject. To account for the variability of the disease dynamics, for each subject two scalar coefficients encoding respectively for the disease onset age and the rate of progression are extracted from longitudinal cognitive evaluations as in [11]. The exp-parallel curve is time-reparametrized accordingly, and finally gives predic- tions for the brain structures. In the proposed experiment, the registrations and geodesic regressions typically feature around 3000 control points in R3 , so that the deformation can be seen as an element of a manifold of dimension 9000. 3.3 Estimating the error associated to a single parallel transport To study the error in this high-dimensional setting, we compute the parallel transport for a varying number of discretization steps N, thus obtaining in- creasingly accurate estimations. We then compute the empirical relative errors, taking the most accurate computation as reference. Arbitrary reference and target subjects being chosen, Figure 4 gives the re- sults for the proposed algorithm and three variations : without enforcing the conservations at step (v) [WEC], using a Runge-Kutta of order 4 at step (ii) [RK4], and using a single perturbed geodesic to compute J at step (iii) [SPG]. 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0 0.5 Length of the step h = 1 N Empirical relative error (%) Fig. 4: Empirical rela- tive error of the paral- lel transport in a high- dimensional setting. In black the proposed algorithm, in green the WEC variant, in red the RK4 variant, and in blue the SPG one. We recover a linear behavior with the length of the step 1 N in all cases. The SPG variant converges much slower, and is excluded from the following considerations. For the other algorithms, the empirical relative error remains below 5% with 15 steps or more, and below 1% with 25 steps or more. The slopes of the asymp- totic linear behaviors, estimated with the last 10 experimental measurements, range from 0.10 for the RK4 method to 0.13 for the WEC one. Finally, an iter- ation takes respectively 4.26, 4.24 and 8.64 seconds for the proposed algorithm, the WEC variant and the RK4 one. Therefore the initially detailed algorithm in section 2.4 seems to achieve the best tradeoff between accuracy and speed in the considered experimental setting. 3.4 Prediction performance Table 1 gathers the predictive performance of the proposed exp-parallelization method. The performance metric is the Dice coefficient, which ranges from 0 for disjoint structures to 1 for a perfect match. A Mann-Witney test is performed to quantify the significance of the results in comparison to a naive methodology, which keeps constant the baseline structures over time. Considering the very high dimension of the manifold, failing to accurately capture the disease progression Method Predicted follow-up visit M12 M24 M36 M48 M60 M72 M96 N=140 N=134 N=123 N=113 N=81 N=62 N=17 [exp] .882 .884 .852 .852 .825 .809  ∗ ∗ .796 .764  ∗ ∗ ∗ .768 .734  ∗ ∗ .756 .706  ∗ ∗ ∗ .730 .636  ∗ ∗ [ref] Table 1: Averaged Dice performance measures. In each cell, the first line gives the average performance of the exp-parallelization-based prediction [exp], and the second line the reference one [ref]. Each column corresponds to an increasingly remote predicted visit from baseline. Significance levels [.05, .01, .001]. trend can quickly translates into unnatural predictions, much worse than the naive approach. The proposed paradigm significantly outperforms the naive prediction three years or later from the baseline, thus demonstrating the relevance of the exp- parallelization concept for disease progression modeling, made computationally tractable thanks to the operational qualities of the fanning scheme for high- dimensional applications. 4 Conclusion We detailed the fanning scheme for parallel transport on a high-dimensional manifold of diffeomorphisms, in the shape analysis context. Our analysis unveiled the operational qualities and computational efficiency of the scheme in high dimensions, with a empirical relative error below 1% for 25 steps only. We then took advantage of the parallel transport for accurately predicting the progression of brain structures in a personalized way, from previously acquired knowledge. References 1. Beg, M.F., Miller, M.I., Trouvé, A., Younes, L.: Computing large deformation metric mappings via geodesic flows of diffeomorphisms. Int J Comput Vis (2005) 2. Durrleman, S., Allassonnière, S., Joshi, S.: Sparse adaptive parameterization of variability in image ensembles. Int J Comput Vis 101(1), 161–183 (2013) 3. Fletcher, T.: Geodesic regression and the theory of least squares on Riemannian manifolds. Int J Comput Vis 105(2), 171–185 (2013) 4. Lenglet, C., Rousson, M., Deriche, R., Faugeras, O.: Statistics on the manifold of multivariate normal distributions: theory and application to diffusion tensor MRI processing. Journal of Mathematical Imaging and Vision 25(3), 423–444 (2006) 5. Lorenzi, M., Ayache, N., Pennec, X.: Schild’s ladder for the parallel transport of deformations in time series of images. IPMI pp. 463–474 (2011) 6. Lorenzi, M., Pennec, X.: Geodesics, parallel transport & one-parameter subgroups for diffeomorphic image registration. Int J Comput Vis 105(2), 111–127 (Nov 2013) 7. Lorenzi, M., Pennec, X.: Parallel transport with pole ladder: Application to defor- mations of time series of images. In: Geometric Science of Information. vol. 8085, pp. 68–75 (2013) 8. Louis, M., Charlier, B., Jusselin, P., Pal, S., Durrleman, S.: A fanning scheme for the parallel transport along geodesics on Riemannian manifolds (Jul 2017), 9. Micheli, M.: The differential geometry of landmark shape manifolds: metrics, geodesics, and curvature. Ph.D. thesis, Providence, RI, USA (2008), aAI3335682 10. Miller, M.I., Trouvé, A., Younes, L.: Geodesic shooting for computational anatomy. Journal of Mathematical Imaging and Vision 24(2), 209–228 (2006) 11. Schiratti, J.B., Allassonnière, S., Colliot, O., Durrleman, S.: Learning spatiotem- poral trajectories from manifold-valued longitudinal data. In: NIPS 2015 12. Younes, L.: Jacobi fields in groups of diffeomorphisms and applications. Quarterly of Applied Mathematics 65(1), 113–134 (2007) 13. Zhang, Miaomiao, Fletcher, P.: Probabilistic principal geodesic analysis. Advances in Neural Information Processing Systems 26 pp. 1178–1186 (2013)