# A Log-Euclidean Framework for Statistics on Diffeomorphisms

Vincent Arsigny, Olivier Commowick, Xavier Pennec, Nicholas Ayache

► To cite this version:

Vincent Arsigny, Olivier Commowick, Xavier Pennec, Nicholas Ayache. A Log-Euclidean Framework for Statistics on Diffeomorphisms. MICCAI 2006 - Medical Image Computing and Computer Assisted Intervention, Part I, 2006, Copenhaguen, Denmark. pp.924-931, 10.1007/11866565_113. inria-00615594

HAL Id: inria-00615594

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

Submitted on 30 Apr 2023

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
# A Log-Euclidean Framework for Statistics on Diffeomorphisms

Vincent Arsigny¹, Olivier Commowick¹,², Xavier Pennec¹, and Nicholas Ayache¹

¹ INRIA Sophia - Epidaure Project, 2004 Route des Lucioles BP 93
06902 Sophia Antipolis Cedex, France
Vincent.Arsigny@Polytechnique.org

² DOSISoft S.A., 45 Avenue Carnot, 94 230 Cachan, France

**Abstract.** In this article, we focus on the computation of statistics of invertible geometrical deformations (i.e., diffeomorphisms), based on the generalization to this type of data of the notion of *principal logarithm*. Remarkably, this logarithm is a simple 3D vector field, and is well-defined for diffeomorphisms close enough to the identity. This allows to perform *vectorial* statistics on diffeomorphisms, while preserving the invertibility constraint, contrary to Euclidean statistics on displacement fields. We also present here two efficient algorithms to compute logarithms of diffeomorphisms and exponentials of vector fields, whose accuracy is studied on synthetic data. Finally, we apply these tools to compute the mean of a set of diffeomorphisms, in the context of a registration experiment between an atlas and a database of 9 T1 MR images of the human brain.

## 1 Introduction

In this article, we focus on the computation of statistics of general *diffeomorphisms*, i.e. of geometrical deformations (non-linear in general) which are both one-to-one and regular (as well as their inverse). To quantitatively compare non-linear registration algorithms, or in order to constrain them, computing statistics on *global* deformations would be very useful as was done in [6] with *local* statistics.

The computation of statistics is closely linked to the issue of the *parameterization* of diffeomorphisms. Many algorithms, as in [5, 9], provide dense transformations which are guaranteed to be diffeomorphic, and parameterize them via their displacement field. However, Euclidean means of displacement fields do not necessarily yield invertible deformations, which makes Euclidean statistics on these parameters (or B-splines as in [10]) problematic for diffeomorphisms. In [8], it was proposed to parameterize arbitrary diffeomorphisms with Geodesic Interpolating Splines control points [4], and then to perform Euclidean operations on these low-dimensional parameters. However, although this guarantees the invertibility of the results, this may not be adequate for the whole variety of invertible transformations used in medical imaging.
To fully take into account the group structure of diffeomorphisms, it has been proposed to parameterize dense deformations with Hilbert spaces of *time-varying* speed vector fields, which yield geometrical deformations via the integration of an Ordinary Differential Equation (ODE) during one unit of time [14, 3]. In [15], it is suggested that the linear space of intial momenta of the geodesics of these spaces could provide an appropriate setting for statistics on diffeomorphisms. However, this is illustrated in [15] only in the case of landmark matching. To our knowledge, this statistical framework has not been used yet in the *general case*, certainly because of the *iterative nature* of the computation of the mean in this setting, which requires very stable numerical algorithms to converge.

In this article, we introduce a novel parameterization of diffeomorphisms, based on the generalization of the *principal logarithm* to non-linear geometrical deformations. Interestingly, this corresponds to parameterizing diffeomorphisms with *stationnary* speed vectors fields. As for matrices, this logarithm is well-defined only for transformations close enough to the identity. However, our preliminary numerical experiments on 3D non-rigid registration suggest that this limitation affects only very large deformations, and may not be problematic for image registration results. This novel setting is the infinite-dimensional analogous of the Log-Euclidean framework proposed in [2] for tensors. In this framework, usual Euclidean statistics can be performed on diffeomorphisms via their logarithms, with excellent mathematical properties like inversion-invariance.

The rest of this paper is organized as follows. In the next Section, we present the Log-Euclidean framework for diffeomorphisms, which is closely linked to the notion of *one-parameter subgroups*. Then, we present two efficient algorithms to compute the exponential of a vector field and the logarithm of a diffeomorphism, which are exemplified on synthetic data. Finally, we apply our framework to non-linear registration results to compute a Log-Euclidean mean deformation between a 3D atlas and a dataset of 9 $T_1$ MR images of human brains.

## 2 A Log-Euclidean Framework for Diffeomorphisms

To generalize the notion of *principal logarithm* to diffeomorphisms, we will rely on its close link with *one-parameter subgroups*. Remarkably, one-parameter subgroups of diffeomorphisms take a simple form, and by analogy with the finite-dimensional case, we will define logarithms of diffeomorphisms as being the *in-finitesimal generators* of these subgroups.

**Exponential and One-Parameter Subgroups.** Let $(G, \cdot)$ be a group, which can typically be a matrix multiplicative group. Then a family of elements $(g(s))_{s \in \mathbb{R}}$ of $G$ is a one-parameter subgroup of $G$ if and only if $g(0)$ is the neutral element $e$ of $G$ and for all $s, t$: $g(s) \cdot g(t) = g(s+t)$. When such a subgroup belongs to a (finite-dimensional) Lie group and is continuous, it is also *differentiable*, and its derivative at 0, $\frac{dg}{ds}(0)$, is called its *infinitesimal generator*.

Examples of such subgroups are given by the exponential in the following way : $(\exp(s \cdot M))_{s \in \mathbb{R}}$ is a one-parameter subgroup, where $M$ belongs to the
Lie algebra of $G$ (i.e., is an element of its tangent space at $e$). Conversely, we have the remarkable result that all continuous one-parameter subgroups are precisely of the form $(\exp(s.M))_{s \in \mathbb{R}}$. For more details, see [12], Chap. V. As a consequence, once one-parameter subgroups are identified in a Lie group, its exponential mapping is immediately known. This also yields the form taken by the logarithm, which is the (local) inverse of the exponential.

**One-Parameter Subgroups of Diffeomorphisms.** In the case of diffeomorphisms, what do continuous one-parameter subgroups look like? Quite intuitively, they are all obtained via the integration of *stationary* ODEs (also called *autonomous*), i.e ODE whose speed vector does not depend on time [13].

Let $\dot{x} = V(x)$ be a stationary ODE, where the vector field $V(x)$ is smooth enough. The flow associated to this ODE is the family of mappings $\Phi_V(\cdot, s) : \mathbb{R}^n \to \mathbb{R}^n$ parameterized by a time parameter $s$, such that for a fixed $x_0$, $s \mapsto \Phi_V(x_0, s)$ is the unique solution of $\dot{x} = V(x)$ with initial condition $x_0$ at time 0.

Intuitively, for a fixed $s$, the flow $x \mapsto \Phi_V(x, s)$ gives the way the ambient space is deformed by the integration of the ODE during $s$ units of time. Remarkably, the flow is a *one-parameter subgroup of diffeomorphisms*. This implies in particular that the deformations of space given at time 1 by $\Phi_V(\cdot, 1)$ are twice that observed at time 0.5 via $\Phi_V(\cdot, 0.5)$. $V(x)$ is the *infinitesimal generator* of this subgroup. Conversely, all continuous one-parameter subgroups of diffeomorphisms are flows of some stationary ODE. See [13] for more details.

**Exponentials of Vector Fields and Logarithms of Diffeomorphisms.**
Here, we propose to define the *exponential* $\exp(V)$ of a (smooth) vector field $V(x)$ as the flow at time 1 of the stationary ODE $\dot{x} = V(x)$. This is the only possible definition generalizing to vector fields the equivalence between one-parameter subgroups and exponential that exists in the finite-dimensional case. In this setting, the *logarithm* $\log(\Phi)$ of a diffeomorphism $\Phi$ close enough to the identity is the unique vector field in a neighborhood of zero such that $\exp(V) = \Phi$.

In the finite-dimensional case, the existence and uniqueness of such a logarithm can be proved by showing that the exponential is continuously differentiable and that its differential mapping at zero is the identity, and is thus locally invertible. Here, in the infinite-dimensional case, let $V$ be a regular vector field, let $t$ be a non-zero scalar, and let $\Phi_{t,V}(x, 1)$ be the flow at time 1 associated to $t.V(x)$. Then, a simple change of variable ($s \mapsto t.s$) shows that:

$$ \lim_{t \to 0} \frac{\Phi_{t,V}(x, 1) - x}{t} = \lim_{t \to 0} \frac{\exp(t.V)(x) - \exp(0)(x)}{t} = V(x), $$

which means intuitively that $\partial_V \exp(0) = V$, which indicates that the differential of the exponential is the identity, and that we have existence and uniqueness of the logarithm locally around the identity. To be entirely rigorous and to complete the proof, it would be necessary to define precisely the smoothness required of $V$, and to show the differentiability of the exponential in this infinite-dimensional setting. We are currently working on these quite technical details.
**Log-Euclidean Statistics on Diffeomorphisms.** On diffeomorphisms whose logarithm is well-defined, one can perform Euclidean operations, since these log-arithms are simple vector fields. In particular, one can define a distance be-tween these diffeomorphisms via a norm $\|\cdot\|$ on vector fields: $\text{dist}(\Phi_1, \Phi_2) = \|\log(\Phi_1) - \log(\Phi_2)\|$. Remarkably, this type of distance is *inversion-invariant*, since $\log(\Phi^{-1}) = -\log(\Phi)$. In the case of Hilbert norms, the point that mini-mizing the weighted sum of squared distances to the data is the *Log-Euclidean mean*, given by $\exp(\sum_i w_i \log(\Phi_i))$. This mean is inversion-invariant, and is also invariant with respect to the taking of square roots, since $\log(\Phi^{1/2}) = \frac{1}{2}\log(\Phi)$. More generally, in this setting, one can perform any kind of statistics on *diffeo-morphisms* via *vectorial statistics* on their logarithms, which allows a straightforward generalization of classical analysis tools like Principal Component Analysis (PCA) on diffeomorphisms.

## 3 Computation of Exponentials and Logarithms

**Fast Computation of Exponentials.** Here, we present an efficient algorithm to compute exponentials of speed vectors on regular grids, which generalizes the Fast Polyaffine Transform we proposed in [1]. It is based on a non-linear generalization of the popular ‘Scaling and Squaring’ method, widely used to compute *matrix exponentials*. The basic idea is that the matrix exponential is much simpler to compute for matrices *close to zero*, for example using Padé approximants. In particular, one can compute very accurately $\exp(\frac{M}{2^N})$ and obtain $\exp(M) = \exp(\frac{M}{2^N})^{2^N}$ by squaring recursively $N$ times the result [7].

In the non-linear case, one can follow exactly the same steps as for matrices and generalize the ‘Scaling and Squaring’ to vector fields in the following way:

1. **Scaling step:** divide $V(x)$ by a factor $2^N$, so that $V(x)/2^N$ is close enough to zero (according to the level of accuracy desired).

2. **Exponentiation step:** $\Phi_V(\frac{1}{2^N}, \cdot)$ is computed with a first-order explicit numerical scheme.

3. **Squaring step:** $N$ recursive squarings of $\Phi_V(\frac{1}{2^N}, \cdot)$ yield an accurate estimation of $\Phi_V(1, \cdot)$ (only $N$ compositions of mappings are used).

Intuitively, this means that we obtain the deformations at time $1$ as a result of $2^N$ times the composition of the (very) small deformations observed at time $\frac{1}{2^N}$. This allows a fast integration of stationary ODEs on regular grids, compared to classical integrations based on fixed time-steps. To perform the composition of (sampled) deformations on regular grids, we use bi- or tri-linear interpolation, which guarantees the continuity of the interpolation in a simple way.

**Fast Computation of Logarithms.** Exactly as for exponentials, we use here a non-linear generalization of the popular ‘Inverse Scaling and Squaring’ method (ISS), widely used to compute *matrix logarithms*. Anew, the idea is that logarithms are much simpler to compute for matrices *close to the identity*, for instance with Padé approximants. To transform a matrix $M$ so that it is closer to
the identity, the ISS algorithm performs recursive computations of square roots. Then the identity $log(M) = 2^N \cdot log(M^{2^{-N}})$ is used to compute $log(M)$.

Let $\Phi$ be a diffeomorphism. To compute its logarithm $V$, we use the following non-linear generalization of the ISS algorithm:

1. **Scaling step:** chose a scaling factor $2^N$ (according to the level of accuracy desired).

2. **Square rooting step:** $\Phi^{2^{-N}}$ is computed by $N$ successive recursive takings of square roots.

3. **Computation of logarithm step:** $log(\Phi)$ is given by $2^N \cdot log(\Phi^{2^{-N}})$, where $log(\Phi^{2^{-N}})$ is simply estimated by $\Phi^{2^{-N}} - Id$ (where $Id$ is the identity).

In order to compute square roots, we perform here a gradient descent on the functional $E_{\text{SQRT}}(T) = \frac{1}{2} \int_{\Omega} \|T \circ T - \Phi\|^2(x)dx$, with $\frac{1}{2} \cdot (\Phi - Id)$ as initialization. The $(L^2)$ gradient of this energy is the following:

$$ \nabla E_{\text{SQRT}}(T) = (DT^t) \circ T \cdot (T \circ T - \Phi) + |\det(D(T^{-1}))|(T - \Phi \circ T^{-1}), $$

where $DT$ is the Jacobian of $T$ and where ‘.’ and ‘$M^t$’ are matrix multiplication and transposition. Exactly as in the matrix case, this requires the inversion of $T$, which we also perform by gradient descent on $E_{\text{INV}}(T) = \frac{1}{2} \int_{\Omega} \|\Phi \circ T - Id\|^2(x)dx$, with an initialization of $-(\Phi - Id)$. The $(L^2)$ gradient is given by:

$$ \nabla E_{\text{INV}}(T) = (D\Phi^t) \circ T \cdot (\Phi \circ T - Id). $$

<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Fast computation of the exponential of a random vector field. From left to right: the two last iterations (scaling of $2^8$) of our fast computation of the exponential and the evolution of the average relative accuracy with $N$. Note how deformations increase exponentially from one iteration to another. The relative accuracy obtained on average converges toward of 0.3%, which is approximately obtained for $N = 7$.</figcaption></figure>

**Experiments on 2D Synthetic Data.** In these experiments, we evaluate the accuracy of our algorithms on 2D random deformations, sampled on a 40 × 40 regular grid. To generate a random vector field, random displacements are computed on a 11 × 11 regular grid (Gaussian white noise of standard deviation
0.2, coordinates of grid points in [-5, 5]) , which are smoothly interpolated on
the finer grid using bilinear interpolation.

We measure the accuracy of the fast computation of the exponential by
computing the relative difference with respect to a very accurate estimation
of the continuous transformation, obtained by a classical integration (i.e., with
fixed time step, here 2⁻⁸) of the stationary ODE associated to the random vector
field, for each voxel of the grid. Fig. 1 presents the results, which show that a
typical accuracy on average of 0.3% can be obtained with 7 squarings.

To evaluate the accuracy of the fast computation of the logarithm, we first
estimate the exponential of the random vector field with our fast algorithm
(with 8 squarings). Then, the logarithm of this diffeomorphism is computed and
is compared to the original vector field (10 iterations are used for each gradient
descent). Fig. 2 presents the results, which show that an accuracy of 3% percent
can be obtained on average with 6 recursive computations of square roots.

<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Fast computation of the logarithm of a random deformation. From left to right: random vector field (logarithm), then difference between original and estimated logarithms (for *N* = 8), evolution of accuracy with *N*. The relative accuracy obtained on average converges toward of 3%.</figcaption></figure>

4 Statistics on 3D Diffeomorphisms

In this Section, we compute Log-Euclidean statistics on the registration results
obtained between a 256 × 256 × 60 artificial T₁ MR image of a human brain
(coming from the BrainWeb <fnref n="3" />), referred to here as the ‘atlas’, and a dataset of 9
T₁ images. The images of the dataset are first globally aligned with the atlas using
a robust block-matching affine registration algorithm. Then, a fine registration
is performed using a dense transformation registration algorithm guaranteeing
the invertibility of the global deformations [11]. Statistics are performed on the
deformations from the atlas to each of the subjects’ geometry obtained in the
non-linear registration step.

To compute logarithms of diffeomorphisms and exponential of vector fields,
we used the fast algorithms presented in the former section. A scaling of 2⁸
is used in both cases, and a maximal number of 10 iterations is used during
each gradient descent.On the 256 × 256 × 60 grid, using a Intel M processor

<footnote n="3">³ Web site: http://www.bic.mni.mcgill.ca/brainweb/</footnote>
2.13GHz with 1 Go of RAM, an exponential was computed in typically 30s and a logarithm in 60 minutes.

Fig. 3 presents the results obtained on mean global deformations. We see in particular that the Log-Euclidean mean deformation tends to increase the size of ventricles, which are smaller in the atlas than in the dataset. We also compared the Log-Euclidean mean to the Euclidean mean of displacements, which is not guaranteed to be invertible. In this example, both means are quite close to each other, although locally, one can observe in the region of large mean deformations relative differences of the order of 30%, for example in the ventricles.

<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Deformation of atlas with Log-Euclidean mean diffeomorphism. From left to right, images of: regular grid deformed by Log-Euclidean mean deformations, Jacobian of Log-Euclidean mean deformations, norm of Log-Euclidean mean deformations, and finally norm of difference between mean Euclidean and Log-Euclidean deformations of atlas. The largest mean deformations are observed in the ventricles, on the cortex and the skull; this is due to the anatomical differences between the atlas and the population. On the Jacobian map, we see in particular that high values are obtained in the ventricles: this comes the fact that the atlas has on average smaller ventricles than in the population. In this example, both means are quite close to each other, although locally, one can observe in the region of large mean deformations relative differences of the order of 30%, for example in the ventricles.</figcaption></figure>

# 5 Conclusions and Perspectives

In this article, we have presented an infinite-dimensional generalization to diffeomorphisms of our Log-Euclidean framework for tensors [2], based on the generalization of the notion of *principal logarithm* to invertible deformations. This logarithm is a simple 3D vector field, and is well-defined for diffeomorphisms close enough to the identity. In the set of diffeomorphisms whose logarithm is well-defined, one can perform Euclidean operations, since these logarithms are simple vector fields. As a consequence, this framework can be used to perform *vectorial* statistics on diffeomorphisms which are always preserve the invertibility constraint (contrary to Euclidean statistics on displacement fields), and are also invariant with respect to inversion.
Furthermore, we have also presented two efficient algorithms which successfully generalize to the non-linear case two popular algorithms used to compute the matrix exponential and logarithm. In practice, this allows for example the computation of Log-Euclidean means of 3D geometrical transformations, as we have shown in the context of a registration experiment between an atlas and a database of 9 subjects. This type of statistics could prove very useful to quantitatively compare registration algorithms, or to constrain them. Also, in the domain of image and shape statistics, the Log-Euclidean framework for diffeomorphisms could provide an interesting setting to build models with a constant topology.

## References

1. V. Arsigny, O. Commowick, X. Pennec, and N. Ayache. A Log-Euclidean polyaffine framework for locally rigid or affine registration. In *Proc. of WBIR'06*, LNCS, Utrecht, The Netherlands, 9 - 11 July 2006. Springer. To appear.

2. V. Arsigny, P. Fillard, X. Pennec, and N. Ayache. Fast and simple calculus on tensors in the Log-Euclidean framework. In *MICCAI (1)*, pages 115–122, 2005.

3. M. F. Beg, M. I. Miller, A. Trouvé, and L. Younes. Computing large deformation metric mappings via geodesic flows of diffeomorphisms. *Int. Jour. Comp. Vis.*, 61(2):139–157, 2005.

4. V. Camion and L. Younes. Geodesic interpolating splines. In M. Figuredo, J. Zerubia, and A. Jain, editors, *Proc. of Energy Minimization Methods in Comp. Vis. and Pat. Rec. (EMMCVPR'01)*, LNCS 2134, pages 513–527, 2001.

5. C. Chefd'hotel, G. Hermosillo, and O. Faugeras. Flows of diffeomorphisms for multimodal image registration. In *Proc. of ISBI*, 2002.

6. O. Commowick, R. Stefanescu, P. Fillard, V. Arsigny, N. Ayache, X. Pennec, and G. Malandain. Incorporating statistical measures of anatomical variability in atlas-to-subject registration for conformal brain radiotherapy. In *Proc. of MICCAI'2005 (II)*, LNCS, pages 927–934, 2005.

7. N. J. Higham. The scaling and squaring method for the matrix exponential revisited. *SIAM J. Matrix Anal. Appl.*, 26(4):1179–1193, 2005.

8. S. Marsland and C. J. Twining. Constructing diffeomorphic representations for the groupwise analysis of nonrigid registrations of medical images. *IEEE Trans. Med. Imaging*, 23(8):1006–1020, 2004.

9. X. Pennec, R. Stefanescu, V. Arsigny, P. Fillard, and N. Ayache. Riemannian elasticity: A statistical regularization framework for non-linear registration. In *Proc. of MICCAI 2005 (II)*, LNCS, pages 943–950, 2005.

10. D. Rueckert, A. F. Frangi, and J. A. Schnabel. Automatic construction of 3D statistical deformation models of the brain using non-rigid registration. *IEEE TMI*, 22(8):1014–1025, 2003.

11. R. Stefanescu, X. Pennec, and N. Ayache. Grid powered nonlinear image registration with locally adaptive regularization. *MedIA*, 8(3):325–342, September 2004.

12. S. Sternberg. *Lectures on Differential Geometry*. Prentice Hall Mathematics Series. Prentice Hall Inc., 1964.

13. M. Tenenbaum and H. Pollard. *Ordinary Differential Equations*. Dover, 1985.

14. A. Trouvé. Diffeomorphisms groups and pattern matching in image analysis. *International Journal of Computer Vision*, 28(3):213–221, 1998.
15. M. Vaillant, M. Miller, L. Younes, and A. Trouvé. Statistics on diffeomorphisms via tangent space representations. *NeuroImage*, 23:S161-S169, 2004.