# Proceedings of the Third International Workshop on Mathematical Foundations of Computational Anatomy (MFCA'11)- Geometrical and Statistical Methods for Modelling Biological Shape Variability

Xavier Pennec, Sarang Joshi, Mads Nielsen

► To cite this version:

Xavier Pennec, Sarang Joshi, Mads Nielsen (Dir.). Proceedings of the Third International Workshop on Mathematical Foundations of Computational Anatomy (MFCA'11)- Geometrical and Statistical Methods for Modelling Biological Shape Variability. Pennec, Xavier and Joshi, Sarang and Nielsen, Mads. HAL, pp.200, 2011. hal-00813806

HAL Id: hal-00813806

https://inria.hal.science/hal-00813806v1

Submitted on 17 Apr 2013

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 CC BY 4.0 - Attribution - International License
<figure><img src="image_2.png" /></figure>

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

Mathematical Foundations of
Computational Anatomy

Geometrical and Statistical Methods
for Biological Shape Variability Modeling

September 22th, 2011, Westin Harbour Castle, Toronto, Canada.

http://www.inria.fr/sophia/asclepios/events/MFCA11

Editors:

Xavier Pennec (Asclepios, INRIA Sophia-Antipolis, France)

Sarang Joshi (SCI, University of Utah, USA)

Mads Nielsen (University of Copenhagen)
<figure><img src="image_4.png" /></figure>
Preface

Computational anatomy is an emerging discipline at the interface of geometry, statistics and image analysis which aims at modeling and analyzing the biological shape of tissues and organs. The goal is to estimate representative organ anatomies across diseases, populations, species or ages, to model the organ development across time (growth or aging), to establish their variability, and to correlate this variability information with other functional, genetic or structural information.

The **Mathematical Foundations of Computational Anatomy (MFCA)** workshop aims at fostering the interactions between the mathematical community around shapes and the MICCAI community in view of computational anatomy applications. It targets more particularly researchers investigating the combination of statistical and geometrical aspects in the modeling of the variability of biological shapes. The workshop is a forum for the exchange of the theoretical ideas and aims at being a source of inspiration for new methodological developments in computational anatomy. A special emphasis is put on theoretical developments, applications and results being welcomed as illustrations. Following the successful first edition of this workshop in 2006¹ and second edition in New-York in 2008², the third edition was held in Toronto on September 22 2011³.

Contributions were solicited in Riemannian and group theoretical methods, geometric measurements of the anatomy, advanced statistics on deformations and shapes, metrics for computational anatomy, statistics of surfaces, modeling of growth and longitudinal shape changes. 22 submissions were reviewed by three members of the program committee. To guaranty a high level program, 11 papers only were selected for oral presentation in 4 sessions. Two of these sessions regroups classical themes of the workshop: statistics on manifolds and diffeomorphisms for surface or longitudinal registration. One session gathers papers exploring new mathematical structures beyond Riemannian geometry while the last oral session deals with the emerging theme of statistics on graphs and trees. Finally, a poster session of 5 papers addresses more application oriented works on computational anatomy.

August 2011

Xavier Pennec
General Chair
MFCA'11

<footnote n="1">¹ http://www.inria.fr/sophia/assembleios/events/MFCA06/</footnote>

<footnote n="2">² http://www.inria.fr/sophia/assembleios/events/MFCA08/</footnote>

<footnote n="3">³ http://www.inria.fr/sophia/assembleios/events/MFCA11/</footnote>
# Organization

## Workshop Chairs

Xavier Pennec

Sarang Joshi

Mads Nielsen

(INRIA Sophia-Antipolis, France)

(SCI, University of Utah, USA)

(University of Copenhagen, Denmark)

## Program Committee

Stéphanie Allassoniere

Rachid Deriche

Stanley Durrleman

Ian L. Dryden

Tom Fletcher

James Gee

Guido Gerig

Polina Golland

Darryl Holm

Natasha Lepore

Yoshitaka Masutani

Stephen Marsland

Michael I. Miller

David Mumford

Marc Niethammer

Salvador Olmos

Sebastien Ourselin

Anand Rangarajan

Daniel Rueckert

Guillermo Sapiro

Martin Styner

Anuj Srivastava

Paul Thompson

Alain Trouvé

Carole Twining

Baba Vemuri

Francois-Xavier Vialard

William M. Wells III

Laurent Younes

(Ecole Polytechnique, France)

(INRIA, France)

(SCI, Univ. Utah, USA)

(University of South Carolina, USA)

(University of Utah, USA)

(Univ. of Pennsylvania, USA)

(University of Utah, USA)

(CSAIL, MIT, USA)

(Imperial College London, UK)

(USC Keek School of Medicine, USA)

(University of Tokyo Hospital, Japan)

(Massey University, New Zealand)

(John Hopkins University, USA)

(Brown University, USA)

(UNC Chapel Hill, USA)

(University of Saragossa, Spain)

(UCL, UK)

(University of Florida, USA)

(Imperial College London, UK)

(University of Minnesota, USA)

(UNC Chapel Hill, USA)

(Florida State University, USA)

(UCLA, Los-Angeles, USA)

(ENS-Cachan, France)

(University of Manchester, UK)

(University of Florida, USA)

(Imperial College London, UK)

(CSAIL, MIY, Boston, USA)

(John Hopkins University, USA)
# Table of Contents

## LDDMM and Longitudinal Deformations

Direct LDDMM of Discrete Currents with Adaptive Finite Elements ... 1
*Andreas Günther, Hans Lamecker, and Martin Weiser*

Frame to Frame Diffeomorphic Motion Analysis from Echocardiographic
Sequences ... 15
*Zhijun Zhang, David J. Sahn, and Xubo Song*

An Improved Estimator of GRID Model for Representing Large
Diffeomorphic Deformations ... 25
*Qian Xie and Anuj Srivastava*

## Beyond Riemannian Geometry

Teichmüller Shape Descriptor and its Application to Alzheimer's
Disease Study ... 38
*Wei Zeng, Rui Shi, Yalin Wang, and Xianfeng David Gu*

On a Geodesic Equation for Planar Conformal Template Matching ... 52
*Stephen Marsland, Robert McLachlan, Klas Modin, and Matthew
Perlmutter*

Geodesics, Parallel Transport & One-parameter Subgroups for
Diffeomorphic Image Registration ... 64
*Marco Lorenzi and Xavier Pennec*

## Statistical Analysis on Riemannian Manifolds

Geodesic Regression on Riemannian Manifolds ... 75
*P. Thomas Fletcher*

Varying Coefficient Models for Modeling Diffusion Tensors Along
White Matter Bundles ... 87
*Ying Yuan, Hongtu Zhu, J. S. Marron, John H. Gilmore, and
Martin Styner*
## Statistical Graphs and Trees

<table><tr><td>Matching and Classification of Images Using The Space of Image Graphs<br><i>Dohyung Seo, Jeffrey Ho, and Baba C. Vemuri</i></td><td>99</td></tr><tr><td>Detecting Long Distance Conditional Correlations Between Anatomical Regions Using Gaussian Graphical Models<br><i>Stéphanie Allassonnière, Pierre Jolivet and Christophe Giraud</i></td><td>111</td></tr><tr><td>An Airway Tree-shape Model for Geodesic Airway Branch Labeling<br><i>Aasa Feragen, Pechin Lo, Vladlena Gorbunova, Mads Nielsen, Asger Dirksen, Joseph M. Reinhardt, François Lauze, and Marleen de Bruijne</i></td><td>123</td></tr></table>

## Poster Session

<table><tr><td>Deformable Organisms and Error Learning for Brain Segmentation<br><i>Gautam Prasad, Anand A. Joshi, Albert Feng, Marina Barysheva, Katie L. McMahon, Greig I. de Zubicaray, Nicholas G. Martin, Margaret J. Wright, Arthur W. Toga, Demetri Terzopoulos, and Paul M. Thompson</i></td><td>135</td></tr><tr><td>Selecting Scales by Multiple Kernel Learning for Shape Diffusion Analysis<br><i>Umberto Castellani, Aydın Ulaş, and Vittorio Murino, Marcella Bellani, Gianluca Rambaldelli, Michele Tansella, Paolo Brambilla</i></td><td>148</td></tr><tr><td>Probabilistic Modeling of Landmark Distances and Structure for Anomaly-proof Landmark Detection<br><i>Shouhei Hanaoka, Yoshitaka Masutani, Mitsutaka Nemoto, Yukihiro Nomura, Takeharu Yoshikawa, Naoto Hayashi, Naoki Yoshioka, and Kuni Ohtomo</i></td><td>159</td></tr><tr><td>Hippocampal Morphometry Study by Automated Surface Fluid Registration and its Application to Alzheimer's Disease<br><i>Jie Shi, Yuting Wang, Paul M. Thompson, Yalin Wang</i></td><td>170</td></tr><tr><td>A New Riemannian Setting for Surface Registration<br><i>Martin Bauer, Martins Bruveris</i></td><td>182</td></tr></table>

<table><tr><td>Author Index</td><td>194</td></tr></table>
Direct LDDMM of Discrete Currents
with Adaptive Finite Elements

Andreas Günther, Hans Lamecker, and Martin Weiser

Zuse Institute Berlin, Takustraße 7, D-14195 Berlin-Dahlem, Germany

guenther@zib.de

**Abstract.** We consider Large Deformation Diffeomorphic Metric Mapping of general *m*-currents. After stating an optimization algorithm in the function space of admissible morph generating velocity fields, two innovative aspects in this framework are presented and numerically investigated: First, we spatially discretize the velocity field with conforming adaptive finite elements and discuss advantages of this new approach. Second, we directly compute the temporal evolution of discrete *m*-current attributes.

# 1 Introduction

The Large Deformation Diffeomorphic Metric Mapping (LDDMM) approach initiated in [4, 18] has attracted considerable attention over the last few years in medical imaging. It allows to match highly deformed objects and as such is capable of performing inter-individual registration. LDDMM constructs a space mapping by evolving a displacement field along a velocity field, we call *wind*. Depending on the regularity of the wind, either *diffeomorphisms* [1, 12] or *homeomorphisms* [21] of the embedded space can be obtained. Thus, it provides a basis for many applications of anatomical shape analysis, where a one-to-one correspondence between different spatial objects is required.

The LDDMM technique is commonly applied for matching *currents* [5]. Currents provide a unified mathematical description of geometrical objects of dimension 0 (points), 1 (curves), 2 (surfaces) or 3 (volumes) [7, 14] which are embedded in $\mathbb{R}^3$. The spaces of *m*-currents are linear and equipped with an inner product and hence are a suitable tool for statistical shape analysis [5]. The induced norm provides a similarity measure for matching of source and target objects.

In [3] a particle-mesh method has been applied to 1-currents in 2D. Therein curves of same topology represented by a parameterization were approximated by a finite point set without any tangential information regardless of some geometrical level of detail for matching. Although attaching a momentum vector at each vertex of the source shape can be proven to be the optimal wind parameterization, Cotter proposes to use cubic B-splines on a fixed grid. In an enclosed efficiency discussion this sub-optimal parameterization is still competitive due to the simpler wind structure. However the numerical scheme there relies on equidistant cartesian grids with constant diffusity in order to apply FFT techniques.
In contrast to [3] we apply in 3D the *Orthogonal Matching Pursuit* (OMP) proposed in [6] to obtain a sparse representation for general *m*-currents at a given geometric resolution (spectral length) in terms of a sum of discrete *Dirac delta m-currents*. Currently, the LDDMM evolution of this representation is only done indirectly via an approximative scheme [5, Rem. 4.13] as depicted in Figure 1 (left) for *m* = 2, because it would otherwise require the computation of the Jacobian of the diffeomorphism, which is a challenge when discretizing the wind using Gaussian kernels [8, 9, 19]. Furthermore, one loses the connectivity of the input mesh structure in this case. However, this is not a significant problem, since the connectivity can be recovered by applying the final displacement field to the input meshes afterwards.

<figure><img src="image_5.png" /><figcaption>**Fig. 1.** Usual discrete 2-current deformation (left) versus the direct approach (right)</figcaption></figure>

In this paper, we study the direct evolution of Dirac delta $m$-currents (right of Figure 1). We show that the direct approach allows to uniformly treat $m$-currents for $m = 0, \dots, 3$ (Sec. 2), which to the best of our knowledge has not been shown before. We show also how to compute the Jacobian in this setting by using finite elements (FEM) to discretize the wind in the LDDMM framework (Sec. 3). Since the compactly supported basis functions are fixed in space the computation is significantly simplified. Although equidistant grids can of course be incorporated, we also in contrast to [3] consider locally refined meshes and exploit the decoupling of the wind and current discretization for an adaptive current representation, giving a significant reduction of degrees of freedom (Sec. 4). Based on our results, the increased spatial flexibility of adaptive FEM may be exploited in the future by implementing hierarchical schemes, error estimators and non-constant spatial anisotropic diffusivity.

# 2 Continuous Matching Problem

For given shapes $S, T \subset \mathbb{R}^3$ we aim at constructing a sufficiently smooth bijection $\phi: S \to T$ such that the distance between $\phi(S)$ and $T$ is minimal. Here we fix the formalisms to describe the matching problem as an optimization task.

## 2.1 Currents

*Currents* are mathematical tools suited for describing geometrical objects such as points, space curves, surfaces and volumes embedded in $\mathbb{R}^3$. Their precise
definition from [7, 14] requires notation for differential forms taken from [15]. Let for $m = 0, 1, 2, 3$ the set $\mathcal{D}^m := C_c^\infty(\mathbb{R}^3, \Lambda^m\mathbb{R}^3)$ denote the vector space of all $C^\infty$ differential $m$-forms on $\mathbb{R}^3$ with compact support. A $m$-current is an element of $\mathcal{D}_m$, the dual space of $\mathcal{D}^m$. The elementary Dirac delta $m$-currents $\delta_x^{u_1 \wedge \dots \wedge u_m} \in \mathcal{D}_m$ act on $\omega \in \mathcal{D}^m$ as

$$ \delta_x^{u_1 \wedge \dots \wedge u_m}(\omega) = \omega(x)(u_1 \wedge \dots \wedge u_m). $$

Following the discussion in [5, Sect. 1.5.1] it turns out that for the purpose of matching currents the testspace of all $C^\infty$ differential $m$-forms is not suited due to a missing bound in variation. Moreover the space $\mathcal{D}^0$ can be identified with scalar $C_c^\infty$ functions on $\mathbb{R}^3$. For $m = 1$ and $m = 2$ the space $\mathcal{D}^m$ is isomorph to the space of vector-valued $C_c^\infty$ functions from $\mathbb{R}^3$ to $\mathbb{R}^3$. An element of $\mathcal{D}^3$ can be written as a scalar $C_c^\infty$ function times the determinant form on $\mathbb{R}^3$.

Both aspects motivate the use of Reproducible Kernel Hilbert Spaces (RKHS) $W^m$ as testspaces.

**Definition 1.** Let $d_m = 1$ for $m \in \{0, 3\}$ and $d_m = 3$ for $m \in \{1, 2\}$. For $m = 0, 1, 2, 3$ let $W^m$ denote the dense span of $d_m$-vectorfields of the form $\omega(x) = k_m(x, y)a$, where $x, y \in \mathbb{R}^3$, $a \in \mathbb{R}^{d_m}$ and $k_m(x, y) = \exp(-\|x - y\|^2/\sigma_m^2)$. The space $W^m$ can be equipped with the inner product $\langle k_m(\cdot, x)a, k_m(\cdot, y)b \rangle_{W^m} = a^*k_m(x, y)b$. Here the symbol $*$ denotes the transpose operation.

An $m$-current in $\mathbb{R}^3$ is a continuous linear functional on $W^m$. $W_m$ denotes the vector space of all $m$-currents in $\mathbb{R}^3$.

For $x \in \mathbb{R}^3$ and attribute $a \in \mathbb{R}^{d_m}$ we define the elementary Dirac delta $m$-currents $\delta_x^a \in W_m$ acting on $\omega \in W^m$ as $\delta_x^a(\omega) = a^*\omega(x)$.

The above inner product induces a norm on $W^m$, which can be computed efficiently via FGT even for a large number of linear combinations of the above basis functions. The chosen Gaussian kernel $k_m$ can be considered as Green's function for some differential operator $L_W$ (see [1, 6, 8]). With the above objects at hand the Riesz representation theorem provides a unique operator $K_W^m : W_m \to W^m$ reflecting the canonical isometry between $W^m$ and $W_m$ defined via

$$ \langle K_W^m f, g \rangle_{W^m} = \langle f, g \rangle_{W_m, W^m} = f(g) $$

for all $f \in W_m$ and $g \in W^m$. It provides for the $m$-current $S^m$ the Riesz representant $K_W^m S^m$ as unique $d_m$-vectorfield on $\mathbb{R}^3$.

## 2.2 Homeomorphisms and Diffeomorphisms

Let $\Omega$ be an open bounded subset of $\mathbb{R}^3$ and consider functions $v_t : \bar{\Omega} \to \mathbb{R}^3$ that vanish on $\partial\Omega$. For given final time $T > 0$ and a time-dependent wind $v = (v_t)_{t \in [0, T]}$ we consider the temporal evolution of the identity map

$$ \frac{\partial \phi_t^v}{\partial t} = v_t(\phi_t^v) \text{ with } \phi_0^v(x) = x. \qquad (1) $$
In what follows it will be useful to define the trajectory $x_t := \phi_t^v(x)$ for some fixed space point $x \in \mathbb{R}^3$ and the map $\phi_{st}^v := \phi_t^v \circ (\phi_s^v)^{-1}$, describing the movement of a particle starting in $x$ at time $s$ towards $\phi_{st}^v(x)$ at time $t$. It is well known (see [21, Thm. C.3]), that (1) is uniquely solvable when for some $x_0 \in \Omega$ the integral $\int_0^T \|v_t(x_0)\|_{\mathbb{R}^3} + \text{Lip}(v_t) dt$ is bounded. Furthermore its solution $\phi_t^v : \mathbb{R}^3 \to \mathbb{R}^3$ is a homeomorphism of $\Omega$ for all times $t \in [0, T]$. Under more restrictive assumptions onto the spatial smoothness of the wind, i.e. $v_t \in C_0^1(\Omega, \mathbb{R}^3) \forall t \in [0, T]$ and $\int_0^T \|v_t\|_{1,\infty} dt < \infty$ the unique solution of (1) is even a diffeomorphism of $\Omega$ for all times $t \in [0, T]$ (see [21, Thm. 8.7]). For convenience we look for the wind $v_t$ in some Hilbert space $V$. Such spaces can be constructed by defining inner products associated to differential operators. Let therefore $L : V \to L^2(\mathbb{R}^3)$ be a differential operator and equip the Hilbert space $V$ with the inner product $\langle v_t, g \rangle_V = \langle Lv_t, Lg \rangle_{L^2} = \langle L^*Lv_t, g \rangle_{V^*,V}$. Here $L^*$ denotes the adjoint operator. For this work we use

$$S := L^*L = (-\operatorname{div}(\sigma_V^2 \nabla) + I)^k = (-\sigma_V^2 \Delta + I)^k \quad (2)$$

and $k=1$ or $k=2$ giving the Sobolev spaces $H^k$ (see [8]). For given $f \in V^*$ we consider solutions $v_t \in V$ of $Sv_t = f$ with homogeneous Dirichlet boundary conditions for $v_t$ (and $v'_t$ if $k=2$). Here the real parameter $\sigma_V > 0$ balances between smoothing and data fitting of the right hand side $f$. For other choices of $L^*L$ and boundary conditions see [12]. Dealing with natural boundary conditions is also possible, but requires a sufficiently large domain to keep all trajectories therein. Analogous to $K_W^m$ we introduce the isometry operator $K_V : V^* \to V$. A mathematically equivalent approach of constructing $V$ consists in defining $K_V$ via the Green's function $k_V(x, y)$ of $L^*L$, see for instance [8, 9, 19, 20].

## 2.3 Diffeomorphic Deformation of Currents

For $m=0, 1, 2, 3$ let currents $S^m \in W_m$ be given. Let $\phi$ denote a diffeomorphism on $\mathbb{R}^3$ and $d_x\phi$ the Jacobian of $\phi$ at $x$. The pushforward $\phi_\sharp(S^m) \in W_m$ of $S^m$ under $\phi$ is rigorously defined in [15] via the pullback of differential forms. For our purpose it is sufficient to mention that if $S^m$ is associated to a submanifold in $\mathbb{R}^3$ its pushforward $\phi_\sharp(S^m)$ under $\phi$ corresponds to the deformed submanifold $\phi(S^m)$. This important property justifies to write also $\phi(S^m) \in W_m$. The explicitly calculated pushforwards for elementary Dirac delta $m$-currents taken from [5, Table 1.2] are given in Table 1.

<table><caption>**Table 1.** Pushforwards of Dirac delta $m$-currents under $\phi$</caption><tr><td>m = 0</td><td>d<sub>0</sub> = 1</td><td>c &isin; &#8477;</td><td>&phi;<sub>&#8455;</sub>(&delta;<sub>x</sub><sup>c</sup>) = &delta;<sub>&phi;<sub>x</sub></sub><sup>c</sup></td></tr><tr><td>m = 1</td><td>d<sub>1</sub> = 3</td><td>&tau; &isin; &#8477;<sup>3</sup></td><td>&phi;<sub>&#8455;</sub>(&delta;<sub>x</sub><sup>&tau;</sup>) = &delta;<sub>&phi;<sub>x</sub></sub><sup>d<sub>x</sub>&phi;(&tau;)</sup></td></tr><tr><td>m = 2</td><td>d<sub>2</sub> = 3</td><td>n &isin; &#8477;<sup>3</sup></td><td>&phi;<sub>&#8455;</sub>(&delta;<sub>x</sub><sup>n</sup>) = &delta;<sub>&phi;<sub>x</sub></sub><sup>|d<sub>x</sub>&phi;|d<sub>x</sub>&phi;<sup>-*(n)</sup></sup></td></tr><tr><td>m = 3</td><td>d<sub>3</sub> = 1</td><td>&rho; &isin; &#8477;</td><td>&phi;<sub>&#8455;</sub>(&delta;<sub>x</sub><sup>&rho;</sup>) = &delta;<sub>&phi;<sub>x</sub></sub><sup>|d<sub>x</sub>&phi;|&rho;</sup></td></tr></table>
Let some wind $v$ be given and consider the family $(\phi_t^v)_t$ of diffeomorphisms generated via (1). The following theorem describes the direct evolution of $m$-current attributes $a_m \in \mathbb{R}^{d_m}$ under $(\phi_t^v)_t$, where $'$ denotes the time derivative.

**Theorem 1.** The pushforwards of $\delta_{x_0}^{c_0}$, $\delta_{x_0}^{\tau_0}$, $\delta_{x_0}^{n_0}$ and $\delta_{x_0}^{\rho_0}$ under $\phi_s^v$ satisfying (1) are $\delta_{x_s}^{c_0}$, $\delta_{x_s}^{\tau_s}$, $\delta_{x_s}^{n_s}$ and $\delta_{x_s}^{\rho_s}$. Their components are given via the ODEs

$$
\begin{align*}
x'_t &= v_t(x_t) & \text{with } x(0) &= x_0 \\
\tau'_t &= (d_{x_t} v_t) \tau_t & \text{with } \tau(0) &= \tau_0 \\
n'_t &= n_t \operatorname{tr}(d_{x_t} v_t) - (d_{x_t} v_t)^* n_t & \text{with } n(0) &= n_0 \\
\rho'_t &= \rho_t \operatorname{tr}(d_{x_t} v_t) & \text{with } \rho(0) &= \rho_0.
\end{align*}
$$

*Proof.* Abbreviating $J_t = d_{x_0}\phi_t^v$ and $A_t = d_{x_t}v_t$ there holds (see [1]) $J'_t = A_t J_t$ with $J(0) = I_3$. Observing the evolution of the Wronskian [13, Thm. 2.14] or via Jacobi's formula one obtains

$$
|J_t'| = |J_t|\operatorname{tr}(J_t^{-1} J_t') = |J_t|\operatorname{tr}(J_t^{-1} A_t J_t) = |J_t|\operatorname{tr}(A_t),
$$

where tr(A) denotes the trace of a matrix A and $A^{-*} = (A^{-1})^*$. Now from Table 1 we read out

$$
x_t = \phi_t^v(x_0), \quad \tau_t = J_t \tau_0, \quad n_t = |J_t| J_t^{-*} n_0, \quad \rho_t = |J_t| \rho_0.
$$

Differentiation of the above equations with respect to $t$ yields

$$
\begin{align*}
x_t' &= \phi_t^v(x_0)' = v_t(\phi_t^v(x_0)) = v_t(x_t) \\
\tau_t' &= J_t' \tau_0 = A_t J_t \tau_0 = A_t \tau_t \\
n_t' &= |J_t'| J_t^{-*} n_0 + |J_t| (J_t^{-*})' n_0 = |J_t| \operatorname{tr}(A_t) J_t^{-*} n_0 - |J_t| A_t^* J_t^{-*} n_0 \\
&= n_t \operatorname{tr}(A_t) - A_t^* n_t \\
\rho_t' &= |J_t'| \rho_0 = |J_t| \operatorname{tr}(A_t) \rho_0 = \rho_t \operatorname{tr}(A_t),
\end{align*}
$$

which proves the assertion. $\square$

*Remark 1.* The authors emphasize the striking advantage that Theorem 1 enables to find the final position and attribute of a Dirac delta $m$-current without computing the Jacobian of the deformation. The appearing ODEs only involve the Jacobian of the velocity fields, which will be given in a closed form in any case.

## 2.4 Optimization Problem in Function Space

Let source $S^m \in W_m$ and target current $T^m \in W_m$ be given for $m = 0, \dots, 3$. For given wind $v$ we define the deformed current $S_t^m := \phi_t^v(S^m)$ at time $t$. Matching means the minimization of the distance of the deformed source current at final time $S_T^m$ with its target current $T^m$, i.e. minimizing the dual norm $\|\phi_T^v(S^m) - T^m\|_{W_m} = \|S_T^m - T^m\|_{W_m}$ in the space of $m$-currents.
Given a regularization parameter $\gamma > 0$ and matching weights $\omega_m \ge 0$ we consider for $v \in L^2([0, T], V)$ the following optimization problem:

$$J(v) := \gamma \int_{0}^{T} \|v_t\|_{V}^{2} dt + \sum_{m=0}^{3} \omega_{m} \|\phi_{T}^{v}(S^{m}) - \mathcal{T}^{m}\|_{W_{m}}^{2} \rightarrow \min . \quad (3)$$

Here the first summand involves the kinetic energy of the wind. The existence of a solution for (3) is proven in [10], however it is generally not unique [2]. Following [9] the gradient of $J$ in $L^2([0, T], V)$ at fixed $v$ is given by $(\nabla J)_t = 2\gamma v_t + 2K_V(f_t)$, where $f_t \in V^*$ is defined by

$$f_t(u) = \sum_{m=0}^{3} \omega_m \langle S_t^m, \nabla(K_W^m(S_T^m - \mathcal{T}^m) \circ \phi_{tT}^v)^* u \rangle_{W_m, W^m} \quad \forall u \in V.$$

For further discussion concerning the choice of the gradients metric we refer the reader to [1]. With the above quantities at hand one is able to state a steepest descent optimization algorithm in the function space of velocity fields $v$.

# 3 Discrete Matching Problem

## 3.1 Discretization of the Wind by Finite Elements (FE)

In the field of optimal current matching mainly wind discretizations of the form

$$v_t(x) = \sum_j k_V(x_{j,t}, x) \alpha_{j,t} \quad (4)$$

have been considered. Here $\alpha_{j,t} \in \mathbb{R}^3$ are the time-dependent momentum vectors and $k_V$ denotes a Gaussian kernel with some global kernel parameter $\sigma_V > 0$, describing the coherent movement of neighboring particles. In order to apply *Fast Gauss Transform* (FGT) for efficient evaluation, $\sigma_V$ is necessarily a constant. Although (4) can be proven to be the optimal wind parameterization, the spatial movement of non-compactly supported basis functions along trajectories $x_{j,t}$ may cause numerical difficulties. Too small distances between them cause a redundant or badly conditioned description of the velocity field while the absence of trajectories in a part of the domain produces almost no wind there for small kernel sizes. The trajectory density varies during optimization and hence is difficult to control. Because the trajectories' starting points are the spatial components of the Dirac delta source currents the number of trajectories is fixed and hence a notion of adaptivity for the velocity field can hardly be introduced. Finally, as mentioned in Sect. 2.2, $C^\infty$ smoothness is not required to solve the evolution equation.

In [16] and [17], some of the above mentioned drawbacks are overcome by incorporating multiple kernel shapes at different scales $\sigma_V$.

Similar to the particle-mesh method proposed in [3], we follow another sub-optimal approach completely decoupling the discretization of the space of $m$-currents $W_m$ from the spatial velocity space $V$. Keeping in mind that fast point
evaluation of the wind is essential for performance, we consider adaptive hexahedral grids for $\Omega$ with hanging nodes saved as an octree. Over such hexahedral grids we construct either $C^1$ conforming Hermite finite elements of third order or simpler $C^0$ conforming Lagrange finite elements of first order. The wind for fixed time $t \in [0, T]$ in the FE basis $\{\varphi_j\}_j$ takes the form

$$v_t(x) = \sum_{j=1}^{n} \varphi_j(x) \alpha_{j,t} . \qquad (5)$$

In contrast to radial basis functions, locally constant functions are contained in the ansatz space and allow to represent local or even global translations with few degrees of freedom (DOF). Due to the compactly supported basis functions there is no need for an approximate evaluation like FGT with further unknown tolerance parameters. Since the basis functions are fixed in space, the underlying mesh provides a natural clustering which can be exploited via a smart parallel octree search algorithm for point evaluation. Furthermore this approach provides a multilevel wind hierarchy with a fraction of DOFs on the coarsest mesh level completely decoupled from the m-current discretization. These advantages also appear in the particle-mesh method with tensor-products of cubic B-splines for instance. But since we do not apply FFT for wind evaluation, we are more flexible with adaptive meshes and do not require a box domain. Moreover non-constant anisotropic diffusity $\sigma_V(x) \in \mathbb{R}^{3 \times 3}$ may be incorporated in future.

A difficulty arises in the computation of the $L^2([0, T], V)$-gradient. It permanently involves the solution of a second ($k=1$) or fourth ($k=2$) order elliptic PDE in every time-step and every iteration. It is clear that one should employ suited preconditioners and / or multigrid solvers. Using existing FE libraries limits implementation overhead. We chose libMesh [11], which provides conforming $C^1$ finite elements on adaptive hexahedral meshes.

The development of adaptive mesh refinement is beyond the scope of this work. Here, we provide a proof of concept that adaptive grids can easily be incorporated. Therefore we simply geometrically refine near $S^m \cup T^m$ considered as subsets of $\mathbb{R}^3$ equally for all times. More sophisticated error indicators suggesting refinements could be the scalar fields $|K_W^m(S_T^m - T^m)| \circ \phi_{tT}^v$, $|v_t|$ or $|L^*L v_t|$. The latter one measures the smoothness of $v_t$. Moreover thinking of hierarchical error estimators one could compute $\|v_t^1 - I_h v_t^1\|_{L^2(Q)}$ or even $\|v_t^1 - v_t^0\|_{L^2(Q)}$ on hexahedrons $Q$, where $v^i$ denote the numerical solutions for the $C^i$ conforming FE discretization and $I_h$ is the usual Lagrange interpolation operator.

All appearing ODEs are numerically integrated via the explicit method of Heun on an equidistant decomposition of the time interval $[0, T]$.

## 3.2 Current Compression and Direct Evolution

For approximating a m-current $S^m \in W_m$ as $\hat{S}^m = \sum_{i=1}^{s_m} \delta_{x_i}^{a_i} \in W_m$ we use the Orthogonal Matching Pursuit (OMP) proposed in [6]. This method iteratively selects the most important points $x_i$ and computes corresponding attributes $a_i$ (i.e. $c_i$, $\tau_i$, $n_i$, $\rho_i$) of a general m-current via a greedy algorithm. It has the advantage of compressing the current information for a characteristical spectral
length $\sigma_m > 0$ towards a fraction. This enables the design of highly efficient numerical solution algorithms. The approximation error in OMP is controlled by a threshold parameter and the grid size of a uniform testgrid.

The obvious drawback of loosing the connectivity between vertices (for $m \ge 1$) can be compensated by applying the obtained optimal diffeomorphism to all connected vertices whenever it is required. This only requires one additional forward flow computation at the end.

In [19, Sec. 3.2] two methods to deform a 2-current $\hat{S}^2$ under a family of diffeomorphisms $(\phi_t^v)_{t \in \mathbb{R}}$ are described. In contrast to all previous work, we will pursue the direct approach motivated by Theorem 1. For 2-currents, only 1 instead of 3 trajectories is needed to evolve the normal $n_0$ (Fig. 1). In general, the direct approach requires only one trajectory per attribute, hence decreasing the number of variables in the computation, whereas in the indirect case an artificial $m$-simplex with $m+1$ vertices is attached.

*Remark 2.* To quote Rem. 4.13 in [5] the direct evolution of current attributes is closer to the analytical concept of currents and is particularly suited for OMP, where no connectivity between the points is provided. But [5] indicates the need of Jacobi matrices (as they arise in Theorem 1) as a disadvantage for numerical implementation. At least for the gradient computation in the next section we benefit from the simpler structure of $v_t$ in (5), which in Lemma 2 enables easy evaluations of $d_{x_t}v_t = \sum_{j=1}^n \alpha_{j,t} \nabla \varphi_j(x_t)^*$ and hence $\text{tr}(d_{x_t}v_t)$, $(d_{x_t}v_t)w$ and $(d_{x_t}v_t)^*w$ for a vector $w \in \mathbb{R}^3$. Note that all sums over $j$ are local sums due to the compact support of the basis functions $\varphi_j$.

### 3.3 Discrete Optimization Problem and its Gradient

Let $a(\cdot, \cdot)$ denote the bilinear form corresponding to the elliptic differential operator $S$ from (2). We define the sparse symmetric, positive definite matrix $S := [a(\varphi_i, \varphi_j)]_{i,j=1}^n$ using the FE basis $\{\varphi_j\}_j$ from Sect. 3.1. Moreover we introduce the block vectors $\alpha_t := [\alpha_{i,t}]_{i=1}^n$, $\mathbf{x}_t := \mathbf{x}_t^m := [x_{i,t}]_{i=1}^{s_m}$ and $\mathbf{a}_t := \mathbf{a}_t^m := [a_{i,t}]_{i=1}^{s_m}$. This notation allows to write the matching terms as

$$E^m = E^m(\mathbf{x}_T, \mathbf{a}_T) = \|\phi_T^v(\hat{S}^m) - \hat{T}^m\|_{W_m}^2 = \|\sum_{i=1}^{s_m} \delta_{x_{i,T}}^{a_{i,T}} - \sum_{j=1}^{r_m} \delta_{y_j}^{b_j}\|_{W_m}^2.$$

Finally the discrete form of the current matching problem (3) is

$$\hat{J}(\alpha_t) := \gamma \int_0^T \|v_t(\alpha_t)\|_V^2 dt + \sum_{m=0}^3 \omega_m \|\phi_T^v(\hat{S}^m) - \hat{T}^m\|_{W_m}^2 \to \min$$

or even shorter via (5) and $\|v_t(\alpha_t)\|_V^2 = a(v_t(\alpha_t), v_t(\alpha_t)) = \alpha_t^* \mathbf{S} \alpha_t$

$$\hat{J}(\alpha_t) = \gamma \int_0^T \alpha_t^* \mathbf{S} \alpha_t dt + \sum_{m=0}^3 \omega_m E^m(\mathbf{x}_T, \mathbf{a}_T) \to \min . \quad (6)$$

The analytical computation of the gradient at given $\alpha_t$ becomes manageable though the simpler wind representation (5). Numerically the computation is
more involved due to presence of Hessians of basis functions. But these are
easily provided via the already mentioned libMesh library.

**Theorem 2.** The gradient of $\hat{J}$ in the $L^2$-metric is

$$ (\nabla \hat{J})_t = 2\gamma \mathbf{S} \alpha_t + \sum_{m=0}^{3} \omega_m ((\varphi_t^m)^* \eta_t^m + (\partial_\alpha \mathbf{g}_t^m)^* \zeta_t^m), \quad (7) $$

$$ \text{with} \quad \varphi_t^m = [\varphi_j(x_{i,t}) I_{d_m}]_{i=1 \dots s_m; j=1 \dots n} $$

$$ \zeta_t^m = \nabla_{\mathbf{a}_T} E^m + \int_t^T (\partial_{\mathbf{a}} \mathbf{g}_s^m)^* \zeta_s^m ds $$

$$ \eta_t^m = \nabla_{\mathbf{x}_T} E^m + \int_t^T (\partial_{\mathbf{x}} \mathbf{g}_s^m)^* \zeta_s^m ds . $$

The proof is postponed to the appendix. The remaining quantities $\nabla_{\mathbf{x}_T} E^m$, $\nabla_{\mathbf{a}_T} E^m$, $\partial_\alpha \mathbf{g}_t^m$, $\partial_{\mathbf{x}} \mathbf{g}_t^m$ and $\partial_{\mathbf{a}} \mathbf{g}_t^m$ from Theorem 2 for each $m$ are specified in the next two lemmas.

**Lemma 1.** Let $f_m(x) = \sum_{i=1}^{s_m} k_m(x_{i,T}, x) a_{i,T} - \sum_{j=1}^{r_m} k_m(y_j, x) b_j$. There hold

$$ \nabla_{\mathbf{x}_T} E^m = [2(d_{x_{i,T}} f_m(x_{i,T}))^* a_{i,T}]_{i=1}^{s_m} \quad \text{and} \quad \nabla_{\mathbf{a}_T} E^m = [2f_m(x_{i,T})]_{i=1}^{s_m} . $$

*Proof.*

$$ (\partial_{x_T} E^m) \eta = 2 [(\partial_{x_T} \sum_{i=1}^{s_m} \delta_{x_{i,T}}^{a_{i,T}}) \eta] (f_m) = 2 a_T^* (d_{x_T} f_m(x_T)) \eta $$

$$ \nabla_{x_T} E^m = 2(d_{x_T} f_m(x_T))^* a_T $$

$$ = 2 \left( \sum_{i=1}^{s_m} (\nabla_2 k_m(x_{i,T}, x_T)) a_{i,T}^* - \sum_{j=1}^{r_m} (\nabla_2 k_m(y_j, x_T)) b_j^* \right) a_T $$

$$ (\partial_{a_T} E^m) \eta = 2 [(\partial_{a_T} \sum_{i=1}^{s_m} \delta_{x_{i,T}}^{a_{i,T}}) \eta] (f_m) = 2 \eta^* f_m(x_T) . $$

□

**Lemma 2.** For $\mathbf{g}_t^m$ in (11) their sparse Jacobians are given via

$$ \mathbf{g}_t^0 = \mathbf{0} $$

$$ \partial_\alpha \mathbf{g}_t^1 = [(\tau_{i,t}^* \nabla \varphi_j(x_{i,t})) I_3]_{i=1 \dots s_1; j=1 \dots n} $$

$$ \partial_\alpha \mathbf{g}_t^2 = [n_{i,t} \nabla \varphi_j(x_{i,t})^* - \nabla \varphi_j(x_{i,t}) n_{i,t}^*]_{i=1 \dots s_2; j=1 \dots n} $$

$$ \partial_\alpha \mathbf{g}_t^3 = [\rho_{i,t} \nabla \varphi_j(x_{i,t})^*]_{i=1 \dots s_3; j=1 \dots n} $$

$$ \partial_{\mathbf{x}} \mathbf{g}_t^1 = \mathrm{diag} \left[ \sum_{j=1}^{n} \alpha_{j,t} \tau_{i,t}^* H_{\varphi_j}(x_{i,t}) \right]_{i=1}^{s_1} $$

$$ \partial_{\mathbf{x}} \mathbf{g}_t^2 = \mathrm{diag} \left[ \sum_{j=1}^{n} n_{i,t} (\alpha_{j,t}^* H_{\varphi_j}(x_{i,t})) - (\alpha_{j,t}^* n_{i,t}) H_{\varphi_j}(x_{i,t}) \right]_{i=1}^{s_2} $$

$$ \partial_{\mathbf{x}} \mathbf{g}_t^3 = \mathrm{diag} \left[ \sum_{j=1}^{n} \rho_{i,t} \alpha_{j,t}^* H_{\varphi_j}(x_{i,t}) \right]_{i=1}^{s_3} $$

$$ \partial_{\tau} \mathbf{g}_t^1 = \mathrm{diag} \left[ \sum_{j=1}^{n} \alpha_{j,t} \nabla \varphi_j(x_{i,t})^* \right]_{i=1}^{s_1} $$
$$ \partial_{\mathbf{n}} \mathbf{g}_t^2 = \mathrm{diag} \left[ \sum_{j=1}^{n} (\alpha_{j,t}^* \nabla \varphi_j(x_{i,t})) I_3 - \nabla \varphi_j(x_{i,t}) \alpha_{j,t}^* \right]_{i=1}^{s_2} $$

$$ \partial_{\rho} \mathbf{g}_t^3 = \mathrm{diag} \left[ \sum_{j=1}^{n} \alpha_{j,t}^* \nabla \varphi_j(x_{i,t}) \right]_{i=1}^{s_3}, $$

where $H_{\varphi_j}(x_{i,t})$ denote the Hessian of $\varphi_j$ at $x_{i,t}$.

*Proof.* The proof for all cases of $m$ can easily be adapted from the case $m = 2$. For this choice the derivatives of $\mathbf{g}_t^2$ follow from direct calculations starting with

$$ \begin{aligned} \mathbf{g}_t^2 &= \mathbf{g}^2(\boldsymbol{\alpha}_t, \mathbf{x}_t, \mathbf{n}_t) = \mathrm{diag} \left[ \mathrm{tr}(d_{x_{i,t}} v_t) I_3 - (d_{x_{i,t}} v_t)^* \right]_{i=1}^{s_2} \mathbf{n}_t \\ &= \mathrm{diag} \left[ \sum_{j=1}^{n} (\alpha_{j,t}^* \nabla \varphi_j(x_{i,t})) I_3 - \nabla \varphi_j(x_{i,t}) \alpha_{j,t}^* \right]_{i=1}^{s_2} \mathbf{n}_t \\ &= \left[ \sum_{j=1}^{n} n_{i,t} (\alpha_{j,t}^* \nabla \varphi_j(x_{i,t})) - \nabla \varphi_j(x_{i,t}) (\alpha_{j,t}^* n_{i,t}) \right]_{i=1}^{s_2}, \end{aligned} $$

where $\mathrm{diag}[\mathbf{v}] = [\delta_{ij}v_i]_{i,j=1}^s$ for $\mathbf{v} \in \mathbb{R}^s$ and $\delta_{ij}$ denotes the Kronecker delta. $\square$

**Corollary 1.** If $\omega_m = 0$ for $m > 0$ Theorem 2 simply provides

$$ (\nabla \hat{J})_t = 2\gamma \mathbf{S} \boldsymbol{\alpha}_t + \omega_0 (\boldsymbol{\varphi}_t^0)^* \nabla_{\mathbf{x}_T} E^0 . $$

*Remark 3.* The $L^2([0, T], V)$-gradient of $\hat{J}$ is immediately obtained by applying $\mathbf{S}^{-1}$ from the left in equation (7).

# 4 Numerical Experiments

Since the numerical implementation is not yet fully tested, we postpone the investigation of the cases $m > 0$ and only consider the case $m = 0$, i.e. $\omega_i = \delta_{i0}$. The surfaces $\mathcal{S}$ and $\mathcal{T}$ are depicted in Fig. 2. To both of them we apply the OMP with $\sigma_0 = 8$ towards $\hat{\mathcal{S}}$ with $s_0 = 1746$ points and $\hat{\mathcal{T}}$ with $r_0 = 2141$ points, which are sketched as set of spheres of diameter 8 in Fig. 3.

<figure><img src="image_6.png" /><figcaption>**Fig. 2.** $\mathcal{S}$ (dark green), $\mathcal{T}$ (light red)</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Fig. 3.** $\hat{\mathcal{S}}$ (green), $\hat{\mathcal{T}}$ (red), grid</figcaption></figure>

We solve the discrete matching problem (6) on $\Omega = (0, 346.56) \times (0, 205.76) \times (0, 256.96)$ with $\gamma = 0$ and tuned $\sigma_V$. For Lagrange FE ($k = 1$) on a hexahedral
adaptive grid from Fig. 3 having 568 nodes (133 of them are hanging nodes) we choose $\sigma_V = 100$. The result $S_T$ is shown in Fig. 4. Secondly we solve problem (6) with Hermite FE ($k = 2$) on a uniform coarse grid having 120 nodes, whose solution is shown in Fig. 5 for $\sigma_V = 15$. Finally in Fig. 6 we compare our results with the software ExoShape<fnref n="1" />, generating $C^\infty$ wind via ansatz (4) with $\sigma_V = 30$. The greyscale highlights the term $\text{dist}_{x \in S_T}(x, T)$. One should keep in mind, that with $\gamma = 0$ the deformation norm vanishes in the example, but different differential operators are still present for velocity field evaluation. Although a proper comparison between all methods should use the same deformation norm $\| \cdot \|_V$, Exoshape does not easily support its change.

<figure><img src="image_8.png" /><figcaption>Fig. 4. $S_T$ for $C^0$ wind</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 5. $S_T$ for $C^1$ wind</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 6. $S_T$ for $C^\infty$ wind</figcaption></figure>

A quantitative comparison between all different wind discretizations is issue of Table 2. Therein the column DOFs denote the number of freely choosable vectors $\alpha_{j,t}$ for fixed $t$. All methods provide acceptable matches with respect to the fixed level of detail $\sigma_0 = 8$. Especially the surface $S_T$ corresponding to $C^1$ wind is also visually closest to $T$ although it is obtained via less wind DOFs compared to the approach from Exoshape. This fact stresses the potential of sub-optimal, albeit simpler wind parameterizations and decoupling the discretization of the spaces $W_m$ from $V$.

<table><caption>**Table 2.** One-sided surface distances between $S_T$ and $T$ for $m = 0$</caption><thead><tr><th>$v_t$</th><th>DOFs</th><th>mean</th><th>stddev</th><th>rms</th><th>max</th></tr></thead><tbody><tr><td>$C^0$</td><td>1·(568-133) = 435</td><td>1.18</td><td>1.14</td><td>1.64</td><td>11.33</td></tr><tr><td>$C^1$</td><td>8·120 = 960</td><td>0.95</td><td>0.91</td><td>1.32</td><td>9.16</td></tr><tr><td>$C^\infty$</td><td>1746</td><td>0.97</td><td>0.93</td><td>1.34</td><td>9.70</td></tr><tr><td>no wind, $\text{dist}_{x \in S}(x, T)$</td><td></td><td>5.08</td><td>3.87</td><td>6.39</td><td>21.05</td></tr></tbody></table>

# Appendix

In the following we prove Theorem 2.

<footnote n="1">¹ http://www-sop.inria.fr/asclepios/projects/Health-e-Child/ShapeAnalysis/</footnote>
Proof. First we consider the variation of the kinetic energy, i.e. $\omega_m = 0$ for all $m$. One directly calculates

$$ (\nabla \hat{J})_t = 2\gamma \mathbf{S} \boldsymbol{\alpha}_t . \qquad (8) $$

Let us now consider the contrary case, i.e. $\gamma = 0$. We aim to compute $\nabla_{\boldsymbol{\alpha}} E^m$ for some fixed $m$. Variation of $E = E^m$ w.r.t. $\boldsymbol{\alpha}_t$ in direction $\tilde{\boldsymbol{\alpha}}_t$ gives

$$ \tilde{E}^m = (\partial_{\mathbf{x}_T} E) \tilde{\mathbf{x}}_T + (\partial_{\mathbf{a}_T} E) \tilde{\mathbf{a}}_T . \qquad (9) $$

There holds

$$ \tilde{\mathbf{x}}_t = \int_0^t \tilde{v}_s(\mathbf{x}_s) \, ds = \int_0^t \varphi_s \tilde{\boldsymbol{\alpha}}_s \, ds . \qquad (10) $$

From Theorem 1 the evolution of $m$-current attributes can be written as

$$ \mathbf{a}'_t = \mathbf{g}(\boldsymbol{\alpha}_t, \mathbf{x}_t, \mathbf{a}_t) = \mathbf{g}_t \quad \text{with} \quad \mathbf{a}(0) = \mathbf{a}_0 . \qquad (11) $$

Its variation in direction $\tilde{\boldsymbol{\alpha}}_t$ satisfies

$$ \tilde{\mathbf{a}}'_t = (\partial_{\boldsymbol{\alpha}} \mathbf{g}_t) \tilde{\boldsymbol{\alpha}}_t + (\partial_{\mathbf{x}} \mathbf{g}_t) \tilde{\mathbf{x}}_t + (\partial_{\mathbf{a}} \mathbf{g}_t) \tilde{\mathbf{a}}_t \quad \text{with} \quad \tilde{\mathbf{a}}(0) = 0 . $$

It remains to express $\tilde{\mathbf{a}}_t$. We therefore introduce the flow $\frac{dF_{st}}{dt} = (\partial_{\mathbf{a}} \mathbf{g}_t) F_{st}$ with $F_{tt} = I$ and get

$$ \begin{aligned} \tilde{\mathbf{a}}_t &= \int_0^t F_{ut}((\partial_{\boldsymbol{\alpha}} \mathbf{g}_u) \tilde{\boldsymbol{\alpha}}_u + (\partial_{\mathbf{x}} \mathbf{g}_u) \tilde{\mathbf{x}}_u) \, du \\ &= \int_0^t F_{ut}(\partial_{\boldsymbol{\alpha}} \mathbf{g}_u) \tilde{\boldsymbol{\alpha}}_u \, du + \int_0^t \int_0^u F_{ut}(\partial_{\mathbf{x}} \mathbf{g}_u) \varphi_s \tilde{\boldsymbol{\alpha}}_s \, ds \, du \\ &= \int_0^t \left( F_{ut}(\partial_{\boldsymbol{\alpha}} \mathbf{g}_u) + \int_u^t F_{st}(\partial_{\mathbf{x}} \mathbf{g}_s) \, ds \, \varphi_u \right) \tilde{\boldsymbol{\alpha}}_u \, du . \end{aligned} $$

In particular there holds

$$ \tilde{\mathbf{a}}_T = \int_0^T \left( F_{tT}(\partial_{\boldsymbol{\alpha}} \mathbf{g}_t) + \int_t^T F_{sT}(\partial_{\mathbf{x}} \mathbf{g}_s) \, ds \, \varphi_t \right) \tilde{\boldsymbol{\alpha}}_t \, dt . \qquad (12) $$

Combining (9), (10) and (12) we have

$$ \begin{aligned} \tilde{E}^m &= \int_0^T (\partial_{\mathbf{x}_T} E) \varphi_t \tilde{\boldsymbol{\alpha}}_t + (\partial_{\mathbf{a}_T} E) \left( F_{tT}(\partial_{\boldsymbol{\alpha}} \mathbf{g}_t) + \int_t^T F_{sT}(\partial_{\mathbf{x}} \mathbf{g}_s) \, ds \, \varphi_t \right) \tilde{\boldsymbol{\alpha}}_t \, dt \\ &= \int_0^T \left( \underbrace{\left( (\partial_{\mathbf{x}_T} E) + \int_t^T (\partial_{\mathbf{a}_T} E) F_{sT}(\partial_{\mathbf{x}} \mathbf{g}_s) \, ds \right)}_{=: \eta_t^*} \varphi_t + \underbrace{(\partial_{\mathbf{a}_T} E) F_{tT}(\partial_{\boldsymbol{\alpha}} \mathbf{g}_t)}_{=: \zeta_t^*} \right) \tilde{\boldsymbol{\alpha}}_t \, dt \\ &= \int_0^T (\eta_t^* \varphi_t + \zeta_t^*(\partial_{\boldsymbol{\alpha}} \mathbf{g}_t)) \tilde{\boldsymbol{\alpha}}_t \, dt . \end{aligned} \qquad (13) $$
Since $F_{st}F_{ts} = I$ and $\frac{dF_{st}^*}{ds} = -(\partial_a g_s)^* F_{st}^*$ we have in particular the integral form
$F_{ts}^* = I + \int_t^s (\partial_a g_u)^* F_{us}^* du$. This helps to simplify

$$
\begin{align}
\zeta_t = F_{tT}^*(\nabla_{a_T} E) &= \left(I + \int_t^T (\partial_a g_s)^* F_{sT}^* ds\right) (\nabla_{a_T} E) \nonumber \\
&= \nabla_{a_T} E + \int_t^T (\partial_a g_s)^* F_{sT}^*(\nabla_{a_T} E) ds = \nabla_{a_T} E + \int_t^T (\partial_a g_s)^* \zeta_s ds \tag{14}
\end{align}
$$

$$
\eta_t = \nabla_{x_T} E + \int_t^T (\partial_{x_T} g_s)^* F_{sT}^*(\nabla_{a_T} E) ds = \nabla_{x_T} E + \int_t^T (\partial_{x_T} g_s)^* \zeta_s ds . \quad (15)
$$

Collecting (8), (13), (14) and (15) yields the assertion. □

Acknowledgement. We thank Stanley Durrleman from the University of Utah for the fruitful discussion and helpful suggestions at the initial phase of this paper. We further acknowledge support by DFG-MATHEON Project F2. Moreover we thank the referees for their comments which helped to improve the quality of this paper.

References

1. Beg, M.F., Miller, M.I., Trouvé, A., Younes, L.: Computing large deformation metric mappings via geodesic flows of diffeomorphisms. IJCV 61(2), 139–157 (2005)

2. Cao, Y., Miller, M.I., Winslow, R.L., Younes, L.: Large deformation diffeomorphic metric mapping of vector fields. IEEE Trans. Med. Imag. 24(9), 1216–1230 (2005)

3. Cotter, C.: The variational particle-mesh method for matching curves. Journal of Physics A: Mathematical and Theoretical 41(34), 344003 (2008)

4. Dupuis, P., Grenander, U., Miller, M.I.: Variational problems on flows of diffeomorphisms for image matching. Q. Appl. Math. 56(3), 587–600 (1998)

5. Durrleman, S.: Statistical models of currents for measuring the variability of anatomical curves, surfaces and their evolution. Phd, Univ. Nice (2010)

6. Durrleman, S., Pennec, X., Trouvé, A., Ayache, N.: Statistical models of sets of curves and surfaces based on currents. MedIA 13(5), 793–808 (2009)

7. Federer, H.: Geometric measure theory. Repr. of the 1969 ed. Classics in Mathematics. Berlin: Springer-Verlag (1996)

8. Glaunès, J., Qiu, A., Miller, M.I., Younes, L.: Large deformation diffeomorphic metric curve mapping. IJCV 80(3), 317–336 (2008)

9. Glaunès, J., Trouvé, A., Younes, L.: Diffeomorphic matching of distributions: a new approach for unlabelled point-sets and sub-manifolds matching. In: CVPR. Los Alamitos: IEEE Comput. Soc. vol. 2, pp. 712–718 (2004)

10. Glaunès, J.: Transport par difféomorphismes de points, de mesures et de courants pour la comparaison de formes et l'anatomie numérique. Ph.D. thesis, Université Paris 13 (2005)

11. Kirk, B.S., Peterson, J.W., Stogner, R.H., Carey, G.F.: libMesh: A C++ Library for Parallel Adaptive Mesh Refinement/Coarsening Simulations. Engineering with Computers 22(3-4), 237–254 (2006)

12. Marsland, S., Twining, C.: Constructing diffeomorphic representations for the groupwise analysis of non-rigid registrations of medical images. IEEE Trans. Med. Imag. 23(8), 1006–1020 (2004)
13. Mattheij, R.M., Molenaar, J.: Ordinary differential equations in theory and practice. Reprint of the 1996 original. Classics in Applied Mathematics. SIAM (2002)

14. Morgan, F.: Geometric measure theory. A beginner's guide. 4th ed. Elsevier (2009)

15. Morita, S.: Geometry of differential forms. American Mathematical Society (2001)

16. Risser, L., Vialard, F.X., Wolz, R., Holm, D.D., Rueckert, D.: Simultaneous fine and coarse diffeomorphic registration: Application to atrophy measurement in alzheimer's disease. In: MICCAI'10. vol. 13, pp. 610–617 (2010)

17. Sommer, S.H., Nielsen, M., Lauze, F.B., Pennec, X.: A multi-scale kernel bundle for LDDMM: Towards sparse deformation description across space and scales. In: IPMI. LNCS, Springer (2011)

18. Trouvé, A.: An infinite dimensional group approach for physics based models in pattern recognition. Tech. rep., Johns Hopkins University (1995)

19. Vaillant, M., Glaunès, J.: Surface matching via currents. In: Christensen, G.E., Sonka, M. (eds.) IPMI, LNCS, vol. 3565, pp. 381–392. Springer (2005)

20. Vaillant, M., Miller, M.I., Younes, L., Trouvé, A.: Statistics on diffeomorphisms via tangent space representations. NeuroImage 23(Supplement 1), 161–169 (2004)

21. Younes, L.: Shapes and diffeomorphisms. App. Math. Sc. 171. Springer (2010)
# Frame to Frame Diffeomorphic Motion Analysis
from Echocardiographic Sequences

Zhijun Zhang¹, David J. Sahn¹,², and Xubo Song¹

¹ Department of Biomedical Engineering

² Department of Pediatric Cardiology

Oregon Health and Science University

20000 NW Walker Road, Beaverton, OR 97006, USA

{zhangzhi,songx,sahn}@ohsu.edu

**Abstract.** Quantitative motion analysis from echocardiography is an important yet challenging problem. We develop a motion estimation algorithm for echocardiographic image sequences based on diffeomorphic image registration in which the velocity field is spatiotemporally smooth. The novelty of this work is that instead of optimizing a functional of velocity field which consists of similarity metrics between a reference image to each of the following images (*first-to-follow*), we optimize a functional which is a sum of similarity metrics of each two consecutive images (*frame-to-frame*). This method can reduce the bias effect of using a single image as reference. It also improves registration accuracy since consecutive frames usually have higher dependency than frames far away. We validate our method by using both simulated images with known ground truth and *in vivo* pig heart images with sonomicrometry. Tests indicate that our frame-to-frame motion estimation method is more accurate than first-to-follow method.

## 1 Introduction

Quantitative analysis of cardiac deformation and motion is important for studying heart function. Many illnesses related to ischemia or infarct can be recognized from the motion and deformation abnormalities [1]. Techniques to discriminate the abnormal motion and accurately locate regions with motion abnormality are critical to identify the disease and to evaluate the treatment. Echocardiography (echo) is the most widely used image modality because it is non-ionizing, real-time, cost-effective and convenient. With the development of the new transducer array technology, 3D echo can now provide real-time images of the whole heart [2]. However, due to the low signal-noise-ratio, general methods for motion estimation do not work well on echo images. In addition, the 4D (3D+t) data is acquired with a compromise that both the spatial and temporal resolutions are reduced comparing to 2D+t sequences. As a result, 3D motion analysis from echo sequences remains a challenging problem.

Cardiac motion analysis algorithms can be classified into three categories: model-based, feature-based and voxel-based methods. Lots of cardiac models
have been proposed for motion analysis and segmentations [3]. Surface models
such as super-quadrics are used for motion analysis by fitting the model with
a sequence of images [4, 5], however, these models only estimate the deforma-
tion on the model surfaces. Comaniciu et al. [6] proposed a Kalman filter based
shape tracking method by using information fusion framework with a probabilis-
tic subspace model constraint. In this work, a shape model needs to be learned
and the motion estimation is limited along the contour points. Volumetric mod-
els such as dynamic finite elements have been used to estimate the deformation
inside the myocardium [7, 8]. Wang et al. [9] tracked myocardial surface points
by maximizing the likelihood of a combined surface and a two-steps motion
prediction model. Both the initial myocardial surface detector and the motion
prediction model need to be learned in advance. Generally speaking, deformable
model based methods needs prior knowledge related to the models and their
generation needs some sort of human interaction. Feature-based methods use
landmarks such as the tagged lines to fit the deformable model such as 4D B-
spline [10]. However, in echocardiography due to lack of stable landmarks in
the myocardium and artificial features such as the tagged line are not available,
feature-based method is difficult to estimate the deformation by using trustable
correspondence. Voxel-based methods require no manual intervention, they esti-
mate spatially dense transformations from all image voxels directly. This method
can be implemented as an automatic method and we focus our work on this ap-
proach. Voxel based image registration methods such as optical-flow [11] and
B-spline based methods [12, 13] have been proposed for cardiac motion analy-
sis from echo images. However, the motion analysis problem is simplified into
a series of independent pairwise image registrations and the temporal motion
smoothness is not considered. To enforce the temporal consistency of particle
motion, many temporal models have been used. Carbayo et al. [14] proposed a
spatiotemporal deformation model for cardiac motion tracking. A 2D+t B-spline
transformation with spatiotemporal smoothness is used with the first frame as
the reference. A 3D+t extension is proposed by Metz et al. [15] and the average
image is used as the reference. Particle trajectory constraint such as polynomial
modeling has been used to regularize the spatiotemporal motion smoothness
[16]. Diffeomorphic image registration is a method which the transformation
is implicitly spatiotemporally smooth. The transformation between two images
is defined as the end point of a velocity field flow which will be obtained by
optimization of an energy functional of it. It has a very useful characteristic
in computational anatomy that the transformation is one-to-one mapping and
topology preserving [17]. This large deformation topology preserving property
is preferred in cardiac motion analysis because of the fact that the deformation
from the reference frame (usually the end of diastolic) to the mostly contracted
frame (the end of systolic) is so large that deformation models without topology
preserving constraint may cause the transformation to fold over or tear apart,
which is not physically plausible. Beg et al. [18] proposed a large deformation
diffeomorphic metric mapping algorithm (LDDMM) in which the smooth veloc-
ity field is estimated by optimizing a sum of squared difference (SSD) energy.
Khan et al. [19] extended the LDDMM method to analyze the anatomical shape evolution in an image sequence. De Craene et al. [20] proposed a method in which the velocity field is defined as temporal piecewise continuous 3D B-spline functions and the B-spline control parameters are estimated by optimization of a parameterized energy function. In a following work [21], the velocity field is defined as a 3D+t spatiotemporal B-spline model to reduce B-spline control points in temporal direction. In both methods, the optimal velocity field minimizes the summed dissimilarity metrics between the first frame and each of the unwarped subsequent frames, which we call them first-to-follow methods. It has been presented theoretically and experimentally that speckle pattern will change under large deformation [22] and that registration of frames further away from the reference is less accurate due to speckle de-correlation [12]. We propose a diffeomorphic registration method with a spatiotemporally smooth velocity field which minimizes the summed SSDs of the unwarped consecutive frames (frame-to-frame method). Our registration method is tested with simulated and in-vivo pig datasets, the results show that the accuracy is improved over first-to-follow method.

# 2 Method

## 2.1 Diffeomorphic Image Sequence Registration

A diffeomorphism flow is a dynamic system with each of the diffeomorphism to be a state in a differentiable manifold [23]. We define a flow $\phi(\mathbf{x}, t), t \in [0, T], \mathbf{x} \in \Omega \subset R^d(d = 2, 3)$ with its smooth velocity field $v(\mathbf{x}, t)$ by using the differential equation of $\frac{d\phi}{dt} = v(\phi(\mathbf{x}, t), t)$. It has been proven in [24] that if $v(\mathbf{x}, t)$ is smooth enough with a differential operator $L$ in a Sobolev space $V$, then the transformation $\phi(\mathbf{x}, t)$ will be a group of diffeomorphisms with $t$ varying from 0 to $T$. The diffeomorphic image registration is stated as a variational problem, that given two images $I_0$ and $I_1$, to find an optimal velocity field $\hat{v}$ which minimizes an energy functional consisting of a sum of squared difference (SSD) and a geodesic distance metric between $\phi(\mathbf{x}, 0)$ and $\phi(\mathbf{x}, T)$ [18]:

$$ \hat{v} = \arg \inf_{v \in V} \lambda \int_{0}^{T} \|v(\mathbf{x}, t)\|_{V}^{2} dt + \int (I_{0}(\mathbf{x}) - I_{1}(\phi(\mathbf{x}, T)))^{2} d\mathbf{x}, \quad (1) $$

with $\lambda$ being the weight to balance these two energies. If we have a sequence of $N_f$ images to be registered, the similarity metric consists of the SSD of the difference between a reference frame $I_0$ and each of the deformed subsequent frames $I_k(\phi(\mathbf{x}, t_k)), k = 1, 2, ..., N_f$. Then the mathematical form can be in a similar form as [19]:

$$ \hat{v} = \arg \inf_{v \in V} \lambda \int_{0}^{T} \|v(\mathbf{x}, t)\|_{V}^{2} dt + \sum_{k=1}^{k=N_f} \int (I_{0}(\mathbf{x}) - I_{k}(\phi(\mathbf{x}, t_k)))^{2} d\mathbf{x}. \quad (2) $$

This scheme has two disadvantages: first, the speckle de-correlation between far away frames are high which may cause correspondence ambiguity between
the two images; second, it take longer time to converge since the difference between reference frame to far away frames is bigger than that of the consecutive frames. Instead of optimizing each deformed subsequent frames to be similar to the reference frame, we propose a variational energy which minimizes the difference between every two deformed consecutive frames $I_{k-1}(\phi(\mathbf{x}, t_{k-1}))$ and $I_k(\phi(\mathbf{x}, t_k))$:

$$ \hat{v} = \arg \inf_{v \in V} \lambda \int_0^T ||v(\mathbf{x}, t)||_V^2 dt + \sum_{k=1}^{k=N_f} \int (I_{k-1}(\phi(\mathbf{x}, t_{k-1})) - I_k(\phi(\mathbf{x}, t_k)))^2 d\mathbf{x}, \quad (3) $$

we denote the two energy terms in Eqn.(3) as $E_{reg}$ and $E_{ssd}$ respectively. This method will find a flow of diffeomorphisms which have the shortest geodesic path in the manifold and simultaneously minimizes summed errors of each two neighboring images. This will improve the accuracy of the transformation since the neighboring frames generally have higher intensity correlation than those which are not neighbors. It can also reduce the chance that the transformation is biased due to noise in the reference images. Our frame-to-frame method is different from the method which simply estimates the diffeomorphisms between consecutive frames and then composites together. In our method, minimizing the difference between two consecutive frames jointly optimizes all the velocity field before and between the time of these two frames.

The direct solution for this variational framework is expensive. Alternatively, a parameterized representation of the velocity field is used [25], where the velocity field is represented as a series of B-spline functions and the displacement field can be expressed as the forward Euler integral of velocity field. We use a spatiotemporally smooth B-spline function to represent the velocity field. It is defined as $v(\mathbf{x}, t_k) = \sum_{i;k} \beta(\mathbf{x} - \mathbf{x}_i)$, with $c_{i;k}$ being the B-spline control vectors at $t_k$ located on a uniform grid of $\mathbf{x}_i$, $\beta(\mathbf{x} - \mathbf{x}_i)$ is the 3D B-spline kernel function which is the tensor product of the 1-D B-spline functions. Define $\phi_k = \phi(\mathbf{x}, t_k)$ the transformation at time step $t_k$, we assume the velocity is piecewise constant within a time step, then we have $\phi_k = \phi_{k-1} + v(\phi_{k-1}, t_{k-1})\Delta t$, with $\phi_0(\mathbf{x}) = \mathbf{x}$, $k = 1, 2, ..., N_t$, with $N_t$ being the total number of time steps of the velocity field. Without loss of generality, we can have $\Delta t = 1$. In our test, we use one time step in velocity field between two neighboring frames since the deformation between them is usually small, that is $N_t = N_f$. However, our method easily generalizes to multiple time steps between frames if the deformation between two consecutive frames is large. The energy functional will be a parameterized function of $c_{i;k}$ and it can be optimized by using a steepest descent method.

## 2.2 Regularization

In order to assure the $\phi(\mathbf{x}, t)$ to be diffeomorphic, we need to define $v(\mathbf{x}, t)$ to be spatiotemporally smooth under a differential operator $L$. The linear operator we choose is: $L = \nabla^2 v + w_t \frac{dv}{dt}$, with $\nabla^2(\cdot)$ being a Laplacian operator and $w_t$ a
constant weight. In the discrete time form of velocity field, the time integral of the norm in V space of Eqn.(3) will be: $E_{reg} = \sum_{k=1}^{N_t} \sum_{\mathbf{x}} (\nabla^2 \mathbf{v}_k)^2 + w_t \sum_{k=2}^{N_t} \sum_{\mathbf{x}} |\mathbf{v}_k(\mathbf{x} + \mathbf{v}_{k-1}\Delta t) - \mathbf{v}_{k-1}|^2$, with $\mathbf{v}_k = \mathbf{v}(\mathbf{x}, k)$. The first term makes the velocity field spatially smooth which is denoted as $E_{sr}$. The second term keeps the particle velocity smooth and it is denoted as $E_{tr}$. The overall effect is to keep the velocity field spatiotemporally smooth.

## 2.3 Optimization

We use a steepest descent method to optimize the parameterized function. The derivative of the total registration energy with respect to the transformation parameters will be calculated analytically. The derivative of the similarity metric with respect to the B-spline parameters $\mathbf{c}_{i;k'}$ is:

$$ \frac{\partial E_{ssd}}{\partial \mathbf{c}_{i;k'}} = \sum_{k=1}^{N_f} (I_k(\phi_k) - I_{k+1}(\phi_{k+1})) (\nabla I_k(\phi_k) \frac{\partial \phi_k}{\partial \mathbf{c}_{i;k'}} - \nabla I_{k+1}(\phi_{k+1}) \frac{\partial \phi_{k+1}}{\partial \mathbf{c}_{i;k'}}), \quad (4) $$

with $\frac{\partial \phi_k}{\partial \mathbf{c}_{i;k'}}$ being the Jacobian matrix of transformation at time step k with respect to the ith B-spline coefficient in k'th frame. It can be calculated with chain rule and it is zero when $k' \ge k$. For detailed computation refer to [21].

For the derivative of the spatial and temporal regularization energies with respect to the mth component of $\mathbf{c}_{i;k}$, we have:

$$ \frac{\partial E_{sr}}{\partial c_{i,m;k}} = \sum_{\mathbf{x} \in \Omega'} \beta''_{m}(\mathbf{x} - \mathbf{x}_{i}), \quad (5) $$

with $\Omega'$ being the local support of the B-spline kernel function, and $\beta''_{m}(\cdot)$ being the second derivative of the B-spline function with respect to mth component. Considering that the displacement between two time step is small, we have:

$$ \frac{\partial E_{tr}}{\partial c_{i,m;k}} \approx w_t \sum_{\mathbf{x} \in \Omega'} (2 * v_{i,m;k} - v_{i,m;k-1} - v_{i,m;k+1}) \beta(\mathbf{x} - \mathbf{x}_i). \quad (6) $$

The registration energy can be optimized by starting from initial position and descending along the negative gradient direction at each iteration until there is no significant decrease.

## 2.4 Implementation

In our implementation, we use a series of B-spline transformations with grid spacing of 10 in each dimension to represent the velocity field. The values of $\lambda$ and $w_t$ are set to be 0.1 and 0.5. The algorithm is implemented with Matlab under a windows XP 64 bit system on a machine with 2.13GHz Xeon 8 cores CPU and 6GB memory. It takes about 1 hour to register a 3D (3 minutes for a 2D) sequence with 20 frames.
## 3 Experiment and Data

We use both simulated and real data to validate our algorithm. In the simulated data experiment, a longitudinal view of a diastolic left ventricle (LV) image with size of 274 × 192 is used as the reference image. This frame is then deformed with a series of continuous displacement field functions. The deformations are symmetrical along the long axis of the LV to simulate the myocardial contraction effect along radial and longitudinal directions. The displacement functions are in form of: $f_x(i) = a_x \sin(\frac{\pi(x-x_c)}{2r_d}) \sin(\frac{i\pi}{N_f})$ and $f_y(i) = a_y \sin(\frac{\pi(y-y_{apex})}{2(y_{base}-y_{apex})}) (\sin(\frac{i\pi}{N_f} + \frac{\pi}{16}) - \sin(\frac{\pi}{16}))$, with $x_c, r_d$ the axis center coordinate and the average axial radius of LV, $y_{apex}$ and $y_{base}$ the height of base and apex planes, $N_f$ and $i$ the number of frames and the frame index, and $a_x, a_y$ are the magnitudes of displacement fields. An image sequence with $N_f + 1$ frames is generated when $i$ varies from 0 to $N_f$ to simulate the cardiac motion in one cycle.

We carry out two experiments for the simulated data. In the first experiment, three sequences with 20 frames each are simulated with multiplicative speckle noise of variance 0.06, 0.08 and 0.10 added. The reference frame and the 10th frame with speckle noise variance 0.10 are shown together with the ground truth displacement field in Fig.1. In the second simulated experiment, we first generate 20 frames without deformation by adding independent speckle noise of variance 0.10 to the reference frame. Then each of the frame $I_i$ will have a percentage $p$ pixels replaced with the intensity at the same position in frame $I_{i-1}$. By updating noisy image one by one we assure that the two consecutive frames to have noise overlap ratio of $p$. Each frames will then be deformed by using the ground truth displacement fields. We simulate two sequences with overlap ratio of 0.2 and 0.4 respectively.

<figure><img src="image_11.png" /><figcaption>**Fig. 1.** The reference frame and the 10th frame in speckle variance 0.1 test and the displacement field (only displacement field inside a bell-shaped mask is displayed).</figcaption></figure>

A real world dataset is acquired from an open-chest pig by using a Philips IE33 system. For validation, we installed six sonomicrometers in the heart wall. The distances between each pair of the sonomicrometers are recorded with the image sequences and are used as ground truth to compare with the tracked distance in the echo images. The images are resampled into volume sequences of 160 × 100 × 128 with voxel size 1mm×1mm×1mm. The crystal coordinates in the reference frame are manually denoted.
# 4 Result

In the first simulated experiment, we compare our frame-to-frame method with the first-to-follow method by tracking the trajectories of the points in the myocardial wall during the motion process. The estimated trajectories of four example points from a small region of the myocardium for speckle noise variance 0.06 test are shown in Fig.2, where the ground truth trajectories are overlaid for comparison. We can see generally coordinates of the points in each time step in our method are closer to the ground truth position than the first-to-follow method.

<figure><img src="image_12.png" /><figcaption>**Fig. 2.** Points trajectories in the first-to-follow method (left) and the frame-to-frame method (right). The ground truth trajectories (blue) are overlaid with the estimated curves (multiple color) for comparison. The arrow shows the velocity at each time step.</figcaption></figure>

In Fig.3 we illustrate the motion estimation errors in both *x* and *y* coordinates in noise level 0.08 dataset. We can see that the motion estimation errors of *x* and *y* coordinates in most of frames in our method are closer to zero than those of the first-to-follow method. The figure also shows that our method has a smaller error variance in both coordinates. The mean of magnitude of errors in the three noise variance levels in frame-to-frame method are 0.23, 0.26 and 0.32, while in first-to-follow method they are 0.38, 0.45 and 0.56 respectively.

In the second test, the motion estimation errors in *y* coordinates for two methods are shown in Fig.4. The results are similar for *x* coordinate errors. We can see the results of frame-to-frame method are better than those in first-to-follow method in means and standard deviations of errors. We can see for the frame-to-frame method, when the intensity correlation between two images is increased, the mean and variance of the registration error are decreased. In the results of first-to-follow methods, the error mean does not change obviously when the correlation between consecutive frames are increased.

In the *in-vivo* open-chest pig test, we compare the performance of the two algorithms by computing the correlations over time between the algorithm-derived pair-wise distances with sonomicrometry, shown in Table.1. Sonomicrometry pro-
<figure><img src="image_13.png" /><figcaption>**Fig. 3.** The motion estimation errors in *x* (left) and *y* (right) coordinates of frame-to-frame method and first-to-follow method. The error bars show the standard deviation of the errors in each frames. The horizontal lines represent the curves of zero mean transformation errors to help comparison.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Fig. 4.** The motion estimation errors in *y* coordinates of frame-to-frame method and first-to-follow method. The left and right figures show the results of 20 and 40 percent noise overlap tests. The horizontal lines show the zeros mean transformation errors.</figcaption></figure>
vide the ground truth distances between each two of the crystals. We can clearly
see the improvement of our proposed method.

<table><caption>**Table 1.** The correlations between the estimated pair-wise distances and those from the sonomicrometry, with frame-to-frame method (numbers to the left) and the first-to-follow method (numbers to the right). Numbers 1-6 index the six sonomicrometry markers.</caption><thead><tr><th></th><th>1</th><th>2</th><th>3</th><th>4</th><th>5</th><th>6</th></tr></thead><tbody><tr><td>1</td><td>1.0/1.0</td><td>0.936/0.907</td><td>0.901/0.885</td><td>0.913/0.902</td><td>0.948/0.923</td><td>0.859/0.831</td></tr><tr><td>2</td><td>0.936/0.907</td><td>1.0/1.0</td><td>0.881/0.856</td><td>0.927/0.904</td><td>0.887/0.838</td><td>0.951/0.916</td></tr><tr><td>3</td><td>0.901/0.885</td><td>0.881/0.856</td><td>1.0/1.0</td><td>0.825/0.786</td><td>0.902/0.905</td><td>0.819/0.788</td></tr><tr><td>4</td><td>0.913/0.902</td><td>0.927/0.904</td><td>0.825/0.786</td><td>1.0/1.0</td><td>0.937/0.919</td><td>0.934/0.902</td></tr><tr><td>5</td><td>0.948/0.923</td><td>0.887/0.838</td><td>0.902/0.905</td><td>0.937/0.919</td><td>1.0/1.0</td><td>0.918/0.873</td></tr><tr><td>6</td><td>0.859/0.831</td><td>0.951/0.916</td><td>0.819/0.788</td><td>0.934/0.902</td><td>0.918/0.873</td><td>1.0/1.0</td></tr></tbody></table>

# 5 Conclusion

We propose a large deformation diffeomorphic registration method by minimizing the difference between every *consecutive images*. Simulation test shows that our frame-to-frame method has a higher estimation accuracy of motion than the first-to-follow method. Validation with sonomicrometry also shows that our motion estimation result has higher consistency with real data.

# Acknowledgement

This paper is supported by a NIH/NHLBI grant 1R01HL102407-01 awarded to Xubo Song and David Sahn.

# References

1. Fukuda, K., Oki, T., Tabata, T., Iuchi, A., Ito, S.: Regional left ventricular wall motion abnormalities in myocardial infarction and mitral annular descent velocities studied with pulsed tissue Doppler imaging. JASE, 11(9), 841–848 (1998)

2. Hung, J., Lang, R., Flachskampf, F., Shernan, S.K., McCulloch, M.L., Adams, D.B., Thomas, J., Vannan, M., Ryan, T.: 3D Echocardiography: A Review of the Current Status and Future Directions. JASE, 20(3), 213–233 (2007)

3. A.F. Frangi, W.J. Niessen, and Max A. Viergever. "Three-Dimensional Modeling for Functional Analysis of Cardiac Images: A Review." IEEE Trans. Med. Imaging, vol.20, no.1, pp.1–25, 2001.

4. Eric Bardinet, Laurent D. Cohen and Nicholas Ayache. "Tracking and motion analysis of the left ventricle with deformable superquadrics. Medical Image analysis, vol.1, no.2, pp:129–149, 1996.

5. Jinah Park, Dimitri Metaxas, Alistair A. Young and Leon Axel. "Deformable Models with Parameter Functions for Cardiac Motion Analysis from Tagged MRI Data." IEEE Transactions on Medical Imaging, vol.15, no.3, pp.278–289, June, 1996.

6. Comaniciu, D, Zhou, X., Krishnan, S.: Robust tracking of myocardial border: An information fusion approach. IEEE Trans. Med. Imag. 23(7), 849–860 (2004)
7. Papademetris, X., Sinusas, A.J., Dione, D.P., Duncan, J.S.: Estimation of 3D left ventricular deformation from echocardiography. Med. Imag. Anal. 5(1), 17–28 (2001)

8. Joel Schaerera, Christopher Casta, Jerome Pousina and Patrick Claryssea. "A dynamic elastic model for segmentation and tracking of the heart in MR image sequences." Medical Image Analysis, vol.14, pp.738–749, 2010.

9. Wang, Y., Georgescu, B., Houle, H., Comaniciu, D.: Volumetric Myocardial Mechanics from 3D+t Ultrasound Data with Multi-model Tracking. In STACOM 2010. LNCS vol. 6364, 184–193. Springer, Heidelberg (2010)

10. J. Huang, D. Abendschein, V. Davila, A.A. Amini, "Four-dimensional LV tissue tracking from tagged MRI with a 4D B-spline model." IEEE Trans Medical Imaging, vol.18, no.10, 1999.

11. Suhling, M., Arigovindan, M., Jansen, C., Hunziker, P., Unser, M.: Myocardial Motion Analysis from B-Mode Echocardiograms. IEEE Trans. Image Process., 14(4) 525–536 (2005)

12. Elen, A., Choi, H.F., Loeckx, D., Gaom, H., Claus, P., Suetens, P., Maes, F., D'hooge, J.: Three-dimensional cardiac strain estimation using spatio-temporal elastic registration of ultrasound images: a feasibility study. IEEE Trans. Med. Imag., 27(11), 1580–1591, (2008)

13. Myronenko, A., Song, X.B., David, J.S: LV Motion Tracking from 3D Echocardiography Using Textural and Structural Information. MICCAI 2007. 428–435 (2007)

14. Ledesma-Carbayo, M.J., Mah-Casado, P., Santos, A., Prez-David, E., GarMA, Desco, M.: Spatio-Temporal Nonrigid Registration for Ultrasound Cardiac Motion Estimation. IEEE Trans. Med. Imag., 24(9), 1113–1126 (2005)

15. Metz, C.T., Klein, S., Schaap, M., Walsum, T., Niessen, W.J.: Nonrigid registration of dynamic medical imaging data using nD+t B-splines and a groupwise optimization approach. Med. Imag. Anal. 15(2), 238–249 (2011)

16. Castillo, E., Castillo, R., Martinez, J., Shenoy, M., Guerrero, T.: Four-dimensional deformable image registration using trajectory modeling. Physics in Medicine and Biology, 55, 305–327 (2010)

17. Miller, M.I; Computational anatomy: shape, growth, and atrophy comparison via diffeomorphisms. Neuroimage, 23(1), 19–33 (2004)

18. Beg, M.F.; Miller, M.I.; Trouve, A.; Younes, L.; Computing Large Deformation Metric Mappings via Geodesic Flows of Diffeomorphisms. INT J COMPUT VISION. 61(2), 139–157 (2005)

19. Khan, A.R.; Beg, M.F.; Representation of time-varying shapes in the large deformation diffeomorphic framework. ISBI 2008, 1521–1524 (2008)

20. Craene, M.; Camara,O.; Bijnens, B.H.; Frangi, A.F.: Large Diffeomorphic FFD Registration for Motion and Strain Quantification from 3D-US Sequences, Functional Imaging and Modeling of the Heart, FIMH’09, LNCS vol.5528, 437–446 (2009)

21. Craene,M.; Piella,G.; Duchateau,N.; Silva,E.; Doltra,A.; Gao,H.; D'hooge,J.; Camara,O.; Brugada,J.; Marta Sitges: Temporal Diffeomorphic Free-Form Deformation for Strain Quantification in 3D-US Images. MICCAI 2010, 1–8 (2010)

22. Meunier, J.: Tissue motion assessment from 3D echographic speckle tracking. Phys. Med. Biol. 43, 1241–1254 (1998)

23. Arrowsmith, D.K.; Place, C. M.: An introduction to dynamical systems. Cambridge Press, UK (1990)

24. Dupuis, P.; Grenander,U.: Variational problems on flows of diffeomorphisms for image matching. Q APPL MATH. 56(3) 587–600 (1998)

25. Ashburner, J.; A fast diffeomorphic image registration algorithm. NeuroImage, 38(1), 95–113 (2007)
# An Improved Estimator of GRID Model for Representing Large Diffeomorphic Deformations

Qian Xie and Anuj Srivastava

Department of Statistics, Florida State University, United States

`{qxie,anuj}@stat.fsu.edu`

**Abstract.** The growth by random iterated diffeomorphisms (GRID) model seeks to decompose large deformations, caused by growth, anomaly, or anatomical differences, into smaller, biologically-meaningful components. These components are spatially local and parametric, and are characterized by radial deformation patterns around randomly-placed seeds. A sequential composition of these components, using the group structure of diffeomorphism group, models the cumulative deformation. The actual decomposition requires estimation of GRID parameters from observations of large growth, typically from 2D or 3D images. While past papers have estimated parameters under certain simplifying assumptions, including that different components are spatially separated and non-interacting, we address the problem of parameter estimation under the original GRID model that advocates sequential composition of arbitrarily interacting components. Using a gradient-based approach, we present an algorithm for estimation of GRID parameters by minimizing an energy function and demonstrate its superiority over the past additive methods.

**Keywords:** Large deformation, GRID model, parameter estimation

## 1 Introduction

The mathematical and statistical modeling of diffeomorphic deformations over time is an important problem with a variety of applications ranging from medical diagnostics to evolutionary biology. The use of medical images, especially the MRI images of human parts, in studying anatomical structures is a growing area of research by itself. Here one uses 2D and/or 3D images taken across time, species, or specimens to compare to extract salient differences in anatomical structures, and to analyze and model their variations both within and across biological classes. These differences may result from standard biological growth, abnormalities, inter-specimen variability, or other reasons. In terms of image-based analysis of anatomical structures, the study of shapes of anatomical parts has become a central idea. For instance, one can use longitudinal image data for tracking biological growth [11, 9, 19, 4, 3] in fetus brains and evaluating tumor growth. A major difficulty in solving such problems is the high dimensionality of image data. The diffeomorphic deformation when estimated from image
sequences can be very high dimensional and not amenable to standard tools
from multivariate statistics. Some current methods simplify this analysis by us-
ing simplistic measures like lengths, sizes, or areas as indicators of overall shape
changes. Some others use relatively simple geometrical models, such as spheres
or ellipsoids, to represent shapes in parametric forms and to study the evolution
of parameters during growth.

We start with the basic question of how to represent large deformations
in a mathematical framework. There is a large body of work on represent-
ing differences in imaged objects using deformations of background space [2,
12, 7, 20, 1, 13]. This approach utilizes diffeomorphisms of the underlying coor-
dinate systems to represent and measure shape and other differences. Let an
image be $I : [0, 1]^d \to \mathbb{R}$, where $d = 1, 2, 3$. A deformation is then a mapping
$\Phi : [0, 1]^d \to [0, 1]^d$, with the resulting deformed image is $I \circ \Phi : [0, 1]^d \to \mathbb{R}$.
Thus, a point on an anatomical landmark is always observed with the same color
intensity; it simply moves to a different location under the deformation. The goal
is to use $\Phi$ to model, understand and analyze large deformations. These deforma-
tions are typically very high dimensional and do not permit standard statistical
analysis directly. Therefore, some tools for reducing dimensionality become im-
portant. One can apply some standard dimension reduction algorithms, such
as PCA, but it is difficult to interpret the resulting representation in biological
terms. Durrleman et al. [5] proposed a parametric way of representing large dif-
feomorphisms by forcing the instantaneous velocity fields to take a parametric
form. One starts with a finite number of so-called control points and for each of
them specifies a vector that defines the deformation at that point. The vector
field over the whole domain is obtained using a Kernel-based interpolation. This
approach provides a data-driven sparse parametric method to estimate the large
diffeomorphic deformation.

Motivated by the need for biologically-interpretable decompositions of large
deformations, Grenander [8] introduced the Growth as Random Iterated Dif-
feomorphisms (GRID) model. It highlights the role of gene control in biolog-
cal growths and uses a combination of local, structured deformations to form
the large composite deformations. This model has been studied extensively, but
mostly from a perspective of synthesis and asymptotics. Some authors proposed
a “thermodynamic limit equation” that approximates the growth pattern in a
macroscopic way [16, 15]. Portman et al. [14] further developed the GRID model
by analyzing the growth patterns at microscopic levels. In addition to synthesis,
one is also interested in the inverse problem where we want to decompose large
biological growth into smaller biologically-interpretable units. Grenander et al.
[18, 8] studied this inverse problem albeit in a limited context. The estimation
of growth components was done in two steps: (1) estimate the full deformation
between a pair of images that represents biological growth, (2) estimate parame-
ters for growth components under the GRID model, with a major simplification
that different components are spatially local and do not interact with each other.
With this assumption, the cumulative growth becomes a simple superimposition
of different components and one can use standard projection procedures to esti-
mate component parameters. In this paper we seek a solution to the problem of parameter estimation under the original GRID model, without assuming spatial independence of components. This model is different from Durrleman et al. [5] in the sense that it is the individual diffeomorphisms that take the parametric forms, rather than the instantaneous velocity fields. This results in different local deformations around the focal points. For example, in our method the diffeomorphism around a seed is restricted to be radial while in Durrleman et al. [5] there is no such structure.

The estimation follows the two steps as Grenander et al. [8]. In the first step, the full deformation $\Phi$ is estimated using the shape matching technique of [10]. As for the second step, since there is a concatenation of deformation associated with different components, the time-ordering of the components becomes important. Due to the nonlinear effect of compositions, it is not possible to solve for GRID parameters using linear methods. This general estimation problem is posed as an optimization problem with a gradient-based minimization of the cost function. The difficulty of getting trapped in local solutions is handled using clever initializations of the gradient algorithm.

## 2 Grenander's GRID Model

We start by describing the general GRID model as introduced by Grenander [6]. In this model the overall large deformation is modeled as a composition of a sequence of local, elementary deformations. At time t, the elementary growth is a diffeomorphism $\phi_t : [0, 1]^d \to [0, 1]^d$ such that the point x moves to $\phi_t(x)$. The full deformation is then expressed in the form of the composition of iterated diffeomorphisms.

$$ \Phi_{t_1}^{t_n} \equiv \phi_{t_1} \circ \dots \circ \phi_{t_2} \circ \phi_{t_1} \quad (1) $$

The next step in the GRID model is to simplify diffemorphic components by expressing each $\phi_t$ in a parametric form. Here the elementary deformation $\phi_t$ is: (1) assumed to be centered at a point of activation called a growth seed $x_{\text{seed}}$, and (2) the growth around the seed is assumed to be radial. Therefore, it is easier to express this local deformation using polar coordinates centered at the seed $x_{\text{seed}}$: $(r, \tau) \mapsto (\rho(r, \tau), \tau)$ with $r = 0$ denoting the seed. Furthermore, the model assumes that the change in radial distance can be decomposed into two independent parts: $\rho(r, \tau) = r + R(r) A(\tau)$. Here, $A: \mathbb{S}^1 \to \mathbb{R}$ is called the angular deformation function (ADF) and $R: \mathbb{R}_+ \to \mathbb{R}_+$ is called the radial deformation function (RDF). These individual deformation functions are allowed to take the following forms:

1. The radial deformation function RDF can be one of the following two types:

$$ R(r) = \begin{cases} re^{-r^2/c^2}, & r \ge 0, c > 0 \\ (r/c)^{p-1} e^{-(r/c)}, & r \ge 0, p, c > 0. \end{cases} \quad (2) $$

In both cases the deformation is zero at the seed ($r = 0$), increases steadily with $r$, reaches a peak, and then decreases for a further increase in $r$. The “zone” of influence of a seed is determined by the parameter $c$.
2. Similarly, the angular deformation function ADF can also take many forms, including:

$$A(\tau) = \begin{cases} a, & \tau \in \mathbb{S}^1, a \in \mathbb{R} \\ ae^{\kappa \cos(\tau - \tau_0)} - \kappa, & \tau \in \mathbb{S}^1, a \in \mathbb{R}, \kappa \ge 0, \text{ and } \tau_0 \in \mathbb{S}^1 \\ \alpha \sin\left(\frac{2\pi\tau}{\sigma}\right), & \tau \in \mathbb{S}^1, \alpha \in \mathbb{R}, \text{ and } \sigma > 0. \end{cases} \quad (3)$$

The first case provides an isotropic deformation, the second provides a uni-modal deformation with a well-defined growth/decay direction, while the last one provides a sinusoidal variation.

In this paper we will use $R(r) = re^{-r^2/c^2}$ and will study two choices for ADFs: (1) $A(\tau) = a$, and (2) $A(\tau) = ae^{\kappa \cos(\tau - \tau_0)}$. It has been shown that the resulting $\phi_t$ is a diffeomorphism as long as $-1 < A(\tau) < 2.2408$ [8]. Each such $\phi_t$ is now characterized by the following set of parameters: $\theta = [\xi, a, c, \kappa, \tau_0] \in \mathbb{R}^6$, where $\xi \in \mathbb{R}^2$ is the seed location.

**Problem Statement:** Having chosen the model, the estimation problem can be described as follows. Let $\Phi$ represent the observed deformation associated with a growth experiment, observed over a time interval $[0, T]$. The goal now is to estimate $n$, the number of diffeomorphism components and the associated parameters $\theta_j \in \mathbb{R}^6$ for each $j = 1, 2, \dots, n$.

This problem has been studied by several papers in the past. However, a common simplifying assumption in the past papers is that different seeds are placed away from each other so that there is no or negligible interaction between the corresponding deformations. In this case, the total displacement field $\Psi(x) = \Phi(x) - x$ can be written as a superposition of the displacements resulting from individual seeds:

$$\Psi_{t_1}^{t_n} \equiv \psi_{t_n} + \dots + \psi_{t_2} + \psi_{t_1}, \quad (4)$$

where $\psi_t(x) = \phi_t(x) - x$. This is a very restrictive assumption and reduces the efficacy of the GRID model. The additive model has several problems, including the fact that the set of diffeomorphisms is not a group under the additive model. In the context of biological growth, it is difficult to interpret growth components under the assumption that there is no spatial interaction between them. Also, as shown in Fig. 1, the results of these different models, composite versus additive models are quite different for the same components. It is also illustrated that, for the composite model, the ordering of components is also important in determining the cumulative deformation. In this paper, we study the problem of parameter estimation under the composite model (Eqn. 1) and compare the results with those obtained under the additive model (Eqn. 4).

## 3 Our Approach

In this section we formulate the problem of parameter estimation as minimization of a certain objective function. The goal is to estimate diffeomorphic components
<figure><img src="image_15.png" /><figcaption>Fig. 1: Cumulative deformation for composite models with: (a) $\Phi = \phi_2 \circ \phi_1$, (b) $\Phi' = \phi_1 \circ \phi_2$, and additive model (c) $\Phi'' = x + \psi_1 + \psi_2$, for the same $\phi_1$ and $\phi_2$. In the remaining panels, we show $\Phi' - \Phi$ (4th panel), and $\Phi'' - \Phi$ (5th panel) as vector fields.</figcaption></figure>

$\phi_1, \phi_2, \dots, \phi_n$ such that their order composition is as closed to the given $\Phi$ as possible. In order words, we can use a distance between $\Phi$ and $(\phi_n \circ \dots \circ \phi_1)$ as the objective function. Although the choice of a geodesic distance in the space of diffeomorphisms (under a suitable Riemannian metric) seems like a natural choice, the use of the $\mathbb{L}^2$ distance simplifies the problem, by forming an energy

$$E = \int_{[0,1]^2} \|\Phi(x) - (\phi_n \circ \dots \circ \phi_1)(x)\|^2 dx . \qquad (5)$$

We justify the use of $\mathbb{L}^2$ distance, over the geodesic distance, with the argument that minimization under one distance often leads to a minimizer under the other.

It is rather difficult to solve for all the parameters (for all the seeds) simultaneously. Indeed, the expression for cumulative deformation with just two local deformations gets complicated. Therefore, we take a sequential approach and add one local deformation to the model at a time. Let $\phi^{(k)} = \phi_k \circ \phi_{k-1} \circ \dots \circ \phi_1$ be the cumulative deformation generated by first $k$ seeds. Define two energy functions associated with this partial inference problem:

$$E^{(k+1)} = \int \| \Phi(x) - (\phi_{k+1} \circ \phi^{(k)})(x) \|^2 dx$$

and

$$\tilde{E}^{(k+1)} = \int \| \Phi(x) - \phi^{(k)}(x) - \psi_{k+1}(x) \|^2 dx .$$

$E^{(k+1)}$ denotes the energy under the composite model for $k+1$ seeds while $\tilde{E}^{(k+1)}$ denotes a similar energy except that the contribution from the last seed is considered additive. (Since this last seed is additive, it is relatively easier to solve for its parameters by minimizing $\tilde{E}^{(k+1)}$.) Our iterative approach is to solve for the parameters of $\phi_{k+1}$ to minimize $E^{(k+1)}$, for $k = 1, 2, \dots, n$, and we will do so using a gradient approach. Similar to any gradient-based solution, the initialization of parameters becomes very important. For the purpose of improving initialization, we will solve for the parameters of $\phi_{k+1}$ under $\tilde{E}^{(k+1)}$ first and use these values as initial conditions in optimization of $E^{(k+1)}$.

We summarize the iterative procedure for estimating GRID parameters.
**Algorithm 1 Set k = 0.**

1. Given the current estimated parameters for the k-seed composite model, {$\theta_j$}_k = {$\langle \xi_j, a_j, c_j, \kappa_j, \tau_{0,j} \rangle \mid j = 1, \dots, k$}, compute the cumulative deformation $\phi^{(k)}$.

2. Find $\theta_{k+1} = (\xi_{k+1}, a_{k+1}, c_{k+1}, \kappa_{k+1}, \tau_{0,k+1})$ for the $(k+1)^{th}$ seed by minimizing $\tilde{E}^{(k+1)}$. Use these values as initial condition for parameters of the $(k+1)^{th}$ seed.

3. For each possible permutation group of the set {1, 2, ..., k + 1}, perform the following. Update each set of parameters {$\theta_j$}_{k+1} using the gradient method to minimize $E^{(k+1)}$. Finally, choose the permutation/parameters that result in the minimum $E^{(k+1)}$.

4. Test the significance of the $(k+1)^{th}$ seed. If it is found significant, set $k = k + 1$ and go to **step 1**; if not, **stop**.

Note that even though the complexity of the composed deformation $\phi^{(k)}$ and thus $E^{(k)}$ increases with the number of seeds (or local elementary deformations), we still have analytical expressions for the gradients using the chain rule. It is important to note that all the previous seeds are re-estimated/updated as new seeds are added to the deformation. Thus, although this process is iterative, it is not incremental.

We have studied two cases for estimating ADFs:

1. **Constant ADF:** In the first case, we simplify the discussion by first assuming that $A(\tau) = a$ for all $\tau \in \mathbb{S}^1$. In this model, the partial derivatives of $E^{(k)}$ with respect to different parameters of $\phi_j$ for $j = 1, \dots, k$ are given by:

$$ \frac{\partial E^{(k)}}{\partial \xi_j^r} = \sum_{i=1}^{n} \frac{\partial E^{(k)}}{\partial \phi^{(k)}} \cdot \frac{\partial \phi^{(k)}}{\partial \phi^{(k-1)}} \cdot \frac{\partial \phi^{(k-1)}}{\partial \phi^{(k-2)}} \cdots \frac{\partial \phi^{(j+1)}}{\partial \phi^{(j)}} \cdot \frac{\partial \phi^{(j)}}{\partial \xi_j^r}, \quad r = 1, 2, $$

$$ \frac{\partial E^{(k)}}{\partial a_j} = \sum_{i=1}^{n} \frac{\partial E^{(k)}}{\partial \phi^{(k)}} \cdot \frac{\partial \phi^{(k)}}{\partial \phi^{(k-1)}} \cdot \frac{\partial \phi^{(k-1)}}{\partial \phi^{(k-2)}} \cdots \frac{\partial \phi^{(j+1)}}{\partial \phi^{(j)}} \cdot \frac{\partial \phi^{(j)}}{\partial a_j}, $$

$$ \frac{\partial E^{(k)}}{\partial c_j} = \sum_{i=1}^{n} \frac{\partial E^{(k)}}{\partial \phi^{(k)}} \cdot \frac{\partial \phi^{(k)}}{\partial \phi^{(k-1)}} \cdot \frac{\partial \phi^{(k-1)}}{\partial \phi^{(k-2)}} \cdots \frac{\partial \phi^{(j+1)}}{\partial \phi^{(j)}} \cdot \frac{\partial \phi^{(j)}}{\partial c_j}. $$

The different terms needed in these expressions, including $\frac{\partial E^{(j)}}{\partial \phi^{(j)}}, \frac{\partial \phi^{(j)}}{\partial \phi^{(j-1)}}, \frac{\partial \phi^{(j)}}{\partial \xi_j}, \frac{\partial \phi^{(j)}}{\partial a_j},$ and $\frac{\partial \phi^{(j)}}{\partial c_j}$, for $j = 1, \dots, k$, are given in the appendix. Similar expressions can also be derived for the gradient of $\tilde{E}^{(k)}$ wrt the GRID parameters and those expressions are, as expected, simpler compared to the gradient given above.

2. **Non-isotropic ADF:** In the general case where the ADF is non-isotropic, we represent it using a scaled von-Mises density (the second term in Eqn. 3) and is parameterized by $(a, \kappa, \tau_0)$. For every deformation $\phi_j$, the two additional parameters $\kappa_j$ and $\tau_{0,j}$ are initialized in two steps. First, the ADF is estimated non-parametrically by integrating the deformation along each direction, as was done in [8]. Then, the parameters $\kappa$ and $\tau_0$ are estimated
from the nonparametric estimate using moment matching. The estimation is
similar to estimating the parameters of von-Mises density from a sample by
treating the non-parametric estimation as a weighted sample[17]. After these
parameters are initialized, they are estimated using the gradient method,
similar to the other parameters, with the gradient expression given in the
appendix.

The last remaining item in Algorithm 1 is the test of significance of an in-
cremental local deformation. In order to test the significance of the model with
one more seed, several general methods for model selection may apply. Possible
methods includes partial F test, AIC or BIC, and adjusted R². In this paper,
a model is selected based on the adjusted R². It is a modification of R², which
denotes the coefficient of determinant, that adjusts for the number of model
parameters. Given any two estimated models, the model with the larger value
of adjusted R² is preferred. In the experiments, we add one more seed if the
improvement of adjusted R² is larger than a small cutoff value.

Since this method is based on a gradient search, it is difficult to claim a global
solution. In principle, the solution obtained in the parameter space is a local one.
However, there are some advantages to using this approach. Firstly, since the
gradients of the energy function are available analytically, the gradient iterations
are computationally fast. Secondly, for relatively small number of seeds in the
model, the search over different orderings is efficient and gets us out of several
local solutions.

# 4 Experimental Illustrations

Here we demonstrate the use of our framework for decomposing large cumulative
deformations into smaller, parametric components using Algorithm 1. We will
use both the simulated and real data to illustrate the estimation process, and
will compare our results with those obtained using the additive model.

## 4.1 Synthetic Data

In order to validate estimation method for the composite model, we perform
two experiments on the following types of synthetic data: (1) a 2D deforma-
tion with constant ADF, and (2) 2D deformation with non-constant ADF. In
these experiments, a cumulative deformation is simulated using Eqn. 1 for an
arbitrary number of local deformations, each with arbitrary parameters, and a
white Gaussian noise is added at the end to form the observed deformation $\Phi$.

**Example 1:** Fig. 2 shows the estimation results for the isotropic model. It
shows the true underlying deformation (a), made up of $n = 2$ seeds, its noisy
observation $\Phi$ (b) and several different GRID estimates from the observed $\Phi$ (c-
e). Firstly, we estimate GRID components as described in Algorithm 1, with the
result shown as $\tilde{\Phi}_c$ in (e). Then, we reverse the order of two estimated seeds and
try to optimize their parameters, with the result shown as $\hat{\Phi}'_c$ in (c). The result of estimated deformation under the additive model is shown as $\hat{\Phi}_a$ in (d). To highlight the differences between different estimated deformations we also show their differences in the remaining panels. The Table 1 provides a quantification of estimation performance. It compares the energy $E$ and estimated parameters for the three models, with the true values. Since the energy for $\hat{\Phi}_c$ is same as that for true underlying parameters, it shows the superiority of that estimation process described in Algorithm 1.

<figure><img src="image_16.png" /><figcaption>Fig. 2: Experiment with isotropic ADF: (a) synthetic deformation; (b) noisy observation $\Phi$, and different estimated deformations (c) $\hat{\Phi}'_c$; (d) $\hat{\Phi}_a$, and (e) $\hat{\Phi}_c$. The differences (f) $\hat{\Phi}_c - \Phi$, (g) $\hat{\Phi}'_c - \Phi$, and (h) $\hat{\Phi}_a - \Phi$</figcaption></figure>

<table><caption>Table 1: Estimation results for Example 1.</caption><thead><tr><th>Model</th><th><i>E</i></th><th><i>&xi;</i><sub>1</sub><sup>(1)</sup></th><th><i>&xi;</i><sub>1</sub><sup>(2)</sup></th><th><i>a</i><sub>1</sub></th><th><i>c</i><sub>1</sub></th><th><i>&xi;</i><sub>2</sub><sup>(1)</sup></th><th><i>&xi;</i><sub>2</sub><sup>(2)</sup></th><th><i>a</i><sub>2</sub></th><th><i>c</i><sub>2</sub></th></tr></thead><tbody><tr><td>True</td><td><b>0.0334</b></td><td>0.6500</td><td>0.4500</td><td>1.0000</td><td>0.1500</td><td>0.5000</td><td>0.6000</td><td>0.8000</td><td>0.1000</td></tr><tr><td>Estimated <i>&Phi;&#x0302;</i><sub>c</sub></td><td><b>0.0334</b></td><td>0.6500</td><td>0.4502</td><td>1.0045</td><td>0.1498</td><td>0.4995</td><td>0.6006</td><td>0.7947</td><td>0.0998</td></tr><tr><td>Reverse order <i>&Phi;&#x0302;</i>'<sub>c</sub></td><td>0.0576</td><td>0.5070</td><td>0.5930</td><td>0.7072</td><td>0.0956</td><td>0.6478</td><td>0.4525</td><td>1.0108</td><td>0.1527</td></tr><tr><td>Additive <i>&Phi;&#x0302;</i><sub>a</sub></td><td>0.0548</td><td>0.6496</td><td>0.4507</td><td>1.0262</td><td>0.1501</td><td>0.5172</td><td>0.5827</td><td>0.7352</td><td>0.1004</td></tr></tbody></table>

**Example 2:** Similarly, Fig. 3(a) shows an experiment involving seeds with non-isotropic ADFs. In this case we show the synthetic deformation (a), its noisy version $\Phi$ (b), and the estimated deformation under the composite model $\hat{\Phi}_c$ (c). For further evaluation of this estimation, we show the true displacement
$\Psi(x) = \Phi(x) - x$ (d) and the estimated displacement $\hat{\Psi}$ (e).

<figure><img src="image_17.png" /><figcaption>Fig. 3: Experiment with non-isotropic ADF: (a)synthetic deformation; (b) noisy observation $\Phi$; (c) estimated deformation map $\hat{\Phi}$; (d) true displacement field $\Psi$; (e) estimated displacement field $\hat{\Psi}$</figcaption></figure>

<table><caption>Table 2: Estimation Results for Example 2.</caption><thead><tr><th>Model</th><th>E</th><th>ξ<sub>1</sub><sup>(1)</sup></th><th>ξ<sub>1</sub><sup>(2)</sup></th><th>a<sub>1</sub></th><th>c<sub>1</sub></th><th>ξ<sub>2</sub><sup>(1)</sup></th><th>ξ<sub>2</sub><sup>(2)</sup></th><th>a<sub>2</sub></th><th>c<sub>2</sub></th></tr></thead><tbody><tr><td>True</td><td>0.0330</td><td>0.6700</td><td>0.4000</td><td>1.0000</td><td>0.2000</td><td>0.3200</td><td>0.6000</td><td>1.0000</td><td>0.2000</td></tr><tr><td>Composite</td><td>0.0333</td><td>0.6677</td><td>0.4014</td><td>1.0004</td><td>0.1975</td><td>0.3233</td><td>0.5980</td><td>1.0248</td><td>0.1948</td></tr><tr><td></td><td></td><td>κ<sub>1</sub></td><td>τ<sub>0,1</sub></td><td></td><td></td><td>κ<sub>2</sub></td><td>τ<sub>0,2</sub></td><td></td><td></td></tr><tr><td>True</td><td></td><td>2.0000</td><td>2.3562</td><td></td><td></td><td>5.0000</td><td>5.4978</td><td></td><td></td></tr><tr><td>Composite</td><td></td><td>1.9200</td><td>2.3379</td><td></td><td></td><td>4.9050</td><td>5.4865</td><td></td><td></td></tr></tbody></table>

**Ensemble Results:** Beyond individual examples, we have exhaustively compared performances of the composite and the additive models using many realizations. Here we use 20 sets of data that are simulated from the 2D composite model with constant ADF and the parameters are then estimated using both the models. The results are summarized in Fig. 4. Panel (a) is a histogram of which model, composite or additive, is closer to the true deformation. Positive number indicates that the composite model outperforms the additive one, and vice-versa. We can see that for most data sets the composite model outperforms the additive model. Panel (b) presents a close up view of the parameter biases and variances. Each line represents the bias for one parameter. If the estimates equals the true parameters, the line will be all horizontal lines with $y = 0$. The top plot shows the biases from the additive model and the bottom one is for the composite model. The lines for composite model are all around zero with the ones for additive model having larger variation.
<figure><img src="image_18.png" /><figcaption>Fig. 4: (a) Relative errors under two models: $\frac{\|\hat{\Phi}_a - \Phi\| - \|\hat{\Phi}_c - \Phi\|}{\|\hat{\Phi}_a - \Phi\|}$; (b) Parameter estimation bias: top: $\hat{\theta}_a - \theta$, bottom: $\hat{\theta}_c - \theta$.</figcaption></figure>

## 4.2 Estimation of GRID Component for Image Data

This method is applied to analysis of differences in human brains observed using MRI scans. In this case we perform the following experiment. We take two image scans as *I* and *J*, which are two slices of MRI of the same *z* coordinates from different subjects, and use a shape-based technique [10] to estimate a cumulative deformation $\Phi$ from *I* to *J*. This $\Phi$ is in the set of diffeomorphisms from $[0, 1]^2$ to itself and minimizes a certain cost function involving $I \circ \Phi$ and *J*. Then, we use Algorithm 1, to estimate components of this deformation under the GRID model. We present two examples of this idea.

**Example 1:** Fig. 5 shows an example of images *I* and *J* and the deformation $\Phi$ that deforms *I* to match with *J*. In addition to plotting the maps $x \mapsto \Phi(x)$ as a surface mesh, one can also plot the displacement vector field $\Psi(x) = \Phi(x) - x$ for better visualization. Shown in the remaining panels are the estimated displacement vector $\hat{\Psi}$, showing the ordered sequence of displacement fields with the GRID components that were found using Algorithm 1, and a couple of ways of comparing the estimated deformations with the true deformations.

**Example 2:** Fig. 6 shows another example of estimating GRID parameters for deformation estimated from MRI images.

## 5 Conclusion

In this paper we have proposed a method to decompose and estimate the parameters in GRID based decomposition of anatomical deformations. The method preserves the iterative structure of the GRID model and gives an analytical form of the gradient for parameter estimation. Experimental results show that impacts from faraway seeds can be approximated by additive seeds model and composite model can not add much to it; however, as for seeds that are close to each other and have interaction, our method improves the estimation for large deformation as well as the model parameters.
<figure><img src="image_19.png" /><figcaption>Fig. 5: Estimation of GRID components for deformation between MRI images.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>Fig. 6: Estimation of GRID components for deformation between MRI images.</figcaption></figure>

References

1. Andresen, P.R., Nielsen, M.: Non-rigid registration by geometry-constraint diffusion. In: MICCAI. vol. 1679, pp. 533–543 (1999)

2. Arsigny, V., Pennec, X., Ayache, N.: Polyaffine transformations: a new class of diffeomorphisms for locally rigid or affine registration. In: Proc. Nat. Acad. Sci. vol. 90 (1993)
3. Chaplain, M.A.J.: Mathematical modeling of tumor growth. Springer, Seires: In-terdisciplinary Applied Mathematics (2006)

4. Clatz, O., Bondiau, P.Y., Delingette, H., Malandain, G., Sermesant, M., Warfield, S.K., Ayache, N.: In silico tumor growth: application to glioblastomas. In: MICCAI. vol. 3217, pp. 337–345 (2004)

5. Durrleman, S., Prastawa, M., Gerig, G., Joshi, S.: Optimal data-driven sparse parameterization of diffeomorphisms for population analysis. In: Inf. Process. Med. Imaging, vol. 22, pp. 123–134 (2011)

6. Grenander, U.: On the mathematics of growth. Quart. of Appl. Math. 65, 205–257 (2007)

7. Grenander, U., Miller, M.I.: Computational anatomy: an emerging discipline. Quart. Appl. Math LVI(4), 617–694 (1998)

8. Grenander, U., Srivastava, A., Saini, S.: A pattern-theoretic characterization of biological growth. IEEE Trans. Med. Imaging 26(5), 648–659 (2007)

9. Kansal, A.R., Torquato, S., IV, G.R.H., Chioccaeb, E.A., Deisboeck, T.S.: Simulated brain tumor growth dynamics using a three-dimensional cellular automaton. J. Theor. Biol. 203, 367–382 (2000)

10. Kurtek, S., Klassen, E., Ding, Z., Jacobso, S., Jacobson, J., Avison, M., Srivastava, A.: Parameterization-invariant shape comparisons of anatomical surfaces. IEEE Trans. Med. Imaging 30(3), 849–858 (2011)

11. Ling, Y., He, B.: Entropic analysis of biological growth models. IEEE Trans. Biomed. Eng. 40(12), 1193–1200 (1993)

12. Miller, M., Christensen, G., Amit, Y., Grenander, U.: Mathematical textbook of deformable neuronatomies. Proc. Nat. Acad. Sci. 90(24), 11944–11948 (1993)

13. Miller, M.I., Younes, L.: Group actions, homeomorphisms, and matching: a general framework. Int. J. Comput. Vision 41(1/2), 61–84 (2001)

14. Portman, N., Grenander, U., Vrscay, E.R.: Direct estimation of biological growth properties from image data using the “grid” model. In: ICIAR ’09. vol. 5627, pp. 832–843 (2009)

15. Portman, N., Grenander, U., Vrscay, E.R.: Grid macroscopic growth law and its application to image inference. Quart. Appl. Math. pp. 227–260 (2011)

16. Portman, N., Vrscay, E.R.: Existence and uniqueness of solutions to the grid macroscopic growth equation. Appl. Math. Comput. 217(21), 8318–8327 (2011)

17. Sra, S.: A short note on parameter approximation for Von Mises-Fisher distributions and a fast implementation of $i_s(x)$. Computational Stat. to appear (2009)

18. Srivastava, A., Saini, S., Ding, Z., Grenander, U.: Maximum-likelihood estimation of biological growth variables. In: EMMCVPR’05. vol. 3757, pp. 107–118 (2005)

19. Swanson, K.R., Bridge, C., Murray, J., Jr, E.C.A.: Virtual and real brain tumors: using mathematical modeling to quantify glioma growth and invasion. J. Neurol. Sci. 216(1), 1–10 (2003)

20. Trouvé, A.: Diffeomorphisms groups and pattern matching in image analysis. Int. J. Comput. Vision 28(3), 213–221 (1998)

## A Gradient Expressions

The terms, $\frac{\partial E^{(j)}}{\partial \phi^{(j)}} \in R^{1 \times d}$, $\frac{\partial \phi^{(j)}}{\partial \phi^{(j-1)}} \in R^{d \times d}$, $\frac{\partial \phi^{(j)}}{\partial \xi_j} \in R^{d \times d}$, $\frac{\partial \phi^{(j)}}{\partial a_j}$, and $\frac{\partial \phi^{(j)}}{\partial c_j} \in R^{d \times 1}$, that are required in gradient computation are given as follows.
$$ \frac{\partial E^{(j)}}{\partial \phi^{(j)}} = 2(\phi^{(j)} - \Phi) , $$

$$ \frac{\partial \phi^{(j)}}{\partial \phi^{(j-1)}} = a_{j-1} \cdot \exp \left( - \left( \phi^{(j-1)} - \xi_{j-1} \right)^2 / c_{j-1}^2 \right) - $$

$$ 2a_{j-1} \cdot \left( \phi^{(j-1)} - \xi_{j-1} \right)^2 / c_{j-1}^2 \cdot \exp \left( - \left( \phi^{(j-1)} - \xi_{j-1} \right)^2 / c_{j-1}^2 \right) + 1 , $$

$$ \frac{\partial \phi^{(j)}}{\partial \xi_j} = a_j \cdot \exp \left( - \left( \phi^{(j-1)} - \xi_j \right)^2 / c_j^2 \right) \cdot \left( 2 \left( \phi^{(j-1)} - \xi_j \right)^2 / c_j^2 - 1 \right) , $$

$$ \frac{\partial \phi^{(j)}}{\partial a_j} = (\phi^{(j-1)} - \xi_j) \cdot \exp \left( - \left( \phi^{(j-1)} - \xi_j \right)^2 / c_j^2 \right) , $$

$$ \frac{\partial \phi^{(j)}}{\partial c_j} = 2a_j \cdot \exp \left( - \left( \phi^{(j-1)} - \xi_j \right)^2 / c_j^2 \right) \cdot \left( \phi^{(j-1)} - \xi_j \right)^3 / c_j^3 . $$

When the ADF takes parametric forms of the second type, the gradient of the energy function with respect to the two parameters $\kappa_j$ and $\tau_{0,j}$ for $j = 1, \dots, k$ is calculated in the similar way as the other parameters. The partial derivatives from the one time deformation is shown as below.

$$ \frac{\partial \phi^{(j)}}{\partial \kappa_j} = a_j \cdot \exp (\kappa_j \cos (\tau - \tau_{0,j}) - \kappa_j) . $$

$$ \left( \phi^{(j-1)} - \xi_j \right) \cdot \exp \left( - \left( \phi^{(j-1)} - \xi_j \right)^2 / c_j^2 \right) \cdot \left( \cos(\tau - \tau_{0,j}) - 1 \right) , $$

$$ \frac{\partial \phi^{(j)}}{\partial \tau_{0,j}} = a_j \cdot \exp (\kappa_j \cos (\tau - \tau_{0,j}) - \kappa_j) . $$

$$ \left( \phi^{(j-1)} - \xi_j \right) \cdot \exp \left( - \left( \phi^{(j-1)} - \xi_j \right)^2 / c_j^2 \right) \cdot \kappa_j \cdot \sin (\tau - \tau_{0,j}) . $$
# Teichmüller Shape Descriptor and its Application to Alzheimer's Disease Study

Wei Zeng¹, Rui Shi¹, Yalin Wang², and Xianfeng David Gu¹

¹Computer Science Department, SUNY at Stony Brook

²Computer Science and Engineering, Arizona State University

**Abstract.** We propose a novel method to apply Teichmüller space theory to study the signature of a family of non-intersecting closed 3D curves on a general genus zero closed surface. Our algorithm provides an efficient method to encode both global surface and local contour shape information. The signature - Teichmüller shape descriptor - is computed by surface Ricci flow method, which is equivalent to solving an elliptic partial differential equation on surfaces and is quite stable. We propose to apply the new signature to analyze abnormalities in brain cortical morphology. Experimental results with 3D MRI data from ADNI dataset (12 healthy controls versus 12 Alzheimer's disease (AD) subjects) demonstrate the effectiveness of our method and illustrate its potential as a novel surface-based cortical morphometry measurement in AD research.

## 1 Introduction

Some neurodegenerative diseases, such as Alzheimer's disease (AD), are characterized by progressive cognitive dysfunction. The underlying disease pathology most probably precedes the onset of cognitive symptoms by many years. Efforts are underway to find early diagnostic biomarkers to evaluate neurodegenerative risk presymptomatically in a sufficiently rapid and rigorous way. Among a number of different brain imaging, biological fluid and other biomarker measurements for use in the early detection and tracking of AD, structural magnetic resonance imaging (MRI) measurements of brain shrinkage are among the best established biomarkers of AD progression and pathology.

In structural MRI studies, early researches [30, 9] have demonstrated that surface-based brain mapping may offer advantages over volume-based brain mapping work [2] to study structural features of the brain, such as cortical gray matter thickness, complexity, and patterns of brain change over time due to disease or developmental processes. In research studies that analyze brain morphology, many surface-based shape analysis methods have been proposed, such as spherical harmonic analysis (SPHARM) [11, 4], minimum description length approaches [7], medial representations (M-reps) [24], cortical gyrification index [32], shape space [21], metamorphosis [33], momentum maps [25] and conformal invariants [34], etc.; these methods may be applied to analyze shape changes or abnormalities in cortical and subcortical brain structures. Even so, a stable method to compute a global intrinsic transformation-invariant shape descriptors would be highly advantageous in this research field.
Here, we propose a novel and intrinsic method to compute the global correla-
tions between various surface region contours in Teichmüller space and apply it
to study brain morphology in AD. The proposed shape signature demonstrates
the global geometric features encoded in the interested regions, as a biomarker
for measurements of AD progression and pathology. It is based on the brain
surface conformal structure [18, 1, 13, 37] and can be accurately computed using
the surface Ricci flow method [35, 20].

1.1 Related work

In brain mapping research, volumetric measures of structures identified on 3D
MRI have been used to study group differences in brain structure and also to
predict diagnosis [2]. Recent work has also used shape-based features [21, 33, 25],
conformal invariants [34], analyzing surface changes using pointwise displace-
ments of surface meshes, local deformation tensors, or surface expansion factors,
such as the Jacobian determinant of a surface based mapping. For closed sur-
faces homotopic to a sphere, spherical harmonics have commonly been used for
shape analysis, as have their generalizations, e.g., eigenfunctions of the Laplace-
Beltrami operator in a system of spherical coordinates. These shape indices are
also rotation invariant, i.e., their values do not depend on the orientation of the
surface in space [30, 11, 28]]. Chung et al.[4] proposed a weighted spherical har-
monic representation. For a specific choice of weights, the weighted SPHARM is
shown to be the least squares approximation to the solution of an anisotropic heat
diffusion on the unit sphere. Davies et al. performed a study of anatomical shape
abnormalities in schizophrenia, using the minimal distance length approach to
statistically align hippocampal parameterizations [7]. For classification, Linear
Discriminant Analysis (LDA) or principal geodesic analysis can be used to find
the discriminant vector in the feature space for distinguishing diseased subjects
from controls. Tosun et al. [32] proposed the use of three different shape mea-
sures to quantify cortical gyrfication and complexity. Gorczowski [12] presented
a framework for discriminant analysis of populations of 3D multi-object sets. In
addition to a sampled medial mesh representation, m-rep [24], they also consid-
ered pose differences as an additional statistical feature to improve the shape
classification results.

For brain surface parameterization research, Schwartz et al. [26] and Tim-
sari and Leahy [31] computed quasi-isometric flat maps of the cerebral cortex.
Hurdal and Stephenson [18] reported a discrete mapping approach that uses cir-
cle packings to produce "flattened" images of cortical surfaces on the sphere,
the Euclidean plane, and the hyperbolic plane. Angenent et al. [1] implemented
a finite element approximation for parameterizing brain surfaces via conformal
mappings. Gu et al. [13] proposed a method to find a unique conformal mapping
between any two genus zero manifolds by minimizing the harmonic energy of
the map. The holomorphic 1-form based conformal parameterization [37] can
conformally parameterize high genus surfaces with boundaries but the result-
ing mappings have singularities. Other brain surface conformal parametrization
methods, the Ricci flow method [35] and slit map method [36] can handle surfaces
with complicated topologies (boundaries and landmarks) without singularities.
Wang et al. [34] applied the Yamabe flow method to study statistical group differences in a group of 40 healthy controls and 40 subjects with Williams syndrome, showing the potential of these surface-based descriptors for localizing cortical shape abnormalities in genetic disorders of brain development.

Conformal mappings have been applied in computer vision for modeling the 2D shape space by Sharon and Mumford [27]. The image plane is separated by a 2D contour, both interior and exterior are conformally mapped to disks, then the contour induces a diffeomorphism of the unit circle, which is the signature of the contour. The signature is invariant under translations and scalings, and able to recover the original contour by conformal welding. Later, this method is generalized to model multiple 2D contours with inner holes in [22]. To the best of our knowledge, our method is the first one to generalized Sharon and Mumford's 2D shape space to 3D surfaces.

## 1.2 Our Approach

For a 3D surface, all the contours represent the 'shape' of the surface. Inspired by the beautiful research work of Sharon and Mumford [27] on 2D shape analysis (recently it has been generalized to model multiple 2D contours [22]), we build a Teichmüller space for 3D shapes by using conformal mappings. In this Teichmüller space, every 3D contour (a simple closed curve) is represented by a point in the space; each point denotes a unique equivalence class of diffeomorphisms up to a Möbius transformation. For a 3D surface, the diffeomorphisms of all the contours form a global shape representation of the surface. By using this signature, the similarities of 3D shapes can be quantitatively analyzed, therefore, the classification and recognition of 3D objects can be performed from their observed contours.

We tested our algorithm in some segmented regions on a set of brain left cortical surfaces extracted from 3D anatomical brain MRI scans. The proposed method can reliably compute signatures on two cortical functional areas by computing the diffeomorphisms of each observed contour. Using the signature as the statistics, our method achieve about 92% accuracy rate to discriminate a set of AD subjects from healthy control subjects.

To the best of our knowledge, it is the first work to apply contour diffeomorphism to brain morphometry research. Our experimental results demonstrated that this novel and simple method may be useful to analyze certain functional areas, and it may shed some lights on understanding detecting abnormality regions in brain surface morphometry. Our major **contributions** in this work include:

1. A new method to compute Teichmüller shape descriptor, in a way that generalized a prior 2D domain conformal mapping work [27].

2. The method is theoretically rigorous and general. It presents a stable way to calculate the diffeomorphisms of contours in general 3D surfaces based on Ricci flow.

3. It involves solving elliptic partial differential equations (PDEs), so it is numerically efficient and computationally stable.
<figure><img src="image_21.png" /><figcaption>**Fig. 1.** Diffeomorphism signature via uniformization mapping for a genus zero surface with 3 simple closed contours $\gamma_1, \gamma_2, \gamma_3$ in (a), which correspond to the boundaries $c_1, c_2, c_3$ of the circle domains $D_1, D_2, D_3$ in (c), respectively. These three contours are also mapped to the boundaries of the base circle domain $D_0$ in (b). The curves in (d) demonstrate the diffeomorphisms for the three contours.</figcaption></figure>

4. The shape descriptors are global and invariant to rigid motion and conformal deformations.

**Pipeline.** Figure 1 shows the pipeline for computing the diffeomorphism signature for a surface with 3 closed contours. Here, we use a human brain hemisphere surface whose functional areas are divided and labeled in different color. The contours (simple closed curves) of functional areas can be used to slice the surface open to connected patches. As shown in frames (a-c), three contours $\gamma_1, \gamma_2, \gamma_3$ are used to divide the whole brain (a genus zero surface $S$) to 4 patches $S_0, S_1, S_2, S_3$; each of them is conformally mapped to a circle domain (e.g., disk or annuli), $D_0, D_1, D_2, D_3$. Note that $\gamma_1$ is the contour of the joint functional areas of precuneus and posterior cingulate. One contour is mapped to two unit circles in two mappings. The representation of the shape according to each contour is a diffeomorphism of the unit circle to itself, defined as the mapping between periodic polar angles ($Angle_1, Angle_2$), $Angle_1, Angle_2 \in [0, 2\pi]$. The proper normalization is employed to remove Möbius ambiguity. The diffeomorphisms induced by the conformal maps of each curve form a diffeomorphism signature, which is the Teichmüller coordinates in Teichmüller space. As shown
in (d), the curves demonstrate the diffeomorphisms for three contours; the area distance is defined as the metric for shape comparison and classification.

# 2 Theoretical Background

In this section, we briefly introduce the theoretical foundations necessary for the current work. For more details, we refer readers to the classical books [10, 16].

## 2.1 Surface Uniformization Mapping

Conformal mapping between two surfaces preserves angles. Suppose $(S_1, g_1)$ and $(S_2, g_2)$ are two surfaces embedded in $\mathbb{R}^3$, $g_1$ and $g_2$ are the Euclidean induced Riemannian metrics. A mapping $\phi: S_1 \to S_2$ is called conformal, if the pull back metric of $g_2$ induced by $\phi$ on $S_1$ differs from $g_1$ by a positive scalar function: $\phi^*g_2 = e^{2\lambda}g_1$, where $\lambda: S_1 \to \mathbb{R}$ is a scalar function, called the conformal factor.

For example, all the conformal automorphisms of the unit disk form the Möbius transformation group of the disk, each mapping is given by

$$z \to e^{i\theta} \frac{z - z_0}{1 - \bar{z}_0 z}.$$

All the conformal automorphism group of the extended complex plane $\mathbb{C} \cup \{\infty\}$ is also called Möbius transformation group, each mapping is given by

$$z \to \frac{az+b}{cz+d}, \quad ad-bc = 1, a, b, d, c \in \mathbb{C}.$$

By stereo-graphic projection, the unit sphere can be conformally mapped to the extended complex plane. Therefore, the Möbius transformation group is also the conformal automorphism group of the unit sphere.

A circle domain on the complex plane is the unit disk with circular holes. A circle domain can be conformally transformed to another circle domain by Möbius transformations, $z \to e^{i\theta} \frac{z-z_0}{1-\bar{z}_0 z}$. All genus zero surfaces with boundaries can be conformally mapped to circle domains:

**Theorem 21 (Uniformization)** Suppose $S$ is a genus zero Riemannian surface with boundaries, then $S$ can be conformally mapped onto a circle domain. All such conformal mappings differ by a Möbius transformation on the unit disk.

This theorem can be proved using Ricci flow straightforwardly. Therefore, the conformal automorphism group of $S Conf(S)$ is given

$$Conf(S) := \{\phi^{-1} \circ \tau \circ \phi|\tau \in \text{Möb}(\mathbb{S}^2)\}.$$

## 2.2 Teichmüller Space

**Definition 22 (Conformal Equivalence)** Suppose $(S_1, g_1)$ and $(S_2, g_2)$ are two Riemannian surfaces. We say $S_1$ and $S_2$ are conformally equivalent if there is a conformal diffeomorphism between them.
All Riemannian surfaces can be classified by the conformal equivalence relation. Each conformal equivalence class shares the same conformal invariants, the so-called conformal module. The conformal module is one of the key component for us to define the unique shape signature.

**Definition 23 (Teichmüller Space)** *Fixing the topology of the surfaces, all the conformal equivalence classes form a manifold, which is called the Teichmüller space.*

For example, all topological disks (genus zero Riemannian surfaces with single boundary) can be conformally mapped to the planar disk. Therefore, the Teichmüller space for topological disks consists of a single point.

Suppose a genus zero Riemannian surface $S$ has $n$ boundary components $\{\gamma_1, \gamma_2, \cdots, \gamma_n\}$, $\partial S = \gamma_1 + \gamma_2 + \cdots + \gamma_n$, $\phi: S \to \mathbb{D}$ is the conformal mapping that maps $S$ to a circle domain $\mathbb{D}$, such that (a). $\phi(\gamma_1)$ is the exterior boundary of the $\mathbb{D}$; (b) $\phi(\gamma_2)$ centers at the origin; (c) The center of $\phi(\gamma_3)$ is on the imaginary axis. Then the conformal module of the surface $S$ (also the circle domain $\mathbb{D}$) is given by $\mathrm{Mod}(S) = \{(c_i, r_i)|i = 1, 2, \cdots, n\}$. This shows the Teichmüller space of genus zero surfaces with $n$ boundaries is of $3n - 6$ dimensional. The Teichmüller space has a so-called Weil-Peterson metric [27], so it is a Riemannian manifold. Furthermore it is with negative sectional curvature, therefore, the geodesic between arbitrary two points is unique.

## 2.3 Surface Ricci Flow

Surface Ricci flow is the powerful tool to compute uniformization. Ricci flow refers to the process of deforming Riemannian metric $g$ proportional to the curvature, such that the curvature $K$ evolves according to a heat diffusion process, eventually the curvature becomes constant everywhere. Suppose the metric $g = (g_{ij})$ in local coordinate. Hamilton [15] introduced the Ricci flow as

$$ \frac{dg_{ij}}{dt} = -Kg_{ij}. $$

Surface Ricci flow conformally deforms the Riemannian metric, and converges to constant curvature metric [3]. Furthermore, Ricci flow can be used to compute the unique conformal Riemannian metric with the prescribed curvature.

**Theorem 24 (Hamilton and Chow [3])** *Suppose $S$ is a closed surface with a Riemannian metric. If the total area is preserved, the surface Ricci flow will converge to a Riemannian metric of constant Gaussian curvature.*

## 2.4 Teichmüller Shape Descriptor

Suppose $\Gamma = \{\gamma_0, \gamma_1, \cdots, \gamma_n\}$ is a set of non-intersecting smooth closed curves on a genus zero closed surface. $\Gamma$ segments the surface to a set of connected components $\{\Omega_0, \Omega_1, \cdots, \Omega_n\}$, each segment $\Omega_i$ is a genus zero surface with boundary components. Construct the uniformization mapping $\phi_k: \Omega_k \to \mathbb{D}_k$
to map each segment $\Omega_k$ to a circle domain $\mathbb{D}_k$, $0 \le k \le n$. Assume $\gamma_i$ is the common boundary between $\Omega_j$ and $\Omega_k$, then $\phi_j(\gamma_i)$ is a circular boundary on the circle domain $\mathbb{D}_j$, $\phi_k(\gamma_i)$ is another circle on $\mathbb{D}_k$. Let $f_i|_{\mathbb{S}^1} := \phi_j \circ \phi_k^{-1}|_{\mathbb{S}^1} : \mathbb{S}^1 \to \mathbb{S}^1$ be the diffeomorphism from the circle to itself. We called the diffeomorphism $f_i$ the signature of $\gamma_i$.

**Definition 25 (Signature of a Family of Loops)** *The signature of a family non-intersecting closed 3D curves $\Gamma = \{\gamma_0, \gamma_1, \dots, \gamma_k\}$ on a genus zero closed surface is defined as: $S(\Gamma) := \{f_0, f_1, \dots, f_k\} \cup \{\mathrm{Mod}(\mathbb{D}_0), \mathrm{Mod}(\mathbb{D}_1), \dots, \mathrm{Mod}(\mathbb{D}_k)\}$.*

The following **main theorem** plays fundamental role for the current work. Note that if a circle domain $\mathbb{D}_k$ is disk, its conformal module can be omitted from the signature.

**Theorem 26 (Main Theorem)** *The family of smooth 3D closed curves $\Gamma$ on a genus zero closed Riemannian surface is determined by its signature $S(\Gamma)$, unique up to a conformal automorphism of the surface $\eta \in \mathrm{Conf}(S)$.*

The proof of Theorem 26 can be found in the appendix section.

The theorem states that the proposed signature determine shapes up to a Möbius transformation. We can further do a normalization that fixes $\infty$ to $\infty$ and that the differential carries the real positive axis at $\infty$ to the real positive axis at $\infty$, as in Sharon and Mumford's paper [27]. The signature can then determine the shapes uniquely up to translation and scaling.

The shape signature $S(\Gamma)$ gives us a complete representation for the space of shapes. It inherits a natural metric. Given two shapes $\Gamma_1$ and $\Gamma_2$. Let $S(\Gamma_i) := \{f_0^i, f_1^i, \dots, f_k^i\} \cup \{\mathrm{Mod}(\mathbb{D}_0^i), \mathrm{Mod}(\mathbb{D}_1^i), \dots, \mathrm{Mod}(\mathbb{D}_k^i)\}$ ($i = 1, 2$). We can define a metric $d(S(\Gamma_1), S(\Gamma_2))$ between the two shape signatures using the natural metric in the Teichmüller space. Our signature is stable under geometric noise. Our algorithm depends on conformal maps from surfaces to circle domains using discrete Ricci flow method.

# 3 Algorithm

In this section, we explain each step of the pipeline in Figure 1 in details.

## 3.1 Circular Uniformization Mapping

We apply discrete Ricci flow method [20] to conformally map the surfaces onto planar circle domains $\phi_k : S_k \to \mathbb{D}$. The surface is represented as a triangle mesh $\Sigma$. A discrete Riemannian metric is represented as the edge length.

We associate each vertex $v_i$ with a circle $(v_i, \gamma_i)$, where $\gamma_i$ is the radius. Let $u_i = \log \gamma_i$ be the discrete conformal factor. The discrete Ricci flow is defined as follows:

$$ \frac{du_i(t)}{dt} = (\bar{K}_i - K_i), \qquad (1) $$

where $\bar{K}_i$ is the user defined target curvature and $K_i$ is the curvature induced by the current metric. The discrete Ricci flow has exactly the same form as the
smooth Ricci flow, which conformally deforms the discrete metric according to
the Gaussian curvature. The computation is based on circle packing metric [20].

Suppose $\Sigma$ is a genus zero mesh with multiple boundary components. The uniformization conformal mapping $\phi: \Sigma \to \mathcal{D}$, where $\mathcal{D}$ is the circle domain, can be computed using Ricci flow by setting the prescribed curvature as follows: (a) The geodesic curvature on the exterior boundary is +1 everywhere; (b) the geodesic curvature on other boundaries are negative constants; (c) the Gaussian curvature on interior points are zeros everywhere. We use this method to compute conformal mapping, and get conformal module and shape descriptor. The main challenge is that the target curvature is dynamically determined by the metric. The metric is evolving, so is the target curvature. The detailed algorithm is reported in [38].

## 3.2 Computing Shape Descriptor

After the computation of the conformal mapping, each connected component is mapped to a circle domain. We define an order for all the loops on the surface, this induces an order for all the boundary components on each segment. Then by the definition for the conformal module of a circle domain, we normalize each circle domain using a Möbius transformation, then compute the conformal modules directly. For those segments, which are simply connected and mapped to the unit disk, we compute its mass center, and use a Möbius transformation to map the center to the origin.

Each loop on the surface becomes the boundary components on two segments, both boundary components are mapped to a circle under the uniformization mapping. Then we compute the signature directly.

# 4 Experimental Results

We demonstrate the efficiency and efficacy of our method by analyzing the human brain cortexes of Alzheimer's disease (AD) and healthy control subjects. The brain surfaces are represented as triangular meshes; a half brain with 100K triangles. We implement the algorithm using generic C++ on windows XP platform, with Intel Xeon CPU 3.39 GHz, 3.98 G RAM. The numerical systems are solved using Matlab C++ library. In general, the signature calculation on each half brain surface with 2 or 3 contours on each half takes less than 1 minute to compute, even on complicated domains.

*Data and preprocessing.* The experimental data include 12 Alzheimer disease patients and 12 healthy control subjects. The structural MRI images were from the AD Neuroimaging Initiative (ADNI [19, 23]. We used Freesurfer's automated processing pipeline [6] for automatic skull stripping, tissue classification, surface extraction, cortical and subcortical parcellations. It calculates volumes of individual grey matter parcellations in *mm*³ and surface area in *mm*². It also provides surface and volume statistics for about 34 different cortical structures, and also computes geometric characteristics such as curvature, curvedness, local foldedness for each of the parcellations [8]. In this work, we studied segmented surface regions for group difference analysis.
<figure><img src="image_22.png" /><figcaption>**Fig. 2.** Diffeomorphism signature ($d_{left}$, $d_{right}$) of a healthy control brain cortex. Each (**left** and **right**) half brain is a genus zero surface with 2 contours.</figcaption></figure>

**Quantitative analysis.** Figure 2 shows an example of diffeomorphism signatures for a brain cortical surface. We selected two contours on the left and the right half brain cortical surfaces, which correspond to superior temporal and the joint areas of precuneus and posterior cingulate. Early researches [17, 14] have indicated that these two areas may have significant atrophy in AD group. These two contours segment a brain hemisphere surface to 3 patches; one topological annulus (called the base domain), two topological disks. The base domain with two boundaries is mapped to an annulus, one boundary to exterior unit circle, the other one to the inner concentric circle. The diffeomorphism signature for each contour is plotted as a monotonic curve within the square $[0, 2\pi] \times [0, 2\pi]$. The area difference between the plotted curves, $d = \int_0^{2\pi} (\text{Angle}_2^2 - \text{Angle}_1^2)d\text{Angle}_1$, is used as the metric to represent the global shape of both contours. So the signature of the whole brain surface is represented as a pair ($d_{left}$, $d_{right}$) for combining the left-hemisphere and right-hemisphere brain shape signatures. The method was tested on 12 AD subjects and 12 healthy subjects, with mean signatures (3.6827, -7.12957) and (5.2752, -5.6036), respectively. Figure 3 shows that with a simple linear discriminant analysis (LDA) model, there were only two subjects that were not correctly classified. It demonstrates that the proposed global diffeomorphism signature of contours is very efficient and may be effective to differentiate the shapes within healthy control and AD subject groups.
<figure><img src="image_23.png" /><figcaption>**Fig. 3.** Distribution of diffeomorphism signature for 12 AD (in red) and 12 healthy control (CTL) (in blue) subjects. Each point denotes the diffeomorphism signature value ($d_{left}, d_{right}$) for a whole brain surface, computed as in Figure 2.</figcaption></figure>

*Discussion.* The proposed work is based on surface Ricci flow research. Computing the conformal module is equivalent to solving an elliptic partial differential equation on surfaces. According to PDE theory, the solution is smoother than its boundary conditions, so the solution process is quite stable.

For surface-based AD research, the state-of-the-art work has used cortical thickness as the measurement [29, 5]. However, recent research [39] indicated that the commonly used cortical thickness and cortical area measurements are genetically and phenotypically independent. The biological meaning of the proposed shape signature is closely related to brain atrophy so it is more related to cortical area changes. Our method provides a unique and intrinsic shape feature to study brain morphometry changes caused by brain atrophy. It studies the sensitivity and reproducibility of shape features computed in the entire brain surface domain. The gained insights help improve our understanding to AD related pathology and discover the precise etiology of the grey matter changes. The preliminary results demonstrated that the shape signature provides a reasonably good discriminant power for AD biomarker research. We currently studied the superior temporal area, which is directly related to medial temporal lobe atrophy. The method can be equally applied to other regions as well. In future, we may study/compare other functional areas in the medial temporal lobe.

# 5 Conclusion

In this paper, we propose a novel method that computes the global shape signatures on specified functional areas on brain cortical surfaces in Teichmüller space. In the future, we will further explore and validate other applications of this global correlation shape signature in neuroimaging and shape analysis research.
References

1. S. Angenent, S. Haker, R. Kikinis, and A. Tannenbaum. Nondistorting flattening maps and the 3D visualization of colon CT images. *IEEE Trans. Med. Imag.*, 19:665–671, 2000.

2. J. Ashburner, C. Hutton, R. Frackowiak, I. Johnsrude, C. Price, and K. Friston. Identifying global anatomical differences: deformation-based morphometry. *Human Brain Mapping*, 6:348–357, 1998.

3. B. Chow, P. Lu, and L. Ni. *Hamilton's Ricci Flow*. American Mathematical Society, 2006.

4. M. K. Chung, K. M. Dalton, and R. J. Davidson. Tensor-based cortical surface morphometry via weighted spherical harmonic representation. *IEEE Trans. Med. Imag.*, 27:1143–1151, 2008.

5. R. Cuingnet, E. Gerardin, J. Tessieras, G. Auzias, S. Lehricy, M.-O. Habert, M. Chupin, H. Benali, and O. Colliot. Automatic classification of patients with Alzheimer's disease from structural MRI: A comparison of ten methods using the ADNI database. *NeuroImage*, 56(2):766–781, 2011.

6. A. M. Dale, B. Fischl, and M. I. Sereno. Cortical surface-based analysis I: segmentation and surface reconstruction. *Neuroimage*, 27:179–194, 1999.

7. R. H. Davies, C. J. Twining, P. D. Allen, T. F. Cootes, and C. J. Taylor. Shape discrimination in the hippocampus using an MDL model. *in Proc. Infor. Proc. Med. Imag. (IPMI)*, 2003.

8. R. S. Desikan, F. Segonne, B. Fischl, B. T. Quinn, B. C. Dickerson, D. Blacker, R. L. Buckner, A. M. Dale, R. P. Maguire, B. T. Hyman, M. S. Albert, and R. J. Killiany. An automated labeling system for subdividing the human cerebral cortex on MRI scans into gyral based regions of interest. *Neuroimage*, 31:968–80, 2006.

9. B. Fischl, M. I. Sereno, and A. M. Dale. Cortical surface-based analysis II: Inflation, flattening, and a surface-based coordinate system. *NeuroImage*, 9:195–207, 1999.

10. F. P. Gardiner and N. Lakic. *Quasiconformal Teichmuller theory*. American Mathematical Society, 2000.

11. G. Gerig, M. Styner, D. Jones, D. Weinberger, and J. Lieberman. Shape analysis of brain ventricles using SPHARM. *in Proc. MMBIA 2001*, pages 171–178, 2001.

12. K. Gorczowski, M. Styner, J.-Y. Jeong, J. S. Marron, J. Piven, H. C. Hazlett, S. M. Pizer, and G. Gerig. Statistical shape analysis of multi-object complexes. *IEEE Conf. Comp. Vis. Patt. Recog. CVPR '07*, pages 1–8, 2007.

13. X. Gu, Y. Wang, T. F. Chan, P. M. Thompson, and S.-T. Yau. Genus zero surface conformal mapping and its application to brain surface mapping. *IEEE Trans. Med. Imag.*, 23:949–958, 2004.

14. X. Guo, Z. Wang, K. Li, Z. Li, Z. Qi, Z. Jin, L. Yao, and K. Chen. Voxel-based assessment of gray and white matter volumes in Alzheimer's disease. *Neurosci Lett*, 468:146–50, 2010.

15. R. S. Hamilton. The Ricci flow on surfaces. *Mathematics and general relativity*, 71:237–262, 1988.

16. P. Henrici. *Applied and Computational Complex Analysis*, volume 3. Wiley-Interscience, 1988.

17. X. Hua, S. Lee, D. P. Hibar, I. Yanovsky, A. D. Leow, A. W. Toga, C. R. J. Jr, M. A. Bernstein, E. M. Reiman, D. J. Harvey, J. Kornak, N. Schuff, G. E. Alexander, M. W. Weiner, and P. M. Thompson. Mapping Alzheimer's disease progression in 1309 MRI scans: Power estimates for different inter-scan intervals. *Neuroimage*, 51:63–75, 2010.

18. M. K. Hurdal and K. Stephenson. Cortical cartography using the discrete conformal approach of circle packings. *NeuroImage*, 23:S119–S128, 2004.
19. C. R. J. Jack, M. A. Bernstein, N. C. Fox, P. M. Thompson, G. Alexander, D. Harvey, B. Borowski, P. J. Britson, J. L. Whitwell, C. Ward, and e. al. The Alzheimer's disease neuroimaging initiative (ADNI): MRI methods. *J. of Mag. Res. Ima.*, 27:685-691, 2007.

20. M. Jin, J. Kim, F. Luo, and X. Gu. Discrete surface Ricci flow. *IEEE Trans. Vis. Comput. Graphics*, 14:1030–1043, September 2008.

21. X. Liu, Y. Shi, I. Dinov, and W. Mio. A computational model of multidimensional shape. *Int J Comput Vis*, 89:69–83, 2010.

22. L. M. Lui, W. Zeng, S.-T. Yau, and X. Gu. Shape analysis of planar objects with arbitrary topologies using conformal geometry. In *ECCV 2010*, 2010.

23. S. G. Mueller, M. W. Weiner, L. J. Thal, R. C. Petersen, C. Jack, W. Jagust, J. Q. Trojanowski, A. W. Toga, and L. Beckett. The Alzheimer's disease neuroimaging initiative. *Neuroimaging clinics of North America*, 15:869–877, 2005.

24. S. Pizer, D. Fritsch, P. Yushkevich, V. Johnson, and E. Chaney. Segmentation, registration, and measurement of shape variation via image object shape. *IEEE Trans. Med. Imag.*, 18:851–865, 1999.

25. A. Qiu and M. I. Miller. Multi-structure network shape analysis via normal surface momentum maps. *NeuroImage*, 42:1430–8, 2008.

26. E. L. Schwartz, A. Shaw, and E. Wolfson. A numerical solution to the generalized Mapmaker's problem: Flattening nonconvex polyhedral surfaces. *IEEE Trans. Patt. Anal. Mach. Intell.*, 11:1005–1008, 1989.

27. E. Sharon and D. Mumford. 2D-shape analysis using conformal mapping. *Int. J. Comput. Vision*, 70:55–75, October 2006.

28. L. Shen, A. J. Saykin, M. K. Chung, and H. Huang. Morphometric analysis of hippocampal shape in mild cognitive impairment: An imaging genetics study. *in IEEE 7th International Conference Bioinformatics and Bioengineering*, 2007.

29. P. M. Thompson, K. M. Hayashi, G. D. Zubicaray, A. L. Janke, S. E. Rose, J. Semple, D. Herman, M. S. Hong, S. S. Dittmer, D. M. Doddrell, and A. W. Toga. Dynamics of gray matter loss in Alzheimer's disease. *J. Neuroscience*, 23:994–1005, 2003.

30. P. M. Thompson and A. W. Toga. A surface-based technique for warping 3-dimensional images of the brain. *IEEE Trans. Med. Imag.*, 15:1–16, 1996.

31. B. Timsari and R. M. Leahy. Optimization method for creating semi-isometric flat maps of the cerebral cortex. *Medical Imaging 2000: Image Processing*, 39:979–987, 2000.

32. D. Tosun, A. Reiss, A. D. Lee, R. A. Dutton, K. M. Hayashi, U. Bellugi, A. M. Galaburda, J. R. Korenberg, D. L. Mills, A. W. Toga, and P. M. Thompson. Use of 3-D cortical morphometry for mapping increased cortical gyrification and complexity in Williams syndrome. *Biom. Imag.: Nano to Macro*, 2006. 3rd IEEE Int. Symp. on, 2006, pages 1172–1175, 2006.

33. A. Trouve and L. Younes. Metamorphoses through Lie group action. *Foundations of Computational Mathematics*, 5:173–198, 2005.

34. Y. Wang, X. Gu, T. F. Chan, and P. M. Thompson. Shape analysis with conformal invariants for multiply connected domains and its application to analyzing brain morphology. *IEEE Conf. Comp. Vis. Patt. Recog. CVPR '09*, pages 202–209, 2009.

35. Y. Wang, X. Gu, T. F. Chan, P. M. Thompson, and S.-T. Yau. Brain surface conformal parameterization with algebraic functions. *Med. Image Comp. Comput.- Assist. Intervention, Proceedings, Part II*, pages 946–954, 2006.

36. Y. Wang, X. Gu, T. F. Chan, P. M. Thompson, and S.-T. Yau. Conformal slit mapping and its applications to brain surface parameterization. *Med. Image Comp. Comput.-Assist. Intervention, Proceedings, Part I*, pages 585–593, 2008.

37. Y. Wang, L. Lui, X. Gu, K. M. Hayashi, T. F. Chan, A. W. Toga, P. M. Thompson, and S.-T. Yau. Brain surface conformal parameterization using Riemann surface structure. *IEEE Trans. Med. Imag.*, 26:853–865, 2007.
<figure><img src="image_24.png" /><figcaption>**Fig. 4.** Proof for the main theorem, the signature uniquely determines the family of closed curves unique up to a Möbius transformation.</figcaption></figure>

38. Y. Wang, J. Shi, X. Yin, X. Gu, T. F. Chan, S.-T. Yau, A. W. Toga, and P. M. Thompson. Brain surface conformal parameterization with the Ricci flow. submitted to *IEEE Trans. Med. Imag.*, 2011.

39. A. M. Winkler, P. Kochunov, J. Blangero, L. Almasy, K. Zilles, P. T. Fox, R. Dugirala, and D. C. Glahn. Cortical thickness or grey matter volume? The importance of selecting the phenotype for imaging genetics studies. *NeuroImage*, 53(3):1135–1146, 2010.

## Appendix: Proof of Theorem 26

*proof* See Figure 4. In the left frame, a family of planar smooth curves $\Gamma = \{\gamma_0, \dots, \gamma_5\}$ divide the plane to segments $\{\Omega_0, \Omega_1, \dots, \Omega_6\}$, where $\Omega_0$ contains the $\infty$ point. We represent the segments and the curves as a tree in the second frame, where each node represents a segment $\Omega_k$, each link represents a curve $\gamma_i$. If $\Omega_j$ is included by $\Omega_i$, and $\Omega_i$ and $\Omega_j$ shares a curve $\gamma_k$, then the link $\gamma_k$ in the tree connects $\Omega_j$ to $\Omega_i$, denoted as $\gamma_k : \Omega_i \to \Omega_j$. In the third frame, each segment $\Omega_k$ is mapped conformally to a circle domain $D_k$ by $\Phi_k$. The signature for each closed curve $\gamma_k$ is computed $f_{ij} = \Phi_i \circ \Phi_j^{-1}|_{\gamma_k}$, where $\gamma_k : \Omega_i \to \Omega_j$ in the tree. In the last frame, we construct a Riemann sphere by gluing circle domains $D_k$'s using $f_{ij}$'s in the following way. The gluing process is of bottom up. We first glue the leaf nodes to their fathers. Let $\gamma_k : D_i \to D_j, D_j$ be a leaf of the tree. For each point $z = re^{i\theta}$ in $D_j$, the extension map: $G_{ij}(re^{i\theta}) = re^{f_{ij}(\theta)}$.

We denote the image of $D_j$ under $G_{ij}$ as $S_j$. Then we glue $S_j$ with $D_i$. By repeating this gluing procedure bottom up, we glue all leafs to their fathers. Then we prune all leaves from the tree. Then we glue all the leaves of the new tree, and prune again. By repeating this procedure, eventually, we get a tree with only the root node, then we get a Riemann sphere, denoted as $S$. Each circle domain $D_k$ is mapped to a segment $S_k$ in the last frame, by a sequence of extension maps. Suppose $D_k$ is a circle domain, a path from the root $D_0$ to $D_k$ is $\{i_0 = 0, i_1, i_2, \dots, i_n = k\}$, then the map from $G_k : D_k \to S_k$ is given by: $G_k = G_{i_0 i_1} \circ G_{i_1 i_2} \circ \dots \circ G_{i_{n-1} i_n}$. Note that, $G_0$ is identity. Then the Beltrami coefficient of $G_k^{-1} : S_k \to D_k$ can be directly computed, denoted as $\mu_k : S_k \to \mathbb{C}$. The composition $\Phi_k \circ G_k^{-1} : S_k \to \Omega_k$ maps $S_k$ to $\Omega_k$, because $\Phi_k$ is conformal, therefore the Beltrami coefficient of $\Phi_k \circ G_k^{-1}$ equals to $\mu_k$.

We want to find a map from the Riemann sphere $S$ to the original Riemann sphere $\Omega$, $\Phi : S \to \Omega$. The Beltrami-coefficient $\mu : S \to \mathbb{C}$ is the union of $\mu_k$'s
each segments: $\mu(z) = \mu_k(z), \forall z \in S_k$. The solution exists and is unique up to a Möbius transformation according to Quasi-conformal Mapping theorem [10].

Note that, the discrete computational method is more direct without explicitly solving the Beltrami equation. From the Beltrami coefficient $\mu$, one can deform the conformal structure of $S_k$ to that of $\Omega_k$, under the conformal structures of $\Omega_k$, $\Phi: S \to \Omega$ becomes a conformal mapping. The conformal structure of $\Omega_k$ is equivalent to that of $D_k$, therefore, one can use the conformal structure of $D_k$ directly. In discrete case, the conformal structure is represented as the angle structure. Therefore in our algorithm, we copy the angle structures of $D_k$'s to $S$, and compute the conformal map $\Phi$ directly.
# On a Geodesic Equation for Planar Conformal
Template Matching

Stephen Marsland¹, Robert McLachlan², Klas Modin¹, and Matthew
Perlmutter²

¹ School of Engineering and Advanced Technology (SEAT)
Massey University, Palmerston North, New Zealand

² Institute of Fundamental Sciences (IFS)
Massey University, Palmerston North, New Zealand

**Abstract.** In this paper we consider planar conformal deformations, motivated by the warps that Wentworth Thompson used to deform images of one species into another. We study an equation for geodesic motion on the infinite dimensional Fréchet manifold Con($D, \mathbb{R}^2$) of conformal embeddings of the disk into the plane. We demonstrate that solutions may be represented as sheets, and use the sheet ansatz to derive a numerical discretization scheme. We also show that the equation admits totally geodesic solutions corresponding to scaling and translation, but not to affine transformations.

## 1 Introduction

The use of diffeomorphic transformations in both image registration and shape analysis is now common and utilised in many machine vision and image analysis tasks. One image or shape is brought into alignment with another by deforming the image until some similarity measure reaches a minimum. The deformation is computed as a geodesic with respect to some metric on the diffeomorphism group. For a general treatment and an overview of the subject see [1] and references therein.

The standard approach to the deformation method is to first perform an affine registration (principally to remove translation and rotation), and then to seek a diffeomorphic warp of the image. However, in what is arguably the most influential demonstration of the application of warping methods – D’Arcy Wentworth Thompson’s seminal book ‘On Growth and Form’ [2] – Thompson transforms images of one species into another using relatively simple warps, so that the gross features of the two match. In a recent review of his work, biologist Arthur Wallace says:

> *This theory cries out for causal explanation, which is something the great man eschewed. [...] His transformations suggest coordinated rather than piece-meal changes to development in the course of evolution, an issue which almost completely disappeared from view in the era of the 'modern synthesis' of evolutionary theory, but which is of central importance again in the era of evo-devo.*
[...] All the tools are now in place to examine the mechanistic basis of transformations. Not only do we have phylogenetic systematics and evo-devo, but, so obvious that it is easy to forget, we have computers, and especially, in this context, advanced computer graphics. We owe it to the great man to put these three things together to investigate the mechanisms that produce the morphological changes that he captured so elegantly with little more than sheets of graph paper and, of course, a brilliant mind. [3]

<table><caption>**Table 1.** Transformation groups used in some transformations in Chapter XII, 'On the Theory of Transformations, or the Comparison of Related Forms', of [2]</caption><thead><tr><th>Figure no. in [2]</th><th>Transformation group</th></tr></thead><tbody><tr><td>515</td><td>x &mapsto; ax, y &mapsto; y</td></tr><tr><td>513.2</td><td>x &mapsto; ax, y &mapsto; by</td></tr><tr><td>509, 510, 518</td><td>x &mapsto; ax, y &mapsto; cx + dy (shears)</td></tr><tr><td>521-22, 513.5</td><td>x &mapsto; ax + by, y &mapsto; cx + dy (affine)</td></tr><tr><td>506,508</td><td>x &mapsto; ax, y &mapsto; g(y)</td></tr><tr><td>511</td><td>x &mapsto; f(x), y &mapsto; g(y)</td></tr><tr><td>517-20, 523, 513.1, 513.3, 513.4, 513.6, 514, 525</td><td>conformal</td></tr><tr><td>524</td><td>'peculiar'</td></tr></tbody></table>

We draw attention to two key aspects of Thompson's examples: (i) the transformations are as simple as possible to achieve what he considers a good enough match (see Table 1); and (ii) the classes of transformations that he considers all forms groups (or pseudo groups), either finite or infinite dimensional. Mostly, he uses conformal transformations, a constraint he is reluctant to give up<fnref n="3" />.

In terms of diffeomorphic image matching, Figs. 506, 508 and 511 of [2] are related to the one-dimensional diffeomorphism group, and hence to the Camassa-Holm family of equations [4]; but the groups usually studied in the literature are the full diffeomorphism group (two functions of two variables) and the volume-preserving group (one function of two variables). Conformal maps are defined by two functions of one variable: a drastically reduced dimensionality.

For applications in image registration we therefore suggest to vary the group from which warps are drawn as well as the metric. If a low-dimensional group gives a close match, then it should be preferred over a similar match from a higher-dimensional group; if necessary, local deformations from the full diffeo-

<footnote n="3">³ "It is true that, in a mathematical sense, it is not a perfectly satisfactory or perfectly regular deformation, for the system is no longer isogonal; but[...] approaches to an isogonal system under certain conditions of friction or constraint." ([2], p. 1064) "[...] it will perhaps be noticed that the correspondence is not always quite accurate in small details. It could easily have been made much more accurate by giving a slightly sinuous curvature to certain of the coordinates. But as they stand, the correspondence indicated is very close, and the simplicity of the figures illustrates all the better the general character of the transformation." (ibid., p. 1074).</footnote>
morphism group can be added later. In this paper we consider the case of conformal transformations.

Although the composition of two conformal maps $\mathbb{R}^2 \to \mathbb{R}^2$ is conformal, it need not be invertible: we need to restrict the domain. The invertible conformal maps $\mathbb{D} \to \mathbb{D}$ do form a group, the disk-preserving Möbius group, but it is only 3 dimensional. We are therefore led to consider the infinite dimensional configuration space $\text{Con}(\mathbb{D}, \mathbb{R}^2)$ of planar conformal embeddings of the closed unit disk $\mathbb{D}$ into the plane. This is not a group, but it is a pseudo group.

## 2 Derivation of the Weak Geodesic Equation

In this section we give a derivation of the weak form of the governing equation.

Let $g$ denote the Euclidean metric on $\mathbb{R}^2$, i.e., in Cartesian coordinates we have $g = dx \otimes dx + dy \otimes dy$. Then $(\mathbb{D}, g)$ is a compact Riemannian manifold with boundary. The linear space $C^\infty(\mathbb{D}, \mathbb{R}^2)$ of smooth maps $\mathbb{D} \to \mathbb{R}^2$ is a Fréchet space (see [5, Sect. I.1] for details on the Fréchet topology used). The set $\text{Emb}(\mathbb{D}, \mathbb{R}^2)$ of smooth embeddings of the disk into the plane is an open subset of $C^\infty(\mathbb{D}, \mathbb{R}^2)$ and has the structure of a Fréchet manifold (see [5, Sect. I.4.1]). The subspace $C_c^\infty(\mathbb{D}, \mathbb{R}^2) = \{\phi \in C^\infty(\mathbb{D}, \mathbb{R}^2); \phi^*g = Fg, F \in \mathcal{F}(\mathbb{D})\}$ of maps preserving the metric up to a scalar function is topologically closed in $C^\infty(\mathbb{D}, \mathbb{R}^2)$. Furthermore, it holds that the set of conformal embeddings $\text{Con}(\mathbb{D}, \mathbb{R}^2) = C_c^\infty(\mathbb{D}, \mathbb{R}^2) \cap \text{Emb}(\mathbb{D}, \mathbb{R}^2)$ is a Fréchet submanifold of $\text{Emb}(\mathbb{D}, \mathbb{R}^2)$.

For standard planar template matching equations, one introduces a positive-definite quadratic Lagrangian function (corresponding to a weak Riemannian structure) on the infinite dimensional manifold $\text{Diff}(\mathbb{D}) = \text{Emb}(\mathbb{D}, \mathbb{D})$. However, as mentioned above, the set $\text{Diff}_c(\mathbb{D})$ of disk preserving conformal maps is small, so we consider instead $\text{Con}(\mathbb{D}, \mathbb{R}^2)$ as configuration manifold. The Riemann mapping theorem asserts that $\text{Con}(\mathbb{D}, \mathbb{R}^2)$ contains a rich set of maps: for any simply connected domain $U \subset \mathbb{R}^2$ we can find $\varphi \in \text{Con}(\mathbb{D}, \mathbb{R}^2)$ such that $\varphi(\mathbb{D}) = U$, and $\varphi$ is then unique up to the disk-preserving Möbius transformations. That is, $\varphi$ must be an element of a unique co-set $[\varphi] \in \text{Con}(\mathbb{D}, \mathbb{R}^2)/\text{Diff}_c(\mathbb{D})$, where $\text{Diff}_c(\mathbb{D})$ acts on $\text{Con}(\mathbb{D}, \mathbb{R}^2)$ by composition from the right. Thus, it holds that $\text{Con}(\mathbb{D}, \mathbb{R}^2)/\text{Diff}_c(\mathbb{D})$ is equivalent to the well known shape space $\text{Emb}(\mathbb{D}, \mathbb{R}^2)/\text{Diff}(\mathbb{D})$, so we expect that the equation studied in this paper will be relevant not only for conformal image matching, but also for planar shape matching$^4$. Also, as developed in [6], a planar shape may, by conformal mappings, be represented by a fingerprint in $\text{Diff}(S^1)$, which suggests that $\text{Diff}(S^1)$ may be used for co-sets.

Since $\text{Con}(\mathbb{D}, \mathbb{R}^2)$ is open in the linear space $C_c^\infty(\mathbb{D}, \mathbb{R}^2)$ its tangent bundle is $T\text{Con}(\mathbb{D}, \mathbb{R}^2) = \text{Con}(\mathbb{D}, \mathbb{R}^2) \times C_c^\infty(\mathbb{D}, \mathbb{R}^2)$. The Lagrangian we are interested in is given by

$$L(\varphi, \dot{\varphi}) = \frac{1}{2} \int_{\varphi(\mathbb{D})} g(\dot{\varphi} \circ \varphi^{-1}, \dot{\varphi} \circ \varphi^{-1}) \, dA \quad (1)$$

<footnote>$^4$ This is only true in $\mathbb{R}^2$: for higher dimensions the set of conformal embeddings of the unit ball $\mathcal{B}$ is very small, and it does not hold that $\text{Con}(\mathcal{B}, \mathbb{R}^n)/\text{Diff}_c(\mathcal{B})$ and $\text{Emb}(\mathcal{B}, \mathbb{R}^n)/\text{Diff}(\mathcal{B})$ are equivalent.</footnote>
where $dA = dx \wedge dy$ is the standard Euclidean area form on $\mathbb{R}^2$ associated with g.

Note that $\varphi^{-1}$ is well-defined as a map $\varphi(\mathbb{D}) \to \mathbb{D}$ and that the Lagrangian is well-defined on the full embedding manifold $\text{Emb}(\mathbb{D}, \mathbb{R}^2)$. It is quadratic in $\dot{\varphi}$ and positive-definite, thus determining a (weak) Riemannian metric. When restricted to the submanifold $\text{Diff}(\mathbb{D})$ of $\text{Emb}(\mathbb{D}, \mathbb{R}^2)$, this metric coincides with the metric used in the $L^2$ template matching equation (TME). Further, $\text{Diff}_c(\mathbb{D})$ is a right symmetry group: if $\phi \in \text{Diff}_c(\mathbb{D})$ then $L(\varphi \circ \phi, \dot{\varphi} \circ \phi) = L(\varphi, \dot{\varphi})$. In turn, this implies that we obtain a reduced variational principle with respect to the reduced variable $\xi = \dot{\varphi} \circ \varphi^{-1} \in \mathcal{X}_c(\varphi(\mathbb{D}))$, where $\mathcal{X}_c(\varphi(\mathbb{D}))$ denotes the set of conformal vector fields (whose flow is conformal). However, we still need to keep track of the correct domain $\varphi(\mathbb{D})$, which corresponds to the “shape space” element in $\text{Con}(\mathbb{D}, \mathbb{R}^2)/\text{Diff}_c(\mathbb{D})$, so the complete set of reduced variables consists of a domain $U \in \mathcal{R}(\mathbb{R}^2)$ (where $\mathcal{R}(\mathbb{R}^2)$ is the Fréchet manifold of simply connected compact planar regions, see [5, Sect. I.4.3]), and a conformal vector field $\xi \in \mathcal{X}_c(U)$ defined on this domain. Thus, the reduced phase space is the Fréchet tangent bundle $\mathcal{X}(\mathcal{R}(\mathbb{R}^2))$ over $\mathcal{R}(\mathbb{R}^2)$, such that the fibre over $U$ is the vector space $\mathcal{X}_c(U)$ of conformal vector fields (see [7] for details). For $\xi, \eta \in \mathcal{X}(U)$ let $\langle \xi, \eta \rangle_U := \int_U g(\xi, \eta)dA$. The variational principle then becomes

$$0 = \left. \frac{d}{d\epsilon} \right|_{\epsilon=0} \int_0^1 L(\varphi_{\epsilon}, \dot{\varphi}_{\epsilon}) \, dt = \left. \frac{d}{d\epsilon} \right|_{\epsilon=0} \int_0^1 \frac{1}{2} \langle \xi_{\epsilon}, \xi_{\epsilon} \rangle_{U_{\epsilon}} \, dt \quad (2)$$

where $\varphi_{\epsilon}$ is a variation of an extremal curve $t \mapsto \varphi(t)$ in $\text{Con}(\mathbb{D}, \mathbb{R}^2)$ and

$$(U_{\epsilon}, \xi_{\epsilon}) := (\varphi_{\epsilon}(\mathbb{D}), \dot{\varphi}_{\epsilon} \circ \varphi_{\epsilon}^{-1}).$$

Since a general variation of $t \mapsto \varphi(t)$ is of the form $\exp(\varepsilon\eta) \circ \varphi$ (where $t \mapsto \eta(t)$ is a path in $\mathcal{X}_c(\varphi(\mathbb{D}))$ with vanishing endpoints) direct calculation yields

$$\left. \frac{d}{d\epsilon} \right|_{\epsilon=0} \xi_{\epsilon} = \dot{\eta} + \mathcal{L}_{\eta} \xi \\ \left. \frac{1}{2} \frac{d}{d\epsilon} \right|_{\epsilon=0} \langle \xi, \xi \rangle_{U_{\epsilon}} = \langle \mathcal{L}_{\eta} \xi, \xi \rangle_{U} + \int_{U} g(\xi, \xi) \, \text{div}(\eta) \, \text{vol}.$$

where $\mathcal{L}_{\eta}$ denotes the Lie derivative along $\eta$. Plugging this into (2) yields the weak form of a planar conformal template matching equation (PCTME) to be studied in this paper:

$$\int_0^1 \langle \xi, \dot{\eta} + 2\mathcal{L}_{\eta}\xi + \text{div}(\eta)\xi \rangle_U \, dt = 0 \quad (3)$$

for all variations $t \mapsto \eta(t) \in \mathcal{X}(U(t))$, where $t \mapsto U(t)$ fulfills the “compatibility” equation $\dot{U} = \xi(\partial U)$ and $U(0) = \mathbb{D}$, i.e., the domain is transported along the flow generated by $\xi$. In practice, however, we do not solve for the domain variable $U$, but instead we immediately solve the reconstruction equation $\dot{\varphi} = \xi \circ \varphi$, and use the fact that $U = \varphi(\mathbb{D})$.
## 2.1 The Complex Form

By identifying $\mathbb{R}^2$ with $\mathbb{C}$ the space $\mathcal{X}_c(U)$ is identified with the space of holomorphic functions in the complex domain $U \subset \mathbb{C}$. Likewise, $C_c^\infty(\mathbb{D}, \mathbb{R}^2)$ is identified with the space of holomorphic function on $\mathbb{D} \subset \mathbb{C}$, and $\text{Con}(\mathbb{D}, \mathbb{R}^2)$ with the space of injective holomorphic functions on $\mathbb{D}$. We therefore identify $\varphi, \dot{\varphi}$ as holomorphic functions on $\mathbb{D}$, and $\varphi^{-1}, \xi$ as holomorphic functions on $\varphi(\mathbb{D})$. The complex derivative of a holomorphic function is denoted by a prime, e.g., $\xi'$ for $z \mapsto \xi'(z)$.

When working with complex entries, it is useful to express the weak form (3) in terms of the Hermitian inner product $\langle\langle\xi, \eta\rangle\rangle_U = \int_U \xi\bar{\eta} dA$ instead of the real-valued inner product.

Direct calculation yields that $g(\xi, \text{div}(\eta)\xi) + ig(\xi, \text{div}(i\eta)\xi) = 2\xi\bar{\eta}'\xi$. Using this and the fact that $\langle\langle\xi, \eta\rangle\rangle_{\varphi(\mathbb{D})} = \langle\xi, \eta\rangle_{\varphi(\mathbb{D})} + i\langle\xi, i\eta\rangle_{\varphi(\mathbb{D})}$ we obtain the complex weak form of the PCTME (3) as (where we write $\varphi(\mathbb{D})$ instead of $U$ to indicate that we choose to simultaneously solve the reconstruction equation $\dot{\varphi} = \xi \circ \varphi$, instead of only solving for $U)$

$$ \int_0^1 \langle\langle \xi, \dot{\eta} + 4\eta'\xi - 2\xi'\eta \rangle\rangle_{\varphi(\mathbb{D})} dt = 0. \quad (4) $$

We investigate a sheet ansatz for this equation in Section 4, leading to a numerical discretization, but first we consider a set of special solutions to this equation, by identifying a totally geodesic submanifold of $\text{Con}(\mathbb{D}, \mathbb{R}^2)$.

## 3 Totally Geodesic Submanifolds

Recall that a submanifold $N \subset M$ of a Riemannian manifold $(M, g)$ is totally geodesic with respect to $(M, g)$ if geodesics in $N$ (with respect to $g$ restricted to $N$) are also geodesics in $M$. For a thorough treatment of totally geodesic subgroups of $\text{Diff}(M)$ (with respect to various metrics), see [8].

Consider the Fréchet submanifold of linear conformal transformations

$$ \text{Lin}(\mathbb{D}, \mathbb{R}^2) = \{ \varphi \in \text{Con}(\mathbb{D}, \mathbb{R}^2); \varphi(z) = cz, c \in \mathbb{C} \}. $$

**Proposition 1.** $\text{Lin}(\mathbb{D}, \mathbb{C})$ is totally geodesic in $\text{Con}(\mathbb{D}, \mathbb{C})$.

*Proof.* If $t \mapsto \varphi(t)$ is a path in $\text{Lin}(\mathbb{D}, \mathbb{C})$, i.e., $\varphi(z) = cz$ with $c \in \mathbb{C}$, then $\xi = \dot{\varphi} \circ \varphi^{-1}$ is of the form $\xi(z) = az$ with $a \in \mathbb{C}$. Now, let $t \mapsto (\varphi, \xi)$ fulfill the variational equation (4) for each variation of the form $\eta(z) = bz$ with $b \in \mathbb{C}$. We need to show that $t \mapsto (\varphi, \xi)$ then fulfills the equation for any variation of the form $\eta(z) = ez^k$ (since the monomials span the space of holomorphic functions). Thus,

$$ \begin{aligned} \langle\langle\xi, \dot{\eta} + 4\eta'\xi - 2\xi'\eta\rangle\rangle_{\varphi(\mathbb{D})} &= \langle\langle\varphi' \cdot \xi \circ \varphi), \varphi' \cdot (\dot{\eta} + 4\eta'\xi - 2\xi'\eta) \circ \varphi\rangle\rangle_{\mathbb{D}} \\ &= |c|^4 \langle\langle az, \dot{b}z^k + 4kbaz^k - 2abz^k\rangle\rangle_{\mathbb{D}}. \end{aligned} \quad (5) $$

where in the first line we use the conformal change of variables formula for integrals. Now, since the monomials are orthogonal with respect to $\langle\langle\cdot, \cdot\rangle\rangle_{\mathbb{D}}$ the expression vanish whenever $k \neq 1$, which concludes the proof. $\square$
From equation (5) we may derive a differential equation for the totally geodesic solutions in Lin($\mathbb{D}$, $\mathbb{C}$) in term of the variables $(c, a)$ corresponding to $\varphi(z) = cz$ and $\xi(z) = az$. Indeed, choosing $k = 1$ and plugging equation (5) into the variational equation (4) we obtain

$$ 
\begin{aligned}
0 &= \int_0^1 |c|^4 (a\bar{b} + 2|a|^2\bar{b}) \, dt = \int_0^1 \left( \frac{d}{dt} a\bar{b}|c|^4 - \bar{b} \left( \frac{d}{dt} a|c|^4 - 2|a|^2|c|^4 \right) \right) \, dt \\
&= - \int_0^1 \bar{b} (\dot{a}|c|^4 + 2a\dot{c}\bar{c}|c|^2 + 2a\bar{c}c|c|^2 - 2|a|^2|c|^4) \, dt \\
&= - \int_0^1 \bar{b}|c|^4 (\dot{a} + 2a)
\end{aligned}
$$

where we have used that $b$ vanish at the endpoints, and in the last line we use that $\dot{c} = ac$, which follows since $\dot{\varphi} = \xi \circ \varphi$. Hence, the governing equations for the totally geodesic solutions are

$$ 
\begin{aligned}
\dot{c} &= ac \\
\dot{a} &= -2a^2
\end{aligned}
\qquad (6)
$$

where $c(0) = 1$ (since $\varphi(z) = z$ initially). We can see that if $a(0)$ is real, then both $a$ and $c$ stay real, so that the smaller submanifold of pure scalings is also totally geodesic. However, pure rotations are not totally geodesic.

By using again the weak form (4) of the governing equation one can further show that the submanifold of translations is not totally geodesic in Con($\mathbb{D}$, $\mathbb{C}$). Nor is the submanifold of affine conformal transformations. The result is somewhat surprising, given that these types of basic transformations are totally geodesic for the $L^2$–TME and $H_\alpha^1$–TME$^5$, if boundary conditions are set up to allow them [8]. In fact, from a template matching point of view, where basic transformations should preferably be totally geodesic, our result indicates that the metric we have picked is not ideal.

## 4 Sheet Ansatz

The $L^2$–TME and $H_\alpha^1$–TME are known to admit solutions corresponding to totally geodesic submanifolds described by momentum particles or sheets [9]. These special solutions can be used for structure preserving numerical discretization [10,11]. Since the ansatz is non-smooth (the momentum is a traveling peakon), it cannot describe solutions to the conformal equation studied here. However, in this section we show that it is possible to express solutions to the PCTME (4) with a different sheet ansatz, using a reproducing kernel.

Let $\mathcal{O} \subset \mathbb{C}$ be an open bounded complex region. Then $A^2(\mathcal{O}) = \{f \in L^2(\mathcal{O}); f \text{ is holomorphic}\}$ is a Hilbert space with the inner product $\langle\cdot, \cdot\rangle_{\mathcal{O}}$ (called

<footnote>$^5$ Also called the *EPDiff equation* or the *averaged template matching equation*.</footnote>
a *Bergman space*, see [12, Ch. 1]). For any $z \in \mathcal{O}$ it holds that point-wise evaluation $E_z : \xi \mapsto \xi(z)$ is continuous with respect to $\langle\cdot, \cdot\rangle_{\mathcal{O}}$. Thus, by the Riesz representation theorem there exists a unique $k_z \in A^2(\mathcal{O})$ such that $f(z) = \langle\langle f, k_z\rangle\rangle_{\mathcal{O}}$ for every $f \in A^2(\mathcal{O})$. Let $\mathcal{U} = \bar{\mathcal{O}}$ be the closure of $\mathcal{O}$. Since $\mathcal{X}_c(\mathcal{U}) \subset A^2(\mathcal{O})$ it clearly holds that $\xi(z) = \langle\langle \xi, k_z\rangle\rangle_{\mathcal{U}}$ for all $\xi \in \mathcal{X}_c(\mathcal{U})$ (however, notice that it does not hold that $k_z \in \mathcal{X}_c(\mathcal{U})$). The function defined by $K_{\mathcal{U}}(z, w) = \bar{k}_z(w)$ is called the *reproducing kernel* (or sometimes the *Bergman kernel function*) for the domain $\mathcal{U}$. On the interior of $\mathcal{U}$ this function is analytic in $z$ and anti-analytic in $w$, with $K(z, w) = \bar{K}(w, z)$. Furthermore, if $\varphi \in \text{Con}(\mathcal{D}, \mathbb{C})$, then

$$ K_{\varphi(\mathcal{D})}(z, w) = \frac{K_{\mathcal{D}}(\varphi^{-1}(z), \varphi^{-1}(w))}{\varphi'(\varphi^{-1}(z))\varphi'(\varphi^{-1}(w))} \quad (7) $$

which follows by a change of variables in the integral over $\varphi(\mathcal{D})$. When working with kernels in practice, e.g., for numerical purposes, this formula is very useful since the kernel function for the disk is known to be

$$ K_{\mathcal{D}}(z, w) = \frac{1}{\pi (1 - z \bar{w})^2} . $$

Now, consider the following ansatz:

$$ \xi(z) = \int_{0}^{2\pi} p_{\varphi}(\rho e^{is}) K_{\varphi(\mathcal{D})}(z, \varphi(\rho e^{is})) ds \quad (8) $$

where $p_{\varphi} \in C_{\mathcal{D}}^{\infty}(\mathcal{D}, \mathbb{C})$ and $0 < \rho < 1$ is a fixed constant. We introduce the periodic functions $\gamma, p_{\gamma} : S^1 \to \mathbb{C}$ defined by $\gamma(s) = \varphi(\rho e^{is})$ and $p_{\gamma}(s) = p_{\varphi}(\rho e^{is})$, so that $\xi(z) = \int_{0}^{2\pi} p_{\gamma} K_{\varphi(\mathcal{D})}(z, \gamma) ds$. One may think of $\gamma(s)$ as just a notational shortcut, but it has other significance: by expanding $\varphi$ as a Taylor series with coefficients $(c_k)_{k=0}^{\infty}$, we see that the periodic curve $\gamma(s)$ only has positive Fourier coefficients, which are given by $(\rho^k c_k)_{k=0}^{\infty}$. An equivalent statement is to say that the Fourier coefficients of $\gamma(s)$ are the Taylor coefficients of $z \mapsto \varphi(\rho z)$. Exactly the same relation holds between $p_{\gamma}$ and $p_{\varphi}$. This observation can be used in the numerical discretization to accelerate the evaluation of the right-hand side using fast Fourier transformations.

Plugging the ansatz (8) into the weak equation (4) yields

$$ \begin{aligned} 0 &= \int_0^1 \int_0^{2\pi} \overline{p_\gamma} \langle\langle \dot{\eta} + 4\eta'\xi - 2\xi'\eta, K_{\varphi(\mathcal{D})}(\cdot, \gamma) \rangle\rangle_{\varphi(\mathcal{D})} \, ds \, dt \\ &= \int_0^1 \int_0^{2\pi} \overline{p_\gamma} (\dot{\eta}(\gamma) + 4\eta'(\gamma)\xi(\gamma) - 2\xi'(\gamma)\eta(\gamma)) \, ds \, dt \\ &= \int_0^1 \int_0^{2\pi} \left( \eta(\gamma) \left(-\overline{p_\gamma} - 2\overline{p_\gamma}\xi'(\gamma)\right) + \overline{p_\gamma}\eta'(\gamma) \left(4\xi(\gamma) - \dot{\eta}\right) \right) \, ds \, dt \end{aligned} $$

where we use the evaluation property of the kernel in the second line and the fact that $t \mapsto \eta(t)$ vanish at the endpoints in the third. Next, using that $\dot{\eta}(s) =$
$\dot{\varphi}(\rho e^{is}) = \xi(\varphi(\rho e^{is})) = \xi(\gamma(s))$, we get

$$ 
\begin{aligned}
0 &= \int_0^1 \int_0^{2\pi} \left( \eta(\gamma) \left( -\overline{\dot{p}_\gamma} - 2\overline{p_\gamma} \xi'(\gamma) \right) + \overline{p_\gamma} \eta'(\gamma) 3\xi(\gamma) \right) \mathrm{d}\sigma \mathrm{d}t \\
&= \int_0^1 \int_0^{2\pi} \left( \eta(\gamma) \left( -\overline{\dot{p}_\gamma} - 2\overline{p_\gamma} \xi'(\gamma) \right) + \eta'(\gamma) \gamma' \left( \frac{3\overline{p_\gamma} \xi(\gamma)}{\gamma'} \right) \right) \mathrm{d}\sigma \mathrm{d}t \\
&= \int_0^1 \int_0^{2\pi} \left( \eta(\gamma) \left( -\overline{\dot{p}_\gamma} - 2\overline{p_\gamma} \xi'(\gamma) \right) - \eta(\gamma) \frac{\partial}{\partial s} \left( \frac{3\overline{p_\gamma} \xi(\gamma)}{\gamma'} \right) \right) \mathrm{d}\sigma \mathrm{d}t \\
&= - \int_0^1 \int_0^{2\pi} \eta(\gamma) \left( \overline{\dot{p}_\gamma} - 2\xi(\gamma) \frac{\partial}{\partial s} \left( \frac{\overline{p_\gamma}}{\gamma'} \right) + 5 \frac{\partial}{\partial s} \left( \frac{\overline{p_\gamma} \xi(\gamma)}{\gamma'} \right) \right) \mathrm{d}\sigma \mathrm{d}t
\end{aligned}
\quad (9) $$

where $\gamma'$ denotes the derivative of $\gamma$ with respect to $s$. Thus, by inserting the ansatz (8) we obtain a weak formulation of the PCTME, equivalent to (4), in terms of the variables $\gamma$ and $p_\gamma$. This form involves an integral operator, since $\xi$ is related to $\gamma$, $p_\gamma$ by the ansatz (8).

In contrast to the $L^2$-TME and $H_\alpha^1$-TME, any solution to the PCTME can be represented by the sheet ansatz, since any $\varphi$ can be reconstructed from $\gamma$ by the Fourier transformation.

A nice feature of the weak formulation (9) is that the test function $\eta$ is isolated, i.e., without derivatives. This allows us to use the fundamental lemma of calculus of variations to get rid of the outermost time integral (so we get a condition at each fixed point in time). However, we do not immediately obtain a strong formulation, because the test function $\eta$ must be holomorphic, so we cannot use the fundamental lemma of calculus of variations for the inner integral without introducing a projection operator. Since the trigonometric monomials $e^{iks}$ are orthogonal with respect to the $s$-integral, the projection is given by neglecting negative Fourier coefficients. Hence, let $\chi_+ : \mathbb{Z} \to \mathbb{R}$ be the sequence

$$ \chi_+(k) = \begin{cases} 1 & \text{if } k \ge 0 \\ 0 & \text{if } k < 0 \end{cases} $$

and let $\mathcal{F}$ denote the Fourier operator taking a $2\pi$-periodic function to its Fourier series. Then the projection operator is given by $\mathcal{P} = \mathcal{F}^{-1}\chi_+\mathcal{F}$ (where $\chi_+$ is acting by element-wise multiplication). We now obtain a strong integral formulation of the PCTME in terms of $\gamma$ and $p_\gamma$ as

$$ 
\begin{aligned}
\dot{\gamma} &= \xi(\gamma) \\
\dot{p}_\gamma &= \mathcal{P} \left( 2\overline{\xi(\gamma)} \frac{\partial}{\partial s} \frac{p_\gamma}{\gamma'} - 5 \frac{\partial}{\partial s} \frac{p_\gamma \overline{\xi(\gamma)}}{\gamma'} \right) .
\end{aligned}
\quad (10) $$

# 5 Numerical Discretization

In this section we describe a method for numerical discretization of equations (10). The idea is to represent the dynamic variables $\gamma(s)$ and $p_\gamma(s)$ by truncated positive Fourier series, or, equivalently, truncated Taylor series of $z \mapsto \varphi(\rho z)$ and
$z \mapsto p_{\varphi}(\rho z)$. Thus, our finite set of dynamic variables are $\mathbf{c} = (c_k)_{k=0}^{n-1} \in \mathbb{C}^n$ and $\mathbf{p}_c = (p_k)_{k=0}^{n-1} \in \mathbb{C}^n$ such that

$$ \gamma(s) = \varphi(\rho e^{is}) = \sum_{k=0}^{n-1} c_k e^{is}, \quad p_{\gamma}(s) = p_{\varphi}(\rho e^{is}) = \sum_{k=0}^{n-1} p_k e^{is}. \quad (11) $$

For the discretization we insert these into the right-hand side of equation (10) and compute the first $n$ positive Fourier coefficients, which then gives an ordinary differential equation for $(\mathbf{c}, \mathbf{p}_c) \in \mathbb{C}^{2n}$. We now describe how this discretized right-hand side is computed efficiently.

1. Compute $\gamma = (\gamma(0), \gamma(\pi/n), \dots, \gamma(\pi(2n-1)/n))$ from $\mathbf{c}$. (Notice that $\gamma$ has length $2n$.) Since $\gamma(s) = \varphi(\rho e^{is})$ and $\mathbf{c}$ are the Taylor coefficients of $z \mapsto \varphi(\rho z)$, we obtain $\gamma$ from the inverse FFT as $\gamma = 2n \text{IFFT}(\mathbf{c}^+)$, where $\mathbf{c}^+ = (\mathbf{c}, 0, \dots, 0)$ is zero-padded to get the same length as $\gamma$ (since in the FFT the last half of the vector corresponds to negative frequencies).

2. Compute $\mathbf{p}_{\gamma} = (p_{\gamma}(0), \dots, p_{\gamma}(\pi(2n-1)/n))$ from $\mathbf{p}_c$. By the same argument we have $\mathbf{p}_{\gamma} = 2n \text{IFFT}(\mathbf{p}_c)$.

3. Compute $\gamma' = (\gamma'(0), \gamma'(2\pi/n), \dots, \gamma'(2\pi(n-1)/n))$. Since the Fourier coefficients of $\gamma(s)$ are given by $\mathbf{c}$, we have $\gamma' = 2n \text{IFFT}(i\mathbf{k}^+ \cdot \mathbf{c}^+)$ where $\mathbf{k} = (0, 1, \dots, n-1)$ and $\mathbf{k}^+ \cdot \mathbf{c}^+$ denotes element-wise multiplication.

4. Compute $\xi = (\xi(\gamma_0), \dots, \xi(\gamma_{2n-1}))$. From the ansatz (8) and the equation for the Bergman kernel (7), we get

$$ \xi(\gamma(s)) = \int_{0}^{2\pi} \frac{p_{\gamma}(\sigma) K_{\mathbb{D}}(\rho e^{is}, \rho e^{i\sigma})}{\varphi'(\rho e^{is}) \overline{\varphi'(\rho e^{i\sigma})}} d\sigma = \frac{1}{\pi \gamma'(s)} \int_{0}^{2\pi} \frac{p_{\gamma}(\sigma) G(\rho^2 e^{i(s-\sigma)})}{\overline{\gamma'(\sigma)}} d\sigma $$

where $G(z) = z/(1-z)^2$. Now, set $f(s) = p_{\gamma}(s)/\overline{\gamma'(s)}$ and $g(s) = G(\rho^2 e^{is})$. Then

$$ \xi(\gamma(s)) = \frac{1}{\pi \gamma'(s)} \int_{0}^{2\pi} f(\sigma) g(s - \sigma) d\sigma = \frac{(f * g)(s)}{\pi \gamma'(s)}. $$

Thus, using that a convolution becomes element-wise multiplication in the Fourier domain, we get $\xi = \frac{2n}{\pi} \text{IFFT}(\hat{f} \cdot \hat{g})/\gamma'$ where $\hat{g} = (0, \rho^2, \dots, \rho^{2(n-1)}(n-1), 0, \dots, 0)$, which follows by computing the Taylor coefficients of $G(z)$, and $\hat{f} = \text{FFT}(\mathbf{p}_{\gamma}/\overline{\gamma'})/2n$ (divisions are carried out element-wise).

5. Compute $\dot{c}$. Using again the correspondence between Taylor and Fourier coefficients, we get $\dot{c} = \text{FFT}(\xi^{-1}/2n)$, where the FFT$(\xi^{-1})$ means that we only keep the first $n$ elements. Thus, we have now computed the first half of the right-hand side.

6. Compute $\dot{p}_c$. From the second equation in (10) we get

$$ \dot{p}_c = \text{FFT}(2\bar{\xi} \cdot \text{IFFT}(i\mathbf{k}^+ \cdot \hat{f}))^{-1} - 5i\mathbf{k} \cdot \text{FFT}(\bar{\xi} \cdot \mathbf{p}_{\gamma}/\overline{\gamma'})^{-1}/2n. $$

Note that the projection (corresponding to the operator $\mathbb{P}$ in equation (10)) occurs in the above computation since we only keep elements corresponding to positive frequencies.
We have now computed the full right-hand side $(\dot{c}, \dot{p}_c)$. By using a time-
stepping method for ODEs (e.g., a Runge-Kutta method), a numerical method
for equation (10) is obtained. Each evaluation of the right-hand side requires
element-wise operations and 9 FFTs on vectors of length 2n, which leads to
complexity $O(n \log(2n))$.

**6 Experimental Results**

In this section we study the dynamical behaviour of some solutions to the PCTME. First, we look at the solutions corresponding to the totally geodesic submanifold derived in Section 3. Thereafter, we use the numerical method derived in Section 5 to study other solutions. In particular, we study the spectrum of small perturbations of a totally geodesic solution.

**6.1 Totally Geodesic Solutions**

Let $\varphi_0(z) = z$ be the identity maps, and let $\varphi_1(z) = c_1z$ where $c_1 \in \mathbb{C}\setminus\{0\}$.
Since both $\varphi_0$ and $\varphi_1$ belong to the totally geodesic submanifold $\text{Lin}(\mathbb{D}, \mathbb{C})$, the
geodesic from $\varphi_0$ to $\varphi_1$ stays in $\text{Lin}(\mathbb{D}, \mathbb{C})$. In Fig. 1 the geodesics for $c_1 =$
0.2 (scaling of the disk) and $c_1 = e^{1.6\pi i}$ (rotation of the disk) are shown. The
solutions confirm what we earlier noticed, that pure scalings are totally geodesic
(the curve $\varphi(t)$ remains a scaling for each $t$) whereas the geodesic corresponding
to the rotated disk does not stay a pure rotation, but contains scaling along its
path.

<figure><img src="image_25.png" /><figcaption>Fig. 1. Geodesic curve from $\varphi_0(z) = z$ to $\varphi_1(z) = c_1z$ for different values of $c_0$. The mesh lines show how the unit circle evolves. Notice that the scaling geodesic stays a scaling (left figure), whereas the rotation geodesic picks up some scaling during its time evolution (right figure).</figcaption></figure>
## 6.2 Development of Cusps in Finite Time

The $L^2$-TME is known to develop singularities in finite time. Originally, we expected solutions to the PCTME be more “well behaved” due to the holomorphic constraint. However, experiments with the numerical method derived in Section 5 indicates that the same phenomenon still occurs despite of the constraint. For example, Fig. 2 shows an experiment with initial data of the form $\dot{\varphi}_0(z) = a_2z^2$. A cusp develops in the solution, which leads to a breakdown of the dynamics.

<figure><img src="image_26.png" /><figcaption>**Fig. 2.** Time evolution of the PCTME with initial conditions of the form $\varphi_0(z) = z$ and $\dot{\varphi}_0(z) = a_2z^2$. From left to right the plots show the development of a cusp. Thus, the geodesic curve ceases to exist in finite time.</figcaption></figure>

## 6.3 Spectral Behaviour

In Fig. 3 the time evolution of the absolute values of $p_c$, computed with the numerical method in Section 5, are shown in two plots. Dark colours correspond to low order coefficients, and light colours to high order coefficients. In the top plot, initial data corresponding to a totally geodesic solution in $\text{Lin}(\mathbb{D}, \mathbb{C})$ are used. In the bottom plot, this initial data is perturbed $\dot{\varphi}_0(z) = a_1z + \delta(z)$, where $\delta(z)$ is a small perturbation of the first five Taylor coefficients.

Notice in both plots that the initially suppressed coefficients grow exponentially fast, which indicates that the totally geodesic solutions are not stable. Also, notice that the growth rate increases with the order, which indicates that the PCTME is ill-conditioned, as is also the case for the $L^2$-TME. Again, our investigation indicates that the dynamics of the PCTME is similar to that of the $L^2$-TME. We plan to carry out a more thorough investigation of the spectral behaviour, by using both analytical and numerical techniques.

## 7 Conclusions

In this paper we have studied a geodesic equation on the manifold of planar conformal embeddings. We showed that the equations have a totally geodesic submanifold corresponding to linear conformal maps. We also showed that the submanifold of affine conformal maps is not totally geodesic. Numerical experiments indicates that the dynamic behaviour of the equation is similar to that of the $L^2$-TME.
<figure><img src="image_27.png" /><figcaption>Fig. 3. Time evolution of the absolute values of $p_c$ for a totally geodesic solution computed with the numerical method. (Top): without perturbation, (bottom): a random perturbation of size $\sim 10^{-4}$ of initial data for the first five coefficients.</figcaption></figure>

In future work we will study the spectral behaviour more thoroughly. Also, as an approach for obtaining more well behaved dynamics, we will consider other metrics, in particular $H_\alpha^1$. We will also look into more advanced numerical techniques for solving the equations.

References

1. Younes, L.: Shapes and Diffeomorphisms. Applied Mathematical Sciences. Springer-Verlag (2010)

2. Thompson, D.: On Growth and Form. Cambridge University Press (1942)

3. Wallace, A.: D'Arcy Thompson and the theory of transformations. Nature Reviews Genetics **7** (2006) 401–406

4. Camassa, R., Holm, D.D.: An integrable shallow water equation with peaked solitons. Phys. Rev. Lett. **71** (1993) 1661–1664

5. Hamilton, R.S.: The inverse function theorem of Nash and Moser. Bull. Amer. Math. Soc. (N.S.) **7** (1982) 65–222

6. Sharon, E., Mumford, D.: 2d-shape analysis using conformal mapping. Int. J. Comput. Vis. **70** (2006) 55–75 10.1007/s11263-006-6121-z.

7. Marsland, S., McLachlan, R.I., Modin, K., Perlmutter, M.: Reduced geodesic equation for planar conformal embeddings. In preparation (2011)

8. Modin, K., Perlmutter, M., Marsland, S., McLachlan, R.I.: On Euler-Arnold equations and totally geodesic subgroups. J. Geom. Phys. **61** (2011) 1446–1461

9. Holm, D.D., Marsden, J.E.: Momentum maps and measure-valued solutions (peakons, filaments, and sheets) for the EPDiff equation. In: The breadth of symplectic and Poisson geometry. Volume 232 of Progr. Math. Birkhäuser Boston, Boston, MA (2005) 203–235

10. McLachlan, R.I., Marsland, S.: N-particle dynamics of the Euler equations for planar diffeomorphisms. Dyn. Syst. **22** (2007) 269–290

11. Chertock, A., Toit, P., Marsden, J.E.: Integration of the EPDiff equation by particle methods. (2009)

12. Duren, P., Schuster, A.: Bergman Spaces. Volume 100 of Mathematical Surveys and Monographs. American Mathematical Society, Providence, RI (2004)
# Geodesics, Parallel Transport & One-parameter Subgroups for Diffeomorphic Image Registration

Marco Lorenzi¹,² and Xavier Pennec¹

¹ Project Team Asclepios, INRIA Sophia Antipolis, France

² LENITEM, IRCCS San Giovanni di Dio, Fatebenefratelli, Brescia, Italy

**Abstract.** The aim of computational anatomy is to develop models for understanding the physiology of organs and tissues. The diffeomorphic non-rigid registration is a validated instrument for the detection of anatomical changes on medical images and is based on a rich mathematical background. For instance, the “large deformation diffeomorphic metric mapping” framework defines a Riemannian setting by providing an opportune right invariant metric on the tangent space, and solves the registration problem by computing geodesics parametrized by time-varying velocity fields. In alternative, stationary velocity fields have been proposed for the diffeomorphic registration based on the one-parameter subgroups from Lie groups theory. In spite of the higher computational efficiency, the geometric setting of the latter method is more vague, especially regarding the relationship between one-parameter subgroups and geodesics. In this study, we present the relevant properties of the Lie groups for the definition of geometrical properties within the one-parameter subgroups parametrization, and we define the geometric structure for computing geodesics and for parallel transporting. The theoretical results are applied to the image registration context, and discussed in light of the practical computational problems.

## 1 Introduction

Main objective of the computational anatomy is to develop suitable statistical models on several subjects for understanding the physiology of organs and tissues. In particular, the longitudinal observations from time series of images are an important source of information for understanding the developmental processes and the dynamics of pathologies. Thus, a reliable method for comparing different longitudinal trajectories is required, in order to develop population-based longitudinal models.

Non-rigid registration is a validated instrument for the detection of anatomical changes on medical images, and it has been widely applied on different clinical contexts for the definition of population-based anatomical atlases ([14],[9],[3]). However, in case of longitudinal data, the optimal method for comparing deformation trajectories across different subjects is still under discussion. In fact, the methods for integrating the subtle inter-subject changes into the group-wise analysis have an important impact on the accuracy and reliability of the subsequent
results. The aim is to preserve as much as possible the biological informations
carried on by the different subjects, while allowing a precise comparison in a
common geometric space.

Among the different techniques proposed for the comparison of longitudinal trajectories ([12],[2],[4]), *parallel transport* represents a promising method which relies on a solid mathematical background. Basically, it consists in transporting the infinitesimal deformation vector across different points by preserving its properties with respect to the space geometry, such as the parallelism.

The parallel transport has been introduced for the first time in medical imaging with the LDDMM setting [16]. LDDMM solves the image registration problem by using a Riemannian framework in which the deformations are parametrized as *diffeomorphisms* living in a suitable space, once provided an opportune right-invariant metric [10]. The registration problem is solved by computing the diffeomorphisms lying on the *geodesics* of the space parametrized by *time-varying* velocity fields under the Riemannian exponential. The setting allows the computation of the parallel transport along geodesics at the cost of a computationally intensive scheme, and this limitation prevents the application on high resolution images or large datasets.

A more efficient solution to the image diffeomorphic registration problem was introduced by the stationary velocity field (SVF) setting [1]. In this case, the diffeomorphisms are parametrized as one-parameter subgroups by *stationary* velocity fields through the Lie group exponential. This restriction allows an efficient numerical scheme for the computation but it does not directly rely on any geometric assumption on the underlying space. This implies that some important mathematical properties are not guaranteed, for instance whether the one-parameter subgroups are still geodesics or if the space is metrically complete. In spite of this lack of knowledge, the framework was found very efficient and reliable in many applications in different contexts ([8],[7],[13]) and, in [6], a framework based on the Schild's Ladder has been proposed for the evaluation of the parallel transport with the SVF.

In this paper, we investigate the relationship between Lie groups and Riemannian geometry and we highlight many interesting properties that might provide the SVF setting with part of the geometrical solidity required. In Section 2 we present the relevant properties of the Lie groups and the relationship with the Riemannian setting for the definition of the geodesics and the parallel transport. In Section 3, the results are introduced and discussed for the image registration context, while in Section 4 we show how these theoretical insights provide a clean, precise, and numerically efficient solution for the parallel transport of deformation trajectories on time series of images.
## 2 Lie Group and covariant differentiation

This section will recall the conceptual basis for the definition of the parallel transport along the one-parameter subgroups.

Let $\mathbb{G}$ an arbitrary finite dimensional Lie group and let $\mathfrak{g}$ the associated Lie algebra defined here with the tangent space at the identity $T_{id}\mathbb{G}$. We define the left translation $L_a$ as the mapping $La : g \mapsto ag$, and we say that a vector field $\mathbf{X} \in T(\mathbb{G})$ is left invariant if $DL_a(\mathbf{X})_b = (\mathbf{X})_{ab}$.

There is a one-to-one correspondence between left-invariant vector fields and elements of the Lie algebra $\mathfrak{g}$, which associates to each $\mathbf{X} \in \mathfrak{g}$ the vector field defined as $\tilde{\mathbf{X}}(g) = DL_g\mathbf{X}$. The left-invariant vector fields are complete and their associated flow $\varphi_t$ is such that $\varphi_t(g) = g\varphi_t(id)$. The association $\mathbf{X} \mapsto \varphi_1(id)$ of $\mathfrak{g}$ into $\mathbb{G}$ is called Lie group exponential and denoted by $\exp$. In particular, the map $\exp$ defines the one-parameter subgroup associated to the vector $\mathbf{X}$ and has the following properties:

$$ - \varphi_t(id) = \exp(t\mathbf{X}), \text{ for each } t \in \mathbb{R} $$

$$ - \exp((t+s)\mathbf{X}) = \exp(t\mathbf{X})\exp(s\mathbf{X}), \text{ for each } t, s \in \mathbb{R} $$

It can be shown that the Lie group exponential is a diffeomorphism from a neighborhood of 0 in $\mathfrak{g}$ to a neighborhood of $id$ in $\mathbb{G}$.

We are going to illustrate the transport of vectors along the exponential path $\exp(t\mathbf{X})$, and in particular the analogies with the classical Riemannian parallel transport defined for geodesics.

An *affine connection* on $\mathbb{G}$ is an operator which assigns to each $X \in T(\mathbb{G})$ a linear mapping $\nabla_X : T(\mathbb{G}) \to T(\mathbb{G})$ such that

$$ \nabla_f \mathbf{X} + g\mathbf{Y} = f\nabla_X + g\nabla_Y \qquad (1) $$

$$ \nabla_X(f\mathbf{Y}) = f\nabla_X(\mathbf{Y}) + (\mathbf{X}f)\mathbf{Y} \qquad (2) $$

A vector field $\mathbf{X}$ is *parallel transported* along a curve $\gamma(t)$ if $\nabla_{\gamma(t)}\mathbf{X} = 0$ for each $t$. In particular, a path $\gamma(t)$ on $\mathbb{G}$ is then said *geodesic* if $\nabla_{\gamma(t)}\gamma = 0$. The definition generalizes the concept of “straight lines”, by requiring to the tangent vector of the path to be covariantly constant.

Given a point $p \in \mathbb{G}$ and a vector $\mathbf{X} \in T_p\mathbb{G}$, there exist a unique geodesic $\gamma(t, p, \mathbf{X})$ such that at the instant $t = 0$ passes through $p$ with velocity $\mathbf{X}$. We define therefore the *Riemannian exponential* as the application $\exp : \mathbb{G} \times T(\mathbb{G}) \to \mathbb{G}$ given by $\exp_p(\mathbf{X}) = \gamma(1, p, \mathbf{X})$.

If, as in the euclidean case, we want to associate to the straight lines the property of minimizing the distance between points, we need to provide the group $\mathbb{G}$ with a Riemannian manifold structure, i.e. with a metric operator $g$ on the tangent space. In this case there is a unique symmetric connection compatible with
the metric in the sense that, for each $X, Y, Z \in T(\mathcal{G})$ the following conditions hold:

$$ \begin{aligned} \mathbf{X} g(\mathbf{Y}, \mathbf{Z}) &= g(\nabla_{\mathbf{X}} \mathbf{Y}, \mathbf{Z}) + g(\mathbf{X}, \nabla_{\mathbf{X}} \mathbf{Z}) && (\text{Compatibility wrt the metric}) \\ \nabla_{\mathbf{X}} \mathbf{Y} - \nabla_{\mathbf{Y}} \mathbf{X} &= [\mathbf{X}, \mathbf{Y}] && (\text{Torsion free}) \end{aligned} $$

With the choice of this special connection, called *Levi-Civita connection*, the corresponding geodesics (Riemannian geodesics) $\gamma(t)$ are the length minimizing path.

## 2.1 Relationship between Riemannian geodesic and one-parameter subgroups

Given a vector $X$ on $T_{id}\mathcal{G}$, we can therefore define two curves on $\mathcal{G}$ passing through $id$ and having $X$ as tangent vector, one given by the Lie group exponential $\exp$ and the other given by the Riemannian exponential $\exp_{id}$. When do they coincide?

The connection $\nabla$ on $\mathcal{G}$ is left-invariant if, for each left translation $L_a$ ($a \in \mathcal{G}$), we have $\nabla_{DL_a} X(DL_a Y) = DL_a \nabla_X Y$.

A left-invariant connection $\nabla$ on a space $\mathcal{G}$ is a *Cartan connection* if, for any element of the Lie algebra $X \in \mathfrak{g}$, the one-parameter subgroups and the Riemannian geodesics coincide, i.e. $\exp(tX) = \exp(t, id, X)$ [11]. For each left-invariant connection $\nabla$ we can univoquely associate a product $\alpha$ (symmetric bilinear operator) on $T_{id}\mathcal{G}$ given by

$$ \alpha(\mathbf{X}, \mathbf{Y}) = (\nabla_{\tilde{\mathbf{X}}} \tilde{\mathbf{Y}})_{id} $$

where $\tilde{\mathbf{X}}, \tilde{\mathbf{Y}}$ are the unique left-invariant vector fields induced by the tangent vectors $\mathbf{X}, \mathbf{Y}$. We note here that a bilinear for can be uniquely decomposed as $\alpha = \alpha' + \alpha''$, where $\alpha' = \frac{1}{2}(\alpha(X, Y) + \alpha(Y, X))$ is commutative, while $\alpha'' = \frac{1}{2}(\alpha(X, Y) - \alpha(Y, X))$ is skew-symmetric.

We deduce that the condition for $\nabla$ to be a Cartan connection is to satisfy $\alpha(\mathbf{X}, \mathbf{X}) = 0$ or, equivalently, to be skew-symmetric, for instance by assigning

$$ \alpha(\mathbf{X}, \mathbf{Y}) = \lambda[\mathbf{X}, \mathbf{Y}] \qquad (3) $$

In this case, the zero curvature connections are given by $\lambda = 0, 1$ (with torsion $T = -[\mathbf{X}, \mathbf{Y}]$ and $T = [\mathbf{X}, \mathbf{Y}]$ respectively) and are called *left* and *right* Cartan connections.

The choice of $\lambda = \frac{1}{2}$ lead to the *symmetric* (or *mean*) Cartan connection $\nabla_X Y = \frac{1}{2}[\mathbf{X}, \mathbf{Y}]$, with curvature $C = -\frac{1}{4}[[\mathbf{X}, \mathbf{Y}], \mathbf{Z}]$ and torsion-free. This connection is the average between left and right Cartan connection. Therefore the Cartan