# Barycentric Subspace Analysis: a new Symmetric Group-wise Paradigm for Cardiac Motion Tracking

Marc-Michel Rohé, Maxime Sermesant, Xavier Pennec

► To cite this version:

Marc-Michel Rohé, Maxime Sermesant, Xavier Pennec. Barycentric Subspace Analysis: a new Symmetric Group-wise Paradigm for Cardiac Motion Tracking. MICCAI 2016 - Medical Image Computing and Computer Assisted Intervention, Oct 2016, Athens, Greece. pp.300-307, 10.1007/978-3-319-46726-9_35. hal-01373706

HAL Id: hal-01373706

https://inria.hal.science/hal-01373706

Submitted on 29 Sep 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Barycentric Subspace Analysis: a new Symmetric Group-wise Paradigm for Cardiac Motion Tracking

Marc-Michel Rohé, Maxime Sermesant and Xavier Pennec

Université Côte d'Azur, Inria, France, Asclepios Research Group

**Abstract.** In this paper, we propose a novel approach to study cardiac motion in 4D image sequences. Whereas traditional approaches rely on the registration of the whole sequence with respect to the first frame usually corresponding to the end-diastole (ED) image, we define a more generic basis using the barycentric subspace spanned by a number of references images of the sequence. These subspaces are implicitly defined as the locus of points which are weighted Karcher means of $k+1$ references images. We build such subspace on the cardiac motion images, to get a Barycentric Template that is no longer defined by a single image but parametrized by coefficients: the barycentric coordinates. We first show that the barycentric coordinates - the coefficients of the projection of the motion during a cardiac sequence - define a meaningful signature for group-wise analysis of dynamics and can efficiently separate two populations. Then, we use the barycentric template as a prior for regularization in cardiac motion tracking, efficiently reducing the error of tracking between end-systole and end-diastole by almost 40% as well as the error of the evaluation of the ejection fraction. Finally, to best exploit the fact that multiple reference images allow to reduce the registration displacement, we derived a symmetric and transitive registration that can be used both for frame-to-frame and frame-to-reference registration and further improves the accuracy of the registration.

## 1 Introduction

Understanding and analyzing the cardiac motion pattern in a patient is an important task in many clinical applications. It can give insight into a pathology, by evaluating for example how the cardiac function is affected by a cardiovascular disease and if a therapy is needed or not. On top of traditional simple parameters such as the ejection fraction (EF), it can also be used to compute more complex parameters - such as strains in different directions - giving deeper insight to the efficiency of the heart motion and function. The cardiac motion is usually studied by finding correspondences - the registration step - between each of the frame of the sequence and the first frame corresponding to the end-diastole (ED) image, yielding a dense displacement field that tracks the motion of the myocardium. Taking the ED image as a reference is natural as it is the starting point of the contraction of the heart which is the most important phase
in evaluating the efficiency of the cardiac function but this specific choice can
lead to important biases in quantifying the motion especialy at end-systole (ES)
where the deformations to be evaluated are large [1].

In this paper, we propose a novel approach to study cardiac motion. Instead of taking an unique image as the reference to evaluate the motion, we build affine subspaces on the manifold of deformations encoding cardiac motion. There are different ways to extend the concept of principal affine spaces from an Euclidian space to something defined on manifolds. The simplest generalization is tangeant PCA, where a covariance matrix is build on the tangeant space of the Karcher or Frechet mean. In Principle Geodesic Analysys (PGA) [2], subspaces are spanned by the geodesics going through a point and the tangent vector is restricted to belong to a linear space of the tangent space. In this paper, we use more general type of family of subspaces on manifolds called Barycentric Subspaces which were first introduced in [3]. With respect to the method previously mentioned, it has the benefit not to be controlled by the central value. This gives a more consistent framework to study data in the case the underlying distribution is either multimodal or simply not sufficiently centered.

In the context of deformation analysis in medical imaging, the points of
the manifold corresponds to 3D images whereas the geodesic are deformations
mapping two images together. Optimal paths (geodesics) are represented by
the initial velocity field at the geodesic path resulting from the registration of
images. In the first part of this article, we define the barycentric subspaces of
manifold and introduce the way to compute the barycentric coefficients and the
projection of an image on a Barycentric Subspace of dimension k based on k + 1
images. Instead of performing registration with respect to a single template, we
build a subspace based on multiple references images and take advantage of the
information of group-wise registration [4], by building a Barycentric Template
of dimension 2 parametrized by the barycentric coefficients. Experiments are
conducted on sequences of healthy and pathological patients and show that the
barycentric coefficients of both populations present significant differences and
two clear clusters appear. Then, we improve the registration of cardiac motion
by relaxing the regularization within the 2-dimensional barycentric template
representing a cardiac sequence. Finally, we further improve the methodology
by deriving a formula leading to symmetric and transitive registration.

## 2 Barycentric Subspaces in Deformation Manifolds

In this section, we introduce barycentric subspaces following the notation de-
scribed in [3]. In order to adapt the framework from Riemanian Manifolds to
the context of computational anatomy (image deformation analysis), we follows
the framework of [5]. Working in the space of images M, we define I as a point
of this space, which can be for example an image of a cardiac sequence and
we identify paths to deformations. In the following, we will use (k + 1) points
R<sub>j</sub>, the *references images*, on this Manifold as well as (k + 1) coefficients λ<sub>j</sub>
the *barycentric coefficients*. The *Barycentric Subspace* of dimension k spanned
<figure><img src="image_2.png" /><figcaption>Fig. 1: Barycentric subspace of dimension 2 built from 3 references images ($R_1, R_2, R_3$). $\hat{I}$ is the projection of the image $I$ within the barycentric subspace such that $\| \hat{v} \|^2$ is minimum under the conditions $\sum_j \lambda_j \hat{v}_j = 0$ and $\hat{v} + \hat{v}_j = v_j$.</figcaption></figure>

by these points is then defined as the set of points (images) $\hat{I}$ in $\mathcal{M}$ such that: $\sum_{j=1}^{k+1} \lambda_j \hat{I} \hat{R}_j = 0$, where $\hat{I} \hat{R}_j = \log_I(R_j)$ is the smallest velocity field that registers $I$ to $R_j$. Contrary to the Riemannian setting where we would have exactly $R_j = \exp_I(\hat{I} \hat{R}_j)$, we obtain through registration an inexact matching that approximates the log vector. This is the tangeant vector of the geodesic shooting $I$ to $R_j$. In the following, we will place ourselves in stationary velocity fields (SVF) framework [6] which gives a simple and yet effective way to parametrize smooth deformations along geodesics using one-parameter sub-group. In this case, the tangent vector $\hat{I} \hat{R}_j$ will be parametrized by the SVF $\hat{v}_j$ and the condition simply becomes $\sum_{j=1}^{k+1} \lambda_j \hat{v}_j = 0$. The notation are summed up in Fig. 1.

## 2.1 Projection on Barycentric Subspace

Having defined the barycentric subspace spanned by a set of $k + 1$ references $R_j$, we are looking to find the projection $\hat{I}$ of any image $I$ in $\mathcal{M}$ on this subspace together with the coefficient $\lambda_j$ representing the coordinates of $\hat{I}$ within the barycentric template. The projection $\hat{I}$ of $I$ is the closest point to $I$ that belongs to the barycentric subspace. We define the SVF $\hat{v}$ which parametrizes the projection of $I$ such that $\hat{I}\hat{I} = \hat{v}$ as well as the SVFs $(v_i)_{i=1,...,k+1}$ such that $\overrightarrow{R_i}\hat{I} = v_i$ as shown in Fig. 1. The distance between $I$ and $\hat{I}$ is represented by the norm of the SVF $\| \hat{v} \|^2$. As seen previously, the constraint that $\hat{I}$ belongs to the
<figure><img src="image_3.png" /><figcaption>Fig. 2: (Left): Optimal references for the two group for each of the 3 references. (Right): Barycentric coefficient curves λ projected on the 2D-plan $\sum_i \lambda_i = 1$.</figcaption></figure>

barycentric subspace can be written as $\sum_j \lambda_j \hat{v}_j = 0$. Using the Baker-Campbell-Hausdorff (BCH) [6] formula, we get a first order development of $v_i = \hat{v} + \hat{v}_i$. The problem can now be written as:

$$ \min_{\hat{v}} \| \hat{v} \|^2, \quad \text{subject to} \quad \sum_i \lambda_i (v_i - \hat{v}) = 0, \quad \sum_i \lambda_i \neq 0 $$

whose set of solutions is:

$$ \lambda = \alpha A^{-1} 1, \quad \alpha \in \mathbb{R}^*, \quad A = \langle v_i | v_j \rangle_{i,j}. $$

With the additional constraint $\sum_j \lambda_j = 1$, we get the unique solution by normalizing the λ.

Finally, given a set of N images $I_n$, we can define the distance of a barycentric subspace spanned by the references $(R_j)_{j=1,...,M}$ to the set of images $(I_n)_{n=1,...,N}$ as the sum of the error over all the residual projection vector $\hat{v}_n$ such that:

$$ \mathcal{E}((R_j)_{j=1,...,M}) = \sum_{n=1,...,N} \| \hat{v}_n \|^2, \quad (1) $$

and find the closest subspace with respect to the cardiac sequence by minimizing this distance over all set of references.

## 2.2 Cardiac Motion Signature for Group-Wise Analysis of Dynamics

We applied the previously defined methodology to compare the cardiac motion signature of two different populations. The first group consists of 15 healthy subjects from the STACOM 2011 cardiac motion tracking challenge dataset [1]: an openly available data-set, and the second group is made of 10 Tetralogy of Fallot (ToF) patients. Short axis cine MRI sequences were acquired with T = 15 to 30 frames. The methodology described was applied by projecting each of
the *T* frame of the cardiac motion to a barycentric subspace of dimension 2 spanned by 3 references. This set of 3 references is chosen by building the optimal barycentric subspace as induced by the distance defined in Eq. 1. Significant differences in the frame for the optimal references can be seen between the two populations (Fig. 2, left). In particular, the second reference - corresponding to the end-systole - is significantly higher for the ToF patients showing that this population has on average longer systolic contraction. Then we project the whole cardiac motion on the barycentric subspace made of these three references to compute the barycentric coefficients (Fig. 2, right). We see significant differences between both groups of curves, especially in the region $\lambda_1 < 0$ showing that this signature of the motion is encoding relevant features and classical machine learning algorithms can separate the two populations using this representation.

# 3 Registration Using Barycentric Subspaces

In this section, we show how the use of barycentric subspaces as a prior on the cardiac motion can improve the registration by relaxing the regularization. Most of the registration methods rely on a trade-off between a fidelity-to-data term - capturing how well the registration matches the intensities of the voxels of the images - and a smoothness regularization term - encoding our prior information about the regularity of the deformations we are looking for. It is standard practice in registration algorithms to consider slowly varying deformation as our prior knowledge of the transformation (either by constraining the deformation to be within a small subspace of all diffeomorphisms or by penalizing large deformations). While this methodology works well to find small deformations, the regularization often leads to an underestimation of the large deformations as the one happening between the ED and ES frame. To overcome this drawback, solutions usually rely on performing the registration in a group-wise manner: a group of images are considered simultaneously and an additional criteria is added to ensure temporal-consistency [7] [8]. In this paper, we propose to use the barycentric template defined by 3 frames of the sequence as an additional prior on the transformations by considering that only the distance to the closest image within this barycentric subspace should be minimized. In the regularization step, we no longer consider the whole velocity field *v* but we run the regularization only with respect to $\hat{v}$ which encodes the distance of the current image to the closest image within the barycentric template representing the cardiac motion.

## 3.1 Barycentric Log-Demons Algorithm

We apply this methodology to the *Symmetric Log-Domain Diffeomorphic Demons* algorithm [6] which successively updates the velocity field to match the data, then smooths the velocity field with a gaussian filter. Instead of performing the regularization on the complete velocity field (the "standard" method), we decompose the velocity field $v_i$ as the sum of $\hat{v}_i$ mapping the reference $R_i$ to the projection $\hat{I}$ inside the template and the residual velocity field $\hat{v}$ of the projection
<figure><img src="image_4.png" /><figcaption>Fig. 3: (Left): Image of heart at end-systole, the contour of the warped mesh of the initial frame 1 (ED) using barycentric (red), standard (blue) registrations are shown together with the ground truth (green). Using barycentric regularization, the registration is less constrained and we manage to get a more accurate contours for the end-systolic endocardium. (Middle): Error of the registration with respect to 3 reference images using the two methods, barycentric (dotted lines) and standard (plain lines). (Right): volume curves induced by the registration and comparison with the ground truth volume.</figcaption></figure>

<figure><img src="image_5.png" /></figure>

(see Fig. 1) and we regularize only the residual $\hat{v}$ with the gaussian filter. The barycentric template is therefore used as a prior on the cardiac motion for which we do not perform regularization. The method was evaluated using a synthetic time serie of $T = 30$ cardiac image frames [9], so that we have ground truth meshes along the sequence allowing us to estimate the accuracy of the registration. First, we find the optimal references by minimizing the energy in Eq. 1 giving us the frames 1, 11 and 21 which will be the three references spanning the barycentric subspace. Then we register each frame $i$ of the sequence using the method described above to get the deformations from each of the three references to the current images. We deform the ground truth meshes at the references frames with these deformations and compare the results with classical registration. As can be seen in Fig. 3, barycentric registration performs better at catching the end-systolic deformation with the contour of the warped mesh at end-systole matching better the ground truth. The estimation of the ejection fraction from the volume curve is also improved, going from 32% with the standard method to 38%, closer to the ground truth (43%), reducing the estimation error by half. Finally, the average point-to-point error for both methods shows (in Fig. 4, right) that, while barycentric registration has largest error for small deformations close to each reference, it has around 30% smaller error for largest deformations as between ED and ES.
<figure><img src="image_6.png" /><figcaption>Fig. 4: (Left): Schematic representation of the symmetric multi-references barycentric registration in the case of a 1-D barycentric subspace spanned by 2 references. $I_t$ and $I_s$ are two frames of the sequence and $\hat{I}_t, \hat{I}_s$ corresponds to the respective projection to the barycentric subspace. (Right): comparison of the error between the standard registration (blue plain), the barycentric method presented in section 3.1 (blue dotted) and the symmetric-barycentric extension presented in section 3.2 (red dotted).</figcaption></figure>

## 3.2 Towards Symmetric Transitive Registration

In this last section, we quickly introduce a way to derive approximately consistent transitive (at the first-order of the BCH approximation) registration from the barycentric SVFs computed in the previous section. Symmetry and transitivity are two important properties for registration methods to improve robustness and reduce the unpredictability of the results [10]. A registration method is said to be symmetric if it associates two points regardless of the order of images that are registered together (in the SVF setting it is equivalent to $v_i^j = -v_j^i$). Transitivity requires that the deformation given by the registration between two images should be equal whether it is done directly or by the composition of the result of the registration with an intermediate image (in the SVF setting it can be stated as $v_i^j = BCH(v_i^k, v_j^k) \simeq v_i^k + v_j^k$ with the BCH at the first order). Most registration methods fail to be transitive due to the accumulation of the registration errors at each step of the registration. Using Barycentric Subspaces as a basis for the registration at each step, we define the symmetric registration using the following formula which is schematically represented in Fig. 4:

$$W_s^t = \hat{v}_t - \hat{v}_s + \frac{1}{2} \sum_i (\lambda_s^i \hat{v}_t^i - \lambda_t^i \hat{v}_s^i). \qquad (2)$$

In this formula, the first two SVFs on the left represent the residual transformations from the barycentric subspace to the two time points, and the sum on the right is a symmetric estimation of the SVF $\hat{W}_t^s$ within the barycentric subspace by going through each reference image forward and backward. This formula defines registrations that are both symmetric and transitive up to higher orders of the BCH in the compositions. It can be used for frame-to-frame as well as for
frame-to-reference registration. In the former case, setting the reference to the first frame ($s = 1$ in the above formula) leads to improved results as shown in Fig. 4: the maximum error over the sequence is reduced by approx. 10% with respect to barycentric registration.

# 4 Conclusion

A new symmetric group-wise paradigm to study cardiac motion was proposed. Our approach relies on building subspaces as the reference for registration instead of choosing a specific arbitrary single image which can introduce bias. These subspaces represent the cardiac motion by meaningful parameters showing different clear patterns between two populations. Using these subspace as a prior, thereby relaxing the regularization on a 2-dimensional subspace, we achieve a better evaluation of the deformation between ED and ES frames and in particular we improve the estimation of the ejection fraction. Finally, the methodology can also be used to perform symmetric transitive registration, for better tracking along the sequence.

*Acknowledgements* The authors acknowledge the partial funding by the EU FP7-funded project MD-Paedigree (Grant Agreement 600932)

# References

1. Tobon-Gomez, C., et al.: Benchmarking framework for myocardial tracking and deformation algorithms: an open access database. Medical Image Analysis (2013)

2. Fletcher, P.T., Lu, C., Pizer, S.M., Joshi, S.: Principal geodesic analysis for the study of nonlinear statistics of shape. IEEE transactions on medical imaging (2004)

3. Pennec, X.: Barycentric Subspaces and Affine Spans in Manifolds. In: Geometric Science of Information GSI'2015, Second International Conference. Lecture Notes in Compuer Science, Palaiseau, France (2015)

4. Yigitsoy, M., Wachinger, C., Navab, N.: Temporal groupwise registration for motion modeling. In: Information Processing in Medical Imaging, Springer (2011)

5. Joshi, S., Davis, B., Jomier, M., Gerig, G.: Unbiased diffeomorphic atlas construction for computational anatomy. NeuroImage (2004)

6. Vercauteren, T., et al.: Symmetric log-domain diffeomorphic registration: A demons-based approach. In: Proc. Medical Image Computing and Computer Assisted Intervention (MICCAI'08), Part I. Lecture Notes in Computer Science (2008)

7. Balci, S.K., Golland, P., Wells, W.: Non-rigid groupwise registration using b-spline deformation model. Open source and open data for MICCAI (2007) 105–121

8. Perperidis, D., Mohiaddin, R.H., Rueckert, D.: Spatio-temporal free-form registration of cardiac mr image sequences. Medical image analysis 9(5) (2005) 441–456

9. Prakosa, A., et al.: Generation of Synthetic but Visually Realistic Time Series of Cardiac Images Combining a Biophysical Model and Clinical Images. IEEE Transactions on Medical Imaging (2013)

10. Škrinjar, O., Bistoquet, A., Tagare, H.: Symmetric and transitive registration of image sequences. Journal of Biomedical Imaging (2008)