# Symmetric Log-Domain Diffeomorphic Registration: A Demons-based Approach

Tom Vercauteren, Xavier Pennec, Aymeric Merchant, Nicholas Ayache

## ► To cite this version:

Tom Vercauteren, Xavier Pennec, Aymeric Merchant, Nicholas Ayache. Symmetric Log-Domain Diffeomorphic Registration: A Demons-based Approach. MICCAI 2008 - Medical Image Computing and Computer Assisted Intervention, Sep 2008, New York, United States. pp.754-761, 10.1007/978-3-540-85988-8_90. inria-00280602

HAL Id: inria-00280602

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

Submitted on 19 May 2008

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>
Symmetric Log-Domain Diffeomorphic
Registration: A Demons-based Approach

Tom Vercauteren¹, Xavier Pennec², Aymeric Perchant¹, and Nicholas Ayache²

¹Mauna Kea Technologies, France

²Asclepios, INRIA Sophia-Antipolis, France

**Abstract.** Modern morphometric studies use non-linear image registration to compare anatomies and perform group analysis. Recently, log-Euclidean approaches have contributed to promote the use of such computational anatomy tools by permitting simple computations of statistics on a rather large class of invertible spatial transformations. In this work, we propose a non-linear registration algorithm perfectly fit for log-Euclidean statistics on diffeomorphisms. Our algorithm works completely in the log-domain, i.e. it uses a stationary velocity field. This implies that we guarantee the invertibility of the deformation and have access to the true inverse transformation. This also means that our output can be directly used for log-Euclidean statistics without relying on the heavy computation of the log of the spatial transformation. As it is often desirable, our algorithm is symmetric with respect to the order of the input images. Furthermore, we use an alternate optimization approach related to Thirion's demons algorithm to provide a fast non-linear registration algorithm. First results show that our algorithm outperforms both the demons algorithm and the recently proposed diffeomorphic demons algorithm in terms of accuracy of the transformation while remaining computationally efficient.

# 1 Introduction

Non-linear image registration has opened the way for computational characterization of morphological evolution and morphological variability. Most computational anatomy tools make use of registration results [1–4] but require that they satisfy some advanced properties such as invertibility and symmetry with respect to the order of the inputs. Image registration schemes can thus only be used if they meet the requirements of these tools. Large deformation diffeomorphic methods have initially been developed for this purpose. Transformations are determined by a time-varying ordinary differential equation (ODE) [5]. Following the seminal work on inverse consistency [6], the large deformations framework has also been extended to enforce the symmetry of the solution [3, 7].

A widely acknowledged issue with the large deformation setting lies in its computational complexity and memory requirements. Recent work has strived towards bridging the gap between these rigorous mathematical tools and very efficient non-linear registration schemes such as Thirion's demons algorithm [8]. On one hand it has been proposed to constrain the large deformation setting by using
transformations that satisfy the initial momentum conservation [9, 10]. Similarly, in [1], the authors proposed to parameterize the diffeomorphisms with *stationary* velocity fields to allow easy computations of statistics on diffeomorphisms. This parameterization was used for registration in [11, 12]. These algorithms are well-fit for further statistical processing [1, 2, 4]. On the other hand, in [13] the authors proposed an efficient diffeomorphic registration scheme based on the demons algorithm that encodes the optimization steps, but not the complete transformation, with such stationary velocity fields. Between these attempts, we believe that there is still a gap to be bridged. Second generation large deformation algorithms still need to solve rather large Euler-Lagrange equations at each iteration and the diffeomorphic demons cannot be directly used by the recent statistical tools we mentioned. Furthermore, the demons algorithm is not symmetric with respect to the order of the images to register.

In this work, we propose an image registration scheme that uses a demons-like alternate optimization approach for efficiency but represents the complete deformation as an exponential of a smooth velocity field. This approach will hereafter be referred to as log-domain one. Thanks to such representation, our results are symmetric and can be directly used for computational anatomy.

The remainder of this paper is organized as follows. Classical and diffeomorphic demons are presented in Section 2. Our log-domain approach is developed in Section 3 and evaluated in Section 4. Finally Section 5 concludes the paper.

## 2 Additive and Diffeomorphic Demons Algorithms

Non-linear image registration aims at finding a well-behaved spatial transformation $s(.)$ that best aligns two given images $I_0(.)$ and $I_1(.)$. Typically, a similarity criterion $Sim(I_0, I_1, s)$ is used to measure the resemblance of the aligned images while a regularization energy $Reg(s)$ estimates the likelihood of the transformation. *Non-parametric* methods need to find the displacement $s(p)$ of each point $p$ in order to optimize the following energy functional:

$$E(s) = \frac{1}{\sigma_i^2} \mathrm{Sim}(I_0, I_1, s) + \frac{1}{\sigma_T^2} \mathrm{Reg}(s),$$

where $\sigma_i$ accounts for the noise on the image intensity, and $\sigma_T$ controls the amount of regularization we need. The desired properties of the final spatial transformation can be encoded within the regularization term or can be enforced by constraining the search space. Instead of looking for a solution in the complete space of non-parametric spatial transformations, it is for example possible to search only within a subspace of diffeomorphisms.

**Additive Demons.** Even if we use a simple transformation space and a classical regularization term, approaching the registration problem directly often leads to computationally expensive iterations that need the solution of an Euler-Lagrange equation. Contrastingly, Thirion's demons algorithm uses an efficient two-step procedure at each iteration [8]. It first looks for an unconstrained update step
with an optical flow computation and then uses a simple Gaussian smoothing on the updated transformation. It has been shown in [14] that the demons algorithm could be cast to the minimization of a well-posed criterion by introducing a hidden variable c for point correspondences. The interest of this auxiliary variable is that it decouples the optimization into easily tractable sub-problems. Each iteration walks towards the optimum of the global energy

$$E(I_0, I_1, c, s) = \frac{1}{\sigma_i^2} \mathrm{Sim}(I_0, I_1, c) + \frac{1}{\sigma_x^2} \mathrm{dist}(s, c)^2 + \frac{1}{\sigma_T^2} \mathrm{Reg}(s), \quad (1)$$

where $\sigma_x$ accounts for a spatial uncertainty on the correspondences. We clas- sically have $\mathrm{Sim}(I_0, I_1, c) = \|I_0 - I_1 \circ c\|^2$, $\mathrm{dist}(s, c) = \|c - s\|$ and $\mathrm{Reg}(s) = \|\nabla s\|^2$ but more advanced criteria can be used [14].

In the additive demons algorithm, the optimization is performed within the complete space of non-parametric transformation using additive updates of the form $s + u$. The optical flow procedure solves for the correspondence energy

$$E_{\mathrm{add}}^{\mathrm{corr}}(I_0, I_1, s, u) = \mathrm{Sim}(I_0, I_1, s + u) + \|u\|^2 \quad (2)$$

with respect to $u$, while the Gaussian smoothing solves for the regularization. Different optimizers lead to different forces that have been justified in [15].

**Diffeomorphic Demons.** In [13], the authors proposed to adapt the demons algorithm to provide diffeomorphisms. The diffeomorphic demons algorithm uses Thirion's alternate optimization approach to maintain the computational efficiency but works in a space of diffeomorphisms to enforce the invertibility.

An efficient computational framework for diffeomorphisms was proposed in [1]. It uses a Lie group structure that defines an exponential mapping from the vector space of smooth stationary velocity fields to diffeomorphisms. The exponential $\exp(u)$ of a velocity field is given by the flow at time one of the stationary ODE: $\partial p(t)/\partial t = u(t)$. The nice property of this framework lies in the low computational requirement needed to compute the exponential.

At each iteration, the diffeomorphic demons takes advantage of this exponential mapping by looking for an update step $u$ in the Lie algebra (the vector space of velocity fields) and then by mapping it in the space of diffeomorphisms through the exponential. The update step is thus of the form $s \circ \exp(u)$. The advantage of this approach is that it can compute $u$ with an unconstrained optimizer that has the same form and complexity as the classical demons forces [13]. The diffeomorphic demons retains the simple Gaussian smoothing of Thirion's algorithm for its efficiency. With this approach the first step optimizes the modified correspondence energy

$$E_{\mathrm{diffeo}}^{\mathrm{corr}}(I_0, I_1, s, u) = \mathrm{Sim}(I_0, I_1, s \circ \exp(u)) + \|u\|^2. \quad (3)$$

## 3 A Log-Domain Approach to Diffeomorphic Demons

The parameterization of diffeomorphic transformations through stationary velocity fields proposed in [1] provides a very efficient means of dealing with dif-
diffeomorphisms. In the diffeomorphic demons of [13], the exponential is used only to encode the adjustment made at each iteration to the current transformation. This leads to a computationally attractive scheme but lacks some of the characteristics of log-domain registration tools [11, 12] that encode the complete transformation with stationary velocity fields: namely the ability to compute the inverse of the transformation at a very low cost, the symmetry of the registration result and the adequacy of the representation for the statistical tools of [1, 2].

Our main contribution in this paper is to show that the diffeomorphic demons can be extended to represent the complete spatial transformation in the log domain. The main idea of the proposed algorithm is to represent the current transformation $s$ as an exponential of a smooth velocity field $v$, i.e. $s = \exp(v)$, and use the diffeomorphic demons to efficiently compute a field $u$ for an update of the form $s \circ \exp(u)$, i.e. $\exp(v) \circ \exp(u)$. With this idea, there are two questions that come to mind. First of all, given that we want to represent everything in the log-domain, we need to know whether for any $v$ and $u$ there exists a velocity field $w$ such that $\exp(w) = \exp(v) \circ \exp(u)$. Then we need to design a regularization scheme that is consistent with the log-domain representation.

**Baker-Campbell-Hausdorff Approximations.** Our first goal is to find a smooth velocity field $Z(v, \varepsilon u)$ such that

$$ \exp(Z(v, \varepsilon u)) \approx \exp(v) \circ \exp(\varepsilon u), \quad (4) $$

where $\varepsilon$ is simply used to emphasize the fact that we look for an approximation valid for small $\varepsilon u$ (to encode the update) but arbitrary $v$ (to encode the complete transformation). Since we deal with an infinite-dimensional space which has a Lie group structure but is not an actual Lie group, the question of the existence of such a velocity field is a tough mathematical one that needs further investigation. In practice though, it has been shown in [2] that the Baker-Campbell-Hausdorff (BCH) formula which is valid for finite-dimensional spaces could be applied successfully on diffeomorphisms. By using the first terms of the BCH formula, the authors of [2] obtained an approximation that seems well-fit for the demanding application of brain atlas construction. In our setting, since only $\varepsilon u$ is assumed to be small, the first order approximation of $Z(v, \varepsilon u)$ is:

$$ Z(v, \varepsilon u) = v + \varepsilon u + \frac{1}{2}[v, \varepsilon u] + \frac{1}{12}[v, [v, \varepsilon u]] + O(\|\varepsilon u\|^2), \quad (5) $$

where the Lie bracket $[v, u]$ provides a velocity field defined at each point $p$ by<fnref n="1" />:

$$ [v, u](p) = \mathrm{Jac}(v)(p) \cdot u(p) - \mathrm{Jac}(u)(p) \cdot v(p). \quad (6) $$

This first-order approximation provides a good candidate for our update rule but is still somewhat complex as it requires three Lie brackets. This might also

<footnote n="1">¹ Most authors define the Lie bracket as the opposite of (6). Numerical simulations, and personal communication with M. Bossa, showed the relevance of this definition. Future research will aim at fully understanding the reason of this discrepancy.</footnote>
lead to an unstable numerical scheme as the imbricated Lie bracket amounts to using second-order derivatives. An ad hoc simplification can be made by simply considering the first terms of the BCH expansion. We chose to evaluate the quality of the following approximations: $Z_A(v, \varepsilon u) \triangleq v + \varepsilon u$, $Z_B(v, \varepsilon u) \triangleq v + \varepsilon u + \frac{1}{2}[v, \varepsilon u]$ and $Z_C(v, \varepsilon u) \triangleq v + \varepsilon u + \frac{1}{2}[v, \varepsilon u] + \frac{1}{12}[v, [v, \varepsilon u]]$.

In order to test the validity of these approximations for our application we set up a small experiment to measure the error between $\exp(Z_X(v, \varepsilon u))$ and $\exp(v) \circ \exp(\varepsilon u)$. We generate a random $v$ at a given energy, a random $u$ at a lower energy and measure $\|\exp(Z_X(v, \varepsilon u)) - \exp(v) \circ \exp(\varepsilon u)\|^2$. Due to space constraints, only the conclusions of these experiments can be presented here. The best results are as expected provided by $Z_C$ with $Z_B$ being only a few percents away from it. $Z_A$ surprisingly still provide decent results but the error is however one order of magnitude away from the error resulting from $Z_C$.

**Log-Domain Diffeomorphic Demons** The BCH approximations allow us to cast the update step $s \leftarrow s \circ \exp(u)$ used in the diffeomorphic demons into a log-domain update $v \leftarrow Z_X(v, u)$ provided that the current transformation $s$ can be expressed as an exponential $s = \exp(v)$. It might however be unclear why one would resort to such a BCH approximation. It could indeed be possible to directly look for an update of the form $v \leftarrow v+u$. The problem with this kind of approach used for example in [11] lies in its computational complexity. Since the exponential is not used around zero, we cannot take advantage of the fact that the derivative of the exponential around zero is the identity. This introduces a coupling between the transformation and the update which doesn't appear here.

Finally, in order to be consistent with the log-domain representation but keep the simplicity of the demons algorithm, we chose to perform a Gaussian smoothing directly in the log-domain. Our framework can simply be linked to (1) by using $\text{dist}(s, c) = \|\log(s^{-1} \circ c)\|$ and $\text{Reg}(s) = \|\nabla \log(s)\|^2$.

**Algorithm 1 (Log-Domain Demons)**

- Choose a starting spatial transformation $s = \exp(v)$

- Iterate until convergence:

* Given the current transformation $s = \exp(v)$, compute a correspondence update field $u$ by minimizing $E_{\text{diff}}^{\text{corr}}(I_0, I_1, s, u)$ with respect to $u$

* For fluid-like regularization let $u \leftarrow K_{\text{fluid}} \star u$

* Let $v \leftarrow Z_X(v, u)$, e.g. $v \leftarrow v + u + \frac{1}{2}[v, u]$

* For diffusion-like regularization let $v \leftarrow K_{\text{diff}} \star v$

**Symmetric Extension** The inverse of a spatial transformation $s$ parameterized in the log-domain $s = \exp(v)$, can be obtained at almost no cost by a backward computation $s^{-1} = \exp(-v)$. A symmetric registration framework can be obtained from the non-symmetric one by symmetrizing the global energy:

$$s_{opt} = \arg \min_{s} (E(I_0, I_1, s) + E(I_1, I_0, s^{-1})), \quad (7)$$

where $c$ has been omitted from the energy (1) for clarity purposes.
Our second main contribution in this work is to provide an efficient scheme for solving this symmetrized system. We formulate it as a constrained optimization using two diffeomorphisms: $[s_{opt}, s_{opt}^{-1}] = \arg \min_{[s,t] | t=s^{-1}} E(I_0, I_1, s) + E(I_1, I_0, t)$. We propose to use an unconstrained optimization step on the pair $[s, t]$ and then to project the new transformations onto the space of symmetric transformations $\{[s, t] | t = s^{-1}\}$. By using a complete log-domain demons iteration starting from $s = \exp(v)$ to optimize the first term $E(F, M, \exp(\zeta))$, we get $\zeta = K_{\text{diff}} \star Z(v, K_{\text{fluid}} \star u^{\text{forw}})$, where $u^{\text{forw}}$ is the demons force. Similarly, the second term $E(M, F, \exp(-\tau))$ is optimized with $\tau = -K_{\text{diff}} \star Z(-v, K_{\text{fluid}} \star u^{\text{back}})$, where $u^{\text{back}}$ is the demons force for reversed inputs.

Thanks to the log-domain representation, we deal with a vector space and can design an easy projection operator that guarantees the symmetry of the results. We simply average, in the log-domain, the forward and backward iterations:

$$v \leftarrow \frac{1}{2} K_{\text{diff}} \star (Z(v, K_{\text{fluid}} \star u^{\text{forw}}) - Z(-v, K_{\text{fluid}} \star u^{\text{back}})). \quad (8)$$

As an example, using $Z_A$ provides $v \leftarrow K_{\text{diff}} \star (v + \frac{1}{2} K_{\text{fluid}} \star (u^{\text{forw}} - u^{\text{back}}))$.

<caption>**Algorithm 2 (Symmetric Iteration using $Z_A$)**</caption>

– Compute the demons forces $u^{\text{forw}}$ to minimize $E_{\text{diffeo}}^{\text{corr}}(I_0, I_1, \exp(v), u^{\text{forw}})$
– Compute the demons forces $u^{\text{back}}$ to minimize $E_{\text{diffeo}}^{\text{corr}}(I_1, I_0, \exp(-v), u^{\text{back}})$
– For fluid-like regularization let $u \leftarrow \frac{1}{2} K_{\text{fluid}} \star (u^{\text{forw}} - u^{\text{back}})$
– For diffusion-like regularization let $v \leftarrow K_{\text{diff}} \star (v + u)$ else let $v \leftarrow v + u$ ■

<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Left column: Reference image, in vivo microscopy of normal colonic mucosa, courtesy of [masked for anonymity], and one example random warp. Other columns: Registration results on 100 random experiments. The log-domain diffeomorphic performs similarly to the diffeomorphic demons while our symmetric approach outperforms it. We see the small impact of the BCH expansion that we use.</figcaption></figure>
## 4 Experiments

The proposed algorithms were evaluated as follows. A reference image $I_{\text{ref}}$ is deformed through a random diffeomorphism. Some random noise is added to both the reference and warped image. The pair of images is registered first using $I_{\text{ref}} = I_0$ and then using $I_{\text{ref}} = I_1$. We compare the additive and diffeomorphic demons, the proposed log-domain demons with two different BCH expansions (using $Z_A$ and $Z_B$), and the proposed symmetric log-domain demons again with two different BCH expansions. To make a fair comparison, each algorithm uses the same expression of the demons forces. Figure 1 shows the evolution of several criterion over the iterations. The choice of the BCH expansion does not significantly change the performance of our schemes. Hence we advocate the use of the simplest one ($Z_A$). Our log-domain schemes perform similarly to the diffeomorphic demons but allow an easy computation of the inverse and are well-fit for statistical analysis. Finally our symmetric extension outperforms the other algorithms in terms of distance to the ground truth transformation. The computational time is only twice the one of the diffeomorphic demons.

Finally Fig. 2 illustrates the adaptation of our algorithms for computational anatomy. A simple atlas is built from the 20 synthetic anatomies [16]. Of course more advanced techniques [1–3] should be used, but this proof of concept opens the way to neat future work. Thanks to our regularization by a simple smoothing, the integration of deformation statistics could be as simple as performing a non-stationary smoothing. Local covariance matrices could indeed be used to replace the norm used in the regularization by a Mahalanobis distance.

## 5 Conclusion

We proposed an efficient diffeomorphic algorithm that combines the efficiency of the demons algorithm and the desirable properties of modern large deformation algorithms such as invertibility with respect to the order of the inputs and memory efficient representations. Since we consider the spatial transformations as exponentials of smooth velocity fields, our results can directly be used by the recent statistical tools of [1, 2]. We focused on a simple similarity criterion but

<figure><img src="image_3.png" /><figcaption>**Fig. 2.** 19 3D synthetic MRs of distinct anatomies were registered to an arbitrary reference. We show the principal direction of variability found by statistical analysis.</figcaption></figure>
our approach can easily be extended to other intensity relationships by borrowing ideas from [7, 14, 17]. The next step will be to integrate deformation statistics within the algorithm by locally adapting the regularization.

## References

1. Arsigny, V., Commowick, O., Pennec, X., Ayache, N.: A Log-Euclidean framework for statistics on diffeomorphisms. In: Proc. MICCAI'06. Volume 4190 of LNCS. (2006) 924-931

2. Bossa, M., Hernandez, M., Olmos, S.: Contributions to 3D diffeomorphic atlas estimation: Application to brain images. In: Proc. MICCAI'07. Volume 4792 of LNCS. (October 2007) 667-674

3. Joshi, S., Davis, B., Jomier, M., Gerig, G.: Unbiased diffeomorphic atlas construction for computational anatomy. Neuroimage 23(S1) (2004) S151-S160

4. Vaillant, M., Miller, M.I., Younes, L., Trouvé, A.: Statistics on diffeomorphisms via tangent space representations. Neuroimage 23(S1) (2004) 161-169

5. Beg, M.F., Miller, M.I., Trouvé, A., Younes, L.: Computing large deformation metric mappings via geodesic flows of diffeomorphisms. Int. J. Comput. Vis. 61(2) (February 2005) 139-157

6. Christensen, G.E., Johnson, H.J.: Consistent image registration. IEEE Trans. Med. Imag. 20(7) (July 2001) 568-582

7. Avants, B.B., Epstein, C.L., Grossman, M., Gee, J.C.: Symmetric diffeomorphic image registration with cross-correlation: Evaluating automated labeling of elderly and neurodegenerative brain. Med. Image Anal. 12(1) (February 2008) 26-41

8. Thirion, J.P.: Image matching as a diffusion process: An analogy with Maxwell's demons. Med. Image Anal. 2(3) (1998) 243-260

9. Marsland, S., McLachlan, R.I.: A Hamiltonian particle method for diffeomorphic image registration. In: Proc. IPMI'07. Volume 4584 of LNCS. (July 2007) 396-407

10. Younes, L.: Jacobi fields in groups of diffeomorphisms and applications. Quart. Appl. Math. 65 (February 2007) 113-134

11. Ashburner, J.: A fast diffeomorphic image registration algorithm. Neuroimage 38(1) (October 2007) 95-113

12. Hernandez, M., Bossa, M.N., Olmos, S.: Registration of anatomical images using geodesic paths of diffeomorphisms parameterized with stationary vector fields. In: Proc. MMBIA Workshop of ICCV'07. (October 2007) 1-8

13. Vercauteren, T., Pennec, X., Perchant, A., Ayache, N.: Non-parametric diffeomorphic image registration with the demons algorithm. In: Proc. MICCAI'07. Volume 4792 of LNCS. (October 2007) 319-326

14. Cachier, P., Bardinet, E., Dormont, D., Pennec, X., Ayache, N.: Iconic feature based nonrigid registration: The PASHA algorithm. Comput. Vis. Image Underst. 89(2-3) (February 2003) 272-298

15. Vercauteren, T., Pennec, X., Malis, E., Perchant, A., Ayache, N.: Insight into efficient image registration techniques and the demons algorithm. In: Proc. IPMI'07. Volume 4584 of LNCS. (July 2007) 495-506

16. Aubert-Broche, B., Griffin, M., Pike, G.B., Evans, A.C., Collins, D.L.: Twenty new digital brain phantoms for creation of validation image data bases. IEEE Trans. Med. Imag. 25(11) (November 2006) 1410-1416

17. Hermosillo, G., Chefd'Hotel, C., Faugeras, O.: Variational methods for multimodal image matching. Int. J. Comput. Vis. 50(3) (December 2002) 329-343