Symmetric Algorithmic Components for Shape Analysis with Diffeomorphisms

TL;DR

This paper introduces a symmetric space approach to shape analysis with diffeomorphisms, reducing numerical errors in registration by leveraging symmetries and involutions, and demonstrates improved consistency on cardiac shape data.

Contribution

It proposes a novel symmetric space framework for diffeomorphic shape analysis, enhancing numerical stability and accuracy in registration tasks.

Findings

Improved numerical consistency in shape registration.

Effective use of symmetries to reduce registration errors.

Successful application to cardiac shape data.

Abstract

In computational anatomy, the statistical analysis of temporal deformations and inter-subject variability relies on shape registration. However, the numerical integration and optimization required in diffeomorphic registration often lead to important numerical errors. In many cases, it is well known that the error can be drastically reduced in the presence of a symmetry. In this work, the leading idea is to approximate the space of deformations and images with a possibly non-metric symmetric space structure using an involution, with the aim to perform parallel transport. Through basic properties of symmetries, we investigate how the implementations of a midpoint and the involution compare with the ones of the Riemannian exponential and logarithm on diffeomorphisms and propose a modification of these maps using registration errors. This leads us to identify transvections, the composition of two symmetries, as a mean to measure how far from symmetric the underlying structure is. We test our method on a set of 138 cardiac shapes and demonstrate improved numerical consistency in the Pole Ladder scheme.

Figures

Tables

Table 1. Table 1 : Mean errors measured on cardiac shapes, in millimeters

	${\alpha }^2=0.01$		${\alpha }^2=1$		${\alpha }^2=1089$
Error Type	Residual	No Residual	Residual	No Residual	Residual	No Residual
Centrality	0.36	0.43	0.10	0.50		5.76
Involutivity	1.42	1.55	0.33	0.80		9.39
Transvectivity	1.98	2.16	0.58	0.62		0.16
Reg. Error	0.23		0.41		5.61
Reg. Norm	42		30		1
Inverse Cons.	0.13		0.14		0.10

Equations

$$<v,v^{\prime }{>}_H={\sum }_i{\sum }_{j}K(c_i,c_j^{\prime }){\mu }_i^T{\mu }_j^{\prime }.$$

$$\{\begin{array}{ll}{\dot{c_k}}^{(t)}& ={\sum }_{j}K(c_k^{(t)},c_j^{(t)}){\mu }_j^{(t)}\\ {\dot{\mu }}^{(t)}& =-{\sum }_j{\nabla }_{1}K(c_k^{(t)},c_j^{(t)}){\mu }_k^{(t{)}^T}{\mu }_j^{(t)}\end{array}$$

$$C(c,\mu )=\Vert S-{\varphi }_1^{c,\mu }(T){\Vert }_2^2+{\alpha }^2\cdot \Vert v_0^{c,\mu }{\Vert }_H^2.$$

$$u^{\prime }-u=\frac{1}{12}\left(({\nabla }_{v}R)(u,v)(5u-2v)+({\nabla }_{u}R)(u,v)(v-2u)\right)+O(\Vert v+u{\Vert }^5).$$

$$S={\varphi }_1(T)+\delta ,$$

$$M=\phi_{\frac{1}{2}}(T)+\frac{1}{2}\delta={\textrm{Exp}}_{T}\Big{(}\frac{1}{2}{\textrm{Log}}_{T}(S)\Big{)}+\frac{1}{2}\delta.$$

$$s_{T}(S)={\textrm{Exp}}_{T}\Big{(}-{\textrm{Log}}_{T}(S)\Big{)}-\delta.$$

$$\forall M,S\in \mathcal{V}{s}_M\circ s_S=s_T\circ s_M,\text{where}T=s_M(S).$$

$$u^{\prime \prime }-u=-\frac{1}{12}\left(({\nabla }_{v}R)(u,v)(5u+v)+({\nabla }_{u}R)(u,v)(v+2u)\right)+O(\Vert v+u{\Vert }^5).$$

$$u^{\prime }-u^{\prime \prime }=\frac{1}{12}\left(({\nabla }_{v}R)(u,v)(10u-v)-4({\nabla }_{u}R)(u,v)u\right)+O(\Vert v+u{\Vert }^5).$$

$$LA{S}_i^{S-S^{\prime }}=\frac{1}{k_i}\sum _{j=1}^{k_i}\frac{(a_j-a_j^{\prime })}{a_j}$$

$$AS{E}^2=\sum _i(LA{S}_i^{S-S^{\prime }}-LA{S}_i^{T-T^{\prime }}{)}^2$$

Full text

1 Introduction

Computational anatomy aims at modeling the temporal evolution and cross-sectional variability of anatomical shapes. The deformations between shapes are obtained by applying non-rigid registration algorithms that seek the smallest transformation - in a sense that will be defined precisely - of the ambient space to match two shapes. In the diffeomorphic registration setting, the deformations are modeled by diffeomorphisms, that provide invertible and folding-free transformations.

In the Large Deformation Diffeomorphic Metric Mapping (LDDMM) framework, the space of diffeomorphisms is endowed with a right invariant metric and deformations of interest are obtained by geodesic flows from the identity transformation. They are parameterized by their initial velocity fields, which are tangent vectors at the identity deformation, defined by initial control points and dual momenta. A transformation is then computed by integration of differential equations. Registration is performed by solving an optimization problem on the initial momenta and control points with a gradient descent. These numerical schemes efficiently implement an exponential and logarithm map on a subspace of diffeomorphisms.

However in practice, the optimization problem is relaxed to enforce smooth deformations and results in inexact matching. Therefore we can think of the exponential map as not going “far enough”. In this work we introduce a modified exponential map that accounts for a residual error due to registration, indifferent to the choice of the regularization parameter.

The choice of the metric and regularization parameter affects the geometric structure of the space of deformations under consideration. Many convergence results depend on the curvature of this space, e. g. [9], and especially on its covariant derivative. This gradient being null in locally symmetric spaces, they form a very convenient setting to perform statistics on shapes. In order to assess how far from symmetric our structure is, it would therefore be valuable to develop a procedure to measure this gradient. In this paper we build on a specific parallel transport scheme.

In the statistical analysis of temporal deformations, parallel transport along geodesics is commonly used to perform inter-subject normalization, that is the vector transport of velocity fields from each subject’s space to a common atlas’ space. An approximation based on Jacobi fields was proposed in [11, 7]. A numerical implementation named Pole Ladder (PL) was proposed in [5] and relies only on the computation of exponential and logarithm maps. Following the Shild’s ladder, it consists in the construction of geodesic parallelograms. The progression between two shapes $S$ and $S^{\prime}$ is transported to $T$ by:

• •

first computing a “midpoint” $M$ on the geodesic between $S$ and $T$. It is seen as the diagonal of the geodesic parallelogram;

• •

then extending the geodesic from $S^{\prime}$ to $M$ by the same length to obtain $T^{\prime\prime}$;

• •

similarly extending the geodesic from $T^{\prime\prime}$ to $T$ to obtain the parallel deformation of the template $T$;

See Fig. 1 for a schematic representation.

For large deformations that typically arise in inter-subject registration, this procedure is usually iterated. We therefore expect numerical errors to grow linearly in the number of steps, and lose crucial numerical accuracy. The accuracy of this scheme was analyzed in [9] and shown to be a third order scheme in general, and exact in symmetric spaces where it is equivalent to a single transvection, that is, a composition of two symmetries, in our case symmetries with respect to $M$ and $T$.

In many cases, it is well known that numerical errors can be drastically reduced in the presence of a symmetry. This is for instance the case to diagonalize a symmetric matrix versus an arbitrary one. Using the Stationary Velocity Fields (SVF) framework for registration, a symmetric variant of Pole Ladder built on a Lie Group intrinsic symmetric structure was proposed in [3]. This procedure is recalled in section 2.2. In this work, we build on this idea, but rely on LDDMM to implement a more general involution that accounts for the registration residual. This elementary algorithmic component constructs the symmetric shape of an original shape with respect to another one. It is presented in section 2.3. In section 2.4, we introduce the basic properties that symmetries must verify in an affine symmetric space and discuss whether these are fulfilled by our implementation. From a theoretical perspective, deviations from these properties are due to a non zero covariant derivative of the curvature tensor with the LDDMM metric. Conversely, we may interpret transvection errors as estimates of the numerical curvature gradient that encompass all the approximations due to the implementation. The numerical experiments of section 3 show that there is an optimal regularisation parameter for which the space of deformations can best be approximated by a symmetric space.

The paper is organised as follows: in section 2.1, we recall the LDDMM framework following [2] and in section 2.2 the Pole Ladder procedure. We then introduce in sections 2.3 and 2.4 the main contributions of the paper, namely accounting for residuals, defining symmetries and their properties: centrality, involutivity and transvectivity. In section 3 we present the numerical experiments and comment on the results.

2 Background and Method

2.1 The LDDMM Framework

In this work we consider shapes represented by 3D meshes. However, the methodology seamlessly applies to images. In order to define a practical finite dimensional parameterization of a subspace of diffeomorphisms $G$ acting on the ambient space ${\mathbb{R}}^{d}$, we consider time-varying velocity fields $v_{t}(x)=\sum_{k=1}^{N_{c}}K(x,c^{(t)}_{k})\mu^{(t)}_{k}$ obtained by convolution of a Gaussian kernel $K(x,y)=\exp(-\frac{\|x-y\|^{2}}{\sigma^{2}})$ over $N_{c}$ control points $c^{(t)}=[c_{k}^{(t)}]_{k}$ and momenta $\mu^{(t)}=[\mu_{k}^{(t)}]_{k}$.

The set of such fields forms a pre-Hilbert space with scalar product between $v=\sum_{k}K(\cdot,c_{k})\mu_{k}$ and $v^{\prime}=\sum_{k}K(\cdot,c_{k}^{\prime})\mu_{k}^{\prime}$ defined by

[TABLE]

Diffeomorphisms are then defined as flows of velocity fields from $\phi_{0}=Id$. This amounts to integrating the ordinary differential equation (ODE) $\partial\phi_{t}(\cdot)=v_{t}[\phi_{t}(\cdot)]$ between [math] and $1$. The scalar product on velocity fields induces a right-invariant Riemannian metric on the obtained subspace of diffeomorphisms, and geodesics of this metric are parameterized by control points and momenta that satisfy the following Hamiltonian equations:

[TABLE]

$\phi_{t}$ and $v_{t}$ are thus uniquely determined by their initial control points $\mu$, this dependence will be explicitly written $\phi_{t}^{c,\mu}$ and $v_{t}^{c,\mu}$. In practice the interval $[0,1]$ is discretized with $n$ time steps and the ODE is solved with an iterative Euler forward or Runge-Kutta 2 method. This implements an exponential map at identity. The registration problem between a template shape $T$ and a target mesh $S$ optimizes the following criterion over initial control points and momenta $c,\mu$:

[TABLE]

For simplicity, we measure the distance between shapes by the $L_{2}$ distance between nodes of the meshes. $\|\cdot\|_{H}$ is the norm defined by the scalar product of eq. 1, which is actually the metric on $G$. The resulting $v_{0}$ is a logarithm at identity of $\phi_{1}$. In fact the metric is scaled by a factor $\alpha$, and this impacts the geometry of the underlying space as will be demonstrated. It allows to smoothly interpolate between solutions that belong to two paradigms:

• •

$\alpha\rightarrow 0$: exact matching between shapes;

• •

$\alpha\rightarrow\infty$: point distribution models (PDM), no deformation.

2.2 Symmetric Pole Ladder for parallel transport

In the context of computational anatomy, the aforementioned registration framework is used to represent a subject-specific temporal deformation between shapes $S$ and $S^{\prime}$, and to transport this deformation to a common atlas or template $T$ along the geodesic segment $[S,T]$. The anatomical shapes are modeled as points in a manifold ${\mathcal{V}}$ under the action of the space of diffeomorphisms $G$ described above. We suppose here that this manifold is equipped with an affine connection, which defines parallel transport and the Exp map. Locally it further defines the Log map.

Algorithm 1 presents a symmetric variant of Pole Ladder [5] introduced in [3] to approximately perform parallel transport in ${\mathcal{V}}$.

A Taylor expansion at the midpoint $M$ of the error between the vector transported by Pole Ladder and exact parallel transport of the vector $u_{S}={\textrm{Log}}_{S}(S^{\prime})$ is derived in [9]. We denote by $u=\Pi_{S}^{M}u_{S}$ the exact transport to $M$, and $u^{\prime}=\Pi_{T}^{M}{\textrm{Log}}_{T}(T^{\prime})$ where $T^{\prime}$ is obtained by Pole Ladder. Let also $v={\textrm{Log}}_{M}(T)$. Then

[TABLE]

In fact this scheme is exact in an affine locally symmetric space, where $\nabla R=0$. In this case, using the local symmetries $s_{X}$ at point X, we have $T^{\prime}=s_{T}\circ s_{M}(S^{\prime})$ meaning that Pole Ladder is equivalent to a transvection. As local symmetries are affine mappings, the following diagrams commute:

{T_{S}{\mathcal{V}}}$${T_{T}{\mathcal{V}}}$${T_{M}{\mathcal{V}}}$${T_{M}{\mathcal{V}}}$$\scriptstyle{\Pi_{S}^{M}}$$\scriptstyle{(ds_{M})_{S}}$$\scriptstyle{\Pi_{T}^{M}}$$\scriptstyle{-Id}

{T_{S}{\mathcal{V}}}$${T_{T}{\mathcal{V}}}$${{\mathcal{V}}}$${{\mathcal{V}}}$$\scriptstyle{{\textrm{Exp}}_{S}}$$\scriptstyle{(ds_{M})_{S}}$$\scriptstyle{{\textrm{Exp}}_{T}}$$\scriptstyle{s_{M}}

Thus, $\Pi_{T}^{M}\circ(ds_{M})_{S}=-\Pi_{S}^{M}$ [10, Prop 4.3]. So with the previous notations, $(ds_{M})_{S}u_{S}=-\Pi_{M}^{T}u$ and $(ds_{M})_{S}u_{S}=Log_{T}(T")=-\Pi_{M}^{T}u^{\prime}$ yielding $u^{\prime}=u$ .

2.3 Accounting for residuals to improve Centrality and Symmetry

In practice the registration is never exact due to the regularization. Thus, we propose to decompose the space of shapes into a deformation part encoding the orbit of the template, and a Euclidean space of residual displacement fields

[TABLE]

where $\delta$ is a displacement field between corresponding points of the two meshes. Of course different possibilities exist to compute the residual, to transport it, and to apply it to different shapes, but our experiments suggest that this very simple formulation may be sufficient, so that we will not detail other approaches in this paper. With this decomposition, a midpoint between $T$ and $S$ is defined by (Fig. 2):

[TABLE]

Unfortunately this formulation is not symmetric in $T$ and $S$, and registering $S$ on $T$ and shooting from $S$ results in a different midpoint in general. We will see however that using residuals decreases the distance between the midpoints obtained with the two initial points.

Similarly, a symmetry is defined by inverting the geodesics and the residuals (Fig. 3):

[TABLE]

The first desired consistency property which we refer to as centrality is the compatibility of the midpoint with the symmetry, namely $s_{M}(S)=T$. Note that this construction of a central midpoint and a local geodesic symmetry is possible in arbitrary affine connection manifolds for close enough points. However, these symmetries are affine mappings if and only if the space is locally symmetric [10, Prop 4.2].

2.4 Involutivity and Transvectivity to measure curvature gradient

By construction our symmetry verifies $s_{T}(T)=T$ for all $T\in{\mathcal{V}}$, and if the $\log$ and $\exp$ maps are exact, $s_{T}\circ s_{T}=Id$. We will evaluate the exactitude of this property, called involutivity. Note that in a Lie Group of transformations, natural symmetries may be defined at $\phi$ by $\psi\mapsto\phi\circ\psi^{-1}\circ\phi$, which, at $\phi=Id$ is in fact the inversion. Involutivity in this case reduces to inverse consistency: $(\psi^{-1})^{-1}=\psi\iff\psi\circ\psi^{-1}=Id$. However, in our framework, both types of errors are different because the metric exponential differs from the one defined by the canonical Cartan-Shouten connection, which is the only one compatible with the group operations.

Finally, in an affine globally symmetric space, symmetries must verify the following property that we will call transvectivity [10, Prop 5.3]:

[TABLE]

We want to evaluate the exactitude of this property, and use the deviation to this ideal case as a proxy to measure the gradient of the curvature of the space in the directions of interest. At this point the role of $\alpha$ becomes clearer, it allows to form a continuum of decompositions of the space of diffeomorphisms under consideration:

• •

$\alpha\rightarrow 0,\;\delta\rightarrow 0$: the Riemannian space of deformations where the metric is not compatible with the Cartan-Shouten connection. This space is not symmetric, which shows in the transvectivity error.

• •

$\alpha\rightarrow\infty,\;v_{0}\rightarrow 0$: the shapes are considered in the ambient space with Euclidean norm, this space is of course symmetric. This is the PDM framework.

We saw in section 2.2 that Pole Ladder was doing the transvection $s_{T}\circ s_{M}$. With the same notations and using a Taylor expansion from [9], the deviation to parallel transport when applying the opposite transvection $s_{M}\circ s_{S}$ is:

[TABLE]

Thus when measuring the transvectivity error $\|s_{M}\circ s_{S}(S^{\prime})-s_{T}\circ s_{M}(S^{\prime})\|_{2}$, we in fact measure $\|{\textrm{Exp}}_{T}(\Pi_{M}^{T}u")-{\textrm{Exp}}_{T}(\Pi_{M}^{T}u^{\prime})\|_{2}$ where

[TABLE]

The transvectivity error thus provides a practical way to measure the gradient of the curvature of the space even in the absence of any closed-form expression. This is noticeable in regards to the complexity of the curvature tensor itself in Mario’s formula [8] and it may lead in the future to new ways of estimating the curvature and its gradient.

3 Experiments and Application to Cardiac Shapes

In this section we assess the consistency of our numerical implementation of the symmetry compared to its theoretical properties. We compare the symmetry with residuals to standard symmetry without residuals ($\delta=0$) for different values of the parameter $\alpha$. We used Deformetrica [1] for all our experiments. We also compare the parallel transport obtained by Pole Ladder with both types of symmetry to the one implemented in Deformetrica [6] using the fanning scheme.

We use a database of cardiac shapes from 138 subjects [3] for which the shapes at two time-points are available: at end-diastole ($S$) and at end-systole ($S^{\prime}$). We use a population atlas as template $T$. Four types of errors are first measured. The first is the distance between midpoints when shooting from $T$ or from $S$. The three others are, where $M$ is the midpoint obtained by shooting from $T$:

• •

$\|s_{M}(T)-S\|_{2}$: the centrality error (Fig. 4(a))

• •

$\|s_{M}\circ s_{M}(S^{\prime})-S^{\prime}\|_{2}$: the involutivity error (Fig. 4(b))

• •

$\|s_{T}\circ s_{M}(S^{\prime})-s_{M}\circ s_{S}(S^{\prime})\|_{2}$: the transvectivity error (Fig. 4(c))

Mean results for extreme values of $\alpha$ are given in table 1. The average registration error of $T$ on each subject’s $S$, as well as the norm of the deformation and inverse consistency (by registering $S$ on $T$) for each regularisation parameter are also given for reference.

These results illustrate the two contributions of this paper. Firstly, using residuals in the symmetry considerably improves the numerical accuracy of the computation of a midpoint. As we can see on Fig. 5(a), the distance between midpoints computed by shooting from $T$ or from $S$ is reduced when using the residuals. This error compares well with the inverse consistency error in general and is even significantly lower for $\alpha\geq 1$.

Moreover, this increase in numerical accuracy is also visible on the centrality and involutivity errors. Indeed, for centrality (Fig. 5(b)), the error when using residuals is consistently smaller than without residuals, and this gain becomes larger as $\alpha$ grows. It is also remarkable that this error is significantly lower than the registration error for $\alpha\geq 1$. The same behavior is observed for the involutivity on Fig. 5(c). This means that for reasonable values of $\alpha$ we obtain reliable implementations of the midpoint and symmetry.

Secondly, the LDDMM space endowed with the right invariant metric is not symmetric, and the transvectivity error reflects the covariant derivative of the curvature. Fig. 5(d) gives further details: pushing registration with a small $\alpha$ generates larger errors, and there is an optimal $\alpha\in[1;2]$ for which the space is closest to being symmetric. Using residuals, this error decreases as the space flattens to a Euclidean space and deformations tend to the identity.

Finally, in order to evaluate the result of symmetric Pole Ladder and compare it to the fanning scheme, we compute the local area strain (LAS) between end-diastole and end-systole at every landmark of the mesh. For two corresponding points $m_{i},m_{i}^{\prime}$ that belong to $k_{i}$ triangular cells, we compute the mean of the difference of area of each of these cells between $S$ and $S^{\prime}$.

[TABLE]

This feature is commonly used by clinicians to characterise the cardiac motion [4]. Here we use it to test the isometric property of the parallel transport scheme: we measure the area strain between the original subject’s meshes $S$ and $S^{\prime}$, and compare it to the one measured between the atlas $T$ and the deformed atlas $T^{\prime}$ obtained by the parallel transport algorithm. We report in Fig. 6(a) the square root of the sum of squared differences over all landmarks:

[TABLE]

As $\alpha$ grows, both the temporal and subject-to-atlas deformations decrease, thus generating less area changes, which explains the growing errors. Furthermore, we can see on Fig. 6(b) that the area strain is dominated by a bending artifact of the valve sections at the borders. Using residuals reduces this effect in this example. Although more suitable methods may exist to directly map scalar functions from one shape to another, these results emphasize the contribution of this paper: using residuals improves the symmetry and parallel transport with Pole Ladder.

4 Conclusion

We introduced residuals from registration errors to compute midpoints and symmetries between shapes. This results in improved numerical consistency for the centrality and involutivity properties. Furthermore the transvectivity error reflects the curvature of the underlying Riemannian space of deformations, and allows to estimate how far from symmetric this space is, depending on the regularisation parameter. Performing the same experiments in the framework of SVF with a similar implementation would yield very interesting comparison as the space is naturally symmetric, and accounting for residuals would thus provide a more consistent method for parallel transport. This gain could be even more interesting when considering registration between images. Indeed images are not in the template’s orbit in practice and residuals would encode intensity bias, which is a key source of error in image registration.

Acknowledgements:

This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (grant agreement G-Statistics No 786854).