Geodesic analysis in Kendall's shape space with epidemiological applications

TL;DR

This paper develops analytical tools for geodesic analysis in Kendall's shape space, enabling efficient Riemannian optimization and applying it to longitudinal epidemiological shape data for early osteoarthritis prediction.

Contribution

It introduces analytical expressions for Jacobi fields and parallel transports in Kendall's shape space, significantly reducing computational costs for geodesic regression.

Findings

Identified shape differences between osteoarthritis and control groups

Demonstrated efficient geodesic regression in shape analysis

Enabled early prediction of knee osteoarthritis using shape data

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.

Equations

$$A=\left(\begin{array}{cc}{A}_1& a\\ -a^t& A_0\end{array}\right)$$

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

$${\left(\frac{{\mathrm{D}}^{2}J}{\mathrm{d}{t}^2}+J\right)}^h$$

$${\left(\frac{{\mathrm{D}}^{2}J}{\mathrm{d}{t}^2}+J\right)}^v$$

$${\mathcal{T}}_{Y}Z$$

$${\mathcal{A}}_{Y}Z$$

$${\left(\frac{{\mathrm{D}}^{2}J}{\mathrm{d}{t}^2}-\mathrm{R}(J,X,X)\right)}^h$$

$${\left(\frac{{\mathrm{D}}^{2}J}{\mathrm{d}{t}^2}-\mathrm{R}(J,X,X)\right)}^v$$

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

$$Ax{x}^t+x{x}^{t}A=v{\xi }^t-\xi v^t.$$

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

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

$$\Gamma (s,t)=\exp (sC)\gamma (t)$$

$${\Gamma }_1(s,t)=\cos (t)\alpha (s)+\sin (t)V(s)$$

$${\Gamma }_2(s,t)=\cos (t)x+\sin (t)(\cos (s)v+\sin (s){\xi }_2).$$

$$\nabla {\rho }_y(x)=-2{\mathrm{Log}}_{x}y=-2({\log }_{x}y{)}^h$$

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

$${\nabla }_{x}f(x,y)=-2P_x\left(\frac{\sin ((1-t)\phi )}{\sin \phi }(W^{\perp }{)}^h+(1-t)W^{\top }\right),$$

$${\nabla }_{x}F(x,y)$$

$${\nabla }_{y}F(x,y)$$

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

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

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

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

$${\mathrm{Hor}}_x=\{u\in M(m,k-1):u{x}^t=x{u}^t\text{ and }⟨x,u⟩=0\}.$$

$$\Phi (t,x,y):={\exp }_x(t{\log }_{x}y)=\frac{\sin ((1-t)\phi )}{\sin \phi }x+\frac{\sin (t\phi )}{\sin \phi }y$$

$$Ax{x}^t+x{x}^{t}A=w{x}^t-x{w}^t.$$

$$\mathrm{rank}(x_1)=j,w=\left(\genfrac{}{}{0ex}{}{w_1}{w_0}\right),$$

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

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

$$\bar{q}$$

$$f(x)=\sum _{i=1}^N{\mathrm{Log}}_x{q}_i,x\in \mathbb{U},$$

Full text

[TABLE]

1. (a)

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

[TABLE] 2. (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 O’Neill (1966), the $\mathcal{A}$ and $\mathcal{T}$-tensor fields are defined as

[TABLE]

Due to (O’Neill, 1967, Theorem 2), $\mathrm{d}\pi(J)$ is a Jacobi field if and only if

[TABLE]

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

[TABLE]

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 $\|X\|=1$. Thus $\mathrm{R}(J,X,X)=-J$ and we arrive at (9) and (10). (b) It follows immediately from $\frac{\mathrm{D}^{2}J^{S}}{\mathrm{d}t^{2}}+J^{S}=0$. $\blacksquare$ In the following, we give a geometric construction for normal Jacobi fields which will be employed for geodesic regression. For $\xi\in T_{x}S$, we set $ver_{x,v}(\xi):=Ax$, where $A$ denotes the solution (cf. Lemma 2.1) of the Sylvester equation

$Axx^{t}+xx^{t}A=v\xi^{t}-\xi v^{t}.$

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

1. (a)

Let $Z_{i}$ and $Y_{1}$ denote the parallel extensions of $\xi_{i}$ resp. $Ax=ver_{x,v}(\xi_{1})$ along $\gamma$. Then

[TABLE] 2. (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

[TABLE]

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 $\|v\|=1$, $\langle x,v\rangle=0$ and $vx^{t}=xv^{t}$. The variation

[TABLE]

is a variation of $\gamma$ through horizontal geodesics since $\dot{\Gamma}\Gamma^{t}$ is symmetric, and defines a vertical Jacobi field $J^{v}=\frac{d\Gamma}{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 {\leavevmode\color[rgb]{0,0,0}\|}\xi_{i}\|=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

[TABLE]

and

[TABLE]

A straightforward computation shows that $\dot{\Gamma_{i}}\Gamma_{i}$ is symmetric, i.e., $\Gamma_{i}$ is a variation of $\gamma$ through horizontal geodesics. Hence $d\pi J_{i}$ is a Jacobi field, where $J_{i}=\frac{d\Gamma_{i}}{ds}(0,.)$. Therefore, $d\pi(J_{1}+J_{2})$ is a Jacobi field. Moreover, $\frac{d\Gamma_{2}}{ds}(0,0)=0$ and $\frac{\mathrm{D}}{dt}\frac{d\Gamma_{2}}{ds}(0,0)=\frac{\mathrm{D}}{ds}\frac{d\Gamma_{2}}{dt}(0,0)=\xi_{2}$. Hence the solution with $J^{h}(0)=0$ and \frac{{\leavevmode\color[rgb]{0,0,0}\mathrm{D}}J^{h}}{dt}(0)=\xi_{2} is given by $J_{2}^{h}$, where $J_{2}(t)=\sin(t)Z_{2}(t)=\frac{d\Gamma_{2}}{ds}(0,t)$. It follows that the solution with $J^{h}(0)=0$ and $\frac{\mathrm{D}J^{h}}{dt}(0)=Ax$ is given by $(t\mapsto\sin(t)Y_{1}(t))^{h}$. Furthermore, $\frac{d\Gamma_{1}}{ds}(0,0)=\xi_{1}$ and $\frac{\mathrm{D}}{dt}\frac{d\Gamma_{2}}{ds}(0,0)=\frac{\mathrm{D}}{ds}\frac{d\Gamma_{2}}{dt}(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_{i}v$, $u_{i}v^{t}=vu_{i}^{t}$ and $B_{i}vv^{t}+vv^{t}B_{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^{t}Q=Q=Qvv^{t}+vv^{t}Q$ (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}=\xi_{i}-w_{i}$ is also horizontal and normal. Therefore its parallel extension is given by $U_{i}(t)=u_{i}$. Using $u_{i}v^{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_{i}v-\langle B_{i}v,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}$. Similarily, for the constant vector field $B_{1}x$, we have $(B_{1}x)^{h}=-\sin(t)W_{1}$. Implying in the expression of $J^{h}$ from part a), we arrive at the desired formula. $\blacksquare$ ; , the so-called slerp (spherical linear interpolation),; Having in mind that (cf. Pennec (2006) and Jost (2017))

[TABLE]

where $\rho_{y}(x)\coloneqq d_{\Sigma}^{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 $\Phi$. Now, for fixed $q$ and $t$, let $\nabla_{x}f$ denote the gradient of $f$ with respect to $x$ where $f(x,y):=\rho_{q}\circ\Phi(t,x,y)$. Then for any $u\in\mathrm{Hor}_{x}$, $d_{x}\mathrm{Exp}_{x}tv\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_{y}x$, $s\in[0,\varphi]$ and $\varphi=\|v\|$. 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.5, $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\Phi(t,x,y)$ parametrizes the reverse geodesic by arc length ($\|\dot{\Phi}(0,x,y)\|=\varphi$), we arrive at $J(t)=\tilde{J}((1-t)\varphi)$. Hence

[TABLE]

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

[TABLE]

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

[TABLE]

where $W_{i}=\mathrm{Log}_{\gamma(t)}q_{i}$. Now, our procedure for geodesic regression can be summarized as presented in Algorithm 1.

; An overview of OA-related dysmorphisms is shown in Figure 1.; (cf. Absil et al. (2007)). In particular, for the intrinsic treatment of the optimization problem underlying the geodesic regression we use the open-source Matlab toolbox manopt (Boumal et al., 2014). In our experiments; routine for nonlinear regression (viz. fitnlm); 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. ; Principal components; Principal components; 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); variability is better caputured by a single variable (time) than the physiological trends in HH;

; 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 Lorenzi and Pennec (2014) 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 vs. HD and HH vs. 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. Figure 1) 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 (Vincent et al., 2012). ; Group-wise analysis of femoral geodesic trends: Magnitude of differences between the group-average trends for HH vs. HD (left column), HH vs. DD (middle column), and HD vs. DD (right column) after transport to common reference shape. Only significantly different displacements ($p<0.01$) are shown (2.0e-4  4.2e-4). ; Group-wise analysis of femoral geodesic trends: Magnitude of differences between the group-average trends for HH vs. HD (left column), HH vs. DD (middle column), and HD vs. DD (right column) after transport to common reference shape. Only significantly different displacements ($p<0.01$) are shown (2.0e-4  4.2e-4). ; 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 $\geq$ 3) as the initial point for the transport path. In line with previous work (Lorenzi and Pennec, 2013), 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. ; First, we would like to use the presented methodology within the mixed-effect framework (see e.g. Bône et al. (2018)), which provides a joined analysis of longitudinal and cross-sectional variability.; 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.)

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 $\bullet$ Shape trajectory $\bullet$ Riemannian metric $\bullet$ Principal geodesic analysis $\bullet$ Geodesic regression $\bullet$ Parallel transport $\bullet$ Jacobi fields

ZIBZuse Institute Berlin, Takustraße 7, 14195 Berlin, Germany. (, , ,

) FUBFreie Universität Berlin, Arnimallee 6, 14195 Berlin, Germany. ()

1 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 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 Huckemann and Hotz (2014)), 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 Fletcher (2013). For the derivation of the corresponding Euler–Lagrange equations for some important manifolds, we refer to Machado and Silva Leite (2007). An extension to intrinsic Riemannian polynomials has been considered in Hinkle et al. (2014). Earlier related results in the framework of large deformation diffeomorphic metric mapping (LDDMM) can be found in Qiu et al. (2008) and Qiu et al. (2009). In Banerjee et al. (2016), 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 Huckemann and Hotz (2014) and Pennec (2006).

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 (Rao et al., 2004; Bossa et al., 2010), constructions based on parallel transport provide the most natural approach and have shown superior sensitivity and stability in the context of diffeomorphic registration (Lorenzi et al., 2011). 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 (Lorenzi et al., 2011) . 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. Huckemann et al. (2010) and Fletcher (2013)). However, for three and more dimensions, because of less restrictive structure, many questions remain open.

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 Kendall et al. (1999). For the relevant tools from Riemannian geometry, we refer to Gallot et al. (2005).

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 .

Let $x\in M(m,k)$, where $M(m,k)$ denotes the space of real $m\times k$ matrices. The result $\mathbb{R}_{m}^{k}\coloneqq\{x\in M(m,k):\>\sum_{i=1}^{k}x_{i}=0\}$, identified with $M(m,k-1)$, will be endowed with its canonical scalar product given by $\langle x,y\rangle=\operatorname{trace}(xy^{t})$. Denoting the Frobenius norm by \|\hbox to5.71527pt{\hss\cdot\hss}\|, we call the sphere $\mathcal{S}_{m}^{k}\coloneqq\{x\in\mathbb{R}_{m}^{k}:\>\|x\|=1\}$ pre-shape space and endow it with the spherical Procrustes metric $d(x,y)\coloneqq\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 $\Sigma_{m}^{k}=\mathcal{S}_{m}^{k}/\mathord{\sim}$. Provided that $k\geq m+1$, the dimension of $\Sigma_{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

[TABLE]

where $\lambda_{1}\geq\cdots\geq|\lambda_{m}|$ denote the pseudo-singular values of $yx^{t}$. Denoting $\mathcal{D}_{j}\coloneqq\{x\in\mathcal{S}_{m}^{k}:\>\operatorname{rank}(x)\leq j\}$, it turns out that $\Sigma_{m,m}^{k}\coloneqq\Sigma_{m}^{k}\setminus\pi(\mathcal{D}_{m-2})$ inherits a differential structure that is compatible with its quotient topology. Following Kendall et al. (1999), we refer to $\pi(\mathcal{D}_{m-2})$ as the singular part of $\Sigma_{m}^{k}$. In particular, $\Sigma_{m}^{k}$ is a strata of manifolds with varying dimensions and $\Sigma_{m,m}^{k}$ is open and dense in $\Sigma_{m}^{k}$. Away from the singular part, the quotient map $\pi$ is a Riemannian submersion . Moreover, for $k\geq 3$, the shape space $\Sigma_{1}^{k}$ (resp. $\Sigma_{2}^{k}$) is isometric to the sphere (resp. projective space). We call $x,y\in\mathcal{S}_{m}^{k}$ well positioned, and write $x\stackrel{{\scriptstyle\omega}}{{\sim}}y$, if and only if $yx^{t}$ is symmetric and $d(x,y)=d_{\Sigma}(x,y)$. For each $x,y\in\mathcal{S}_{m}^{k}$, there exists an optimal rotation $R\in\mathrm{SO}_{m}$ such that $x\stackrel{{\scriptstyle\omega}}{{\sim}}Ry$. Note that $R$ does not need to be unique. Let denote a neighborhood in $\mathcal{S}_{m}^{k}$ with radius smaller then $\pi/4$ (the diameter of $\Sigma_{m}^{k}$ is $\pi/2$) such that

[TABLE]

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

[TABLE]

is well-defined.

Due to Kendall et al. (1999) the vertical space at $x\in\mathcal{S}_{m}^{k}$ is given by

[TABLE]

and the horizontal space is given by

[TABLE]

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

Now, let $\exp$ and $\log$ denote the exponential and logarithm map of the pre-shape space. For $x\stackrel{{\scriptstyle\omega}}{{\sim}}y$ the geodesic from $x$ to $y$ given by

[TABLE]

with $\varphi=\arccos(\langle x,y\rangle),\>0\leq t\leq 1$, is horizontal. Hence $\Phi$ 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 2.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}$.

1. (a)

$\mathrm{ver}_{x}(w)=Ax$* if and only if $A$ solves the Sylvester equation*

[TABLE]

. 2. (b)

Fix $R\in\mathrm{SO}_{m}$. Then $\mathrm{ver}_{Rx}(Rw)=R\mathrm{ver}_{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\coloneqq\operatorname{rank}(x)$. Suppose without loss of generality that $j>1$ and write $x=\binom{x_{1}}{0}$ with

[TABLE]

where $w_{1}$ is $j\times k$. .

For (b) note that $\langle Rw,Rx\rangle=\langle w,x\rangle=0$, i.e., $w\in\mathrm{T}_{x}S$ implies $Rw\in\mathrm{T}_{Rx}S$. Now, $\mathrm{ver}_{Rx}(Rw)=BRx$ where $B$ is the solution of . Hence $B=RAR^{t}$, which implies that $\mathrm{ver}_{Rx}(Rw)=R.\mathrm{ver}_{x}(w)$ . $\blacksquare$

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\coloneqq\{x\in\mathcal{S}_{m}^{k}:\>\operatorname{rank}(x)\geq 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)}\Sigma_{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

[TABLE]

In the following [\hbox to5.71527pt{\hss\cdot\hss},\hbox to5.71527pt{\hss\cdot\hss}] 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 $\|\gamma\|=1$, and $\frac{\mathrm{D}^{2}W}{\mathrm{d}t^{2}}=\ddot{W}$, etc. We set222Note that the Riemannian exponential map of the shape space denoted by $\widetilde{\exp}$ satisfies $\pi(\exp_{x}u)=\widetilde{\exp}_{\pi(x)}(\mathrm{d}_{x}\pi(u))=\widetilde{\exp}_{\pi(x)}(\mathrm{d}_{x}\pi(u^{h}))$.

[TABLE]

For the computation of the Fréchet mean (cf. Huckemann et al. (2010) and Pennec (2006)) $\pi(\bar{q})$ of the shapes $\pi(q_{1}),\cdots,\pi(q_{N})$ with $q_{i}\in$, i.e.,

[TABLE]

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

[TABLE]

and set

[TABLE]

A suitable initial value is the normalized Euclidean mean

[TABLE]

The total variance of $q=(q_{1},\cdots,q_{N})$ reads

[TABLE]

2.2 Parallel Transport

Next, we parallel transport in the shape space and its relation to parallel transport in the pre-shape space.333Essentially, part (a) of Proposition 2.2 was recently also obtained by Kim et al. (2018).

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

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

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

[TABLE] 2. (b)

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

[TABLE] 3. (c)

Let $Cv$ denote the orthogonal projection of $u$ on $\mathrm{Skew}_{m}\cdot v$, i.e. $Cvv^{t}+vv^{t}C=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

[TABLE]

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

[TABLE]

Proof. (a) $\mathrm{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 2.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 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\coloneqq\|W\gamma^{t}-\gamma W^{t}\|^{2}+\langle W,\gamma\rangle^{2}$, since $f(0)=0$. If equation (4) holds, then

[TABLE]

and

[TABLE]

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^{t}A=\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$. Now, (4) implies

[TABLE]

(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 . To prove (7), insert $v=\frac{1}{\varphi}\log_{x}y=\frac{y-x\cos(\varphi)}{\sin(\varphi)}$ and $\dot{\gamma}=\frac{-1}{\varphi}\log_{y}x$ into (6). $\blacksquare$

Note that, , the differential equation for the h-parallel transport can also be written as

[TABLE]

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

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\coloneqq\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

Hence, for a horizontal geodesic $\gamma$, any geodesic variation of $\pi\circ\gamma$ in the latter space reads $\pi\circ\Gamma$ with $\Gamma$ a variation of $\gamma$ through horizontal geodesic. Thus the variation field \frac{\mathrm{d}}{\mathrm{d}s}(\pi\circ\Gamma(s,\hbox to5.71527pt{\hss\cdot\hss}))|_{s=0}=\mathrm{d}\pi(\frac{\mathrm{d}}{\mathrm{d}s}\Gamma(s,\hbox to5.71527pt{\hss\cdot\hss})|_{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.

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. Fletcher (2013) and Jost (2017)).

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_{\Sigma}$, i.e.

[TABLE]

We can assume that $t_{1}=0$ and $t_{N}=1$. While the authors of Fletcher (2013) and Machado and Silva Leite (2007) 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

[TABLE]

The reason is that geodesic computations in terms of the function $\Phi$ defined in equation (1) are more efficient. Model estimation is then formulated as the least-squares problem

[TABLE]

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 (Boumal et al., 2014) with a Hessian approximation based on finite differences and use $(q_{1},\omega(q_{1},q_{N}))$ as initial guess.

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.

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 (Kellgren and Lawrence, 1957) 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 mm matrix, 0.7 mm slice thickness) using a state-of-the-art automatic segmentation approach (Ambellan et al., 2018). 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.

4.2 Geodesic Modeling of Femoral Trajectories

We apply the geodesic regression approach detailed in Section 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 , 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 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 Figure 2.

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 Fletcher (2013), a generalization to manifolds is defined as

[TABLE]

with $F(\gamma)$ and $G(x)$ as defined in equations (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 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 . Based on the coefficient of determination we also test for the significance of the estimated trends employing permutation tests as suggested in Fletcher (2013). 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

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.

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 Muralidharan and Fletcher (2012). 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.

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

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 OAI444The 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 (Boumal et al., 2014).