Geodesic analysis in Kendall's shape space with epidemiological applications
Esfandiar Nava-Yazdani Hans-Christian Hege T. J. Sullivan Christoph von Tycowicz



Geodesic analysis in Kendall's shape space with epidemiological applications


E. Nava-Yazdani, H.-C. Hege, T. J. Sullivan, and C. von Tycowicz





Geodesic analysis in Kendall's shape space
with epidemiological applications[This is a preprint version of an article to appear in the Journal of Mathematical Imaging and Vision  Springer with the DOI https://doi.org/10.1007/s10851-020-00945-w10.1007/s10851-020-00945-w and differs from the final published version in layout and  typographical detail.]
    
	Esfandiar Nava-YazdaniZIB
	Hans-Christian HegeZIB
	T. J. SullivanZIB,FUB
	Christoph von TycowiczZIB

    December 30, 2023
================================================================================================================================================================================================================================================================================================================================================================



1ZIBZuse Institute Berlin, Takustraße 7, 14195 Berlin, Germany. (navayazdani@zib.de, hege@zib.de, sullivan@sib.de, vontycowicz@zib.de)
2FUBFreie Universität Berlin, Arnimallee 6, 14195 Berlin, Germany. (t.j.sullivan@fu-berlin.de)





  
Abstract:
We analytically determine Jacobi fields and parallel transports and compute geodesic regression in Kendall's shape space.
Using the derived expressions, we can fully leverage the geometry via Riemannian optimization and thereby reduce the computational expense by several orders of magnitude over common, nonlinear constrained approaches.
The methodology is demonstrated by performing a longitudinal statistical analysis of epidemiological shape data.
As an example application we have chosen 3D shapes of knee bones, reconstructed from image data of the Osteoarthritis Initiative (OAI).
Comparing subject groups with incident and developing osteoarthritis versus normal controls, we find clear differences in the temporal development of femur shapes.
This paves the way for early prediction of incident knee osteoarthritis, using geometry data alone.



  
Keywords:
Longitudinal modeling ∙ Shape trajectory ∙ Riemannian metric ∙ Principal geodesic analysis ∙ Geodesic regression ∙ Parallel transport ∙ Jacobi fields




§ INTRODUCTION



In recent years, there has been an increased interest in statistical analysis of geometric shapes.
Such analyses are especially often performed in the field of morphometry, but mostly for static forms.
A frequently-encountered situation, however, is that instead of a set of discrete shapes, series of shapes are given, often together with co-varying parameters.
For example, longitudinal imaging studies track biological shape changes over time within and across individuals to gain insight into dynamical processes such as ageing or disease progression.
Statistical modeling and analysis of shapes is of critical importance for a better understanding of such temporal shape data.




The main challenge is that shape variability is inherently nonlinear and high-dimensional, so that classical statistical approaches are not always appropriate.
One way to address this is linearization.
The quality of the resulting statistical model, however, then depends strongly on the validity of the linearity assumption, i.e. that the observed data points lie to a good approximation in a flat Euclidean subspace.
Since the natural variability in populations often leads to a large spread in shape space and the observed data may lie in highly-curved regions (see <cit.>), linearity often cannot be assumed in practical applications.

In the context of longitudinal studies, an important task is to estimate continuous trajectories from sparse and potentially noisy samples.
For smooth individual biological changes, subject-specific spatiotemporal regression models are adequate.
They also provide a way to describe the data at unobserved times (i.e. shape changes between observation times and — within certain limits — also at future times) and to compare trends across subjects in the presence of unbalanced data (e.g. due to drop-outs).
One approach in use is to approximate the observed temporal shape data by geodesics in shape space and, based on these, to estimate overall trends within groups.
Geodesic models are attractive as they feature a compact representation (similar to the slope and intercept term in linear regression) and therefore allow for computationally efficient inference.

The intrinsic theory of least squares and geodesic regression in shape spaces has been introduced in <cit.>.
For the derivation of the corresponding Euler–Lagrange equations for some important manifolds, we refer to <cit.>.
An extension to intrinsic Riemannian polynomials has been considered in <cit.>.
Earlier related results in the framework of large deformation diffeomorphic metric mapping (LDDMM) can be found in <cit.> and <cit.>. 
In <cit.>, the authors present a kernel-based generalization of geodesic regression to manifold-valued longitudinal parameters.
For an overview of statistical analysis on Riemannian manifolds see <cit.> and <cit.>.

An additional challenge in the analysis of shape trajectories is to distinguish between morphological differences due to (i) temporal shape evolutions of a single individual and (ii) the geometric variability in a population of an object class under study.
To obtain a statistically significant localization of structural changes at the population level (group-wise statistics), the subject-specific trajectories need to be transferred in a standard reference frame.
Among the different techniques proposed for normalizing longitudinal deformations <cit.>, constructions based on parallel transport provide the most natural approach 
and have shown superior sensitivity and stability in the context of diffeomorphic registration <cit.>.
Note also that, for general trajectories, the simple transport of each shape is not suitable because the distances between the shapes are not preserved.
However, if the shapes belong to the same geodesic, this problem does not arise, which is another advantage of geodesic regression.

As parallel transport in curved shape spaces is rarely given in closed form, in general it has to be approximated numerically, e.g. employing Schild's ladder <cit.> for fanning <cit.>.

For shapes in 2D, Kendall's shape space is isomorphic to the projective space, which is a symmetric space, so that the essential geometric quantities are well known (cf. <cit.> and <cit.>).
However, for three and more dimensions, because of less restrictive structure, many questions remain open.
Utilizing closed form expressions of the pre-shape sphere, we reduce parallel transport to the solution of a homogeneous first-order differential equation that allows for highly efficient computations.
Moreover, we reduce the important case of parallel transport along a geodesic path to the solution of a low-dimensional equation that only depends on the dimension of the ambient space and not on the spatial resolution of the discrete representation.















The paper is organized as follows.
In <Ref>, after a short overview of Kendall's shape space, we provide a computationally efficient approach (via the so-called Sylvester equation) for the canonical decomposition of tangent vectors into horizontal and vertical components, which is essential for the geometry and analysis of shapes and trajectories.
Moreover, we determine parallel transport and Jacobi fields, which will be employed for geodesic regression.
Parallel transport is essential for statistical normalization, alignment of trajectories and also computation of Jacobi fields.
The latter describes the variability of trajectories that will be modeled as best-fitting geodesics in <Ref>, where we also present our algorithm for the computation of geodesic regression.
In <Ref> we apply this algorithm to yield longitudinal statistical analysis of femur data from an epidemiological study dealing with osteoarthritis and discuss the numerical results.





§ GEODESIC ANALYSIS IN SHAPE SPACE



A pre-shape is a k-ad of landmarks (i.e. particular points) in ^m after removing translations and similarity transformations.
A shape is a pre-shape with rotations removed.
For a comprehensive introduction to Kendall's shape space and details on the subjects of this section, we refer to <cit.>.
For the relevant tools from Riemannian geometry, we refer to <cit.>.



 §.§ Shape Space


In the following we present a brief overview of Kendall's shape space, provide a computationally efficient method to determine horizontal and vertical 
components of tangent vectors of the pre-shape space, and also prove the corresponding equivariance under rotations.

Let x∈ M(m,k), where M(m,k) denotes the space of real m× k matrices. Denoting the columns of x by x_i and their Euclidean mean by x̅, in order to remove translations, 
we replace x_i by x_i-x̅. The result {x∈ M(m,k): ∑_i=1^kx_i=0}, identified with M(m,k-1), 
will be endowed with its canonical scalar product given by x,y=(xy^t).
Denoting the Frobenius norm by , we call the sphere {x∈: x=1} pre-shape space and endow it with the spherical Procrustes metric d(x,y) arccos ( x,y).
Now, the left action of _m on  given by (R,x)↦ Rx defines an equivalence relation given by x∼ y if and only if y=Rx for some R∈_m.
Kendall's shape space is defined as =/∼.
Provided that k≥ m+1, the dimension of  is m(k-1)-1/2m(m-1)-1.
Now, denoting the canonical projection of ∼ by π, the induced distance between any two shapes π(x) and π (y) is given by

    d_Σ(x,y) min_R∈_m d(x,Ry)=arccos∑_i=1^mλ_i

where λ_1≥⋯≥λ_m-1≥  |λ_m| denote the pseudo-singular values of yx^t.
Denoting _j {x∈: (x)≤ j}, it turns out that ∖π(_m-2) inherits a differential structure that is compatible with its quotient topology.
Following <cit.>, we refer to π(_m-2) as the singular part of .
In particular,  is a strata of manifolds with varying dimensions and  is open and dense in .
Away from the singular part, the quotient map π is a Riemannian submersion with respect to the metric induced by the ambient Euclidean space.
Moreover, for k≥ 3, the shape space Σ_1^k (resp. Σ_2^k) is isometric to the sphere (resp. projective space).
We call x,y∈ well positioned, and write x y, if and only if yx^t is symmetric and d(x,y)=d_Σ(x,y).
For each x,y∈, there exists an optimal rotation R∈_m such that x Ry.
Note that R does not need to be unique.
Let 𝕌 denote a neighborhood in  with radius smaller then π/4 (the diameter of  is π/2) such that

    λ_m-1+λ_m> 0  for all x,y∈𝕌.

For x,y∈𝕌 the optimal rotation R is unique and the function

    ∋ y↦ω(x,y)  Ry

is well-defined.

We recall that, for a Riemannian submersion f  M→ N and y∈ N, f^-1(y) is a submanifold of M. 
For any x∈ M, denoting the kernel of _x f by _x, the tangent space _x M to M at x admits an orthogonal decomposition 
_x M=_x⊕_x where _x and the _x are the so-called horizontal and vertical subspaces. Due to <cit.> 
the vertical space at x∈ is given by

    _x={Ax: A+A^t=0} ,

and the horizontal space is given by

    _x={u∈ M(m,k-1): ux^t=xu^t and x,u=0}.

We denote the vector space of m× m skew-symmetric real matrices by _m.
Thus _x=_m· x.
Furthermore, a smooth curve is called horizontal if and only if its tangent field is horizontal.
Geodesics in the shape space are equivalence classes of horizontal geodesics.
Now, let exp and log denote the exponential and logarithm map of the pre-shape space.
For x y the geodesic from x to y given by

    Φ(t,x,y) exp_x(tlog_x y)=sin((1-t))/sinx+sin(t)/siny

with =arccos(x,y), 0≤ t≤ 1, is horizontal.
Hence Φ realizes the minimizing geodesic from π(x) to π(y).
The following result concerns determination and _m-equivariance for horizontal and vertical projection.

	
	Fix x∈ and w∈_x.
	Let _x resp. _x denote the restriction of vertical resp. horizontal projection to _x.
	
		
		
  * _x(w)=Ax if and only if A solves the Sylvester equation
		
    Axx^t+xx^tA=wx^t-xw^t.

		Moreover, the above equation has a unique skew-symmetric solution if (x)≥ m-1.
		
  * Fix R∈_m.
		Then _Rx(Rw)=R_x(w) and _Rx(Rw)=R_x(w).
	



	For (a), let _x(w)=Ax, i.e., w=u+Ax with ux^t symmetric and A∈_m.
	A straightforward computation eliminating ux^t implies that (<ref>) holds.
	To prove the converse, let j (x).
	Suppose without loss of generality that j>1 and write x=x_10 with
	
    (x_1)=j, w=w_1w_0 ,

	where w_1 is j× k.
	We observe that both equations A_1x_1x_1^t+x_1x_1^tA_1=w_1x_1^t-x_1w_1^t and a^tx_1x_1^t=-w_0x_1^t are uniquely solvable, since x_1x_1^t is invertible.
	Furthermore, the solution of the first equation is skew-symmetric, since its right-hand side is skew-symmetric. It follows that
	
    A=[  A_1    a; -a^t  A_0 ]

	with A_0∈_m-j arbitrary, is skew-symmetric and solves the Sylvester equation (<ref>) which also implies that
	(w-Ax)x^t is symmetric. Hence Ax is the vertical component of w.
	If (x)=m-1, then A_0=0. If x has full rank, then A=A_1.
	
	For (b) note that Rw,Rx=w,x=0, i.e., w∈_xS implies Rw∈_RxS.
	Now, _Rx (Rw)=BRx where B is the solution of BRxx^tR^t+Rxx^tR^tB=R(wx^t-xw^t)R^t.
	Hence B=RAR^t, which implies that _Rx(Rw)=R._x (w) and _Rx(Rw)=R_x(w).



Henceforth the superscript v (resp. h) denotes the vertical (resp. horizontal) component, i.e., for any w∈ we have the orthogonal decomposition w=w,xx+w^h+w^v.
Due to the explicit computation above, (R.w)^v=R.w^v and (R.w)^h=R.w^h, i.e., horizontal and vertical projections are _m-equivariant.
Note that this property holds even if π (x) belongs to the singular part of the shape space.
As appropriate for our applications and for brevity, unless otherwise specified, we restrict our data to the
open and dense set S {x∈: (x)≥ m-1} on which π is a Riemannian submersion, thus the geometry of the shape space is mainly described by its horizontal lift in the pre-shape space.
In particular, for x∈ S the Sylvester equation (<ref>) has a unique solution determining horizontal and vertical projections and the restriction of _xπ to _x 
is an isometry of Euclidean vector spaces _x and _π(x).
Denoting the covariant derivatives in the pre-shape and shape space by ∇ resp. ∇, for horizontal vector fields X and Y we have

    (∇_π Xπ Y) ∘π = π(∇_X Y).

In the following [  ,  ] denotes the Lie bracket in , i.e., [U,V]=  V(U)-  U(V) ( Euclidean).
For the Euclidean derivative of a vector field W along a curve  in  we
use / t and also for simplicity of notation a dot, i.e., ∇_W=Ẇ-Ẇ, if =1, and
^2W/ t^2=Ẅ, etc.


We set[Note that the Riemannian exponential map of the shape space denoted by exp satisfies
π(exp_x u)=exp_π(x)(_xπ(u))=exp_π(x)(_xπ(u^h)).]

    _xy log_xω(x,y), _xu exp_x u^h, u∈_x.

For the computation of the Fréchet mean (cf. <cit.> and <cit.>) π(q̅) of the
shapes π(q_1),⋯,π(q_N) with q_i∈𝕌, i.e.,

    q̅   _x G(x),  G(x)  ∑_i=1^N d_Σ^2(x,q_i) ,

we apply Newton's method to Karcher's equation ∑_i=1^N_xq_i=0 as follows.
We search for the unique zero q̅ of the function f defined by

    f(x)=∑_i=1^N_xq_i, x∈𝕌 ,

and set

    x_k+1=_x_k(-(_x_kf)^-1f(x_k)).

A suitable initial value is the normalized Euclidean mean

    x_0=1/∑_i=1^Nq_i∑_i=1^Nq_i.

The total variance of q=(q_1,⋯,q_N) reads

    (q)=1/NG(q̅)=1/N∑_i=1^N_q̅q_i^2.




 §.§ Parallel Transport

Next, we derive formulas for parallel transport in the shape space and its relation to parallel transport in the pre-shape space.[Essentially, part (a) of <Ref> was recently also obtained by <cit.>.]

We call a vector field W along a horizontal curve  horizontally parallel (for brevity h-parallel) if and only if W is horizontal
and π W is parallel along π∘. In the following, we derive the differential equation for the h-parallelism of W and a corresponding
constructive approach using a Sylvester equation in certain cases.









	
	Let :[0,τ]→ S be a smooth horizontal curve with initial velocity v, u a horizontal vector at x (0) and W a vector field along  with W(0)=u.
	
		
		
  * The vector field W is h-parallel transport of u along  if and only if Ẇ=A-W, where A is the unique solution of
		
    A^t+^tA=W^t-W^t.

		
  * Suppose that  is a unit-speed geodesic.
		Then equation (<ref>) reduces to
		
    Ȧ^t+^tȦ+3(A^t+^tA)=0.

		
  * Let Cv denote the orthogonal projection of u on _m· v, i.e. Cvv^t+vv^tC=uv^t-vu^t.
		Suppose that C is horizontal. If  is a unit-speed geodesic, then the h-parallel transport of u is given by
		
    W=U+(u,v+C)(-v)

		where U denotes the Euclidean parallel extension of u along , i.e., U(t)=u f.a. t. If y =  (φ) with φ=d (x,y), then
		the h-parallel transport W_y of u along  to y reads
		
    W_y = U- 2u,y+Csin(φ)/x+y^2(x+y)

	








	(a) π W is parallel along π∘ if and only if π(∇_W)=0, i.e., infinitesimal variation of W must be vertical.
	Hence ∇_W=(∇_W)^v, which due to <Ref> equals
	A with A^t+^t A=(∇_W)^t-(∇_W)^t=Ẇ^t-Ẇ^t.
	Moreover, _m-equivariance of vertical projection implies the well-definedness, i.e., if π W is parallel, then π(Rw)
	is parallel for all R∈_m.
	Note that existence and uniqueness of the solution for (<ref>) with W(0)=u is immediate from the existence and uniqueness of parallel transport and vertical projection.
	Now, W is horizontal if and only if ḟ=0 where f W^t- W^t^2+W,^2, since f(0)=0.
	If equation (<ref>) holds,
	then
	
    Ẇ^t-Ẇ^t = (∇_W)^t-(∇_W)^t = A^t+^tA=W^t-W^t

	and
	
    Ẇ,+W, = A-W,,+W, = A,=0.

	The last equation follows from the fact that A is skew-symmetric and ^t is symmetric, hence their product is trace-free. Now, we arrive at
	f=0, i.e., W remains horizontal. To prove the converse, note that if W is horizontal, then f and therefore ḟ vanishes. Hence
	Ẇ^t-Ẇ^t=W^t-W^t and W,+Ẇ,=0 and the Sylvester equation for the
	vertical component of Ẇ reads A^t+^tA=W^t-W^t.
	Thus (<ref>) follows.
	
	(b) Note that W^t and ^t are symmetric and +=0. Thus W^t=-W^t is also symmetric.
	Now, (<ref>) implies
	
    Ȧ^t+^tȦ+2(A^t+^tA)
    		   =W^t-W^t+Ẇ^t-Ẇ^t
       =(A-W,)^t -(A-W,)^t
       =-(A^t+^tA).

	
	(c) Obviously W given by (<ref>) satisfies the initial condition W(0)=u. Moreover, it satisfies Ẇ=-C-W,, i.e.,
	(<ref>) holds with A(t)=-C. To prove (<ref>), insert
	v=1/φlog_xy=y-xcos (φ )/sin(φ) and =-1/φlog_yx into (<ref>).


Note that, due to skew-symmetry of ^t(∇_W)^t, the differential equation for the h-parallel transport can also be written as

    (∇_W)^t^t+^t(∇_W)^t=(γ̇W^t-Wγ̇^t)^t.

Hence, a vector field along a curve in π(S) is parallel if and only if it has a horizontal lift satisfying the above equation.



	We mention two cases such that (<ref>) and (<ref>) apply. First, W coincides with the spherical parallel transport of u if and only if uv^t=vu^t or, equivalently, C=0. Secondly, for planar shapes. To see this, let χ_i and η_i denote the rows of a shape χ and η a horizontal vector
	at χ. Fix μ∈ and let C μ[  0  1; -1  0 ].
	Then Cηχ^t=μ[  η_2χ_1^t  η_2χ_2^t; -η_1χ_1^t -η_1χ_2^t ] is symmetric since η,χ=0.
	Hence, C is horizontal for arbitrary C∈_2.





 §.§ Jacobi Fields

Next, we derive the differential equation for Jacobi fields and
provide a constructive approach to its solution utilizing parallel
transport.

We recall that a smooth horizontal curve  in S is a
geodesic if and only if π∘ is a geodesic in π(S).
Hence, for a horizontal geodesic , any geodesic variation of
π∘ in the latter space reads π∘ with  a
variation of  through horizontal geodesic.
Thus the variation field / s(π∘(s, 
))|_s=0=π ( / s(s,  )|_s=0) is a Jacobi
field of the shape space.
Recall that a vector field J along  is called normal if and
only if J,=0 and the tangential component of any Jacobi
field is
just given by (a+bt)(t) with a,b∈, which is obviously
horizontal.
Thus the challenge is to find those normal vector fields that project to a
Jacobi field in the shape space.


	Let J be a normal vector field along  and denote
	
    K = ( J^v/ t)^v+2( 
    J^h/ t)^v.

	
		
		
  * π (J) is a Jacobi field if and only if
		
    ( ^2J/ t^2 + J)^h    =
    2( K/ t)^h
    ( ^2J/ t^2 + J)^v    = ( K/ t)^v

		
  * A normal Jacobi field J^S of the pre-shape sphere projects to a
Jacobi field if and only if ( K/ t)^h=0.
	



	(a) Obviously, solutions of the equations are, due to
_m-equivariance of horizontal and vertical projection, invariant
under _m action.
	
    Let Y and Z be vector fields on . Following <cit.>,
the  and 𝒯-tensor fields are defined as
	
    _YZ    =(∇_Y^vZ^v)^h+(∇_Y^vZ^h)^v,
    _YZ    =(∇_Y^hY^v)^h+(∇_Y^hZ^h)^v.

	Due to <cit.>, π(J) is a Jacobi field if and
only if
	
    (^2J/ t^2  - (J,X,X) )^h    =2_X K,
    (^2J/ t^2- (J,X,X) )^v    =(
    K/ t)^v+_K X,

	where X=, /dt stands for ∇_X and the
vector field K is given by
	
    K(J)=( J^v/ t)^v-_J^vX+2_X J^h.

	Note that as J is normal, its covariant and Euclidean derivative
coincide. In our setting, the fibers are totally geodesic, hence
≡ 0. Therefore
    K = ( J^v/ t)^v+2(  J^h/
t)^v. Moreover, we may suppose X=1. Thus
(J,X,X)=-J and we arrive at (<ref>) and (<ref>).

	(b) It follows immediately from ^2J^S/ t^2+J^S=0.



In the following, we give a geometric construction for normal Jacobi
fields which will be employed for geodesic regression.  
For ξ∈ T_xS, we set ver_x,v(ξ):=Ax, where A denotes the solution (cf.
<Ref>) of the Sylvester equation 
    Axx^t+xx^tA=vξ^t-ξ v^t.


In the sequel x= (0), v=(0), ξ_1,ξ_2∈ Hor_x, ⟨ξ_1,v⟩ =⟨ξ_2,v⟩ =0 and i=1,2.

	
	Suppose that  has unit speed v. Let J be the solution of the
	differential equations (<ref>), (<ref>) with J^v(0)=Cx, 
	J^h(0)=ξ_1 and J^h /d t(0)=ξ_2. Then J^v=C and the
	following hold.
	
		
		
  * Let Z_i and Y_1 denote the parallel extensions of ξ_i resp. Ax=ver_x,v(ξ_1) along . Then 
    J^h(t)=(cos (t)Z_1(t)-sin (t)(Y_1(t)+Z_2(t)))^h.

		
  * Suppose that m=2. Let w_i denote the orthogonal projection of ξ_i tangent to _2· v and u_i its orthogonal complement. Then   
		
    J^h=cos (t) U_1(t)+cos (2t) W_1(t)+sin(t)U_2(t)+1/2sin(2t)W_2(t),

		where W_i and U_i are parallel extensions of w_i and u_i.
	



We may write (t)=cos(t)x+sin(t)v with v=1, x,v=0 and
vx^t=xv^t. The variation

    (s,t)=exp(sC) (t)

is a variation of  through horizontal geodesics since ^t is symmetric, and defines a vertical Jacobi field J^v=d/ds(0,.)=C with J(0)=Cx.

(a) As any Jacobi field is a linear combination of parallel vector fields, and those vector fields conserve orthogonality
and length, we may assume ξ_i=1. Furthermore, due to ⟨ξ_i,v⟩=0 and ⟨ Ax,v⟩ =0, Z_i(t)=ξ and Y(t)=Ax. Now, let V denote the h-parallel transport of v along the geodesic α given by α(s)=cos (s) x+sin (s) ξ_1. Due to (<ref>), we have V^' (0)=Ax. Now, consider the variations of  given by

    _1(s,t)=cos (t)α (s) +sin (t) V(s)

and

    _2(s,t)=cos (t)x +sin (t)(cos (s)v+sin (s)ξ_2).

A straightforward computation shows that _̇i̇_i is symmetric, i.e., _i is a variation of  through horizontal
geodesics. Hence dπ J_i is a Jacobi field, where J_i=d _i /d s(0,.). Therefore, dπ (J_1+J_2) is a Jacobi field. Moreover, d _2 /d s(0,0)=0 and /d td _2 /d s(0,0)=/d sd _2 /d t(0,0)=ξ_2. 
Hence the solution with J^h(0)=0 and J^h /d t(0)=ξ_2 is given by J_2^h, where J_2(t)=sin (t)Z_2(t)=d _2 /d s(0,t). It follows that the solution with J^h(0)=0 and J^h /d t(0)=Ax is given by (t↦sin (t)Y_1(t))^h. 
Furthermore, d _1 /d s(0,0)=ξ_1 and /d td _2 /d s(0,0)=/d sd _2 /d t(0,0)=Ax. The fact that the space of horizontal vector fields along  is linear, completes the proof.

(b) Let Q=[  0  1; -1  0 ]. Then ξ_i enjoys the orthogonal decomposition ξ_i=u_i+w_i, where w_i=B_iv, u_iv^t=vu_i^t and B_ivv^t+vv^tB_i=ξ_i v^t-vξ_i^t. Moreover, A=μ Q and B_i=λ_i Q for some μ, λ_i∈.
A straightforward computation shows Qxx^t+xx^tQ=Q=Qvv^t+vv^tQ (note that x=v=1). Hence A=-B_1. Now, the vector w_i is horizontal (cf. Remark 1) and normal. 
Hence the vector u_i=ξ_i-w_i is also horizontal and normal. Therefore its parallel extension is given by U_i(t)=u_i. Using u_iv^t=vu_i^t, we arrive at U_i ^t= U_i^t, i.e. U_i^h=U_i. Moreover, utilizing the fact that Q^2 is minus indentity, we have W_i^h=W_i-W_i^v=B_iv-⟨ B_iv,Q⟩ Q =λ_i(Qv-⟨ Qv,Q⟩ Q)=cos (t) B_i=cos (t)W_i. Similarily, for the constant vector field B_1x, we have (B_1x)^h=-sin (t)W_1. Implying in the expression of J^h from part a), we arrive at the desired formula.




We recall that for m=2, the shape space is isometric to the complex projective space endowed with 
its standard (Fubini–Study) metric. The given formula for J^h in this case is well-known (cf. <cit.> and <cit.>).


§ GEODESIC REGRESSION



In the following, we employ the results of the previous section to derive an efficient and robust approach for finding the relation between an independent scalar variable, i.e. time, and a dependent shape-valued random variable.

Regression analysis is a fundamental tool for the spatiotemporal modeling of longitudinal observations.
Given scalars t_1 < t_2 < ⋯ < t_N and distinct pre-shapes q_1,⋯,q_N, the goal of geodesic regression is to find a geodesic curve in shape space that best fits the data in a least-squares sense.
In particular for a horizontal geodesic  from x to y with v=(0), we define the misfit between the data and the geodesic as a sum of squared distances with respect to d_Σ, i.e.

    F() ∑_i=1^Nd_Σ^2(q_i,(t_i)).

We can assume that t_1=0 and t_N=1.
While the authors of <cit.> and <cit.> identify geodesics by their initial point and velocity — and hence they consider F(x,v) — we use for the identification their endpoints, i.e., we consider

    F(x,y) = ∑_i=1^Nd_Σ^2(q_i,(t_i)) = ∑_i=1^Nd_Σ^2(q_i,Φ(t_i,x,y)).

The reason is that geodesic computations in terms of the function Φ defined in equation (<ref>), the so-called slerp (spherical linear interpolation), are more efficient. 
Model estimation is then formulated as the least-squares problem

    (x^∗,y^∗) = _(x,y) F(x,y),  x y.


In the absence of an analytic solution, the regression problem has to be solved numerically.
To this end, we employ a Riemannian trust-regions solver <cit.> with a Hessian approximation based on finite differences and use
(q_1,ω(q_1,q_N)) as initial guess.
Having in mind that (cf. <cit.> and <cit.>)

    ∇ρ_y(x) = -2_x y=-2(log_xy)^h

where ρ_y(x)  d_Σ^2(x,y), the gradient of the cost function F can be computed using Jacobi fields,
since they express the derivatives of the exponential map and therefore those of Φ. 
Now, for fixed q and t, let ∇_xf denote the gradient of f with respect to x where f(x,y):=ρ_q∘Φ(t,x,y). 
Then for any u∈_x, d_x_xtv· u=J(t) where dπ J is the horizontal Jacobi field along  with J(0)=u and J(1)=0. 

In the following, ⊤ and ⊥ denote tangent resp. orthogonal components of vectors. Now, let α denote the unit speed horizontal geodesic from y to x, i.e., α (s)=cos (s) y+sin (s)v with v=log_yx, s∈ [0,φ] and φ=v. Denoting the horizontal component of the parallel extension of u along  by U, and J̃(s)=sin s/sinφ (U^⊥)^h, due to Theorem <ref>, dπJ̃ is a Jacobi field with J̃(0)=0 and 1/sinφJ̇̃̇(0)=(U^⊥)^h. Reparametrization only changes the tangent component of the Jacobi field. Moreover, horizontal projection does not depend on the parametrization. Due to the fact that t↦Φ(t,x,y) parametrizes the reverse geodesic by arc length (Φ̇(0,x,y)=φ), we arrive at J(t)=J̃((1-t)φ). Hence

    J(t)=sin((1-t)φ)/sinφ(U^⊥)^h+(1-t)U^⊤.

Let P_x denote the h-parallel transport to x along . 
In view of (<ref>), -1/2P_(t)∇_xf(x,y) is the adjoint of the mapping u↦ J(t). 
As the latter is self-adjoint it follows that

    ∇_xf(x,y)=-2P_x ( sin((1-t)φ)/sinφ(W^⊥)^h+(1-t)W^⊤) ,

where W=_ (t)q.
To get the gradient of f with respect to y, we simply replace 1-t by t (another advantage of employing the parametrization (<ref>)) and arrive at

    ∇_xF(x,y)   =-2P_x∑_i=1^N ( sin((1-t_i)φ)/sinφ(W_i^⊥)^h+(1-t_i)W_i^⊤)
    ∇_yF(x,y)   =-2P_y∑_i=1^N ( sin(t_iφ)/sinφ(W_i^⊥)^h+t_iW_i^⊤) ,

where W_i=_ (t)q_i. 

Now, our procedure for geodesic regression can be summarized as presented in Algorithm <ref>.





§ APPLICATION TO EPIDEMIOLOGICAL DATA



In this section, we analyze the morphological variability in longitudinal data of human distal femora in order to quantify shape changes that are associated with femoral osteoarthritis.



 §.§ Data Description


We apply the derived scheme to the analysis of group differences in longitudinal femur shapes of subjects with incident and developing osteoarthritis (OA) versus normal controls.
An overview of OA-related dysmorphisms is shown in <Ref>.
The dataset is derived from the Osteoarthritis Initiative (OAI), which is a longitudinal study of knee osteoarthritis maintaining (among others) clinical evaluation data and radiological images from 4,796 men and women of age 45–79.
The data are available for public access at <http://www.oai.ucsf.edu/>.

From the OAI database, we determined three groups of shapes trajectories: HH (healthy, i.e. no OA), HD (healthy to diseased, i.e. onset and progression to severe OA), and DD (diseased, i.e. OA at baseline) according to the Kellgren–Lawrence 
score <cit.> of grade 0 for all visits, an increase of at least 3 grades over the course of the study, and grade 3 or 4 for all visits, respectively.
We extracted surfaces of the distal femora from the respective 3D weDESS MR images (0.37×0.37 mm matrix, 0.7 mm slice thickness) using a state-of-the-art automatic segmentation approach <cit.>.
For each group, we collected 22 trajectories (all available data for group DD minus a record that exhibited inconsistencies, and the same number for groups HD and HH, randomly selected), each of which comprises shapes of all acquired MR images, 
i.e. at baseline, the 12-, 24-, 36-, 48- and 72-month visits.
In a supervised post-process, the quality of segmentations as well as the correspondence of the resulting meshes (8,988 vertices) were ensured.



 §.§ Geodesic Modeling of Femoral Trajectories


We apply the geodesic regression approach detailed in <Ref> to the femoral shape trajectories described above and represented in Kendall's shape space.
Due to the expressions derived for the parallel transport and Jacobi fields, we can fully leverage the geometry using 
Riemannian optimization procedures (cf. <cit.>).
In particular, for the intrinsic treatment of the optimization problem underlying the geodesic regression we use the open-source Matlab toolbox manopt <cit.>.
In our experiments, we observed a superlinear convergence of the intrinsic trust-region solver for most of the shape trajectories.
Solving the high-dimensional (54k degrees of freedom) regression problem on a laptop computer with Intel Core i7-7500U (2×2.70GHz) CPU took about 0.3s on average.
In contrast, the generic Matlab routine for nonlinear regression (viz. fitnlm) required about 25s to determine a solution, thus being two orders of magnitude slower.

The resulting estimated geodesics along with the original trajectories are visualized in <Ref>.

The geodesic representation provides a less cluttered visualization of the trajectory population
making it easier to identify trends within as well as across groups.
For 2-dimensional visualization we perform dimension reduction for the trajectories X^1,⋯,X^k with X^j=(x_1^j,⋯,x_n^j), i.e. we apply tangent PCA to (x_i^j)_i=1,⋯,n^j=1,⋯,k at the mean of all baseline shapes in HH. In the remainder, the latter is referred to as reference shape Ref.




Next we would like to answer the question of how well the observed data is replicated by the estimated geodesic trends.
A common approach to test this is to compute the coefficient of determination, denoted as R^2, that is the proportion of the total variance in the data explained by the model.
Following <cit.>, a generalization to manifolds is defined as

    R^2 = 1 - unexplained variance/total variance = 1 - min_ F()/min_x G(x),

with F() and G(x) as defined in equations (<ref>) and (<ref>), respectively.
As the unexplained variance cannot exceed the total variance (since the Fréchet mean lies in the search space of the regression problem) and both variances are nonnegative, R^2 must 
lie in the interval [0,1] (with larger values indicating a higher proportion of the variance being explained by the model).

The coefficients of determination were computed for all estimated trends amounting to
group-wise medians (95% confidence intervals) of 0.40 (0.33–0.46), 0.55 (0.48–0.63), and 0.51 (0.40–0.72) for group HH, DD, and HD, respectively.
While for all groups the geodesic model is able to describe a relatively large portion of the shape variability, there is a clear difference between the control group HH and the groups DD and HD associated to osteoarthritis.
In particular, pairwise Mann–Whitney U tests confirm that the differences are highly unlikely due to random chance (with p-values of <10^-3, and 0.005 for HH vs. DD, and HH vs. HD, respectively).
These findings indicate that the OA-related shape variability is better caputured by a single variable (time) than the physiological trends in HH.

Based on the coefficient of determination we also test for the significance of the estimated trends employing permutation tests as suggested in <cit.>.
For each of the trajectories we performed 1,000 permutations and considered the results as statistically significant for p-values less than 0.01.
In almost all cases (63 out of 66) the trends were significant, such that we can expect them to be highly unlikely due to random chance.



 §.§ Group-wise Analysis of Longitudinal Trends







In order to perform group-wise analysis of longitudinal shape changes we compare the estimated geodesic trends of the femoral trajectories.
This requires the consistent integration of intra- and inter-subject variability in order to obtain statistically significant localization of changes at the population level.
In fact, the comparison of longitudinal shape changes is usually performed after normalizing (i.e. transporting) them into a common system of coordinates (see <cit.> and the references therein).
Such a normalization can be realized by adapting parallel transport as presented in Algorithm <ref>.
In particular, for geodesic trends this scheme reduces to parallel transport of the subject-specific velocity along the baseline-to-reference shape geodesic.
The group-wise longitudinal progression was modeled as the mean of the transported velocities.
The areas of significant differences between longitudinal changes were investigated by two-sample Hotelling's T^2 tests on the vertex-wise displacements corresponding to the transported velocity-fields.
While the displacements differ significantly (p<0.01,  after Benjamini-Hochberg false discovery correction) between normal controls and the OA groups (for 55% and 19% of the vertices for HH vs. HD and HH vs. DD, respectively), there are no differences in the longitudinal changes in-between both OA groups.
<Ref> shows a qualitative visualization of the group tests in terms of the magnitude of the difference between the group-wise means.
Visible are changes along the ridge of the cartilage plate (characteristic regions for osteophytic growth, cf. <Ref>) in comparison of both HH vs. HD as well as HH vs. DD, albeit the latter are less pronounced suggesting a saturation of morphological developments.
Additionally, the changes are more developed on the medial compartment, which is in line with previous findings <cit.>.


While velocities are constant for subject-specific geodesics, their parallel transport depends on the path (an effect called holonomy).
To investigate this path dependency, we repeated the above experiment using different paths for the HD group.
In particular, we chose the shape at the onset time (transition time to severe OA, viz. Kellgren–Lawrence score ≥ 3) as the initial point for the transport path.
In line with previous work <cit.>, we found that the results are not sensitive to the path.
More precisely, the results of the group tests agreed for 99.70% and 100.00% of the vertices for HH vs. HD and HD vs. DD, respectively.





§ CONCLUDING REMARKS


This work presented characterizations of and computationally efficient methods for the determination of parallel transport, Jacobi fields and geodesic regression of data represented as shapes in Kendall's space.
Furthermore, an application to longitudinal statistical analysis of epidemiological data (femur data for analysis of knee osteoarthritis) has been shown.
An advantage of modeling trajectories by geodesics is the following: A main task
in longitudinal analysis is to translate trajectories to start at a reference shape.
The intermediate distances between the shapes of a geodesic are preserved by parallel transport, which is not the case for general transports.
Moreover, data inconsistencies are minimized by considering the best-fitting geodesics, and Jacobi fields can be employed to analyze the variability of the geodesics, hence providing a canonical descriptor of trends and differences for the trajectories.

There are many potential avenues for future work.
First, we would like to use the presented methodology within the mixed-effect framework (see e.g. <cit.>), which provides a joined analysis of longitudinal and cross-sectional variability.
In particular, group-wise means of the geodesics can be computed with respect to a natural metric in the tangent bundle (e.g. the Sasaki metric) to determine the group parameters as described in <cit.>.
Second, an extension of the method to higher-dimensional longitudinal parameters instead of just time can be examined, to achieve even more differentiated results.
Third, spline regression poses a natural generalization providing more degrees of freedom.

On the application side, based on the results found, it can be said in summary that the shape trajectories
of the healthy subjects expose significantly different temporal changes than those found in groups with incident and developing OA.
Our analysis delivered detailed insights into the complex morphological changes that fit medical knowledge.
It seems possible to make a correct assignment to one of the three groups based on just two measurements.
The aim of further investigations must be to substantiate this statement, by determining with what reliability a prediction can be made about the onset of knee osteoarthritis depending on
the baseline shape and trend as well as the sensitivity of the latter with respect to the number of observations made and the time intervals between them.



§ ACKNOWLEDGEMENTS

tocsectionAcknowledgements
This research was carried out in the framework of Matheon supported by the Einstein Foundation Berlin (E. Nava-Yazdani),
DFG project In-vivo soft tissue kinematics from dynamic MRI (C. von Tycowicz), and the Freie Universität Berlin through the Excellence Initiative of the Deutsche Forschungsgemeinschaft (T. J. Sullivan).
Furthermore we are grateful for the open-access dataset provided by the OAI[The Osteoarthritis Initiative is a public-private partnership comprised of five contracts
(N01-AR-2-2258; N01-AR-2-2259; N01-AR-2-2260; N01-AR-2-2261; N01-AR-2-2262)
funded by the National Institutes of Health, a branch of the Department of Health and
Human Services, and conducted by the OAI Study Investigators.
Private funding partners
include Merck Research Laboratories; Novartis Pharmaceuticals Corporation,
GlaxoSmithKline; and Pfizer, Inc.
Private sector funding for the OAI is managed by the
Foundation for the National Institutes of Health.
This manuscript was prepared using an OAI
public use data set and does not necessarily reflect the opinions or views of the OAI investigators, the NIH, or the private funding partners.]
as well as for the open-source software Manopt <cit.>.

abbrvnat
   
tocsectionReferences

