Geodesic Analysis in Kendall’s Shape Space with Epidemiological Applications

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. 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.

1 Introduction

In recent years, there has been an increasing 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 aging 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 [10]), 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 [7]. For the derivation of the corresponding Euler–Lagrange equations for some important manifolds, we refer to [21]. An extension to intrinsic Riemannian polynomials has been considered in [9]. Earlier related results in the framework of large deformation diffeomorphic metric mapping (LDDMM) can be found in [25] and [26]. In [3], authors present a kernel-based generalization of geodesic regression, to manifold-valued longitudinal parameters. For an overview of statistical analysis on Riemannian manifolds see [10] and [24].

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 [5, 27], constructions based on parallel transport provide the most natural approach and have shown superior sensitivity and stability in the context of diffeomorphic registration [16]. 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 [16] for fanning [19]. 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. [11] and [7]). 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 Sect. 2, 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 Sect. 3, where we also present our algorithm for the computation of geodesic regression. In Sect. 4, we apply this algorithm to yield longitudinal statistical analysis of femur data from an epidemiological study dealing with osteoarthritis and discuss the numerical results.

2 Geodesic Analysis in Shape Space

A pre-shape is a k-ad of landmarks (i.e., particular points) in ${\mathbb {R}}^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 [14]. For the relevant tools from Riemannian geometry, we refer to [8].

2.1 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\in M(m,k)$, where M(m, k) denotes the space of real $m\times k$ matrices. Denoting the columns of x by $x_i$ and their Euclidean mean by ${\bar{x}}$, in order to remove translations, we replace $x_i$ by $x_i-{\bar{x}}$. The result ${\mathbb {R}}_m^k{:}{=}\{x\in M(m,k):\, \sum _{i=1}^kx_i=0\}$, identified with $M(m,k-1)$, will be endowed with its canonical scalar product given by $\langle x,y \rangle ={{\,\mathrm{trace}\,}}(xy^t)$. Denoting the Frobenius norm by $\Vert \cdot \Vert$, we call the sphere ${\mathcal {S}}_m^k{:}{=}\{x\in {\mathbb {R}}_m^k:\, \Vert x\Vert =1\}$pre-shape space and endow it with the spherical Procrustes metric $d(x,y) {:}{=}\arccos (\langle x,y \rangle )$. Now, the left action of $\mathrm {SO}_m$ on ${\mathcal {S}}_m^k$ given by $(R,x)\mapsto Rx$ defines an equivalence relation given by $x\sim y$ if and only if $y=Rx$ for some $R\in \mathrm {SO}_m$. Kendall’s shape space is defined as $\varSigma _m^k={\mathcal {S}}_m^k/\mathord {\sim }$. Provided that $k\ge m+1$, the dimension of $\varSigma _m^k$ is $m(k-1)-\frac{1}{2}m(m-1)-1$. Now, denoting the canonical projection of $\sim$ by $\pi$, the induced distance between any two shapes $\pi (x)$ and $\pi (y)$ is given by

$$\begin{aligned} d_{\varSigma }(x,y) {:}{=}\min _{R\in \mathrm {SO}_m} d(x,Ry)=\arccos \sum _{i=1}^m\lambda _i \end{aligned}$$

where $\lambda _1\ge \cdots \ge {\lambda _{m-1}\ge } |\lambda _m|$ denote the pseudo-singular values of $yx^t$. Denoting ${\mathcal {D}}_j {:}{=}\{x\in {\mathcal {S}}_m^k:\,{{\,\mathrm{rank}\,}}(x)\le j\}$, it turns out that $\varSigma _{m,m}^k{:}{=}\varSigma _m^k\setminus \pi ({\mathcal {D}}_{m-2})$ inherits a differential structure that is compatible with its quotient topology. Following [14], we refer to $\pi ({\mathcal {D}}_{m-2})$ as the singular part of $\varSigma _m^k$. In particular, $\varSigma _m^k$ is a strata of manifolds with varying dimensions and $\varSigma _{m,m}^k$ is open and dense in $\varSigma _m^k$. Away from the singular part, the quotient map $\pi$ is a Riemannian submersion with respect to the metric induced by the ambient Euclidean space. Moreover, for $k\ge 3$, the shape space $\varSigma _1^k$ (resp. $\varSigma _2^k$) is isometric to the sphere (resp. projective space). We call $x,y\in {\mathcal {S}}_m^k$well positioned, and write $x{\mathop {\sim }\limits ^{\omega }}y$, if and only if $yx^t$ is symmetric and $d(x,y)=d_{\varSigma }(x,y)$. For each $x,y\in {\mathcal {S}}_m^k$, there exists an optimal rotation $R\in \mathrm {SO}_m$ such that $x{\mathop {\sim }\limits ^{\omega }}Ry$. Note that R does not need to be unique. Let ${\mathbb {U}}$ denote a neighborhood in ${\mathcal {S}}_m^k$ with radius smaller then $\pi /4$ (the diameter of $\varSigma _m^k$ is $\pi /2$) such that

$$\begin{aligned} \lambda _{m-1}+\lambda _m> 0\,\text { for all }x,y\in {\mathbb {U}}. \end{aligned}$$

For $x,y\in {\mathbb {U}}$ the optimal rotation R is unique and the function

$$\begin{aligned} {\mathcal {S}}_m^k\ni y\mapsto \omega (x,y) {:}{=}Ry \end{aligned}$$

is well-defined.

We recall that, for a Riemannian submersion $f :M\rightarrow N$ and $y\in N$, $f^{-1}(y)$ is a submanifold of M. For any $x\in M$, denoting the kernel of $\mathrm {d}_x f$ by $\mathrm {Ver}_x$, the tangent space $\mathrm {T}_x M$ to M at x admits an orthogonal decomposition $\mathrm {T}_x M=\mathrm {Hor}_x\oplus \mathrm {Ver}_x$ where $\mathrm {Hor}_x$ and the $\mathrm {Ver}_x$ are the so-called horizontal and vertical subspaces. Due to [14], the vertical space at $x\in {\mathcal {S}}_m^k$ is given by

$$\begin{aligned} \mathrm {Ver}_x=\{Ax:\,A+A^t=0\} , \end{aligned}$$

and the horizontal space is given by

$$\begin{aligned} \mathrm {Hor}_x=\{u\in M(m,k-1):\,ux^t=xu^t\text { and } \langle x,u \rangle =0\}. \end{aligned}$$

We denote the vector space of $m\times m$ skew-symmetric real matrices by $\mathrm {Skew}_m$. Thus $\mathrm {Ver}_x=\mathrm {Skew}_m\cdot 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{\mathop {\sim }\limits ^{\omega }}y$ the geodesic from x to y given by

$$\begin{aligned} \varPhi (t,x,y) {:}{=}\exp _x(t\log _x y)=\frac{\sin ((1-t)\varphi )}{\sin \varphi }x+\frac{\sin (t\varphi )}{\sin \varphi }y\nonumber \\ \end{aligned}$$

(1)

with $\varphi =\arccos (\langle x,y \rangle ),\,0\le t\le 1$, is horizontal. Hence $\varPhi$ realizes the minimizing geodesic from $\pi (x)$ to $\pi (y)$. The following result concerns determination and $\mathrm {SO}_m$-equivariance for horizontal and vertical projection.

Lemma 1

Fix $x\in {\mathcal {S}}_m^k$ and $w\in \mathrm {T}_x{\mathcal {S}}_m^k$. Let $\mathrm {ver}_x$ resp. $\mathrm {hor}_x$ denote the restriction of vertical resp. horizontal projection to $\mathrm {T}_x{\mathcal {S}}_m^k$.

• (a) $\mathrm {ver}_x(w)=Ax$ if and only if A solves the Sylvester equation

  $$\begin{aligned} Axx^t+xx^tA=wx^t-xw^t. \end{aligned}$$

  (2)

  Moreover, the above equation has a unique skew-symmetric solution if ${{\,\mathrm{rank}\,}}(x)\ge m-1$.

• (b) Fix $R\in \mathrm {SO}_m$. Then $\mathrm {ver}_{Rx}(Rw)=R\mathrm {ver}_x(w)$ and $\mathrm {hor}_{Rx}(Rw)=R\mathrm {hor}_x(w)$.

Proof

For (a), let $\mathrm {ver}_x(w)=Ax$, i.e., $w=u+Ax$ with $ux^t$ symmetric and $A\in \mathrm {Skew}_m$. A straightforward computation eliminating $ux^t$ implies that (2) holds. To prove the converse, let $j {:}{=}{{\,\mathrm{rank}\,}}(x)$. Suppose without loss of generality that $j>1$ and write $x=\left( {\begin{array}{c}x_1\\ 0\end{array}}\right)$ with

$$\begin{aligned} {{\,\mathrm{rank}\,}}(x_1)=j,\,w=\left( {\begin{array}{c}w_1\\ w_0\end{array}}\right) , \end{aligned}$$

where $w_1$ is $j\times 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

$$\begin{aligned} A=\begin{pmatrix}A_1 \quad &{} a\\ -a^t &{} \quad A_0\end{pmatrix} \end{aligned}$$

with $A_0\in \mathrm {Skew}_{m-j}$ arbitrary is skew-symmetric and solves the Sylvester equation (2) which also implies that $(w-Ax)x^t$ is symmetric. Hence Ax is the vertical component of w. If ${{\,\mathrm{rank}\,}}(x)=m-1$, then $A_0=0$. If x has full rank, then $A=A_1$.

For (b) note that $\langle Rw,Rx \rangle =\langle w,x \rangle =0$, i.e., $w\in \mathrm {T}_xS$ implies $Rw\in \mathrm {T}_{Rx}S$. Now, $\mathrm {ver}_{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 $\mathrm {ver}_{Rx}(Rw)=R.\mathrm {ver}_x (w)$ and $\mathrm {hor}_{Rx}(Rw)=R\mathrm {hor}_x(w)$. $\square$

Henceforth the superscript v (resp. h) denotes the vertical (resp. horizontal) component, i.e., for any $w\in {\mathbb {R}}_m^k$ we have the orthogonal decomposition $w=\langle w,x \rangle x+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 $\mathrm {SO}_m$-equivariant. Note that this property holds even if $\pi (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\in {\mathcal {S}}_m^k:\,{{\,\mathrm{rank}\,}}(x)\ge m-1\}$ on which $\pi$ 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\in S$ the Sylvester equation (2) has a unique solution determining horizontal and vertical projections and the restriction of $\mathrm {d}_x\pi$ to $\mathrm {Hor}_x$ is an isometry of Euclidean vector spaces $\mathrm {Hor}_x$ and $\mathrm {T}_{\pi (x)}\varSigma _{m,m}^k$. Denoting the covariant derivatives in the pre-shape and shape space by $\nabla$ resp. $\tilde{\nabla }$, for horizontal vector fields X and Y we have

$$\begin{aligned} (\tilde{\nabla }_{\mathrm {d}\pi X} \mathrm {d}\pi Y) \circ \pi = \mathrm {d}\pi (\nabla _X Y). \end{aligned}$$

In the following $[ \cdot , \cdot ]$ denotes the Lie bracket in ${\mathbb {R}}_m^k$, i.e., $[U,V]= \mathrm {D}V(U)- \mathrm {D}U(V)$ ($\mathrm {D}$ Euclidean). For the Euclidean derivative of a vector field W along a curve $\gamma$ in ${\mathbb {R}}_m^k$ we use $\frac{\mathrm {D}}{\mathrm {d}t}$ and also for simplicity of notation a dot, i.e., $\nabla _{{\dot{\gamma }}}W={\dot{W}}-\langle {\dot{W}},\gamma \rangle \gamma$ if $\Vert \gamma \Vert =1$, and $\frac{\mathrm {D}^2W}{\mathrm {d}t^2}={\ddot{W}}$, etc. We set¹

$$\begin{aligned} \mathrm {Log}_xy {:}{=}\log _x\omega (x,y),\,\mathrm {Exp}_xu {:}{=}\exp _x u^h,\,u\in \mathrm {T}_x{\mathcal {S}}_m^k. \end{aligned}$$

For the computation of the Fréchet mean (cf. [11] and [24]) $\pi ({\bar{q}})$ of the shapes $\pi (q_1),\cdots ,\pi (q_N)$ with $q_i\in {{\mathbb {U}}}$, i.e.,

$$\begin{aligned} {\bar{q}}&{:}{=}{{\,\mathrm{arg\,min}\,}}_{x} G(x),\, G(x) {:}{=}\sum _{i=1}^N d_{\varSigma }^2(x,q_i) , \end{aligned}$$

(3)

we apply Newton’s method to Karcher’s equation $\sum _{i=1}^N\mathrm {Log}_xq_i=0$ as follows. We search for the unique zero ${\bar{q}}$ of the function f defined by

$$\begin{aligned} f(x)=\sum _{i=1}^N\mathrm {Log}_xq_i,\,x\in {\mathbb {U}} , \end{aligned}$$

and set

$$\begin{aligned} x_{k+1}=\mathrm {Exp}_{x_k}(-(\mathrm {d}_{x_k}f)^{-1}f(x_k)). \end{aligned}$$

A suitable initial value is the normalized Euclidean mean

$$\begin{aligned} x_0=\frac{1}{\Vert \sum _{i=1}^Nq_i\Vert }\sum _{i=1}^Nq_i. \end{aligned}$$

The total variance of $q=(q_1,\cdots ,q_N)$ reads

$$\begin{aligned} {{\,\mathrm{var}\,}}(q)=\frac{1}{N}G({\bar{q}})=\frac{1}{N}\sum _{i=1}^N\Vert \mathrm {Log}_{{\bar{q}}}q_i\Vert ^2. \end{aligned}$$

2.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.²

We call a vector field W along a horizontal curve $\gamma$horizontally parallel (for brevity h-parallel) if and only if W is horizontal and $\mathrm {d}\pi W$ is parallel along $\pi \circ \gamma$. 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.

Proposition 1

Let $\gamma :[0,\tau ]\rightarrow S$ be a smooth horizontal curve with initial velocity v, u a horizontal vector at $x {:}{=}\gamma (0)$ and W a vector field along $\gamma$ with $W(0)=u$.

• (a) The vector field W is h-parallel transport of u along $\gamma$ if and only if ${\dot{W}}=A\gamma -\langle W,{\dot{\gamma }} \rangle \gamma$ where A is the unique solution of

  $$\begin{aligned} A\gamma \gamma ^t+\gamma \gamma ^tA={\dot{\gamma }}W^t-W{\dot{\gamma }}^t. \end{aligned}$$

  (4)

• (b) Suppose that $\gamma$ is a unit-speed geodesic. Then Eq. (4) reduces to

  $$\begin{aligned} {\dot{A}}\gamma \gamma ^t+\gamma \gamma ^t{\dot{A}}+3(A{\dot{\gamma }}\gamma ^t+\gamma {\dot{\gamma }}^tA)=0. \end{aligned}$$

  (5)

• (c) Let Cv denote the orthogonal projection of u on $\mathrm {Skew}_m\cdot v$, i.e., $Cvv^t+vv^tC=uv^t-vu^t$. Suppose that $C{\dot{\gamma }}$ is horizontal. If $\gamma$ is a unit-speed geodesic, then the h-parallel transport of u is given by

  $$\begin{aligned} W=U+(\langle u,v \rangle +C)({\dot{\gamma }}-v) \end{aligned}$$

  (6)

  where U denotes the Euclidean parallel extension of u along $\gamma$, i.e., $U(t)=u$ f.a. t. If $y = \gamma (\varphi )$ with $\varphi =d (x,y)$, then the h-parallel transport $W_y$ of u along $\gamma$ to y reads

  $$\begin{aligned} W_y = U- 2\frac{\langle u,y \rangle +C\sin (\varphi )}{\Vert x+y\Vert ^2}(x+y) \end{aligned}$$

  (7)

Proof

(a) $d\pi W$ is parallel along $\pi \circ \gamma$ if and only if $\mathrm {d}\pi (\nabla _{{\dot{\gamma }}}W)=0$, i.e., infinitesimal variation of W must be vertical. Hence $\nabla _{{\dot{\gamma }}}W=(\nabla _{{\dot{\gamma }}}W)^v$, which due to Lemma 1 equals $A\gamma$ with $A\gamma \gamma ^t+\gamma \gamma ^t A=(\nabla _{{\dot{\gamma }}}W)\gamma ^t-\gamma (\nabla _{{\dot{\gamma }}}W)^t={\dot{W}}\gamma ^t-\gamma {\dot{W}}^t$. Moreover, $\mathrm {SO}_m$-equivariance of vertical projection implies the well-definedness, i.e., if $\mathrm {d}\pi W$ is parallel, then $\mathrm {d}\pi (Rw)$ is parallel for all $R\in \mathrm {SO}_m$. Note that existence and uniqueness of the solution for (4) 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 ${\dot{f}}=0$ where $f {:}{=}\Vert W\gamma ^t-\gamma W^t\Vert ^2+\langle W,\gamma \rangle ^2$, since $f(0)=0$. If Eq. (4) holds, then

$$\begin{aligned} {\dot{W}}\gamma ^t-\gamma {\dot{W}}^{t}= & {} (\nabla _{{\dot{\gamma }}}W)\gamma ^t-\gamma (\nabla _{{\dot{\gamma }}}W)^t \\= & {} A\gamma \gamma ^t+\gamma \gamma ^tA={\dot{\gamma }}W^t-W{\dot{\gamma }}^t \end{aligned}$$

and

$$\begin{aligned} \langle {\dot{W}},\gamma \rangle +\langle W,{\dot{\gamma }} \rangle= & {} \langle A\gamma -\langle W,\dot{\gamma } \rangle \gamma ,\gamma \rangle +\langle W,\dot{\gamma } \rangle \\ {}= & {} \langle A\gamma ,\gamma \rangle =0. \end{aligned}$$

The last equation follows from the fact that A is skew-symmetric and $\gamma \gamma ^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 ${\dot{f}}$ vanishes. Hence ${\dot{W}}\gamma ^t-\gamma {\dot{W}}^t=\dot{\gamma }W^t-W\dot{\gamma }^t$ and $\langle W,\dot{\gamma } \rangle +\langle {\dot{W}},\gamma \rangle =0$ and the Sylvester equation for the vertical component of ${\dot{W}}$ reads $A\gamma \gamma ^t+\gamma \gamma ^tA=\dot{\gamma }W^t-W\dot{\gamma }^t$. Thus (4) follows.

(b) Note that $W\gamma ^t$ and $\dot{\gamma }\gamma ^t$ are symmetric and $\ddot{\gamma }+\gamma =0$. Thus $W\ddot{\gamma }^t=-W\gamma ^t$ is also symmetric. Now, (4) implies

$$\begin{aligned}&{\dot{A}}\gamma \gamma ^t+\gamma \gamma ^t{\dot{A}}+2(A\dot{\gamma }\gamma ^t+\gamma \dot{\gamma }^tA)\\&\quad =\ddot{\gamma }W^t-W\ddot{\gamma }^t+\dot{\gamma }{\dot{W}}^t-{\dot{W}}\dot{\gamma }^t\\&\quad =\dot{\gamma }(A\gamma -\langle W,\dot{\gamma } \rangle \gamma )^t -(A\gamma -\langle W,\dot{\gamma } \rangle \gamma )\dot{\gamma }^t\\&\quad =-(A\gamma \dot{\gamma }^t+\dot{\gamma }\gamma ^tA). \end{aligned}$$

(c) Obviously W given by (6) satisfies the initial condition $W(0)=u$. Moreover, it satisfies ${\dot{W}}=-C\gamma -\langle W,\dot{\gamma } \rangle \gamma$, i.e., (5) holds with $A(t)=-C$. To prove (7), insert $v=\frac{1}{\varphi }\log _xy=\frac{y-x\cos (\varphi )}{\sin (\varphi )}$ and $\dot{\gamma }=\frac{-1}{\varphi }\log _yx$ into (6).$\square$

Note that, due to skew-symmetry of $\gamma \gamma ^t(\nabla _{\dot{\gamma }}W)\gamma ^t$, the differential equation for the h-parallel transport can also be written as

$$\begin{aligned} (\nabla _{\dot{\gamma }}W)\gamma ^t\gamma \gamma ^t+\gamma \gamma ^t(\nabla _{\dot{\gamma }}W)\gamma ^t=(\dot{\gamma }W^t-W\dot{\gamma }^t)\gamma \gamma ^t. \end{aligned}$$

(8)

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

Remark 1

We mention two cases such that (6) and (7) 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 $\chi _i$ and $\eta _i$ denote the rows of a shape $\chi$ and $\eta$ a horizontal vector at $\chi$. Fix $\mu \in {\mathbb {R}}$ and let $C {:}{=}\mu \begin{pmatrix} 0&{} 1\\ -1&{} 0\end{pmatrix}$. Then $C\eta \chi ^t=\mu \begin{pmatrix} \eta _2\chi _1^t&{} \eta _2\chi _2^t\\ -\eta _1\chi _1^t&{} -\eta _1\chi _2^t\end{pmatrix}$ is symmetric since $\langle \eta ,\chi \rangle =0$. Hence, $C\dot{\gamma }$ is horizontal for arbitrary $C\in \mathrm {Skew}_2$.

2.3 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 $\gamma$ in S is a geodesic if and only if $\pi \circ \gamma$ is a geodesic in $\pi (S)$. Hence, for a horizontal geodesic $\gamma$, any geodesic variation of $\pi \circ \gamma$ in the latter space reads $\pi \circ \varGamma$ with $\varGamma$ a variation of $\gamma$ through horizontal geodesic. Thus the variation field $\frac{\mathrm {d}}{\mathrm {d}s}(\pi \circ \varGamma (s, \cdot ))|_{s=0}=\mathrm {d}\pi ( \frac{\mathrm {d}}{\mathrm {d}s}\varGamma (s, \cdot )|_{s=0})$ is a Jacobi field of the shape space. Recall that a vector field J along $\gamma$ is called normal if and only if $\langle J,\dot{\gamma } \rangle =0$ and the tangential component of any Jacobi field is just given by $(a+bt)\dot{\gamma }(t)$ with $a,b\in {\mathbb {R}}$, which is obviously horizontal. Thus the challenge is to find those normal vector fields that project to a Jacobi field in the shape space.

Theorem 1

Let J be a normal vector field along $\gamma$ and denote

$$\begin{aligned} K = \left( \frac{\mathrm {D}J^v}{\mathrm {d}t}\right) ^v+2\left( \frac{\mathrm {D}J^h}{\mathrm {d}t} \right) ^v. \end{aligned}$$

• (a) $\mathrm {d}\pi (J)$ is a Jacobi field if and only if

  $$\begin{aligned} \left( \frac{\mathrm {D}^2J}{\mathrm {d}t^2} + J\right) ^h&= 2\left( \frac{\mathrm {D}K}{\mathrm {d}t} \right) ^h \end{aligned}$$

  (9)
  $$\begin{aligned} \left( \frac{\mathrm {D}^2J}{\mathrm {d}t^2} + J\right) ^v&= \left( \frac{\mathrm {D}K}{\mathrm {d}t} \right) ^v \end{aligned}$$

  (10)

• (b) A normal Jacobi field $J^S$ of the pre-shape sphere projects to a Jacobi field if and only if $\left( \frac{\mathrm {D}K}{\mathrm {d}t} \right) ^h=0$.

Proof

(a) Obviously, solutions of the equations are, due to $\mathrm {SO}_m$-equivariance of horizontal and vertical projection, invariant under $\mathrm {SO}_m$ action. Let Y and Z be vector fields on ${\mathcal {S}}_m^k$. Following [22], the ${\mathcal {A}}$ and ${\mathcal {T}}$-tensor fields are defined as

$$\begin{aligned} {\mathcal {T}}_YZ&=(\nabla _{Y^v}Z^v)^h+(\nabla _{Y^v}Z^h)^v,\\ {\mathcal {A}}_YZ&=(\nabla _{Y^h}Y^v)^h+(\nabla _{Y^h}Z^h)^v. \end{aligned}$$

Due to [23, Theorem 2], $\mathrm {d}\pi (J)$ is a Jacobi field if and only if

$$\begin{aligned} \left( \frac{\mathrm {D}^2J}{\mathrm {d}t^2} -\mathrm {R}(J,X,X) \right) ^h&=2{\mathcal {A}}_X K,\\ \left( \frac{\mathrm {D}^2J}{\mathrm {d}t^2}-\mathrm {R}(J,X,X) \right) ^v&=\left( \frac{\mathrm {D}K}{\mathrm {d}t}\right) ^v+{\mathcal {T}}_K X, \end{aligned}$$

where $X=\dot{\gamma }$, $\frac{\mathrm {D}}{dt}$ stands for $\nabla _X$ and the vector field K is given by

$$\begin{aligned} K(J)=\left( \frac{\mathrm {D}J^v}{\mathrm {d}t}\right) ^v-{\mathcal {T}}_{J^v}X+2{\mathcal {A}}_X J^h. \end{aligned}$$

Note that as J is normal, its covariant and Euclidean derivative coincide. In our setting, the fibers are totally geodesic, hence ${\mathcal {T}}\equiv 0$. Therefore $K = \left( \frac{\mathrm {D}J^v}{\mathrm {d}t}\right) ^v+2\left( \frac{\mathrm {D}J^h}{\mathrm {d}t} \right) ^v$. Moreover, we may suppose $\Vert X\Vert =1$. Thus $\mathrm {R}(J,X,X)=-J$ and we arrive at (9) and (10).

(b) It follows immediately from $\frac{\mathrm {D}^2J^S}{\mathrm {d}t^2}+J^S=0$.$\square$

In the following, we give a geometric construction for normal Jacobi fields which will be employed for geodesic regression. For $\xi \in T_xS$, we set $ver_{x,v}(\xi ):=Ax$, where A denotes the solution (cf. Lemma 1) of the Sylvester equation

$$\begin{aligned} Axx^t+xx^tA=v\xi ^t-\xi v^t. \end{aligned}$$

In the sequel $x=\gamma (0)$, $v=\dot{\gamma }(0)$, $\xi _1,\xi _2\in Hor_x$, $\langle \xi _1,v\rangle =\langle \xi _2,v\rangle =0$ and $i=1,2$.

Theorem 2

Suppose that $\gamma$ has unit speed v. Let J be the solution of the differential equations (9), (10) with $J^v(0)=Cx$, $J^h(0)=\xi _1$ and $\frac{{\mathrm {D}} J^h }{d t}(0)=\xi _2$. Then $J^v=C\gamma$ and the following hold.

(a) Let $Z_i$ and $Y_1$ denote the parallel extensions of $\xi _i$ resp. $Ax=ver_{x,v}(\xi _1)$ along $\gamma$. Then

$$\begin{aligned} J^h(t)=(\cos (t)Z_1(t)-\sin (t)(Y_1(t)+Z_2(t)))^h. \end{aligned}$$

(b) Suppose that $m=2$. Let $w_i$ denote the orthogonal projection of $\xi _i$ tangent to $\mathrm {Skew}_2\cdot v$ and $u_i$ its orthogonal complement. Then

$$\begin{aligned} J^h&=\cos (t) U_1(t)+\cos (2t) W_1(t)+\sin (t)U_2(t)\\&\qquad +\frac{1}{2}\sin (2t)W_2(t), \end{aligned}$$

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

Proof

We may write $\gamma (t)=\cos (t)x+\sin (t)v$ with $\Vert v\Vert =1$, $\langle x,v \rangle =0$ and $vx^t=xv^t$. The variation

$$\begin{aligned} \varGamma (s,t)=\exp (sC)\gamma (t) \end{aligned}$$

is a variation of $\gamma$ through horizontal geodesics since $\dot{\varGamma }\varGamma ^t$ is symmetric, and defines a vertical Jacobi field $J^v=\frac{d\varGamma }{ds}(0,.)=C\gamma$ 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 ${\Vert }\xi _i\Vert =1$. Furthermore, due to $\langle \xi _i,v\rangle =0$ and $\langle Ax,v\rangle =0$, $Z_i(t)=\xi$ and $Y(t)=Ax$. Now, let V denote the h-parallel transport of v along the geodesic $\alpha$ given by $\alpha (s)=\cos (s) x+\sin (s) \xi _1$. Due to (5), we have $V^\prime (0)=Ax$. Now, consider the variations of $\gamma$ given by

$$\begin{aligned} \varGamma _1(s,t)=\cos (t)\alpha (s) +\sin (t) V(s) \end{aligned}$$

and

$$\begin{aligned} \varGamma _2(s,t)=\cos (t)x +\sin (t)(\cos (s)v+\sin (s)\xi _2). \end{aligned}$$

A straightforward computation shows that $\dot{\varGamma _i}\varGamma _i$ is symmetric, i.e., $\varGamma _i$ is a variation of $\gamma$ through horizontal geodesics. Hence $d\pi J_i$ is a Jacobi field, where $J_i=\frac{d \varGamma _i }{d s}(0,.)$. Therefore, $d\pi (J_1+J_2)$ is a Jacobi field. Moreover, $\frac{d \varGamma _2 }{d s}(0,0)=0$ and $\frac{\mathrm {D}}{d t}\frac{d \varGamma _2 }{d s}(0,0)=\frac{\mathrm {D}}{d s}\frac{d \varGamma _2 }{d t}(0,0)=\xi _2$. Hence the solution with $J^h(0)=0$ and $\frac{{\mathrm {D}} J^h }{d t}(0)=\xi _2$ is given by $J_2^h$, where $J_2(t)=\sin (t)Z_2(t)=\frac{d \varGamma _2 }{d s}(0,t)$. It follows that the solution with $J^h(0)=0$ and $\frac{\mathrm {D}J^h }{d t}(0)=Ax$ is given by $(t\mapsto \sin (t)Y_1(t))^h$. Furthermore, $\frac{d \varGamma _1 }{d s}(0,0)=\xi _1$ and $\frac{\mathrm {D}}{d t}\frac{d \varGamma _2 }{d s}(0,0)=\frac{\mathrm {D}}{d s}\frac{d \varGamma _2 }{d t}(0,0)=Ax$. The fact that the space of horizontal vector fields along $\gamma$ is linear, completes the proof.

(b) Let $Q=\begin{pmatrix}0&{} 1\\ -1 &{} 0\end{pmatrix}$. Then $\xi _i$ enjoys the orthogonal decomposition $\xi _i=u_i+w_i$, where $w_i=B_iv$, $u_iv^t=vu_i^t$ and $B_ivv^t+vv^tB_i=\xi _i v^t-v\xi _i^t$. Moreover, $A=\mu Q$ and $B_i=\lambda _i Q$ for some $\mu , \lambda _i\in {\mathbb {R}}$. A straightforward computation shows $Qxx^t+xx^tQ=Q=Qvv^t+vv^tQ$ (note that $\Vert x\Vert =\Vert v\Vert =1$). Hence $A=-B_1$. Now, the vector $w_i$ is horizontal (cf. Remark 1) and normal. Hence the vector $u_i=\xi _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 \gamma ^t=\gamma 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-\langle B_iv,Q\gamma \rangle Q\gamma =\lambda _i(Qv-\langle Qv,Q\gamma \rangle Q\gamma )=\cos (t) B_i\dot{\gamma }=\cos (t)W_i$. Similarly, 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. $\square$

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. [7] and [12]).

3 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< \cdots < t_N$ and distinct pre-shapes $q_1,\cdots ,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 $\gamma$ from x to y with $v=\dot{\gamma }(0)$, we define the misfit between the data and the geodesic as a sum of squared distances with respect to $d_\varSigma$, i.e.,

$$\begin{aligned} F(\gamma ) {:}{=}\sum _{i=1}^Nd_\varSigma ^2(q_i,\gamma (t_i)). \end{aligned}$$

(11)

We can assume that $t_1=0$ and $t_N=1$. While the authors of [7] and [21] 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

$$\begin{aligned} F(x,y) = \sum _{i=1}^Nd_\varSigma ^2(q_i,\gamma (t_i)) = \sum _{i=1}^Nd_\varSigma ^2(q_i,\varPhi (t_i,x,y)). \end{aligned}$$

The reason is that geodesic computations in terms of the function $\varPhi$ defined in Eq. (1), the so-called slerp (spherical linear interpolation), are more efficient. Model estimation is then formulated as the least squares problem

$$\begin{aligned} (x^*,y^*) = {{\,\mathrm{arg\,min}\,}}_{(x,y)} F(x,y),\, x{\mathop {\sim }\limits ^{\omega }}y. \end{aligned}$$

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 [6] with a Hessian approximation based on finite differences and use $(q_1,\omega (q_1,q_N))$ as initial guess. Having in mind that (cf. [24] and [12])

$$\begin{aligned} \nabla \rho _y(x) = -2\mathrm {Log}_x y=-2(\log _xy)^h \end{aligned}$$

(12)

where $\rho _y(x) {:}{=}d_\varSigma ^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 $\varPhi$. Now, for fixed q and t, let $\nabla _xf$ denote the gradient of f with respect to x where $f(x,y):=\rho _q\circ \varPhi (t,x,y)$. Then for any $u\in \mathrm {Hor}_x$, $d_x\mathrm {Exp}_xtv\cdot u=J(t)$ where $d\pi J$ is the horizontal Jacobi field along $\gamma$ with $J(0)=u$ and $J(1)=0$.

In the following, $\top$ and $\perp$ denote tangent resp. orthogonal components of vectors. Now, let $\alpha$ denote the unit speed horizontal geodesic from y to x, i.e., $\alpha (s)=\cos (s) y+\sin (s)v$ with $v=\log _yx$, $s\in [0,\varphi ]$ and $\varphi =\Vert v\Vert$. Denoting the horizontal component of the parallel extension of u along $\gamma$ by U, and ${\tilde{J}}(s)=\frac{\sin s}{\sin \varphi } (U^\perp )^h$, due to Theorem 2, $d\pi {\tilde{J}}$ is a Jacobi field with ${\tilde{J}}(0)=0$ and $\frac{1}{\sin \varphi }\dot{{\tilde{J}}}(0)=(U^\perp )^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\mapsto \varPhi (t,x,y)$ parametrizes the reverse geodesic by arc length ($\Vert \dot{\varPhi }(0,x,y)\Vert =\varphi$), we arrive at $J(t)={\tilde{J}}((1-t)\varphi )$. Hence

$$\begin{aligned} J(t)=\frac{\sin ((1-t)\varphi )}{\sin \varphi }(U^\perp )^h+(1-t)U^\top . \end{aligned}$$

Let $P_x$ denote the h-parallel transport to x along $\gamma$. In view of (12), $-\frac{1}{2}P_{\gamma (t)}\nabla _xf(x,y)$ is the adjoint of the mapping $u\mapsto J(t)$. As the latter is self-adjoint it follows that

$$\begin{aligned} \nabla _xf(x,y)=-2P_x \left( \frac{\sin ((1-t)\varphi )}{\sin \varphi }(W^\perp )^h+(1-t)W^\top \right) , \end{aligned}$$

where $W=\mathrm {Log}_{\gamma (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 (1)) and arrive at

Fig. 1

Healthy (left) and osteoarthritic (right) distal femur with delineated pathological changes in shape

$$\begin{aligned} \nabla _xF(x,y)&=-2P_x\sum _{i=1}^N \left( \frac{\sin ((1-t_i)\varphi )}{\sin \varphi }(W_i^\perp )^h+(1-t_i)W_i^\top \right) \\ \nabla _yF(x,y)&=-2P_y\sum _{i=1}^N \left( \frac{\sin (t_i\varphi )}{\sin \varphi }(W_i^\perp )^h+t_iW_i^\top \right) , \end{aligned}$$

where $W_i=\mathrm {Log}_{\gamma (t)}q_i$.

Now, our procedure for geodesic regression can be summarized as presented in Algorithm 1.

4 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.

4.1 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 Fig. 1. 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 [13] 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 \times 0.37\hbox { mm}$ matrix, 0.7 mm slice thickness) using a state-of-the-art automatic segmentation approach [2]. 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 postprocess, the quality of segmentations as well as the correspondence of the resulting meshes (8,988 vertices) were ensured.

4.2 Geodesic Modeling of Femoral Trajectories

We apply the geodesic regression approach detailed in Sect. 3 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. [1]). In particular, for the intrinsic treatment of the optimization problem underlying the geodesic regression we use the open-source Matlab toolbox manopt [6]. 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\times {}2.70$ GHz) 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 Fig. 2. 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,\cdots ,X^k$ with $X^j=(x_1^j,\cdots ,x_n^j)$, i.e., we apply tangent PCA to $(x_i^j)_{i=1,\cdots ,n}^{j=1,\cdots ,k}$ at the mean of all baseline shapes in HH. In the remainder, the latter is referred to as reference shape Ref.

Fig. 2

Principal components for femoral shape trajectories of subjects with no (HH), progressing (HD), and severe (DD) osteoarthritis (left) and their qualitatively estimated shape trajectories via geodesic regression (right). Note that points on the left show the observed shapes, while those on the right show the corresponding points on the fitted geodesic

Fig. 3

Group-wise analysis of femoral geodesic trends: Magnitude of differences between the group-average trends for HH versus HD (left column), HH versus DD (middle column), and HD versus DD (right column) after transport to common reference shape. Only significantly different displacements ($p<0.01$) are shown (2.0e−4   4.2e−4)

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 [7], a generalization to manifolds is defined as

$$\begin{aligned} R^2 = 1 - \frac{\text {unexplained variance}}{\text {total variance}} = 1 - \frac{\min _{\gamma } F(\gamma )}{\min _x G(x)}, \end{aligned}$$

with $F(\gamma )$ and G(x) as defined in Eqs. (11) and (3), 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 $\left[ 0,1\right]$ (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 versus DD, and HH versus HD, respectively). These findings indicate that the OA-related shape variability is better captured 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 [7]. 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.

4.3 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 [18] and the references therein). Such a normalization can be realized by adapting parallel transport as presented in Algorithm 2. 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 versus HD and HH versus DD, respectively), there are no differences in the longitudinal changes in-between both OA groups. Figure 3 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. Fig. 1) in comparison with both HH versus HD as well as HH versus 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 [28]. 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 $\ge$ 3) as the initial point for the transport path. In line with previous work [17], 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 versus HD and HD versus DD, respectively.

5 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., [4]), 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 [20]. 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.