Geodesics in the TPS Space

Abstract

In shape analysis, the interpolation of shapes’ trajectories is often performed by means of geodesics in an appropriate Riemannian Shape Space. Over the past several decades, different metrics and shape spaces have been proposed, including Kendall shape space, LDDMM based approaches, and elastic contour, among others. Once a Riemannian space is chosen, geodesics and parallel transports can be used to build splines or piecewise geodesics paths. In a recent paper, we introduced a new Riemannian shape space named TPS Space based on the Thin Plate Spline interpolant and characterized by an appropriate metric and parallel transport rule. In the present paper, we further explore the geometry of the TPS Space by characterizing the properties of its geodesics. Several applications show the capability of the proposed formulation to conserve important physical properties of deformation, such as local strains and global elastic energy.

1. Introduction

In shape analysis, the interpolation of shapes’ trajectories is often performed by means of geodesics in appropriate Riemannian Shape Spaces. Different proposed metrics and shape spaces include Kendall shape space [1,2], LDDMM-based approaches [3,4,5,6], and elastic contour [7]. Once a Riemannian space is chosen, geodesics and parallel transports can be used to build splines or piecewise geodesics paths [1,8,9]. If the torsion of the connection defined on the Riemannian space does not vanish, then a difference can appear between geodesics and autoparallel lines [10]. In recent papers [11,12,13,14], the present authors introduced and developed a new Riemannian shape space named TPS Space based on the Thin Plate Spline interpolant and characterized by an appropriate metric and a parallel transport rule, and the efficiency of this TPS parallel transport was compared with other methods in [14]. One of the main features of the TPS connection is its independence from the path. In fact, TPS parallel transport is named Direct Transport. This feature leads to a flat space (vanishing Riemannian Curvature) with torsion. In previous contributions, the geodesics of the TPS Space have never been studied. The present paper aims at studying the characterization of geodesics and parallel lines in the TPS space together with numerical techniques to compute them. Looking for geodesics’ computation is of particular importance in a variety of morphological analyses, such as in spline regression, which have applications in a wide range of shape analysis tasks spanning fields from medical/clinical investigations to biological research. As will be clarified in the following pages, in a connection with torsion, geodesics and autoparallel lines can coincide; however, this is not a rule [10,15]. In particular, it happens only if the torsion is completely skew-symmetric. On the other hand, the TPS connection is defined by directly assigning a parallel transport rule without analytically defining the corresponding covariant derivative, ∇. For this reason, it is not possible to calculate the Christoffel symbols of the connection or the components of the torsion to establish whether it is completely skew-symmetric [15]. In the following, we exploit the qualitative definitions of Direct Transport and TPS metric in 2D to build the geodesics and the autoparallel lines of the TPS space numerically. In particular, geodesics are built by minimising the length of the path connecting two given points (the initial and final shapes), while autoparallel lines are calculated via shooting from a point (the initial configuration) and a vector (the initial deformation velocity). Our set of examples allows the geometry of the TPS Space to be explored by characterizing the properties of its geodesics and parallel lines. The paper is organised as follows:

• In Section 2, the general definitions of the two families of curves in Riemannian spaces are summarized.
• In Section 3, the main concept defining the TPS Space are recalled.
• In Section 4, the novel contribution of this paper is presented, that is, the construction of and comparison between autoparallel and geodesic lines in TPS Space.
• In Section 6, the numerical results are shown in order to discuss and compare the main features of autoparallel and geodesic lines in TPS Space.

2. Geodesics in Riemannian Manifolds

In this section, we sketch several concepts in differential geometry, referring to [10,15,16] for details.

In differential geometry, a Riemannian manifold $(\mathcal{M},g)$ is a smooth manifold $\mathcal{M}$ equipped with a positive-definite inner product g on the tangent space $T_p\mathcal{M}$ at each point p.

A connection on the manifold is a rule allowing for parallel transport vectors along a smooth curve $\gamma \left(t\right)\in \mathcal{M}$. This rule allows for the comparison of vectors belonging to different tangent spaces.

To be more precise, any parallel transport rule ${\tau }_{b,a}$ along a path from a to b has to fulfill the following properties:

$$\begin{array}{ccc}\hfill {\tau }_{b,a}:T_a\mathcal{M}& \to & T_b\mathcal{M},\mathrm{is}\mathrm{linear},\mathrm{and}\text{non-singular}.\hfill \\ \hfill V_a& \mapsto & V_b;\hfill \end{array}$$

moreover, for any point c on the path

$${\tau }_{b,c}\circ {\tau }_{c,a}={\tau }_{b,a}$$

(1)

It follows from this that ${\tau }_{a,a}$ is the identity on $T_a\mathcal{M}$, and ${\tau }_{a,b}={\left({\tau }_{b,a}\right)}^{-1}$.

This procedure defines the connection, inducing a covariant derivative, ∇, that can be then calculated by the following limit:

$${\nabla }_{V_p}U=\underset{h\to 0}{lim}\frac{{\tau }_{h,0}^{-1}{U}_{\gamma \left(h\right)}-U_{\gamma \left(0\right)}}{h}$$

(2)

A vector field V is said to be parallel along the curve $\gamma$ if the following holds:

$${\nabla }_{\dot{\gamma }}V=0.$$

(3)

A connection is considered compatible with a metric g if the parallel transport is an isometry, i.e., if $g_a(V_a,W_a)=g_b({\tau }_{b,a}\left(V_a\right),{\tau }_{b,a}\left(W_a\right))$ for each pair of vectors $V_a,W_a$ along each path.

In terms of connection, this means that:

$$g({\nabla }_{X}Y,Z)+g({\nabla }_{X}Z,Y)=X·g(Y,Z)$$

(4)

The torsion of the connection ∇ is a tensor field defined as

$${\nabla }_{V}W-{\nabla }_{W}V-\left[V,W\right],$$

(5)

with $\left[·,·\right]$ the Lie bracket.

A connection is called symmetric when the torsion is null for all $V,W$. A fundamental result of Riemannian Geometry is the existence of a unique symmetric connection compatible with the metric g, named the Levi–Civita (LC) connection, which we call ${\nabla }^g$.

In its more general meaning, a geodesic is a curve $\gamma \left(t\right)\in \mathcal{M}$ the tangent field $\dot{\gamma }\left(t\right)$ of which is parallel along the curve itself:

$${\nabla }_{\dot{\gamma }}\dot{\gamma }=0.$$

(6)

In other terms, a geodesic is the curve on the manifold where one walks, maintaining the same direction, according to the given parallel transport rule. In this general sense, the geodesic can be called an autoparallel line [10].

On the other hand, if the connection is the Levi–Civita one, the geodesics acquire an additional property. Given two points, $p,q$, belonging to a convex neighbourhood of $(\mathcal{M},g)$, a geodesic from p to q is the shortest path joining p and q according to the metric g. The distance $d(p,q)$ between p and q can be calculated as ${\int }_p^q\sqrt{g(\dot{\gamma }\left(t\right),\dot{\gamma }\left(t\right))}dt$.

In this paper, following [10], by geodesics we mean only those geodesics of the Levi–Civita connection, ${\nabla }^g$; we refer to the geodesics of a different connection, ∇, as autoparallel lines.

In general, the difference between a Levi–Civita connection, ${\nabla }^g$, and any connection ∇ is a $(2,1)$ tensor field D [15]:

$${\nabla }_{X}Y={\nabla }_X^{g}Y+D(X,Y)X,Y\in T\mathcal{M}$$

(7)

The symmetric and the antisymmetric part of D have direct geometric meanings:

• A connection ∇ is torsion-free if and only if D is symmetric.
• A connection ∇ has the same geodesics as the Levi–Civita connection ${\nabla }^g$ if and only if D is skew-symmetric.

As a consequence of (4), the connection ∇ is compatible with the metric g if and only if D belongs to the space

$${\mathcal{D}}^g:=T\mathcal{M}\otimes \left({\mathsf{\Lambda }}^{2}T\mathcal{M}\right)=\left\{D\in {\otimes }^{3}T\mathcal{M}|D(X,V,W)+D(X,W,V)=0\right\}$$

(8)

Finally, by means of (6), (7), and (8), we find that a connection ∇ on $(\mathcal{M},g)$ is both metric and geodesic-preserving if and only if its torsion lies in ${\mathsf{\Lambda }}^{3}T\mathcal{M}$, i.e., if it is completely skew-symmetric. In this case, $2D=T$ and

$${\nabla }_{X}Y={\nabla }_X^{g}Y+\frac{1}{2}T(X,Y,-)$$

(9)

3. Geometry of the TPS Space

In the present section, we summarize the formulation of the TPS Space presented and developed in [11,12,13]. Both the metric g and parallel transport are based on the interpolation function called Thin Plate Spline (TPS). Then, in order to introduce the Riemannian structure of the TPS Space, we need to first summarize the formulation of the TPS. Furthermore, in order to introduce the metric g we need to explain strain energy, bending energy, and body bending energy and how they are used to build g.

3.1. Thin Plate Spline

Let ${\mathcal{E}}^m$ be the m-dimensional Euclidean space and ${\mathsf{\Omega }}_X,{\mathsf{\Omega }}_{X^{\prime }}\subset {\mathcal{E}}^m$ be two regular regions representing the undeformed (source) and deformed (target) configurations of a body (Figure 1), respectively. We label as x the points in ${\mathsf{\Omega }}_X$ and as $x^{\prime }$ the points in ${\mathsf{\Omega }}_{X^{\prime }}$. The displacement field is represented by the following difference vectors:

$$u\left(x\right)=x^{\prime }-x$$

(10)

If the configurations are sampled in k points (landmarks), then ${\mathsf{\Omega }}_X$ and ${\mathsf{\Omega }}_{X^{\prime }}$ are named k-configurations and represented by the $k\times m$ matrices X and $X^{\prime }$, respectively. The displacements experienced by the k landmarks can be collected in the $k\times m$ matrix $U=X^{\prime }-X$. Because translations do not affect the shape, they are filtered out by centring configurations in the origin of ${\mathcal{E}}^m$. The set of all centred k-configurations is named Centred Configuration Space $\mathcal{C}{\mathcal{C}}_m^k$. Given a configuration X, the corresponding centred configuration can be represented in two different ways: by a $k\times m$ matrix obtained as $X_C=CX$ or, alternatively, by a $(k-1)\times m$ matrix $X_H=HX$, where $C=I_k-\frac{1}{k}{1}_k{1}_k^T$, $I_k$ is the $k\times k$ identity matrix and $1_k$ is a $k\times 1$ column of ones, while H is the Helmert sub-matrix. The $jth$ row of the Helmert sub-matrix H is obtained by

$$(h_j,...,h_j,-j{h}_j,0,...,0),h_j=-{\left(j(j+1)\right)}^{-1/2}$$

and thus, the $jth$ row consists of $h_j$ repeated j times followed by $j{h}_j$ and then $k-j-1$ zeros, $j=1,...,k-1$. The notable property $H^{T}H=C$ can be used to switch from one parametrization to the other. In the following we use only centred configurations; we therefore remove the subscript, C or H, by specifying, when necessary, which parametrization we are using. Let $X,X^{\prime }\in \mathcal{C}{\mathcal{C}}_m^k$ be a pair of centred configurations. The Thin Plate Spline (TPS) $\mathsf{\Phi }$ is a function that interpolates, in ${\mathcal{E}}^m$, the deformation from X to $X^{\prime }$. The TPS is parametrized by the pair $(A_X,W_X)$ where $A_X\in GL\left(m\right)$ is a linear transformation of ${\mathcal{E}}^m$ represented by a $m\times m$ matrix and $W_X$ is a $k\times m$ matrix. Given a point $x\in {\mathcal{E}}^m$ and a centred configuration $X\in \mathcal{C}{\mathcal{C}}_m^k$, we obtain

$$x^{\prime }=\mathsf{\Phi }\left(x\right)=A_{X}x+W_X^{T}s\left(x\right)$$

(11)

where $s\left(x\right)={(\sigma (x-x_1),...,\sigma (x-x_k))}^{T}C$ is a $(k\times 1)$ matrix, $x_i\in X$ is the position of the i-th landmark, and

$$\sigma \left(h\right)=\left\{\begin{array}{cc}{\left|\right|h\left|\right|}^2{log\left(\right|\left|h\right||}^2)\hfill & if\left|\right|h\left|\right|>0;\hfill \\ 0\hfill & if\left|\right|h\left|\right|=0.\hfill \end{array}\right.\mathrm{for}m=2$$

$$\sigma \left(h\right)=\left\{\begin{array}{cc}-\left|\right|h\left|\right|\hfill & if\left|\right|h\left|\right|>0;\hfill \\ 0\hfill & if\left|\right|h\left|\right|=0.\hfill \end{array}\right.\mathrm{for}m=3$$

Equation (11), applied landmark-wise to X, reads:

$$X^{\prime }=X{A}_X^T+S_X{W}_X,\mathrm{with}{\left(S_X\right)}_{ij}=C_{ik}\sigma (x_k-x_l)C_{lj}.$$

(12)

Because X and $X^{\prime }$ are centred, Equation (12) represents $m\times (k-1)$ interpolation constraints, while the matrices $(A_X,W_X)$ consist of $m\times m+k\times m=(m+k)m$ parameters. In order to solve the interpolation problem we need to introduce $m\times (m+1)$ more constraints on $W_X$, uncoupling the affine and non-affine parts:

$$1_k^T{W}_X=0,X^T{W}_X=0.$$

(13)

For a given pair $(X,X^{\prime })$ there exists a unique set of parameters for the pair $\left(A_X,W_X\right)$ that solve the problem (12), constrained with (13):

$$A_X^T={\mathsf{\Gamma }}_{21X}{X}^{\prime },W_X={\mathsf{\Gamma }}_{11X}{X}^{\prime },$$

(14)

where

$$\begin{array}{ccc}\hfill {\mathsf{\Gamma }}_{21X}& =& {\left(X^T{S}_X^{-1}X\right)}^{-1}{X}^T{S}_X^{-1}\hfill \\ \hfill {\mathsf{\Gamma }}_{11X}& =& S_X^{-1}-S_X^{-1}X{\mathsf{\Gamma }}_{21X}\hfill \end{array}$$

are a $m\times k$ and a $k\times k$ matrices, respectively, which only depend on the source configuration X. We note that the inverse of the singular matrix $S_X$ is obtained by means of the Helmert matrix as

$$S_X^{-1}=H^T{\left(H{S}_X{H}^T\right)}^{-1}H$$

(15)

Finally the target $X^{\prime }$ can be represented as the deformation of X:

$$X^{\prime }=X{\mathsf{\Gamma }}_{21X}{X}^{\prime }+S_X{\mathsf{\Gamma }}_{11X}{X}^{\prime }=X{A}_X^T+S_X{W}_X.$$

(16)

The $k\times k$ matrix ${\mathsf{\Gamma }}_{11X}$ is called the Bending Energy matrix and is used to extract the non linear part of the deformation. Because it vanishes on affine deformations, its eigenvectors associated with non vanishing eigenvalues are only $(k-1-m)$. These are called principal warp eigenvectors, and represent the principal modes of deformation of the shape X [11].

3.2. Energies

It has been proven [17] that TPS is the only interpolating function that minimizes the bending energy J, which gauges the second derivative of the displacement field:

$$J={\int }_{{\mathbb{R}}^m}{\nabla }^{2}u·{\nabla }^{2}u$$

(17)

Note that the bending energy is defined as an integral on the whole, ${\mathbb{R}}^m$. In [12], in order to provide a mechanical interpretation of bending energy, the concept of body-bending energy was introduced, allowing the integration to be performed entirely inside the body:

$$J_{\mathsf{\Omega }}={\int }_{\mathsf{\Omega }}{\nabla }^{2}u·{\nabla }^{2}u$$

(18)

The body-bending energy is slightly smaller than the bending energy, and the decay $\rho =J/J_{\mathsf{\Omega }}$ can be used to quantify the difference. Both J and $J_{\mathsf{\Omega }}$ can be used as pseudo-metrics on $T\mathcal{C}{\mathcal{C}}_m^k$, as they measure the difference between two configurations and vanish in affine deformations. In order to endow $X\in \mathcal{C}{\mathcal{C}}_m^k$ with a Riemannian metric we need a non-singular distance. Here, we propose slightly modifying the Dirichlet Energy used in the deformable templates method [18] to obtain the following expression for the strain energy, $\phi$, stored by the body:

$$\phi =\frac{1}{2}{\int }_{\mathsf{\Omega }}E·E$$

(19)

where

$$E=\frac{\nabla u+{(\nabla u)}^T}{2}$$

(20)

is the strain tensor (in the case of small displacements). We note that $\phi$ vanishes on the rotational part of the local deformation. On the other hand we stress that (19) is a global measure that should be calculated via integration, starting from local measures, by means of a discretization of the whole domain $\mathsf{\Omega }$. In [12] it is shown that (19) can be calculated directly as a global quantity starting from landmark displacements, at least for bilinear deformations in 2D. In fact, in that case, when a rectangle bends in a generic trapezoid the deformation can be parametrized as follows:

$$u=(A-I)x+\left[\left(\chi \otimes e_1\otimes e_2\right)x\right]x$$

(21)

where $\chi ={\chi }_1{e}_1+{\chi }_1{e}_2$ are the bending with respect to the two axes, A is a linear transformation, and I is the $(m\times m)$ identity matrix. Then, the strain energy can be calculated as the sum of one contribution depending on the norm of A and a second contribution proportional to the bending energy, J. The proportionality coefficient depends only on geometrical quantities of $\mathsf{\Omega }$:

$$\phi =\frac{1}{2}\mathcal{A}(A-I)·(A-I)+\frac{{\mathcal{I}}_p}{4\rho \mathcal{A}}J$$

(22)

where ${\mathcal{I}}_p$ and $\mathcal{A}$ are the polar inertia and the area of $\mathsf{\Omega }$, respectively. In the next section, we show that both J and $(A-I)·(A-I)$ can be calculated by means of the TPS parameters ${\mathsf{\Gamma }}_{21X}$ and ${\mathsf{\Gamma }}_{11X}$. Then, the expression (22) can be globally calculated starting from landmark displacements by means of TPS.

3.3. TPS Metric

In [11], it is shown that the value of the bending energy J associated with a displacement, U, of a configuration, X, obtained using the quadratic form

$$J\left(U\right)=\nu \pi \mathrm{Tr}\left(U^{T}BU\right).$$

(23)

where $\nu =16$ for $m=2$ and is $\nu =8$ for $m=3$. For this reason the matrix $B:={\mathsf{\Gamma }}_{11X}$ is called the Bending Energy Matrix of X. This fact allows us to evaluate the bending energy directly by means of a closed form expression, avoiding the need to discretize the configuration in a huge number of triangles. Furthermore we note that the Bending Energy Matrix depends only on X and can be used as a pseudo-metric on $T\mathcal{C}{\mathcal{C}}_m^k$. Let two given configurations, X and $X^{\prime }=X+U$, related landmark-wise by a bilinear deformation as (21); then, the strain energy (22) can be calculated as

$$\phi \left(U\right)=\frac{1}{2}\mathcal{A}\mathrm{Tr}\left(U^T{\mathsf{\Gamma }}_{21X}^T{\mathsf{\Gamma }}_{21X}U\right)+\frac{4\pi {\mathcal{I}}_p}{\rho \mathcal{A}}\mathrm{Tr}\left(U^{T}BU\right)$$

(24)

and the average strain energy on the body can be obtained as

$$\begin{array}{cc}\hfill \overline{\phi }\left(U\right)& =\frac{\phi \left(U\right)}{\mathcal{A}}=\frac{1}{2}\mathrm{Tr}\left(U^T{\mathsf{\Gamma }}_{21X}^T{\mathsf{\Gamma }}_{21X}U\right)+\frac{4\pi {\mathcal{I}}_p}{\rho {\mathcal{A}}^2}\mathrm{Tr}\left(U^{T}BU\right)\end{array}$$

(25)

While this expression is valid only for bilinear deformations, it can be generalized by assuming certain approximations concerning, in particular, the decay, $\rho$. In [13], it is further shown that for the body-bending energy calculation it is possible to define a symmetric matrix, $B_{\mathsf{\Omega }}$, such that the following holds:

$$J_{\mathsf{\Omega }}\left(U\right)=\mathrm{Tr}\left(U^T{B}_{\mathsf{\Omega }}U\right)$$

(26)

Then, in general, the decay, $\rho \left(U\right)$, is not isotropic and can be calculated by means of the Rayleigh quotient:

$$\rho \left(U\right)=\frac{J\left(U\right)}{J_{\mathsf{\Omega }}\left(U\right)}=\frac{16\pi \mathrm{Tr}\left(U^{T}BU\right)}{\mathrm{Tr}\left(U^T{B}_{\mathsf{\Omega }}U\right)}$$

(27)

In the two-dimensional case, the BEB matrix $B_{\mathsf{\Omega }}$ is defined as follows:

$$B_{\mathsf{\Omega }}={\mathsf{\Gamma }}_{11X}{C}_{\mathsf{\Omega }}{\mathsf{\Gamma }}_{11X}+{\mathsf{\Gamma }}_{11X}{C}_{\partial \mathsf{\Omega }}{\mathsf{\Gamma }}_{11X}$$

(28)

$$\begin{array}{ccc}\hfill {\left(C_{\mathsf{\Omega }}\right)}_{ij}& =& 8{\alpha }_i\sigma (x_i-x_j)\hfill \\ \hfill {\left(C_{\partial \mathsf{\Omega }}\right)}_{ij}& =& \sum _{p=1}^{k-q}{\int }_0^1\left[{\nabla }^2{s}_i(x_p+\zeta {\ell }_p)\nabla s_j(x_p+\zeta {\ell }_p)\right.\hfill \\ \hfill & -& \left.s_i(x_p+\zeta {\ell }_p)\mathrm{div}\left({\nabla }^2{s}_j(x_p+\zeta {\ell }_p)\right)\right]·*\left({\ell }_p\right)d\zeta \hfill \end{array}$$

(29)

where $s\left(x\right)={(\sigma (x-x_1),...,\sigma (x-x_k))}^T$ and the values of the angles ${\alpha }_i$ are

$${\alpha }_i=\left\{\begin{array}{cc}2\pi \hfill & \mathrm{if}{x}_i\in \mathsf{\Omega };\hfill \\ arccos(\frac{{\ell }_i}{\parallel {\ell }_i\parallel }·\frac{{\ell }_{i+1}}{\parallel {\ell }_{i+1}\parallel })\hfill & \mathrm{if}{x}_i\in \partial \mathsf{\Omega }.\hfill \end{array}\right.$$

(30)

while ${\ell }_i=(x_{i+1}-x_i)$, $*\left({\ell }_p\right)$ is the vector ${\ell }_p$ rotated clockwise by $\pi /2$, q is the number of landmarks that does not lie on the boundary $\partial \mathsf{\Omega }$, and $(k-p)$ is the number of landmarks lying on $\partial \mathsf{\Omega }$ (see Figure 1). As $\rho$ can assume different values depending on the direction of the deformation, while we need a metric depending only on X, in the following we assume an isotropic decay $\overline{\rho }$:

$$\overline{\rho }=\frac{16\pi }{(k-z)}Tr\left(B{B}_{\mathsf{\Omega }}^{-1}\right)$$

(31)

where $z\ge (m+1)$ is the number of vanishing eigenvalues of $\left(B{B}_{\mathsf{\Omega }}^{-1}\right)$. Then, in the following, we approximate the calculus of the BEB by assuming

$$B_{\mathsf{\Omega }}\simeq \overline{\rho }B$$

(32)

Finally, we generalize (25), defining the distance between two generic configurations X and $X^{\prime }$ (called $\mathsf{\Gamma }$-Energy), as follows:

$$\mathsf{\Gamma }(X,X^{\prime }):=\mathrm{Tr}\left({(X^{\prime }-X)}^{T}G(X^{\prime }-X)\right)$$

(33)

where

$$G:={\mu }_1{\mathsf{\Gamma }}_{21X}^T{\mathsf{\Gamma }}_{21X}+{\mu }_2{\mathsf{\Gamma }}_{11X}{\mu }_1=\frac{1}{2}\mathrm{and}{\mu }_2=\frac{4\pi {\mathcal{I}}_p}{\overline{\rho }{\mathcal{A}}^2}$$

(34)

We note that ${\mathsf{\Gamma }}_{21X},{\mathsf{\Gamma }}_{11X},{\mu }_2$ depend only on the source configuration, X. From a mechanical point of view, we can define the $\mathsf{\Gamma }$-energy $\mathsf{\Gamma }\left(U\right)$ as the average strain energy $\overline{\phi }\left(\tilde{U}\right)$, evaluated on a more simple deformation $\tilde{U}$ characterized by the same uniform component ${\tilde{U}}_u=U_u$ of U and a bilinear deformation ${\tilde{U}}_{nu}$ storing the same body bending energy of $U_{nu}$, i.e., such that $J_{\mathsf{\Omega }}\left({\tilde{U}}_{nu}\right)=J_{\mathsf{\Omega }}\left(U_{nu}\right)$. In [12], it is shown that the $\mathsf{\Gamma }$-energy is a good approximation of the strain energy for more general deformations as well as bilinear ones. Then, the TPS Space [11,12,13] can be defined as the $\mathcal{C}{\mathcal{C}}_m^k$ equipped with the TPS metric tensor:

$$g(U,V):=\mathrm{Tr}\left(U^{T}GV\right)$$

(35)

In particular, the affine and non-affine components of the metric are defined by the sub-metrics

$$g_u(U,V):=\mathrm{Tr}\left(U^T{G}_{u}V\right)g_{nu}(U,V):=\mathrm{Tr}\left(U^T{G}_{b}V\right)$$

(36)

where

$$G_u:={\mu }_1{\mathsf{\Gamma }}_{21X}^T{\mathsf{\Gamma }}_{21X}{G}_{nu}={\mu }_2{\mathsf{\Gamma }}_{11X}$$

(37)

Alignments: OPA and MOPA Techniques

After the TPS metric tensor G is introduced, we introduce a technique for managing rotations in order to align two configurations, X and Y, based on minimization of the distance, defined as

$$d(X,Y)=\underset{Q\in S{O}_m}{inf}\sqrt{\mathrm{Tr}\left({(YQ-X)}^T{G}^{\alpha }(YQ-X)\right)}.$$

The aligned configuration, $\widehat{Y}$, is obtained by means of an optimal rotation, $\widehat{Q}$, minimizing d.

$$\widehat{Y}=Y\widehat{Q}$$

where $\widehat{Q}=\mathrm{argmin}g\left((YQ-X),(YQ-X)\right)$. According to this definition, $\widehat{Q}$ turns out to be the rotational component of the polar decomposition of $Y^{T}GX$. When $\alpha =0$, we obtain the rotational component of $Y^{T}X$, the classical Ordinary Procrustes Analysis (OPA). When $\alpha =1$, then $Y^{T}GX=A_X$. In the latter case, we define the alignment Modified OPA, or MOPA [11].

3.4. TPS Direct Transport

The connection called the TPS connection was introduced in [11] and developed in [13]. It has the following properties:

• It is compatible with the TPS metric;
• It is compatible with the decomposition provided in (12);
• It is independent of the path.

We assume all the configurations to be centred and represented by Helmertized landmarks; then, if not otherwise specified, each matrix is a $(k-1)\times m$ matrix. Furthermore, deformation vectors have a subscript denoting the starting point, that is, the source configuration; for details, see [11].

Let X and Y be two source configurations and let $V_X$ and $V_Y$ be the two associated deformation vectors, provided by

$$\begin{array}{ccc}\hfill & V_X& =X^{\prime }-X=X(A_X^T-I)+S_X{W}_X,\hfill \\ & V_Y& =Y^{\prime }-Y=Y(A_Y^T-I)+S_Y{W}_Y.\hfill \end{array}$$

(38)

We can then say that $V_Y$ is the parallel transport of a given $V_X$, that is, $V_Y={\tau }_{Y,X}\left(V_X\right)$, if and only if the uniform part of $V_Y$ equals that of $V_X$

$$A_Y=A_X;$$

(39)

and the non uniform part $W_Y$ of $V_Y$ solves the linear systems

$$Y^T{W}_Y=X^T{W}_X=0Q_Y{E}_Y^T{W}_Y=Q_X{E}_X^T{W}_X,$$

(40)

where the $(k-1)\times (k-1-m)$ body principal warps matrix $E_X$ collects all the body principal warps of X and $Q_X$ is a suitable $(k-1-m)\times (k-1-m)$ orthogonal matrix (i.e., $Q_X^T{Q}_X=I$) defined on each configuration X and representing a rotation or reflection of the principal warps. After being chosen, a configuration P as a Pole for the space $Q_X$ is estimated, minimising the Euclidean distance $\parallel E_X{Q}_X^{\mathcal{B}}-E_P\parallel$ between the rotated principal warps of X and the corresponding basis on the pole, P.

The principal warps matrix can be built as follows:

• Perform a TPS analysis on X and find the $S_X$ and ${\mathsf{\Gamma }}_{11X}$;
• Perform an eigenvalue analysis on ${\mathsf{\Gamma }}_{11X}$ and obtain ${\mathsf{\Gamma }}_{11X}=\mathsf{\Gamma }\mathsf{\Lambda }{\mathsf{\Gamma }}^T$, where $\mathsf{\Gamma }$ is the $(k-1)\times (k-1)$ matrix containing the eigenvectors ${\gamma }_i$ in column and $\mathsf{\Lambda }$ is the diagonal $(k-1)\times (k-1)$ matrix of the eigenvalues ${\lambda }_1,\dots ,{\lambda }_{k-1}$ ordered by increasing magnitude (the first m eigenvalues will be equal to 0);
• Drop the first m columns from $\mathsf{\Gamma }$ by obtaining the $(k-1)\times (k-1-m)$ matrix $\overline{\mathsf{\Gamma }}$, containing the principal warp eigenvectors by column;
• Drop the first m rows and the first m columns from $\mathsf{\Lambda }$ by obtaining the $(k-1-m)\times (k-1-m)$ matrix $\overline{\mathsf{\Lambda }}$;
• Define the $(k-1)\times (k-1-m)$ matrix $E_X=S_X\overline{\mathsf{\Gamma }}{\overline{\mathsf{\Lambda }}}^{1/2}$.

The same steps must be used to build the principal warps matrix $E_Y$ on the target configuration.

The first equation of (40) constrains $W_b$ to be orthogonal to the affine part, while the second defines the isometry in the subspace of the non-affine deformations. This last requirement implies the conservation of the total bending energy. The system (40) can be written as

$$\left[\begin{array}{c}{Y}^T\\ Q_Y{E}_Y^T\end{array}\right]\left[\begin{array}{c}{W}_Y\end{array}\right]=\left[\begin{array}{c}{X}^T\\ Q_X{E}_X^T\end{array}\right]\left[\begin{array}{c}{W}_X\end{array}\right].$$

This can be re-written as

$$W_Y=M_Y^{-1}{M}_X{W}_X$$

And so:

$$V_Y=\left(Y{\mathsf{\Gamma }}_{21X}+\sqrt{\frac{{\mu }_2\left(X\right)}{{\mu }_2\left(Y\right)}}{S}_Y{M}_Y^{-1}{M}_X{\mathsf{\Gamma }}_{11X}\right)V_X,$$

(41)

It is worth noting that Equation (41), characterizing $V_Y$ as the parallel transport of $V_X$, depends only on quantities related to the startpoint, X, and endpoint, Y, of the transport, and does not depend on the path. For this reason, the TPS connection is characterized by vanishing curvature and non-vanishing torsion. Moreover, it is easy to check that (41) is compatible with the TPS metric G and with the decomposition provided in (12).

4. Geodesics and Autoparallel Lines in TPS Space

In the present section, we introduce the main contribution of the present paper, that is, to show the most important features of the geodesics and autoparallel lines in the TPS space and compare the two families of lines.

In the previous section, the TPS connection has been defined by directly assigning the parallel transport rule (41) without analytically defining the corresponding covariant derivative ∇. For this reason, it is not immediately necessary to calculate the Christoffel symbols of the connection or the Christoffel symbols of the corresponding Levi–Civita connection ${\nabla }^g$, nor the components of the torsion for the purpose of establishing whether it is completely skew-symmetric [15] and thus whether or not the autoparallel lines and geodesic lines coincide.

After $V_X,V_Y$ in (41) is substituted with $\dot{X}0,\dot{X}\left(t\right)$ and $X,Y$ with $X0,X\left(t\right)$, we obtain

$$\dot{X}\left(t\right)=\left(X\left(t\right){\mathsf{\Gamma }}_{21X0}+\sqrt{\frac{{\mu }_2\left(X0\right)}{{\mu }_2\left(X\left(t\right)\right)}}{S}_X\left(t\right)M_{X\left(t\right)}^{-1}{M}_{X0}{\mathsf{\Gamma }}_{11X0}\right)\dot{X}0,$$

(42)

Equation (42) can be integrated to shoot the autoparallel line starting from $X0$ with initial velocity $\dot{X}0$. In the general case, this integration is not simple because it involves the construction and inversion of the matrix $M_Y$ and alignment with the pole, P. Fortunately, the integration of (42) is not as complicated in the case of purely affine deformations. In the next Section 4.1, we show the analytical solution of geodesic and autoparallel lines for the case of purely affine deformations. Then, in Section 4.2 and Section 4.3, we exploit the qualitative definitions of Direct Transport and TPS metrics in 2D to numerically build the geodesics and autoparallel lines of the TPS space, respectively, in the general case. In particular, geodesics are built by minimising the length of the path connecting two given points (initial and final shapes), while autoparallel lines are calculated via shooting from a point (initial configuration) and a vector (initial deformation velocity).

4.1. Analytical Solution for the Affine Subspace

A trajectory of affine transformations of $X_0$ can be represented as

$$X\left(t\right)=X_0{A}_{X_0}^T\left(t\right)$$

(43)

with $A_{X_0}\left(t\right)\in SL\left(m\right)\forall t$ and $A\left(0\right)=I$

$$X_0=X\left(t\right)A_{X_0}^{-T}\left(t\right)$$

(44)

$$\dot{X}\left(t\right)=X_0{\dot{A}}_{X_0}^T\left(t\right)=X\left(t\right)A_{X_0}^{-T}\left(t\right){\dot{A}}_{X_0}^T\left(t\right)$$

(45)

$$\dot{X}\left(0\right)=X_0{A}_{X_0}^{-T}\left(0\right){\dot{A}}_{X_0}^T\left(0\right)=X_0{\dot{A}}_{X_0}^T\left(0\right)$$

(46)

The trajectory is an autoparallel line if and only if

$$\dot{X}\left(t\right)={\tau }_{X\left(t\right),X_0}\dot{X}\left(0\right)\forall t\in [0,1]$$

(47)

that is, by means of the (39)

$$\dot{X}\left(t\right)=X\left(t\right)A_{X_0}^{-T}\left(t\right){\dot{A}}_{X_0}^T\left(t\right)=X\left(t\right){\dot{A}}_{X_0}^T\left(0\right)$$

(48)

that is,

$${\dot{A}}_{X_0}\left(t\right)A_{X_0}^{-1}\left(t\right)={\dot{A}}_{X_0}\left(0\right)$$

(49)

we note that this equation is the same as that characterizing the autoparallel lines in $GL\left(m\right)$. The solution is as follows:

$$X\left(t\right)=X_{0}exp{\left[t{\dot{A}}_{X_0}\left(0\right)\right]}^{T}t\in [0,1]$$

(50)

as is well known, in $GL\left(m\right)$ geodesics and autoparallel lines coincide [15,19]; then, this property holds for TPS-geodesics and TPS-autoparallel lines as well in the case of affine deformations. The geodesic from $X_0$ to $X_1$ can be calculated as

$$X\left(t\right)=X_{0}exp\left[tlog\left({\mathsf{\Gamma }}_{21X_0}{X}_1\right)\right]t\in [0,1]$$

(51)

4.2. Geodesics Calculation: Objective Function Optimisation and Equality Constraints

Given two configurations $X_0,X_1$, we calculate the geodesics from $X_0$ to $X_1$ by exploiting the property of minimizing the Riemannian distance. The geodesic trajectory $\left\{X\left(t\right)\right|t\in [0,1],X\left(0\right)=X_0,X\left(1\right)=X_1\}$ can be calculated by minimizing the functional:

$$d(X_0,X_f)={\int }_0^1\sqrt{g(\dot{X}\left(t\right),\dot{X}\left(t\right))}dt$$

(52)

In addition, we require that the curve has a constant speed, ${\left(g(\dot{X},\dot{X})\right)}^{·}=0$, and we enforce this constraint by adding a Lagrange multiplier, k, to the objective function:

$$f\left(X\left(t\right)\right)=d(X_0,X_f)+k{\int }_0^1{\left(g(\dot{X},\dot{X})\right)}^{·}dt=0$$

(53)

The optimization problem is solved numerically by discretizing $X\left(t\right)$ in a finite number n of steps $X_i$ using Algorithm 1, sketched below.

Algorithm 1: geodesic algorithm.

Result: Geodesic path with initial configuration X0 and final configuration Xf using n discretization steps.

Then, the objective function is minimized by the R optimizer Solnp, R package version 1.16 [20]. The solver is an indirect solver implementing the augmented Lagrange multiplier method with an SQP interior algorithm.

4.3. Autoparallel Lines Calculation Algorithm via Shooting

Autoparallel lines do not minimize any distance, and are built directly by means of the parallel transport rule (41) via shooting. Let $Xo\in \mathcal{M}$ be the initial configuration and $V\in T\mathcal{M}$ be the initial deformation velocity; V can be called shooting vector and the path $X\left(t\right)$ such that $X\left(0\right)=X_0,\dot{X}\left(0\right)=V,{\nabla }_{\dot{X}}\dot{X}=0\forall t\in [0,1]$ is called shooting path of $Xo$ and V. In order to interpolate between $X_0$ and $X_1$ with an autoparallel line, shooting should be used iteratively to find a shooting vector $\overline{V}$ such that a shooting path starting at $X_0$ and with $\overline{V}$ as shooting vector reaches $X_1$ in a unit time (see [7]). In the present work, we are interested in comparing the behaviour of geodesics and autoparallel lines; thus, we avoid the iterative procedure by limiting ourselves to implementing a single shooting procedure starting from a configuration $X_o$ and a deformation $V_o$ and then using Algorithm 2, sketched below.

Algorithm 2: shooting algorithm.

Result: Shooting path with initial position X0 and initial velocity V0 using n discretization steps.

5. Examples

We propose five experiments aimed at finding TPS geodesics and comparing them with original shapes, shapes inputted in the optimizer, and shapes from geodesic shooting. For each experiment, eight shapes were generated.

• Affine case, spherical: in this simple case, the starting rectangle experiences only a size increase.
• Affine case, general: in this case, the rectangle undergoes only a pure affine transformation.
• Non-affine case bending: in this case, the rectangle experiences pure bending parameterized according to the parameters specified below.
• Non-affine case bending+size: in this case, the rectangle experiences pure bending parameterized according to the parameters specified below, with the addition of a size increase.
• Non-affine case bending+general affine component: in this case, the rectangle experiences pure bending parameterized according to the parameters specified below, with the addition of the same parameters of affine transformation as in case 2.

5.1. Dataset

A set of parametric shape paths were generated, starting from a rectangle $1\times 3$, by means of the following formula:

$$\left(\begin{array}{c}x\left(t\right)\\ y\left(t\right)\end{array}\right)=\left(\begin{array}{cc}{F}_{11}\left(t\right)& F_{12}\left(t\right)\\ F_{21}\left(t\right)& F_{22}\left(t\right)\end{array}\right)\left[{\displaystyle \frac{1+\chi \left(t\right)x_o}{\chi \left(t\right)}}\left(\begin{array}{c}sin\left(\chi \left(t\right)y_o\right)\\ cos\left(\chi \left(t\right)y_o\right)-1\end{array}\right)\right].$$

(54)

where $F_{11}\left(t\right),F_{22}\left(t\right),F_{12}\left(t\right),F_{21}\left(t\right)$ parametrize the affine transformations and $\chi \left(t\right)$ is the amount of bending. Each experiment is articulated in the following steps:

• A parametric trajectory $X\left(t\right)$ is generated by the mean of the (54). In this way, the initial and final points, $X\left(0\right)$ and $X\left(1\right)$, are identified.
• Linear interpolation between $X\left(0\right)$ and $X\left(1\right)$.
• The geodesic $Y\left(t\right)$, such that $Y\left(0\right)=X\left(0\right)$ and $Y\left(1\right)=X\left(1\right)$ are calculated following the procedure sketched in Section 4.2. The distance, $d\left(X\right(0),X(1\left)\right)$, and the initial tangent, $Y^{\prime }\left(0\right)$, are calculated.
• The autoparallel line $Z\left(t\right)$, starting from $X\left(0\right)$, is built by shooting $Z^{\prime }\left(0\right)=Y^{\prime }\left(0\right)$ by means of Direct Transport for a distance of $\ell =d\left(X\right(0),X(1\left)\right)$, following the procedure sketched in Section 4.3.

For all cases, we computed the linear interpolation between the first shape (the undeformed rectangle in all cases) and the last shape of any experiment. These linearized shapes (eight shapes) were the input for the optimizer, imposing as equality constraints the maintenance of the first and last shapes in order to force the geodesic to pass between these shapes. Finally, we performed a common Principal Component Analysis (PCA) for each experiment except for the first (which was trivial in terms of pure shape change), including all of the four sets of shapes: original, linearized, optimized, and shooted.

In each one of the experiments, we checked:

• The trend of the $\mathsf{\Gamma }$-energy.
• The trend of the components of the $\mathsf{\Gamma }$-energy.

5.2. Affine Case: Spherical

A parametric trajectory of eight configurations was generated by applying (54) to the initial rectangular shape by selecting eight values of the parameters in the following ranges: $F_{11}\left(t\right)=F_{22}\left(t\right)\in \left[1,2\right];F_{12}\left(t\right)=F_{21}\left(t\right)=0;\chi \left(t\right)=0.$ Figure 2 and

Table 1. Values of the energies for the affine spherical case; ${\mathbf{d}}_{\mathbf{u}}$ represents the affine component of the $\mathsf{\Gamma }$-energy, ${\mathbf{d}}_{\mathbf{b}}$ the non-affine component and ${\mathbf{d}}_{\mathbf{tot}}$ the total $\mathsf{\Gamma }$-energy.

	LINEAR			GEODESIC			A-PARALLEL
${\mathbf{d}}_{\mathbf{u}}$	${\mathbf{d}}_{\mathbf{b}}$	${\mathbf{d}}_{\mathbf{tot}}$	${\mathbf{d}}_{\mathbf{u}}$	${\mathbf{d}}_{\mathbf{b}}$	${\mathbf{d}}_{\mathbf{tot}}$	${\mathbf{d}}_{\mathbf{u}}$	${\mathbf{d}}_{\mathbf{b}}$	${\mathbf{d}}_{\mathbf{tot}}$
0.082	0.000	0.082	0.029	0.000	0.029	0.029	0.000	0.029
0.049	0.000	0.049	0.029	0.000	0.029	0.029	0.000	0.029
0.033	0.000	0.033	0.029	0.000	0.029	0.029	0.000	0.029
0.024	0.000	0.024	0.029	0.000	0.029	0.029	0.000	0.029
0.018	0.000	0.018	0.029	0.000	0.029	0.029	0.000	0.029
0.014	0.000	0.014	0.029	0.000	0.029	0.029	0.000	0.029
0.011	0.000	0.011	0.030	0.000	0.030	0.029	0.000	0.029

show the results for the size-only case. Geodesic searching via optimization satisfactorily recovers both the size change and the equally spaced $\mathsf{\Gamma }$-energy steps between consecutive shapes. The non-affine component of the $\mathsf{\Gamma }$-energy, $d_b$, is of course equal to zero. Optimized geodesics and shooting quietly coincide. Figure 3 shows scatterplots of Table 1 values.

5.3. Affine Case: General Case

A parametric trajectory of eight configurations is generated by applying (54) to the initial rectangular shape by selecting eight values of the parameters in the following ranges:

$$\begin{array}{ccc}\hfill F_{11}\left(t\right)& \in & \left[-0.2,1.2\right];F_{12}\left(t\right)\in \left[-0.2,1.2\right];\hfill \\ \hfill F_{21}\left(t\right)& \in & \left[-0.2,1.2\right];F_{22}\left(t\right)\in \left[0,0.6\right];\hfill \\ \hfill \chi \left(t\right)& =& 0.\hfill \end{array}$$

Figure 4 and

Table 2. Values of the energies for the affine general case; ${\mathbf{d}}_{\mathbf{u}}$ represents the affine component of the $\mathsf{\Gamma }$-energy, ${\mathbf{d}}_{\mathbf{b}}$ the non-affine component, and ${\mathbf{d}}_{\mathbf{tot}}$ the total $\mathsf{\Gamma }$-energy.

	LINEAR			GEODESIC			A-PARALLEL
${\mathbf{d}}_{\mathbf{u}}$	${\mathbf{d}}_{\mathbf{b}}$	${\mathbf{d}}_{\mathbf{tot}}$	${\mathbf{d}}_{\mathbf{u}}$	${\mathbf{d}}_{\mathbf{b}}$	${\mathbf{d}}_{\mathbf{tot}}$	${\mathbf{d}}_{\mathbf{u}}$	${\mathbf{d}}_{\mathbf{b}}$	${\mathbf{d}}_{\mathbf{tot}}$
0.020	0.000	0.020	0.010	0.000	0.010	0.010	0.000	0.010
0.015	0.000	0.015	0.010	0.000	0.010	0.010	0.000	0.010
0.011	0.000	0.011	0.010	0.000	0.010	0.010	0.000	0.010
0.009	0.000	0.009	0.010	0.000	0.010	0.010	0.000	0.010
0.007	0.000	0.007	0.010	0.000	0.010	0.010	0.000	0.010
0.006	0.000	0.006	0.010	0.000	0.010	0.010	0.000	0.010
0.005	0.000	0.005	0.010	0.000	0.010	0.010	0.000	0.010

show the results of the general affine-only case. Optimized geodesics correctly recover the original parameterized deformation, with equally spaced gamma-energy steps between consecutive shapes and db equal to zero. Figure 5 shows scatterplots of Table 2 values. Figure 6 shows the first two PCs resulting from PCA performed on all four types of datasets of the general affine case (parameterized, linearized, optimized geodesic, shooting).

5.4. Non Affine Case: Bending

A parametric trajectory of eight configurations is generated by applying (54) to the initial rectangular shape by selecting eight values of the parameters in the following ranges:

$$\begin{array}{ccc}\hfill F_{11}\left(t\right)& =& 0;F_{12}\left(t\right)=0;\hfill \\ \hfill F_{21}\left(t\right)& =& 0;F_{22}\left(t\right)=0;\hfill \\ \hfill \chi \left(t\right)& \in & \left[0,1.5\right].\hfill \end{array}$$

The results of the general non-affine case are shown in Figure 7, Figure 8 and Figure 9 and

Table 3. Values of the energies for the non-affine case of bending: ${\mathbf{d}}_{\mathbf{u}}$ represents the affine component of the $\mathsf{\Gamma }$-energy, ${\mathbf{d}}_{\mathbf{b}}$ the non-affine component, and ${\mathbf{d}}_{\mathbf{tot}}$ the total $\mathsf{\Gamma }$-energy.

	LINEAR			GEODESIC			A-PARALLEL
${\mathbf{d}}_{\mathbf{u}}$	${\mathbf{d}}_{\mathbf{b}}$	${\mathbf{d}}_{\mathbf{tot}}$	${\mathbf{d}}_{\mathbf{u}}$	${\mathbf{d}}_{\mathbf{b}}$	${\mathbf{d}}_{\mathbf{tot}}$	${\mathbf{d}}_{\mathbf{u}}$	${\mathbf{d}}_{\mathbf{b}}$	${\mathbf{d}}_{\mathbf{tot}}$
0.033	0.042	0.075	0.005	0.072	0.077	0.005	0.106	0.111
0.034	0.070	0.104	0.003	0.072	0.075	0.005	0.106	0.111
0.033	0.215	0.248	0.002	0.071	0.073	0.005	0.107	0.112
0.001	0.084	0.085	0.001	0.068	0.069	0.005	0.110	0.115
0.006	0.084	0.090	0.001	0.062	0.063	0.005	0.114	0.119
0.020	0.195	0.214	0.002	0.059	0.061	0.005	0.115	0.120
0.008	0.129	0.137	0.003	0.061	0.064	0.005	0.114	0.119

. This simulation is particularly challenging due to the particular non-affine transformation experienced by the rectangle. Despite this, the geodesic optimization finds shapes characterized by approximately equally spaced steps in terms of both gamma energy and its two components, du and db. The PCA scatterplot in Figure 9 shows coherent behaviour for all datasets except, as expected, for the linearized shapes.

5.5. Non-Affine Case: Bending and Scaling

A parametric trajectory of eight configurations is generated by applying (54) to the initial rectangular shape by selecting eight values of the parameters in the following ranges:

$$\begin{array}{ccc}\hfill F_{11}\left(t\right)& \in & \left[1,2\right];F_{12}\left(t\right)=0;\hfill \\ \hfill F_{21}\left(t\right)& =& 0;F_{22}\left(t\right)\in \left[1,2\right];\hfill \\ \hfill \chi \left(t\right)& \in & \left[0,1.5\right].\hfill \end{array}$$

Adding a significant size change to the previous experiment led to the results shown in Figure 10, Figure 11 and Figure 12 and

Table 4. Values of the energies for the non-affine case of bending and scaling; ${\mathbf{d}}_{\mathbf{u}}$ represents the affine component of the $\mathsf{\Gamma }$-energy, ${\mathbf{d}}_{\mathbf{b}}$ the non-affine component, and ${\mathbf{d}}_{\mathbf{tot}}$ the total $\mathsf{\Gamma }$-energy.

	LINEAR			GEODESIC			A-PARALLEL
${\mathbf{d}}_{\mathbf{u}}$	${\mathbf{d}}_{\mathbf{b}}$	${\mathbf{d}}_{\mathbf{tot}}$	${\mathbf{d}}_{\mathbf{u}}$	${\mathbf{d}}_{\mathbf{b}}$	${\mathbf{d}}_{\mathbf{tot}}$	${\mathbf{d}}_{\mathbf{u}}$	${\mathbf{d}}_{\mathbf{b}}$	${\mathbf{d}}_{\mathbf{tot}}$
0.060	0.37	0.431	0.017	0.120	0.137	0.017	0.150	0.167
0.003	0.41	0.412	0.003	0.128	0.131	0.017	0.150	0.167
0.147	0.42	0.571	0.007	0.115	0.122	0.017	0.150	0.168
0.028	0.04	0.073	0.021	0.101	0.122	0.017	0.159	0.177
0.021	0.03	0.056	0.029	0.087	0.116	0.017	0.155	0.172
0.036	0.07	0.104	0.025	0.091	0.115	0.017	0.176	0.193
0.012	0.03	0.040	0.044	0.065	0.110	0.017	0.175	0.193

. The equal spacing of the $\mathsf{\Gamma }$-energy of the optimized geodesics is rather acceptable, while its behaviour in the PCA space behaves more coherently than that of the parametrized or shooted shapes.

5.6. Non-Affine Case: Bending and General Affine Component

A parametric trajectory of eight configurations is generated by applying (54) to the initial rectangular shape by selecting eight values of the parameters in the following ranges:

$$\begin{array}{ccc}\hfill F_{11}\left(t\right)& \in & \left[-0.2,1.2\right];F_{12}\left(t\right)\in \left[-0.2,1.2\right];\hfill \\ \hfill F_{21}\left(t\right)& \in & \left[-0.2,1.2\right];F_{22}\left(t\right)\in \left[0,0.6\right];\hfill \\ \hfill \chi \left(t\right)& \in & \left[0,1.5\right].\hfill \end{array}$$

The last experiment is represented by a combination of non-affine and affine components. Results relative to this deformation are shown in Figure 13, Figure 14 and Figure 15 and

Table 5. Values of the energies for the non-affine case of bending and general affine deformation; ${\mathbf{d}}_{\mathbf{u}}$ represents the affine component of the $\mathsf{\Gamma }$-energy, ${\mathbf{d}}_{\mathbf{b}}$ the non-affine component, and ${\mathbf{d}}_{\mathbf{tot}}$ the total $\mathsf{\Gamma }$-energy.

	LINEAR			GEODESIC			A-PARALLEL
${\mathbf{d}}_{\mathbf{u}}$	${\mathbf{d}}_{\mathbf{b}}$	${\mathbf{d}}_{\mathbf{tot}}$	${\mathbf{d}}_{\mathbf{u}}$	${\mathbf{d}}_{\mathbf{b}}$	${\mathbf{d}}_{\mathbf{tot}}$	${\mathbf{d}}_{\mathbf{u}}$	${\mathbf{d}}_{\mathbf{b}}$	${\mathbf{d}}_{\mathbf{tot}}$
0.037	0.196	0.233	0.004	0.096	0.100	0.004	0.110	0.114
0.017	0.453	0.470	0.005	0.097	0.101	0.004	0.110	0.114
0.039	0.689	0.729	0.005	0.100	0.104	0.004	0.113	0.117
0.010	0.158	0.168	0.001	0.108	0.110	0.004	0.120	0.123
0.012	0.152	0.164	0.013	0.101	0.114	0.004	0.114	0.117
0.025	0.289	0.315	0.028	0.089	0.117	0.004	0.115	0.119
0.009	0.124	0.132	0.017	0.103	0.120	0.004	0.128	0.132

. Optimized geodesics struggle to find proper shape at fourth and fifth step of the deformation series; in the end, however, the final series results behave consistently in terms of both equal gamma energy spacing and general morphology.

6. Discussion

The set of performed numerical analyses, together with the analytical results of Section 4.1, showed that: The TPS geodesics are able to catch important qualitative behaviours of the parametric deformations. In particular, in each example the horizontal sides of the initial rectangle remain straight for the whole path while the vertical sides bends, exactly as happens in the parametric path. For affine deformations, TPS geodesics coincide with TPS autoparallel lines and both coincide with the LC geodesics of the group GL(m). For non-affine deformations TPS geodesics and TPS autoparallel lines do not coincide. In particular, the autoparallel lines are not very similar, expecially in the last steps, to the parametric path. This was expected for the procedure of shooting used here. On the other hand, certain qualitative behaviours are lost, as was expected. TPS Autoparallel lines conserve the percentage of affine and non-affine energies, while this does not happen in optimized geodesic. Future directions of the present work will involve the use of 3D data. This will certainly depend on the possibility of its being independent of the body-bending energy matrix computation. Three-dimensional data offer many practical applications, spanning a wide range of scientific disciplines such as cardiology [11,12], vertebrate paleontology [21,22], and paleoanthropology.