Indirect Image Registration with Large Diffeomorphic Deformations∗

arXiv:1706.04048v3 [math.NA] 11 Oct 2017

Chong Chen† and Ozan Öktem‡

Abstract. The paper adapts the large deformation diffeomorphic metric mapping framework for image registration to the indirect setting where a template is registered against a target that is given through
indirect noisy observations. The registration uses diffeomorphisms that transform the template
through a (group) action. These diffeomorphisms are generated by solving a flow equation that is
defined by a velocity field with certain regularity. The theoretical analysis includes a proof that
indirect image registration has solutions (existence) that are stable and that converge as the data
error tends so zero, so it becomes a well-defined regularization method. The paper concludes with
examples of indirect image registration in 2D tomography with very sparse and/or highly noisy data.
Key words. indirect image registration, shape theory, large deformations, diffeomorphisms, shape regularization, image reconstruction, tomography
AMS subject classifications. 65F22, 65R32, 65R30, 65D18, 94A12, 94A08, 92C55, 54C56, 57N25, 47A52

1. Introduction. Image registration (matching) is the task of transforming a template
image so that it matches a target image. This arises in many fields, such as quality control
in industrial manufacturing [27], various applications in remote sensing [25], face recognition
[67, 51], robotic navigation [14], and medical imaging [57, 53, 60], etc. The variant that is considered here is indirect image registration, i.e., when the template image is registered against a
target that is known only through indirect noisy observations, such as in tomographic imaging.
It makes use of the large deformation diffeomorphic metric mapping (LDDMM) framework for
diffeomorphic image registration and thereby extends the indirect image registration scheme
in [46], which uses (small) linearized deformations.
An important development in regularization theory is ongoing where reconstruction and
feature extraction steps are pursued simultaneously, such as joint segmentation and reconstruction. Indirect image registration can be seen as part of this development where the
feature in question is the “shape” of the structures in the image and their temporal variability. This is highly relevant for spatiotemporal imaging, see [46] for a survey and other use
cases of indirect image registration.
There is an extensive literature on image registration where a template is registered by
means of a diffeomorphism against a target so that their “shapes” agree, see [57] for a nice survey. Image registration is then recast as the problem of finding a suitable diffeomorphism that
deforms the template into the target image [21]. The underlying assumption is that the target
∗

Version 1.9 revised 2017-10-09.
Funding: The work by Chong Chen and Ozan Öktem has been supported by the Swedish Foundation for
Strategic Research grant AM13-0049. Chen was also supported in part by the National Natural Science Foundation
of China (NSFC) under the grant 11301520.
†
Department of Mathematics, KTH–Royal Institute of Technology, 100 44 Stockholm, Sweden; LSEC, ICMSEC,
Academy of Mathematics and Systems Science, Chinese Academy of Sciences, Beijing 100190, People’s Republic of
China (chench@lsec.cc.ac.cn).
‡
Department of Mathematics, KTH–Royal Institute of Technology, 100 44 Stockholm, Sweden (ozan@kth.se).
1

2

C. CHEN AND O. ÖKTEM

image is contained in the orbit of the template under the group action of diffeomorphisms.
This principle can be stated in a very general setting where diffeomorphism acts on various
image features, like landmark points, curves, surfaces, scalar images, or even vector/tensor
valued images [65]. It also adapts readily to the indirect image registration setting and this
is worked out in detail for the case when one seeks to register scalar images.
2. Overview of paper and specific contributions. The main contribution is to develop
a variational framework for indirect image registration (section 6) where a template is registered against a target that is known only through indirect noisy observations (data). This is
done by appropriately adapting the LDDMM framework where the template is deformed by
diffeomorphisms that act on images through a group action as explained in section 5. This
section also contains an overview of the LDDMM theory.
An important theoretical topic is to investigate to what extent this indirect registration
scheme is a regularization. Existence of solutions along with stability and convergence results are stated in section 7. In conclusion, the LDDMM based indirect registration scheme
is formally a well-defined regularization method. On the other hand, as most of the other
variational schemes for image registration, one cannot expect to have uniqueness due to lack
of convexity. This theoretical investigation is followed by explicit calculations of the derivative
and gradient of the objective functional in the variational formulation of the indirect registration problem (section 8). Numerical implementation, which relies on theory of reproducing
kernel Hilbert space (section 4) is outlined in section 9.
Section 10 contains some numerical experiments from tomography that show performance
of indirect image registration. The results support the claim that shape information contained
in the template has a strong regularizing effect, and especially so for noisy highly undersampled inverse problems. Furthermore, the experiments also suggest that the prior shape
information does not have to be that accurate, which is important for the used cases mentioned
in [46]. None of these claims are however proved formally by mathematical theorems.
Finally, section 11 discusses possible extensions of indirect image registration to the case
when the template also needs to be recovered. It also outlines how indirect image registration
relates to joint reconstruction and motion estimation in spatiotemporal imaging.
3. Inverse problems, ill-posedness and variational regularization. The purpose of this
section is to set some notations and concepts used throughout the paper.
Image reconstruction. The goal in image reconstruction is to estimate some spatially distributed quantity (image) from indirect noisy observations (measured data). Stated mathematically, the aim is to reconstruct an image ftrue ∈ X from data g ∈ Y where
(1)

g = T (ftrue ) + ∆g.

Here, X (reconstruction space) is the vector space of all possible images, so it is a suitable
Hilbert space of functions defined on a fixed domain Ω ⊂ Rn . Next, Y (data space) is the
vector space of all possible data, which for digitized data is a subset of Rm . Furthermore,
T : X → Y (forward operator ) models how a given image gives rise to data in absence of
noise and measurement errors, and ∆g ∈ Y (data noise component) is a sample of a Y –valued
random element whose probability distribution is (data noise model ) assumed to be explicitly
expressible in terms of T (ftrue ).

INDIRECT IMAGE REGISTRATION WITH LARGE DIFFEOMORPHIC DEFORMATIONS

3

Ill-posedness. A naive approach for reconstructing the true (unknown) image ftrue is to
solve the equation T (f ) = g. Often there are no solutions to this equation since measured
data is not in the range of T for a variety of reasons (noise, modeling errors, etc.). This is
commonly addressed by relaxing the notion of a solution by considering

min D T (f ), g .

(2)

f ∈X

The mapping D : Y × Y → R (data discrepancy functional ) quantifies the data misfit and a
natural candidate is to choose it as a suitable affine transformation of the negative log likelihood of data. In such case, solving (2) amounts to finding maximum likelihood (ML) solutions,
which works well when (2) has a unique solution (uniqueness) that depends continuously on
data (stability). This is however not the case when (1) is ill-posed, in which case one needs to
use regularization that introduces stability, and preferably uniqueness, by making use of prior
knowledge about ftrue .
Variational regularization. The idea here is to add a penalization term to the objective
functional in (2) resulting in a variational problem of the type
h
i
(3)
min µR(f ) + D T (f ), g
for some given µ ≥ 0.
f ∈X

Such regularization methods have gained much attention lately, and especially so in imaging
[52]. The functional R : X → R introduces stability, and preferably also uniqueness, often
by encoding some priori known regularity property of ftrue , e.g., assuming X ⊂ L 2 (Ω, R)
and taking the L 2 -norm of the gradient magnitude (Dirichlet energy) is known to produce
smooth solutions whereas taking the L 1 -norm of the gradient magnitude (total variation)
yields solutions that preserve edges while smooth variations may be suppressed [18].
4. Reproducing kernel Hilbert spaces. The diffeomorphic deformations constructed in
the LDDMM framework (subsection 5.3.1) will make use of velocity fields that at each time
point are contained in a reproducing kernel Hilbert space (RKHS). The gradient computations
in section 8 also rely on this assumption. The short introduction to the theory of RKHS
provided here gives the necessary background needed in subsequent sections.
The theory of reproducing kernel Hilbert spaces was initialized in the 1940s for spaces
of scalar-valued functions [4], which was later extended to spaces of functions with values in
locally convex topological spaces [55]. It has lately gained significant interest due to applications in machine learning [54, 12, 44]. The starting point is an abstract Hilbert space H
whose elements are functions defined on a fixed domain Ω ⊂ Rn that take values in a real
Hilbert space U . Such a space is a RKHS if evaluation functionals δxa : H → R defined as
δxa (ν) := ν(x), a U are bounded for every x ∈ Ω and a ∈ U . One way to construct a RKHS is
to specify a reproducing kernel function. It maps a pair of points in Ω into the Banach space
L (U, U ) of bounded linear operators on U , i.e., an operator K : Ω × Ω → L (U, U ) such that
(i) K( · , x)(a) ∈ H for all x ∈ Ω and a ∈ U .
(ii) The reproducing property holds for K, i.e., if x ∈ Ω then
(4)

ν(x), a

U

= ν, K( · , x)(a)

H

for any ν ∈ H and a ∈ U .

4

C. CHEN AND O. ÖKTEM

An important characterization is that a Hilbert space H of U -valued functions that is
continuously embedded in C (Ω, U ) is a RKHS if and only if it has a continuous reproducing
kernel K : Ω × Ω → L (U, U ).
Square integrable vector fields. If H is a Hilbert space of vector fields on Rn that is admissible (see subsection 5.3.1), then it is in particular continuously embedded in C (Ω, Rn ).
Hence, it is a RKHS if and only if it has a continuous positive definite reproducing kernel
K : Ω × Ω → L (Rn , Rn ), which in turn can be represented by a continuous positive definite
n,m denotes the vector space of all (n × m) matrices
function K : Ω × Ω → Mn×n
+ . Here, M
n,m
and M+ denotes those matrices that are positive definite. Assuming in addition that H is a
RKHS that is continuously embedded in L2 (Ω, Rn ), then there exists an continuous imbedding
ι : H → L 2 (Ω, Rn ). Hence, by Riesz theorem there exist ι∗ : L 2 (Ω, Rn ) → H such that
(5)

ι(ν), η

L 2 (Ω,Rn )

= ν, ι∗ (η)

for all ν, η ∈ H.

H

A natural task is to further examine the relation between the L 2 -inner product and the RKHS
inner product on H. As we next show, one can prove that
D Z
 E
(6)
hη, νiL 2 = ν, K( · , x) η(x) dx .
H

Ω

To prove (6), observe first that by (5)
D
E
ι(ν), ι(η) L 2 (Ω,Rn ) = ν, ι ι∗ (η)

H

for any ν, η ∈ H.

Hence, (6) follows directly from
∗


ι ι (η) =

(7)

Z


K( · , x) η(x) dx.

Ω

To prove (7), observe first that it is sufficient to prove it point wise, i.e., to show that
Z D
D
E

 E
∗
(8)
ι ι (η) (y), a n =
K(y, x) η(x) , a n dx holds for any y ∈ Ω and a ∈ Rn .
R

Ω

R

The equality in (8) follows from the next calculation:
D
E
E


 D


 D
E
ι ι∗ (η) (y), a n = By (4) = ι∗ ι(η) , K( · , y)(a) = By (5) = ι(η), ι K( · , y)(a)
R
H
L2
Z D
Z D
E
 E
=
ι(η)(x), ι K( · , y)(a) (x) n dx =
η(x), K(x, y)(a) n dx.
Ω

R

Ω

R

5. Shape theory. The overall aim is to develop a quantitative framework for studying
shapes and their variability. The approach considered here is based on deformable templates
and it can be traced back to work by D’Arcy Thompson in beginning of 1900’s [24]. Shape
theory based on deformable templates is now an active field of research [31, 65, 58, 43].
The starting point is to specify a shape space, which here refers to a set whose elements
are image features with shapes that are to be analysed. There is no formal definition of

INDIRECT IMAGE REGISTRATION WITH LARGE DIFFEOMORPHIC DEFORMATIONS

5

what constitutes an image feature, but intuitively it is a representation of the image that
retains information relevant for its interpretation. Hence, the choice of shape space is highly
task dependent and examples are landmark points, curves, surfaces, and scalar/vector/tensor
valued images. This is followed by specifying a set of deformations whose elements map
the shape space into itself and thereby model shape changes. Clearly, the identity mapping
preserves all aspects of shape. And in this context, it is therefore the “smallest” deformation.
A key step will be to introduce a metric on the set of deformations that induces a shape
similarity measure in shape space.
5.1. Shape space. There are many image features important for interpretation [7], see
also [38, 39] for an axiomatic characterization of features natural for visual perception. For
(indirect) image registration, it is important that such features are deformable. Namely, it
should be feasible to act on it by means of a deformation.
T
This paper considers grey-scale images, more precisely X = SBV (Ω, R) L ∞ (Ω, R). The
requirement that elements in X are in L ∞ (Ω, R) (essentially bounded) is reasonable for images
since grey-scale values are bounded. Furthermore, SBV (Ω, R) denotes the set of real-valued
functions of special bounded variation over Ω ⊂ Rn . This space was first introduced in [3]
(see also [26]) and its formal definition is somewhat involved, see [59, Definition 22 on p. 141].
Intuitively, these are functions of bounded variation whose singular part of the distributional
derivative is supported by an (n − 1)-dimensional rectifiable set in Ω. In particular, such
functions are not necessarily smooth. The latter is important since images may contain edges.
Occasionally X will be equipped with an L 2 –inner product.
5.2. Set of deformations. Deformations are operations that transform elements in a
shape space. In our setting, deformations is represented by mappings from the image domain
Ω ⊂ Rn into Rn along with an action describing how they deform deformable objects. The
set of deformations should be rich enough to capture the shape variability arising in the
application, yet restricted enough to allow for a coherent mathematical theory.
Considering only rigid body motions is often too limited since many applications involve
non-rigid shape deformations. Furthermore, composing two deformations should yield a valid
deformation, i.e., the set G of deformations closed under composition. Here, the identity
mapping becomes the natural “zero” deformation. Furthermore, it should also be possible
to reverse a deformation, i.e., G closed under inversion. Taken together, this implies that G
forms a group under the group law given by composition of functions. The group structure
also implies that a deformation transforms an element in the shape space by means of a group
action. Finally, often it also makes sense to assume that deformable objects do not tear, fold
or collapse the image, i.e., the deformation preserving the topology.
C p -diffeomorphisms (see below) form a group of non-rigid transformations satisfying the
above requirements. These mappings can also act on real valued functions as in subsection 5.4.
n
o
Diff p (Rn ) := φ ∈ C p (Rn , Rn ) : φ is bijective with φ−1 ∈ C p (Rn , Rn ) .
However, formulating and proving mathematical results relevant for (indirect) image registration, and the need to perform computations with such objects, require one to further restrict
the set of diffeomorphisms, i.e., considering suitable sub-groups of Diff p (Rn ).

6

C. CHEN AND O. ÖKTEM

5.3. Large diffeomorphic deformations. A natural approach is to consider deformations
given by additively perturbing the identity map with a vector field that is sufficiently regular
[46]. The resulting set of deformations is closed under composition, i.e., it forms a semi-group.
Furthermore, such a deformation can be seen as linear approximation of a diffeomorphism,
but it is not necessarily invertible unless the aforementioned vector field is sufficiently small
and regular [65, Proposition 8.6]. Hence, this is not a suitable framework for diffeomorphic
image registration where the template and target have large differences.
One approach to address the above issue of invertibility was presented in [21]. The idea
is to consider transformations given as a composition of infinitesimally small linearised deformations. More precisely, diffeomorphic large deformations are obtained by integrating the
identity map along a velocity field (curve of vector fields). If properly designed, the velocity field induces an isotopy (a curve of diffeomorphisms). Thereby, the set of velocity fields
parametrises the set of deformations and a natural requirement is that velocity fields under
consideration have integrable trajectories. This provides a framework for generating a rich
sets of deformations, which then can be adapted to specific image registration problems by
appropriately choosing the velocity fields.
The indirect image registration method in this paper is based on large deformation diffeomorphic metric mapping (LDDMM), see subsection 11.2 for a brief survey of alternative
frameworks for large diffeomorphic deformations. In LDDMM, the velocity field is square integrable and its corresponding vector field at each time is contained in a fixed Hilbert space of
vector fields that is an admissible RKHS. The deformations are generated by integrating the
identity map along trajectories of such velocity fields and these form a sub-group of Diff 1 (Rn ),
see subsection 5.3.1 for more details. As shown there, this group can furthermore be equipped
with a metric, which by a group action (subsection 5.4) induces a metric on the shape space
that can be used to quantify shape similarity (section 6).
5.3.1. The LDDMM framework. The starting point in the LDDMM framework for large
diffeomorphic deformations is to define the vector space of velocity fields that generate the
diffeomorphisms. An important point is to show that these diffeomorphisms for a group.
Another is to define a right invariant metric so that the group is a complete metric space
w.r.t. this metric. The metric is later used to induce a metric on the shape space, which in
turn is used as a regularizer in indirect image registration (section 6). A key part of showing
that indirect image registration is a well-defined regularization method (section 7) is to relate
convergence of velocity fields to convergence of corresponding diffeomorphisms.
Admissible vector fields and the flow equation. The idea in LDDMM is to generate large
diffeomorphic deformations by solving a flow equation. This is only possible when the velocity
field enjoys certain regularity that ensures the flow equation has a unique solution. More
precisely, fix a domain Ω ⊂ Rn and let V ⊂ C0p (Ω, Rn ) denote a p–admissible vector space of
vector fields, i.e., V can be continuously embedded into C0p (Ω, Rn ) under the topology induced
by the k · kp,∞ –norm. Furthermore, let L q ([0, 1], V ) denote the vector space of velocity fields
that are L q -integrable in time and belong to V at any time point, i.e.,
(9)

n
o
L q ([0, 1], V ) := ν : [0, 1] × Ω → Rn : ν(t, · ) ∈ V and kνkL q ([0,1],V ) < ∞

INDIRECT IMAGE REGISTRATION WITH LARGE DIFFEOMORPHIC DEFORMATIONS

with
1

Z
kνkL q ([0,1],V ) :=

ν(t, · )
0

q
V

1/q
dt

7

for an integer q ≥ 1.

This is a normed vector space and if V is a Hilbert space, then L 2 ([0, 1], V ) can also be
equipped with a Hilbert space structure. Focus will henceforth be on the Hilbert space
L 2 ([0, 1], V ), but some statements will involve L 1 ([0, 1], V ). The following theorem forms
the basis of the LDDMM approach, see [65, Appendix C] for its proof.
Theorem 5.1 (Existence of flows). Given ν ∈ L 1 ([0, 1], V ), there exists a unique continuous curve [0, 1] 3 t 7→ φν0,t ∈ Diff p (Rn ) that satisfies the following flow equation:



 d φν (x) = ν t, φν (x)
0,t
0,t
dt
(10)
for all x ∈ Ω.

φν (x) = x.
0,0
By Theorem 5.1, φνs,t := φν0,t ◦ (φν0,s )−1 becomes a well defined element in Diff p (Rn ) for any
s, t ∈ [0, 1], and by [65, Proposition C.6] one gets
(11)

φνs,t = φντ,t ◦ φνs,τ

and

(φνs,t )−1 = φνt,s

for any 0 ≤ s, τ, t ≤ 1.

Diffeomorphisms generated by flows. Assume V ⊂ C0p (Ω, Rn ) is a p-admissible Banach
space and define
n
o
(12)
GV := φν0,1 : Ω → Rn : φν0,1 solves (10) with ν ∈ L 1 ([0, 1], V ) .
As already concluded, GV ⊂ Diff p (Ω) is a sub-group. Next, the norm on L 1 ([0, 1], V ) can be
used to define a (right invariant) metric dG : GV × GV → R+ on GV as
(13)

dG (φ, ψ) :=

inf

ν∈L 1 ([0,1],V )
ψ=φ◦φν
0,1

kνkL 1 ([0,1],V )

for φ, ψ ∈ GV .

With this construction, (GV , dG ) is a complete metric space whenever V is p-admissible [65,
Theorem 8.15]. The right invariance of the metric is important since
dG (φ, ψ) = dG (φ ◦ ϕ, ψ ◦ ϕ)

for any ϕ ∈ GV .

Furthermore, with this metric one can prove existence of a minimizing velcocity field between
two diffeomorphisms in GV . More precisely, if V is p-admissible, then by [65, Theorems 8.18
and 8.20], for any φ, ψ ∈ GV there exists a velocity field ν ∈ L 2 ([0, 1], V ) such that
(14)

dG (φ, ψ) = kνkL 2 ([0,1],V )

and ψ = φ ◦ φν0,1 .

Hence, when computing the metric in (13), the infimum over L 1 ([0, 1], V ) can be replaced by
a minimum over the Hilbert space L 2 ([0, 1], V ).
Developing a theory for image registration based on arbitrary admissible spaces is difficult,
mainly due to issues related to assigning topological and smooth structures on GV . See
subsection 11.1 for a discussion. This paper considers admissible Hilbert space of vector fields
V constructed from RKHS theory (section 4) that are continuously embedded in L 2 (Ω, Rn ).
The specific application to tomography in section 10 makes use of a diagonal Gaussian kernel.

8

C. CHEN AND O. ÖKTEM

Remark 5.2. It is possible to use other RKHS kernels than diagonal ones as long as the
kernel is represented by a function continuously differentiable in both variables up to order p, in
which case V ⊂ C p (Rn , Rn ) [59, Theorem 9]. See [42] for admissible vector fields constructed
using other non-scalar matrix-valued kernels. An alternative way to construct admissible
vector spaces is by introducing a norm defined through a differential (or pseudo-differential)
operator L:
Z
kνk2V :=

L(ν)(x), ν(x)

Rn
s
(Id − 4) with

Rn

dx.

As an example, choosing L :=
4 as the Laplacian operator yields the Sobolev
s
n
n
space V = H (R , R ). See [65, Section 9.2] for more on building admissible vector spaces
from operators, and in particular [65, Theorem 9.12] for the relation to RKHS theory.
Strong and weak convergence results. The aim here is to relate convergence of vector fields
to convergence of corresponding diffeomorphisms. The next result shows that weak convergence of velocity fields implies uniform convergence of corresponding diffeomorphism on every
compact set. This is an important part for proving the results in section 7. The precise
statement is given below, see [59, Theorem 3 on p. 12] for the proof.
Theorem 5.3 (Uniform convergence of compact sets). Let V ⊂ C p (Rn , Rn ) be p–admissible
with p ≥ 1. If ν k → ν in the L 1 ([0, 1], V )-topology (strong convergence of velocity fields),
νk
) → Dj (φν0,t ) uniformly on any compact subset of Ω for any t ∈ [0, 1] and j ≤ p.
then Dj (φ0,t
Likewise, if ν k * ν in the L 2 ([0, 1], V )-topology (weak convergence of velocity fields), then
νk
Dj (φ0,t
) → Dj (φν0,t ) uniformly on any compact subset of Ω for any t ∈ [0, 1] and j ≤ p − 1.
The case p = 1 is of specific interest:
1. ν k → ν in L 1 ([0, 1], V ), then D(φν0,tk ) → D(φν0,t ) uniformly on every compact set.
2. ν k * ν in L 2 ([0, 1], V ), then φν0,tk → φν0,t uniformly on every compact set.
Remark 5.4. The flow map is defined as the mapping L 1 ([0, 1], V ) 3 ν 7→ φν0,1 ∈ GV and
by Theorem 5.3, it is well-defined and continuous. This mapping is also called the “exponential
map” (not to be confused with the exponential map of a Riemannian metric). This exponential
map is the analogue of the exponential map of a finite dimensional Lie group. One can argue
that the group Diff ∞ (Ω) is an infinite dimensional Lie group with Lie algebra the Lie algebra
V provided that Ω is compact. However, the nice relations between a finite dimensional Lie
group and its Lie algebra, cease to exist in the infinite dimensional setting. For instance the
exponential map fails to be one-to-one or surjective near the identity, see [8, Section 1.3.6]
and subsection 11.1 for further discussion on these matters.
5.4. Group actions. Elements in the group GV of diffeomorphisms in (12) can act on images X ⊂ L 2 (Rn , R), which are real-valued functions defined on some fixed bounded domain
Ω ⊂ Rn . The group action defines the operator
W : GV × X → X.

(15)

In imaging, two group actions are natural:
Geometric deformation: Deforms images without changing grey-scale values. This choice
is suitable for inverse problems where shape, not texture, is the main image feature.
(16)

W(φ, I) := I ◦ φ−1

for I ∈ X and φ ∈ GV .

INDIRECT IMAGE REGISTRATION WITH LARGE DIFFEOMORPHIC DEFORMATIONS

9

Mass-preserving deformation: Deformation changes intensities but ensures mass is preserved. This choice is suitable for inverse problems where intensities are allowed to
change while preserving mass.
(17)

W(φ, I) := D(φ−1 ) (I ◦ φ−1 )

for I ∈ X and φ ∈ GV .

In the above, D(φ−1 ) denotes the Jacobian determinant of φ−1 .
6. Indirect image registration. Consider the inverse problem in (1) where the shape space
X has elements representing grey-scale images defined over some fixed image domain Ω ⊂ Rn .
Next, assume a priori that the true (unknown) target image ftrue ∈ X in (1) can be written
as an admissible deformation of a given shape template I : Ω → R, i.e.,
(18)

ftrue = W(φ∗ , I)

on Ω for some φ∗ ∈ GV .

The metric (13) can, through the action of the group GV on the shape space X, be used to
define a shape similarity measure between objects in X. To define the latter, fix an element
in shape space, the template I ∈ X, and introduce the shape functional S( · , I) : X → R+ as
(19)

S(f, I) :=

inf dG (φ, Id)2

for f ∈ X.

φ∈GV
f =W(φ,I)

It is now possible to register I onto an indirectly observed target ftrue by solving the following
variational problem:

h
i

 inf γS(f, I) + µR(f ) + D T (f ), g
f ∈X
(20)

f = W(φ, I) on Ω for some φ ∈ G .
V
In the above, D : Y × Y → R is the data discrepancy functional introduced in (2), R : X → R
is the regularity functional that encodes further regularity properties of the target ftrue that
are known before hand, and S( · , I) : X → R+ is the shape functional defined in (19). Finally,
γ, µ ≥ 0 are regularization parameters in which γ regulates the influence of the a priori shape
information and µ regulates the a priori regularity information.
The constraint in (20) simply states that the solution must be contained in the orbit of I
under the group action. Furthermore,
(21)

S( · , I) ◦ W( · , I) = dG (Id, · )2

on GV .

Hence, (20) can be reformulated as
(22)

inf

φ∈GV

h
i
γdG (Id, φ)2 + L ◦ W(φ, I)

where L : X → R is given by
(23)

L(f ) := µR(f ) + D T (f ), g



for f ∈ X.

10

C. CHEN AND O. ÖKTEM

6.1. Reformulating the variational problem. If GV is a Riemannian manifold, then one
could solve (22) using a Riemannian gradient method [50]. It is however often more convenient
to work with operators defined on vector spaces and GV lacks a natural vector space structure,
e.g., the point-wise sum of two diffeomorphisms is not necessarily a diffeomorphism. This poses
both mathematical and numerical difficulties.
On the other hand, diffeomorphisms in GV are generated by velocity fields in L 1 ([0, 1], V )
through (10) and by (14), it is sufficient to compute dG (Id, · )2 for diffeomorphisms generated
by velocity fields in L 2 ([0, 1], V ). Hence, there is a natural vector space associated to GV ,
namely L 2 ([0, 1], V ), and a straightforward re-formulation of (22) as an optimization over
L 2 ([0, 1], V ) yields
h
i
(24)
inf
γkνk2L 2 ([0,1],V ) + L ◦ V(ν, I)
ν∈L 2 ([0,1],V )

where V : L 2 ([0, 1], V ) × X → X is the deformation operator that is defined as
(25)

V(ν, I) := W(φν0,1 , I)

where φν0,1 ∈ GV solves (10).

If V is admissible, then (22) and (24) are equivalent as shown in [65, Theorem 11.2 and
Lemma 11.3], i.e.,
(26)

ν ∗ ∈ L 2 ([0, 1], V ) solves (24) ⇐⇒ f ∗ := V(ν ∗ , I) solves (20).

Note also that V is often infinite dimensional Hilbert space, in which case (24) is a minimization over the infinite dimensional Hilbert space L 2 ([0, 1], V ).
6.1.1. PDE constrained formulation. Note that evaluating ν 7→ V(ν, I) requires solving
the ordinary differential equation (ODE) in (10), so the variational problem in (24) is an ODE
constrained optimization problem:





2

min
γkνkL 2 ([0,1],V ) + L ◦ W φ(1, · ), I



ν∈L 2 ([0,1],V )

d
(27)

φ(t, · ) = ν t, φ(t, · ) on Ω and t ∈ [0, 1],


dt



φ(0, · ) = Id on Ω.
For image registration the above can also be formulated as solving a PDE constrained optimization problem with a time dependent image, see, e.g., [34, eq. (1)]. As to be shown next,
such a reformulation easily adapts to the indirect image registration setting.
Re-formulating (27) as a PDE constrained minimization is based on observing that (10)
is strictly related (via the method of characteristics) to a continuity equation. More precisely,
start by considering the time derivative of

(28)
f (t, x) := W φ(t, · ), I (x) for x ∈ Ω and t ∈ [0, 1].
Since φ(t, · ) is given by ν ∈ L 2 ([0, 1], V ) through (10), it should be possible to express the
time derivative of the right-hand-side of (28) entirely in terms of f and ν. Furthermore,

f (0, · ) = I and f (1, · ) = W φ(1, · ), I
on Ω.

INDIRECT IMAGE REGISTRATION WITH LARGE DIFFEOMORPHIC DEFORMATIONS

11

Hence, (27) can be re-stated as a PDE constrained minimization with objective functional

(29)
ν 7→ γkνk2L 2 ([0,1],V ) + L f (1, · ) .
Note that f (1, · ) above depends on ν since it depends on φ(1, · ) through (28) and φ(1, · )
depends on ν by the ODE constraint in (27). The precise form for the PDE depends on the
choice of group action, see subsection 5.4 for a list of some natural group actions.
Geometric group action. Let f : [0, 1] × Ω → R be given by (28) and consider the group
action is given by (16). Then,

(30)
f t, φ(t, · ) = I on Ω for t ∈ [0, 1] where φ solves (10).
Differentiating (30) w.r.t. time t and using (10) yields


 dφ(t, · )
∂t f t, φ(t, · ) + ∇f t, φ(t, · ) ,
dt



Rn

 D

E
= 0.
= ∂t f t, φ(t, · ) + ∇f t, φ(t, · ) , ν t, φ(t, · )
n
R

Since the above holds on Ω, it is equivalent to
D
E
(31)
∂t f (t, · ) + ∇f (t, · ), ν(t, · ) n= 0

on Ω for t ∈ [0, 1].

R

Hence, using the geometric group action in (24) yields the following PDE constrained formulation:
h

i
2

min
γkνk
+
L
f
(1,
·
)

L 2 ([0,1],V )
ν∈L 2 ([0,1],V )


D
E
(32)
∂t f (t, · ) + ∇f (t, · ), ν(t, · ) n= 0 on Ω and t ∈ [0, 1],


R



f (0, · ) = I on Ω.
Mass-preserving group action. Let f : [0, 1] × Ω → R be given by (28) and consider the
group action is given by (17). Then,
(33)



D φ(t, · ) f t, φ(t, · ) = I

on Ω and t ∈ [0, 1] with φ solving (10).

Differentiating (33) w.r.t. time t yields






dφ(t, · )
∂t f t, φ(t, · ) + ∇ · f (t, · )
= ∂t f t, φ(t, · ) + div f (t, · ) ν t, φ(t, · ) = 0.
dt
The last equality above makes use of (10) and the definition of the divergence operator. Since
the above holds on Ω, it is equivalent to

(34)
∂t f (t, · ) + div f (t, · ) ν(t, · ) = 0 on Ω for t ∈ [0, 1].

12

C. CHEN AND O. ÖKTEM

Hence, using the mass-preserving group action in (24) yields the following PDE constrained
formulation:

h
i
2

min
γkνk
+
L
f
(1,
·
)

L 2 ([0,1],V )


ν∈L 2 ([0,1],V )

(35)
∂t f (t, · ) + div f (t, · ) ν(t, · ) = 0 on Ω and t ∈ [0, 1],




f (0, · ) = I on Ω.
Remark 6.1. PDEs are frequently used for defining similarity measures between images.
As an example, registration by optical flow handles shapes as boundaries of objects which are
then treated as fluids [10]. See also [57] for further examples from medical imaging.
In our setting, the PDE constraint encodes both the ODE in (10) and the diffeomorphic
group action. An advantage with the ODE constrained formulation is that these two components, the generative model for the diffeomorphisms and their action on images, are explicit
whereas in the PDE constrained formulation they are “hidden” in the PDE. Another potential
advantage relates to computational feasibility. Within the LDDMM theory, there are several
numerical methods based on various characterizations of minimizers to (24), see e.g., [65].
On the other hand, the PDE constrained formulation may be more suitable for considering
a setting where the velocity fields are non-smooth. As outlined in subsection 5.3.1, admissibility was a key assumption for ensuring that (10) is uniquely solvable. In such case, the
velocity field is Lipschitz w.r.t. to space uniformly in time. An alternative to the ODE based
arguments in subsection 5.3.1 is to use Cauchy-Lipschitz theory and classical PDE arguments
to prove existence and uniqueness for (10). These PDE based techniques extend to certain
non-smooth cases, e.g., [28, 2] proves existence and uniqueness of a certain solution to (10)
(regular Lagrangian flow) even in the case for Sobolev and BV vector fields, see [22] for further
details. Finally, the PDE constrained formulation may also be better suited in applications
where the actual image trajectory is of interest alongside the final deformed image.
7. Regularizing properties. The goal here is to investigate the regularizing properties of
the variational reconstruction scheme in (20) for solving (1) under the simplification where
µ = 0. This involves determining wether (20) has solutions, if they are unique, and assessing
stability and convergence properties of these solutions.
Following regularization theory, a reconstruction scheme is a regularization if existence
holds along with stability and convergence [52, Chapter 3] (see also the notion of “well-defined
regularization method” in [30]). Hence, even though formal uniqueness is a desirable property,
it is not required from a regularization scheme. The starting point is thus to state continuity
of the deformation operator w.r.t. to the velocity field (subsection 7.1). This is followed by
arguments showing that regularizing properties for indirect registration can be transferred
to regularizing properties of (20) (subsection 7.2). The analysis of indirect registration is
concluded by studying existence (subsection 7.3), uniqueness (subsection 7.4), and stability
& convergence (subsection 7.5).
7.1. Convergence of deformed templates. The aim here is to couple weak convergence
of velocity fields T
in L 2 ([0, 1], V ) to L 2 –convergence of corresponding deformed templates in
X = SBV (Ω, R) L ∞ (Ω, R).

INDIRECT IMAGE REGISTRATION WITH LARGE DIFFEOMORPHIC DEFORMATIONS

13

Theorem 7.1. Assume Ω ⊂ Rn isTa bounded domain, V ⊂ C01 (Ω, Rn ) is an admissible
Hilbert space, and X = SBV (Ω, R) L ∞ (Ω, R). Furthermore, let I ∈ X and consider a
sequence {ν k }k ⊂ L 2 ([0, 1], V ) that converges weakly to some ν ∈ L 2 ([0, 1], V ), i.e., ν k * ν
in L 2 ([0, 1], V ). Then,
(36)

W(φν0,1k , I) → W(φν0,1 , I)

in L 2 (Ω, R).

In the above, W : GV × X → X is given by (16) for geometric group action and by (17) for
mass preserving group action.
Proof. The claim (36) for the case when W is given by (16) follows directly from [59, Theorem 20 on p. 59]. Consider now the case when W is given
by (17). From Theorem 5.3 we know

ν k −1 
2
that ν k * ν in L ([0, 1], V ) implies that D (φ0,1 )
converges uniformly to D (φν0,1 )−1 on
compact subsets of Ω. Combined with previous argument proves the claim in (36) also when
W is given by (17). This concludes the proof.
Remark 7.2. The above theorem handles the case of a non-differentiable template. If
the template I is differentiable, then one does not have to introduce SBV -functions and the
corresponding result follows from the more widely know result [65, Theorem 8.11].
It is also possible to state Theorem 7.1 in terms of the deformation operator given in (25). The
corresponding claim, which is given in Corollary 7.3, follows directly from Theorem 5.3, which
states that weak convergence of velocity fields in L 2 ([0, 1], V ) implies uniform convergence of
corresponding diffeomorphisms on compact subsets.
Corollary 7.3. Consider the deformation operator V : L 2 ([0, 1], V )×X → X defined by (25)
and assume the assumptions in Theorem 7.1 hold. Then V(ν k , I) → V(ν, I) in L 2 (Ω, R).
7.2. Reformulation. With the a priori assumption in (18), the inverse problem in (1)
can be re-phrased as the indirect image registration problem of recovering the velocity field
ν ∗ ∈ L 2 ([0, 1], V ) given data g ∈ Y and template I ∈ X such that

(37)
g = T V(ν ∗ , I) + ∆g.
In the above, V( · , I) : L 2 ([0, 1], V ) → X is given by (25) whereas T : X → Y and ∆g ∈ Y
come from (1). The variational scheme for solving (37) is given in (24). Hence, a natural
question is whether regularizing properties for (24) carry over to (20).
By (26), any solution for (22) yields a solution for (20) whenever L : X → R is given as in
(23). Next, by [65, Lemma 11.3], (22) is equivalent to (24). Hence, existence
and uniqueness

of solutions for (24) implies the same for (20). Next, if f → D T (f ), g is continuous, then
by Corollary 7.3, stability and convergence for (24) implies the same for (20). Hence, under
the priori assumption in (18), regulaizing properties of (24) carry over to (20).
Finally, we consider the somewhat simplified case when µ = 0 in (20), so the corresponding
image registration scheme in (24) that will be analyzed is to minimize the objective functional
Eg,λ : L 2 ([0, 1], V ) → R given by

 
(38)
Eg,γ (ν) := γkνk2L 2 ([0,1],V ) + D T V(ν, I) , g
for ν ∈ L 2 ([0, 1], V )
with V : L 2 ([0, 1], V ) × X → X given by (25) and D : Y × Y → R+ is the data discrepancy
functional introduced in (2).

14

C. CHEN AND O. ÖKTEM

7.3. Existence. Existence of solutions is a basic property, since otherwise the reconstruction scheme is of limited usefulness for reconstruction.
Theorem 7.4 (Existence). Assume Ω T
⊂ Rn is a bounded domain, V ⊂ C01 (Ω, Rn ) is an
admissible Hilbert space, I ∈ SBV (Ω, R) L ∞ (Ω, R), and both D T ( · ), g , R( · ) : X → R
are lower semi-continuous. Then, Eg,λ : L 2 ([0, 1], V ) → R in (38) has a minimizer, i.e.,
the variational problem in (24) has a solution expressible as V(ν, I) ∈ X for some ν ∈
L 2 ([0, 1], V ).
Proof. First, V is separable so L 2 ([0, 1], V ) is separable, which in turn implies that
bounded balls are weakly compact. Next, let {ν k }k ⊂ L 2 ([0, 1], V ) be a minimizing sequence.
Then, there exists a subsequence that weakly converges to an element ν ∈ L 2 ([0, 1], V ), i.e.,
ν k * ν in L 2 ([0, 1], V ). Now, if the objective functional in (38) is lower semi-continuous in
L 2 ([0, 1], V ) w.r.t. the weak topology and if, then ν is also minimizer. This proves existence
of solutions for (24).
From the above, existence of solutions for (24) holds whenever the functional in (38) is
lower semi-continuous on L 2 ([0, 1], V ) w.r.t. the weak topology.
T This requires making use of
the admissibility of V and the assumption that I ∈ SBV (Ω, R) L ∞ (Ω, R). By Theorem 7.1
ν k * ν in L 2 ([0, 1], V )

=⇒

I ◦ φν1,0k → I ◦ φν1,0 in L 2 (Ω, R).

Furthermore, φν1,0k and D(φν1,0k ) converge uniformly on Ω to φν1,0 and D(φν1,0 ), respectively.
In particular, ν 7→ V(ν, I) is continuous in the weaktopology on L 2 ([0, 1], V ) for the group
actions listed in subsection 5.4. Next, f 7→ D T (f ), g is by assumption lower semi-continuous
on X, so ν 7→ L ◦ V(ν, I) is lower semi-continuous in the weak topology of L 2 ([0, 1], V ), i.e.,
the objective functional in (24) lower semi-continuous. This concludes the proof.
7.4. Uniqueness. Another desirable property is uniqueness, which in this context means
(24) has a unique solution for given data. Unfortunately, ν 7→ V(ν, I) is not necessarily convex,
so the objective functional in (24) is not necessarily convex even when both f → R(f ) and
f → D T (f ), g are strictly convex. Hence, (24) may not have a unique solution.
As with all reconstruction methods that involve solving non-convex optimization problems,
there is always the issue of getting stuck in local extrema. One option to address this is to
further restrict the set V of velocity fields, but it is highly non-trivial to work out conditions
on V that would guarantee uniqueness.
7.5. Stability and convergence. Stability essentially means that reconstructions are continuous w.r.t. data and convergence refers to the case when data in (1) tends to ideal data,
i.e., as g → T (ftrue ) in Y . Both are essential parts of a regularization scheme and this section proves stability and convergence results. Preferably such results are complemented with
convergence rates and stability estimates, but they are left for the future.
The starting point is to show that minimizing Eg,γ : L 2 ([0, 1], V ) → R in (38) defines a
reconstruction scheme for (37) that is stable in the weak sense. The precise formulation reads
as follows:
Theorem 7.5 (Stability). In addition to the assumptions in Theorem 7.4, let g ∈ Y be fixed,
X = L 2 (Ω, R), the forward operator T : X → Y is continuous, and the data discrepancy

INDIRECT IMAGE REGISTRATION WITH LARGE DIFFEOMORPHIC DEFORMATIONS

15


D T ( · ), g : X → R+ is continuous. Finally, assume furthermore that {gk }k ⊂ Y where
gk → g. If ν k ∈ L 2 ([0, 1], V ) as a minimizer of ν 7→ Egk ,γ (ν) for some k and γ ≥ 0, then
there exists a sub-sequence of {ν k }k that converges weakly towards a minimizer of Eg,γ in (38).
Proof. First, since the data discrepancy is positive and ν k ∈ L 2 ([0, 1], V ) as a minimizer
of ν 7→ Egk ,γ (ν) in (38), one has
kν k k2L 2 ([0,1],V ) ≤

1
1
Egk ,γ (ν k ) ≤ Egk ,γ (0)
γ
γ

for each k and γ ≥ 0.

In the above, 0 ∈ L 2 ([0, 1], V ) denotes the zero-velocity field. This velocity field has zero
L 2 ([0, 1], V )-norm. Next, both the forward operator and the data discrepancy are continuous,
so


 
 
Egk ,γ (0) = D T V(0, I) , gk → D T V(0, I) , g
as gk → g in Y .
Hence, Egk ,γ (0) is a bounded sequence and therefore {ν k }k ⊂ L 2 ([0, 1], V ) is bounded. Hence,
b ∈ L 2 ([0, 1], V ). The continuity of the
it has a subsequence that converges weakly to some ν
forward operator along with Corollary 7.3 implies


b in L 2 ([0, 1], V ).
T V(ν k , I) → T V(b
ν , I) in Y as ν k * ν
Since the data discrepancy is continuous,


 
 
D T V(ν k , I) , gk → D T V(b
ν , I) , g
as gk → g in Y .
Finally, kb
ν kL 2 ([0,1],V ) ≤ lim inf k kν k kL 2 ([0,1],V ) and ν k is a minimizer of Egk ,γ , so
Eg,γ (b
ν ) ≤ lim inf Egk ,γ (ν k ) ≤ lim inf Egk ,γ (ν)
k

k

for any ν ∈ L 2 ([0, 1], V ).

b is a
Hence, Eg,γ (b
ν ) ≤ Eg,γ (ν) as gk → g, but this holds for any ν ∈ L 2 ([0, 1], V ), so ν
minimizer of Eg,γ .
In conclusion, the reconstructed velocity fields depend continuously on data in the weak
sense, which in turn concludes the proof.
Remark 7.6. A natural open question is to study whether the same result holds in the
strong sense, i.e., whether there exists a sub-sequence of {ν k }k that converges strongly to a
minimizer of Eg,γ .
Next is to consider convergence, again in the weak sense.
Theorem 7.7 (Convergence). Assume the conditions
in Theorem 7.5 holds and there exists

some ν ∗ ∈ L 2 ([0, 1], V ) such that T V(ν ∗ , I) = g. Furthermore, let {gk }k ⊂ Y be sequence
of data and δk is a sequence of positive real numbers converging to 0 such that
D(gk , g) ≤ δk2

for each k.

Finally, assume that there is a parameter choice rule γ : R → R such that
(39)

γ(δ) → 0

and

δ2
→0
γ(δ)

when δ → 0,

16

C. CHEN AND O. ÖKTEM

and let ν k ∈ L 2 ([0, 1], V ) be a minimizer of Egk ,γ(δk ) : L 2 ([0, 1], V ) → R given by (38). Then,
there exists a sub-sequence of {ν k }k that converges weakly towards an element in the pre-image
of g under the mapping T V( · , I) .

Proof. By assumption, there is a ν ∗ ∈ L 2 ([0, 1], V ) such that T V(ν ∗ , I) = g. Then,
kν k k2L 2 ([0,1],V ) ≤

1
1
Egk ,γk (ν k ) ≤ Egk ,γk (ν ∗ )
γk
γk

where γk := γ(δk ).

The last inequality above follows from the assumption that ν k is a minimizer of ν 7→ Egk ,γk (ν).
Moreover, (39) gives
i
1
1h
Egk ,γk (ν ∗ ) =
D(g, gk ) + γk kν ∗ k2L 2 ([0,1],V )
γk
γk
δk 2
≤
+ kν ∗ k2L 2 ([0,1],V ) → kν ∗ k2L 2 ([0,1],V )
γk

as δk → 0.

Hence, {ν k }k ⊂ L 2 ([0, 1], V ) is bounded, so it has a sub-sequence that converges weakly to
e ∈ L 2 ([0, 1], V ). Since both the forward operator and data discrepancy are continuous, so
ν


 
 
ν , I) , g
as δk → 0.
(40)
D T V(ν k , I) , gk → D T V(e
Furthermore,

 
(41) D T V(ν k , I) , gk ≤ Egk ,γ (ν k ) ≤ Egk ,γ (ν ∗ ) = D(g, gk ) + γk ν ∗

2
L 2 ([0,1],V )

for any k.

However, D(g, gk ) → 0 as δk → 0 and γk → 0 by (39), so (40) and (41) implies that
T V(e
ν , I) = g, which implies that {ν k }k ⊂ L 2 ([0, 1], V ) has a sub-sequence
that converges

weakly as data error goes to zero to ν ∗ ∈ L 2 ([0, 1], V ) where T V(ν ∗ , I) = g. This concludes
the proof of convergence.
8. Derivative and gradient computations. To goal here is to compute the derivative and
gradient of the objective functional E : L 2 ([0, 1], V ) → R in (24), i.e., so
(42)

E(ν) := γkνk2L 2 ([0,1],V ) + L ◦ V(ν, I)

for ν ∈ L 2 ([0, 1], V )

where L : X → R given by (23) and V : L 2 ([0, 1], V ) × X → X given by (25). Since E is
defined on the Hilbert space L 2 ([0, 1], V ), the gradient associated with its Gâteaux derivative
will be based on the Hilbert structure of L 2 ([0, 1], V ), which is the natural inner-product
space for the optimization problem (24). The derivation will
a differentiable template,
T assume
∞
but one can in fact work with templates in SBV (Ω, R) L (Ω, R that are non-smooth as
shown in [59, Chapter 2].
The starting point is to compute the derivative of the deformation operator V (subsection 8.1). Next step is to use this in computing the gradient of the registration functional
L ◦ V( · , I) : L 2 ([0, 1], V ) → R. This is done in subsection 8.2 for an abstract registration
functional L, the case of the 2-norm is worked out in subsection 8.3. The gradient of the

INDIRECT IMAGE REGISTRATION WITH LARGE DIFFEOMORPHIC DEFORMATIONS

17

k · k2L 2 ([0,1],V ) is straightforward (subsection 8.4) and all these pieces are put together in subsection 8.5 to obtain expressions for the gradient of the objective functional in (42). The explicit
expressions for the gradient assume V is an admissible RKHS with a continuous positive defin×n
nite reproducing kernel K : Ω × Ω → L (Rn , Rn ) represented by a function K : Ω × Ω → M+
(see section 4).
8.1. The deformation operator. The aim here is to compute the derivative of the deformation operator V : L 2 ([0, 1], V ) × X → X in (25). The expressions will depend on the group
action, see subsection 8.1.2 for the case with geometric group action and subsection 8.1.3
for the case with mass-preserving group action. The starting point however is to derive the
derivative of diffeomorphisms in GV w.r.t. the underlying velocity field.
8.1.1. Derivative of diffeomorphisms w.r.t. the velocity field. Each diffeomorphism in
GV depends on an underlying velocity field in L 2 ([0, 1], V ). The derivative of such a diffeomorphism w.r.t. the underlying velocity field can be explicitly computed and is given by the
following theorem.
Proposition 8.1. Let ν, η ∈ L 2 ([0, 1], V ) and let φνs,t ∈ GV be a a solution to (10). Then,
(43)
(44)

d ν+η 
φ
d s,t

d
(φν+η
)−1
s,t
d

Z

t




D φντ,t φνs,τ (x) η τ, φνs,τ (x) dτ
=0
s
Z th


i−1 
η τ, φνt,τ (x) dτ
D φνs,τ φνt,s (x)
(x) = −
(x) =

=0

s

for x ∈ Ω and 0 ≤ s, t ≤ 1.
Proof. The first equality (43) follows directly from [31, Lemma 11.5] (see also [65, Theorem 8.10]) and the second, (44), follows directly from [31, Lemma 12.8] (see also [65,
eq. (10.16)]) and (11).
8.1.2. Geometric group action. For the (left) geometric group action in (16), the deformation operator becomes
V(ν, I) = I ◦ (φν0,1 )−1 = I ◦ φν1,0

(45)

where the last equality follows from (11). Its Gâteaux derivative is given by the next result.
Theorem 8.2. Let V be admissible and I ∈ X is differentiable. The deformation operator
V( · , I) : L 2 ([0, 1], V ) → X in (45) is then Gâteaux differentiable at ν ∈ L 2 ([0, 1], V ) and its
Gâteaux derivative is
Z 1


ν
ν
ν
(46)
∂V(ν, I)(η)(x) = −
∇(I ◦ φt,0 ) φ1,t (x) , η t, φ1,t (x)
dt
0

Rn

for x ∈ Ω and η ∈ L 2 ([0, 1], V ).
Proof. The Gâteaux derivative is a linear mapping ∂V(ν, I) : L 2 ([0, 1], V ) → X given as
∂V(ν, I)(η)(x) :=

d
V(ν + η, I)(x)
d

=0

=


d
ν+η −1
I ◦ (φ0,1
)
(x)
d

=0

18

C. CHEN AND O. ÖKTEM

for η ∈ L 2 ([0, 1], V ) and x ∈ Ω. Since I ∈ X ⊂ L 2 (Ω, R) is differentiable, the chain rule
yields

(47)

∇I

∂V(ν, I)(η)(x) =

φν1,0 (x)

 d
ν+η −1 
(φ0,1
)
,
d


=0

for x ∈ Ω.

(x)
Rn

The second term in the scalar product on the right hand side of (47) is the derivative of a flow
with respect to variations in the associated field. The following equation now follows from
(44) in Proposition 8.1.
Z 1
(48)

∂V(ν, I)(η)(x) = −

∇I


 h
 ν
i−1 

ν
ν
, D φ0,t φ1,0 (x)
η t, φ1,t (x)
dt.

φν1,0 (x)

0

Rn

To prove (46), consider first the chain rule:

D(φν0,t ◦ φν1,0 )(x) = D(φν0,t ) φν1,0 (x) ◦ D(φν1,0 )(x).
Next, φν0,t ◦ φν1,0 = φν1,t , so


−1
D(φν0,t ) φν1,0 (x) = D(φν1,t )(x)◦ D(φν1,0 )(x) ,
which in turn implies that,
(49)

h

D(φν0,t ) φν1,0 (x)

i−1

h
i−1
= D(φν1,0 )(x) ◦ D(φν1,t )(x)
.

Inserting (49) into (48) yields
Z 1


h
i−1 

ν
ν
dt
∂V(ν, I)(η)(x) = −
∇I
◦ D(φ1,t )(x)
η t, φ1,t (x)
0
Rn

Z 1
i−1 


∗ 
 h
ν
ν
ν
ν
η t, φ1,t (x)
dt
=−
D(φ1,0 )(x) ∇I φ1,0 (x) , D(φ1,t )(x)
0
Rn

Z 1

−1 

ν
ν
ν
=−
∇(I ◦ φ1,0 )(x), D(φ1,t )(x)
η t, φ1,t (x)
dt
0
Rn
Z 1


ν
ν
ν
=−
∇(I ◦ φt,0 ) φ1,t (x) , η t, φ1,t (x)
dt.
φν1,0 (x)



, D(φν1,0 )(x)

0

Rn

The last equality above follows from

∗ 

∇(I ◦ φν1,0 )(x) = ∇(I ◦ φνt,0 ◦ φν1,t )(x) = D(φν1,t )(x) ∇(I ◦ φνt,0 ) φν1,t (x) .
This concludes the proof of Theorem 8.2.

INDIRECT IMAGE REGISTRATION WITH LARGE DIFFEOMORPHIC DEFORMATIONS

19

8.1.3. Mass-preserving group action. Under the (left) mass-preserving group action in
(17), the deformation operator becomes



(50)
V(ν, I) = D (φν0,1 )−1 I ◦ (φν0,1 )−1 = D(φν1,0 ) I ◦ φν1,0
where the last equality follows from (11). Its Gâteaux derivative is given in the next result.
Theorem 8.3. If I ∈ X is differentiable, then V( · , I) : L 2 ([0, 1], V ) → X in (50) is
Gâteaux differentiable at ν ∈ L 2 ([0, 1], V ) and its Gâteaux derivative is
(51)


∂V(ν, I)(η)(x) = D φν1,0 (x) div I( · )h( · ) ◦ φν1,0 (x) for x ∈ Ω and η ∈ L 2 ([0, 1], V ),
where
Z

1

h(x) := −

(52)




D φνt,0 φν0,t (x) η t, φν0,t (x) dt.

0

Proof. The Gâteaux derivative is the linear mapping ∂V(ν, I) : L 2 ([0, 1], V ) → X given
as
∂V(ν, I)(η)(x) :=
for η ∈

L 2 ([0, 1], V

d ν+η
φ
d 1,0

=0

=

d
ν+η 
ν+η
D φ1,0
(x) I ◦ φ1,0
(x)
d

=0

) and x ∈ Ω. Now, by (43)
Z

=0

d
V(ν + η, I)(x)
d

(x) =

0




D φνt,0 φν1,t (x) η t, φν1,t (x) dt

1

Z 1



D φνt,0 φν0,t ◦ φν1,0 (x) η t, φν0,t ◦ φν1,0 (x) dt
=−
0
 Z

1
 ν


ν
ν
= −
D φt,0 φ0,t ( · ) η t, φ0,t ( · ) dt ◦ φν1,0 (x) = h( · ) ◦ φν1,0 (x),
0

where the last equation uses the definition of h( · ) by (52).
ν+η
(x) at  = 0 gives
Taking a Taylor expansion of x 7→ φ1,0
(53)

ν
ν
φν+η
1,0 (x) = φ1,0 (x) +  (h ◦ φ1,0 )(x) + o()

which in turn implies




D φν+η
(x) = D φν1,0 (x) +  D φν1,0 (x) div(h) ◦ φν1,0 (x) + o().
1,0
Taking the derivative w.r.t.  and evaluating it at  = 0 yields
(54)


d
D φν+η
(x)
1,0
d

=0



= D φν1,0 (x) div(h) ◦ φν1,0 (x).

The chain rule now implies that
(55)


d
I ◦ φν+η
(x)
1,0
d

=0

D
E

= ∇I φν1,0 (x) , (h ◦ φν1,0 )(x)

Rn

.

20

C. CHEN AND O. ÖKTEM

Hence


∂V(ν, I)(η)(x) = D φν1,0 (x) div(h) ◦ φν1,0 ( · )I( · ) ◦ φν1,0 (x)
D
E


+ D φν1,0 (x) ∇I φν1,0 (x) , (h ◦ φν1,0 )(x) n
R


= D φν1,0 (x) div I( · )h( · ) ◦ φν1,0 (x).
This proves (51).
8.2. Registration functionals in the abstract setting. The aim is to compute the derivative and gradient of a registration functional JI : L 2 ([0, 1], V ) → R of the form in (24), i.e.,
JI (ν) := L ◦ V(ν, I)

(56)

for ν ∈ L 2 ([0, 1], V )

where L : X → R is Gâteaux differentiable, and V is given by (25).
By the chain rule one immediately obtains an expression for the Gâteaux derivative of the
registration functional JI in (56):


(57)
∂JI (ν)(η) = ∂L V(ν, I) ∂V(ν, I)(η) for ν, η ∈ L 2 ([0, 1], V ).
Furthermore, the expression for its L 2 ([0, 1], V )–gradient is derived as
D
E

(58)
∇JI (ν), η L 2 ([0,1],V ) = ∇L V(ν, I) , ∂V(ν, I)(η) .
X

L 2 ([0, 1], V

In the above, ∂V(ν, I) :
∂V(ν, I)∗ : X → L 2 ([0, 1], V ), so
(59)

∇JI (ν), η

L 2 ([0,1],V

) → X is a linear map that has a L 2 ([0, 1], V )–adjoint

)

D

 E
= ∂V(ν, I)∗ ∇L V(ν, I) , η

L 2 ([0,1],V )

.

Hence, the L 2 ([0, 1], V )–gradient of JI in (56) is given by


(60)
∇JI (ν) = ∂V(ν, I)∗ ∇L V(ν, I) .
More explicit expressions requires choosing a specific group action and Hilbert space structures on X and V (for gradient expressions). This is done in the following sections.
8.2.1. Geometric group action. The deformation operator is here given by the geometric
group action in (45). Furthermore, consider the case when X has the L 2 –Hilbert space structure and V is a RKHS. Under these conditions, it is possible to provide explicit expressions
for the derivative and gradient of JI .
Theorem 8.4. Let the registration functional JI : L 2 ([0, 1], V ) → R be given as in (56)
where I ∈ X is differentiable and L : X → R is Gâteaux differentiable on X. Furthermore,
assume GV acts on X by means of the (left) geometric group action (16). Then, the Gâteaux
derivative of JI is given as (57), where ∂V(ν, I)(η) : Ω → R is given as (46).
Furthermore, if X ⊂ L 2 (Ω, R) and V is a RKHS with a reproducing kernel represented
n×n
by a symmetric and positive definite function K : Ω × Ω → M+
(see section 4), then the
2
L ([0, 1], V )–gradient is
Z


(61)
∇JI (ν)(t, x) = −
D(φνt,1 )(y) ∇L V(ν, I) φνt,1 (y) K(x, y) · ∇(I ◦ φνt,0 )(y) dy
Ω

for x ∈ Ω and 0 ≤ t ≤ 1.

INDIRECT IMAGE REGISTRATION WITH LARGE DIFFEOMORPHIC DEFORMATIONS

21

Proof. The deformation operator V( · , I) is given by (45) for geometric group action and its
derivative is given by (46) in Theorem 8.2. The gives the statement. To derive the expression
in (61) for the gradient, one needs to use the assumption that X ⊂ L 2 (Ω, R). An expression
for ∂V(ν, I)(η) is given by (46), which then can be inserted into (58):


Z 1D


E
ν
ν
ν
∂JI (ν)(η) = ∇L V(ν, I) ( · ), −
∇(I ◦ φt,0 ) φ1,t ( · ) , η t, φ1,t ( · )
dt
n
=−

L 2 (Ω,R)

R

0

Z 1




∇L V(ν, I) ( · )∇(I ◦ φνt,0 ) φν1,t ( · ) , η t, φν1,t ( · )

0


dt.
L 2 (Ω,Rn )

The last equality makes use of the fact that the inner product in L 2 (Ω, Rn ) (square integrable
Rn -valued functions) is expressible as the inner product in Rn followed by the inner product
in L 2 (Ω, R) (square integrable real valued functions). Note also that


x 7→ ∇L V(ν, I) (x)∇(I ◦ φνt,0 ) φν1,t (x) ,

x 7→ η t, φν1,t (x)
are both mappings from Ω into Rn , so the integrand in the last expression, which is given by
the L 2 (Ω, Rn ) inner product, is well defined.
Now, introduce the variable y := φν1,t (x), so x = φνt,1 (y) by (11), and use the fact that the
inner product is symmetric. This gives
Z 1D
E


dt.
∂JI (ν)(η) = −
η(t, · ), D(φνt,1 )( · ) ∇L V(ν, I) φνt,1 ( · ) ∇(I ◦ φνt,0 )( · ) 2
n
L (Ω,R )

0

Finally, V is a RKHS with reproducing kernel K : Ω × Ω → L (Rn , Rn ), so (6) implies that

Z 1
Z

∂JI (ν)(η) = −
η(t, · ), K( · , x) ν̃(t, x) dx dt
0

Ω

V

where ν̃(t, · ) : Ω → Rn is defined as


ν̃(t, x) := D(φνt,1 )(x) ∇L V(ν, I) φνt,1 (x) ∇(I ◦ φνt,0 )(x).
It is now possible to read off the expression for the L 2 ([0, 1], V )–gradient of JI :
Z

∇JI (ν)(t, x) = − K(x, y) ν̃(t, y) dy,
Ω
n×n
and inserting the matrix valued function K : Ω × Ω → M+
representing the reproducing
kernel yields (61). This concludes the proof of Theorem 8.4.

8.2.2. Mass-preserving group action. The deformation operator is here given by the
mass-preserving group action in (50). Furthermore, consider the case when X has the L 2 –
Hilbert space structure and V is a RKHS. Under these conditions, it is possible to provide
explicit expressions for the derivative and gradient of JI .

22

C. CHEN AND O. ÖKTEM

Theorem 8.5. Let the assumptions in Theorem 8.4 hold with the (left) mass-preserving
group action (17) instead of the geometric one. Then, the Gâteaux derivative of JI is
given as in (57) with ∂V(ν, I)(η) : Ω → R is given as (51). Furthermore, the corresponding L 2 ([0, 1], V )–gradient is
Z
(62)




D(φνt,0 )(y) (I ◦ φνt,0 )(y)K(x, y) · ∇ ∇L V(ν, I) ◦ φνt,1 (y) dy

∇JI (ν)(t, x) =
Ω

for x ∈ Ω and 0 ≤ t ≤ 1.
Proof. Under the mass-preserving group action, the deformation operator V( · , I) is given
by (50) and its derivative is calculated by (51) in Theorem 8.3. By chain rule, the first
statement is immediately proved.
To prove (62), start by inserting the expression for ∂V(ν, I)(η) into (58):
Z 1
∂JI (ν)(η) =






 ν 

ν
ν
ν
ν
∇L V(ν, I) , − D φ1,0 div D φt,0 φ0,t η t, φ0,t I ◦ φ1,0

dt

L 2 (Ω,R)

0

Z

1

D

=

φνt,0



I◦

φνt,0 ∇






∇L V(ν, I) ◦ φνt,1 dt, η t, ·

0


.
L 2 (Ω,Rn )

Now, (62) follows from combining the above with (6) and reading off the gradient term.
8.3. L 2 -based registration functionals. Consider (1) where X ⊂ L 2 (Ω, R) and Y ⊂
L 2 (M, R) are Hilbert spaces endowed with the L 2 –Hilbert space structure. Here, M denotes a smooth manifold providing coordinates for data. Finally, assume that T is Gâteaux
differentiable and ∆g is independent of T (ftrue ) and Gaussian. Given these assumptions, the
natural registration functional is
(63)

L(f ) := T (f ) − g

2
Y

Z
T (f )(y) − g(y)

=

2

dy

for f ∈ X,

M

for given data g ∈ Y . The corresponding registration functional JI : L 2 ([0, 1], V ) → R is then
(64)


JI (ν) := T V(ν, I) − g

2
L 2 (M,R)

.

Note also that the Gâteaux derivative of L and its corresponding gradient are given by

(65)

D
 E
∂L(f )(h) = 2 ∂T (f )∗ T (f ) − g , h
X

∇L(f ) = 2∂T (f )∗ T (f ) − g

for f, h ∈ X.

Here, “∗” denotes the Hilbert space adjoint and ∂T is the Gâteaux derivative of T . When
T is linear, then ∂T (f ) = T in (65). Next, compute the Gâteaux derivative of JI and its
corresponding gradient under the two group actions.

INDIRECT IMAGE REGISTRATION WITH LARGE DIFFEOMORPHIC DEFORMATIONS

23

8.3.1. Geometric group action. Here, the registration functional is given by (64) with
the deformation operator given as in (45), i.e.,
JI (ν) =

(66)

Z 

2

T I ◦ φν1,0 (x) − g(x) dx.

M

Corollary 8.6. Let the assumptions in Theorem 8.4 hold and let JI : L 2 ([0, 1], V ) → R be
given as in (66) with T : X → Y Gâteaux differentiable. Then
(67)
D
E
∗

for ν, η ∈ L 2 ([0, 1], V ),
∂JI (ν)(η) = 2 ∂T V(ν, I) T V(ν, I) − g , ∂V(ν, I)(η) 2
L (Ω,R)

where ∂V(ν, I)(η) : Ω → R is given as (46).
Furthermore, if X ⊂ L 2 (Ω, R) and V is a RKHS with a reproducing kernel represented
by the symmetric and positive definite function K : Ω × Ω → Mn×n
(see section 4), then the
+
2
corresponding L ([0, 1], V )–gradient is
(68) ∇JI (ν)(t, x)
Z
∗


D(φνt,1 )(y) ∂T V(ν, I) T V(ν, I) − g φνt,1 (y) K(x, y) · ∇(I ◦ φνt,0 )(y) dy
=−
Ω

for x ∈ Ω and 0 ≤ t ≤ 1.
Proof. The proof of (67) follows directly from inserting the expression for ∂L in (65) into
(57). Similarly, (68) follows directly from inserting the expression for ∇L in (65) into (61).
8.3.2. Mass-preserving group action. The registration functional is here given by (64)
with the deformation operator given as in (50), i.e.,
(69)

Z  

2
JI (ν) =
T D(φν1,0 ) I ◦ φν1,0 (x) − g(x) dx.
M

Corollary 8.7. Let the assumptions in Corollary 8.6 hold but with mass-preserving group
action instead of the geometric one. Then, η 7→ ∂JI (ν)(η) is given as in (67) but with
∂V(ν, I)(η) : Ω → R is given as (51). Furthermore, the corresponding L 2 ([0, 1], V )–gradient
is
Z
(70) ∇JI (ν)(t, x) = 2
Ω


D φνt,0 (y) (I ◦ φνt,0 )(y)K(x, y)·
h
i

∗ 

∇ ∂T V(ν, I) T V(ν, I) − g ◦ φνt,1 (y) dy

for x ∈ Ω and 0 ≤ t ≤ 1.
The proof is analogous to the proof of Corollary 8.6.

24

C. CHEN AND O. ÖKTEM

8.4. The shape functional. Define Λ : L 2 ([0, 1], V ) → R+ as Λ(ν) := γ dG (Id, φν0,1 )2 .
Then, (19), (21), and (25), yields
(71)

Λ(ν) = γdG (Id, φν0,1 )2 = γkνk2L 2 ([0,1],V )

for ν ∈ L 2 ([0, 1], V ).

The Gâteaux derivative and associated L 2 ([0, 1], V )–gradient of Λ is then
(72)

∂Λ(ν)(η) = 2γhν, ηiL 2 ([0,1],V )
∇Λ(ν) = 2γν

for ν, η ∈ L 2 ([0, 1], V ).

8.5. Gradient of the objective functional. The goal here is to provide expressions for
the gradient of objective functional in (42) for the two group actions, geometric and masspreserving.
8.5.1. Geometric group action. The following generalization of [31, Theorem 16.2] provides an explicit expression for the gradient of the objective functional in (42).
Corollary 8.8. Let the assumptions in Theorem 8.4 hold. Then, the L 2 ([0, 1], V )–gradient
of ν 7→ E(ν) in (42) is
Z


D(φνt,1 )(y) ∇L V(ν, I) φνt,1 (y) K(x, y) · ∇(I ◦ φνt,0 )(y) dy
(73) ∇E(ν)(t, x) = 2γν(t, x) −
Ω

for x ∈ Ω and 0 ≤ t ≤ 1.
This is an immediate consequence of Theorem 8.4 and (72).
8.5.2. Mass-preserving group action. The following result is the version of Corollary 8.8
under the mass-preserving group action where the deformation operator is given as in (50).
Corollary 8.9. Let the assumptions in Theorem 8.5 hold. Then, the L 2 ([0, 1], V )–gradient
of ν 7→ E(ν) in (42) is
Z




D φνt,0 (y) (I ◦φνt,0 )(y)K(x, y)·∇ ∇L V(ν, I) ◦φνt,1 (y) dy
(74) ∇E(ν)(t, x) = 2γν(t, x)+
Ω

for x ∈ Ω and 0 ≤ t ≤ 1.
This is an immediate consequence of Theorem 8.5 and (72).
9. Numerical implementation. The focus here is to describe the numerical methods used
for solving (24). The corresponding implementation is available used in this paper is available from https://github.com/chongchenmath/odl lddmm and it makes use of the Operator
Discretization Library (http://github.com/odlgroup/odl).
9.1. Optimization strategies. Section 8 gives explicit expressions for the L 2 ([0, 1], V )–
gradient of the objective functional in (24). Hence, one can use any optimization algorithm
that makes use of gradient information, like gradient descent. The numerical implementation
of such a gradient descent scheme is outlined in subsection 9.2.
An alternative approach is to consider the gradient in the Hamiltonian form, which directly
incorporates the dimension reduction that arises from the projection onto a finite dimensional

INDIRECT IMAGE REGISTRATION WITH LARGE DIFFEOMORPHIC DEFORMATIONS

25

subset of the RKHS V , see [65, section 11.6.2] (in particular algorithm 4 on [65, p. 275]) for
further details. Yet another approach is the shooting method outlined in [65, section 11.6.4].
Here, one makes use of the characterization of a minimizer based on the momentum geodesic
equations. The numerical implementation can make use of the natural finite dimensional
counterpart obtained by projecting onto a finite dimensional subset of the RKHS V .
Even though the above approaches may require less iterates than gradient descent, they
have an important drawback. They rely on discretizing the vector fields in the RKHS V
by control points that do not remain stationary over iterates, so the control points do not
necessarily remain on a fixed regular grid. Thus, evaluating the velocity field at some point
in time by convolving against a time-dependent kernel function in the RKHS cannot easily
make use of efficient Fast Fourier transform (FFT) based schemes. This becomes in particular
troublesome when the shape space is scalar valued functions defined in 2D/3D. In contrast, the
gradient descent scheme can be formulated so that the RKHS kernel used to evaluate the vector
fields is given by a fixed set of control points. For this reason, the gradient descent scheme
competes favourably regarding computational complexity against the shooting method.
9.2. Gradient descent. The gradient descent scheme for solving (24) is a first-order iterative optimization algorithm given as
ν k+1 = ν k − α∇E(ν k ).

(75)

Here, α > 0 is the step-size and ∇E is the L 2 ([0, 1], V )–gradient of the objective functional
E in (42).
Corollary 8.8 gives an expression for ∇E in (73) where the reproducing kernel function
K : Ω × Ω → Mn×n
is evaluated on points that do not move as iterates proceed. Choosing
+
a translation invariant kernel and points on a regular grid in Ω, allows one to use computationally efficient FFT-based schemes for computing the velocity field at each iterate (see
[65, section 11.7.2]). This is necessary for numerically evaluating ∇E(x) for x ∈ Ω in the
aforementioned grid and it is computationally more feasible than letting the kernel depend on
points that move in time as in the shooting method.
Computing diffeomorphic deformations. Evaluating ν 7→ ∇E(ν) in (75) requires computing
the diffeomorphic deformations φνt,0 and φνt,1 . Since this is done repeatedly within an iterative
scheme, one needs to have an efficient implementation for computing these diffeomorphisms.
First recall that φνs,t solves the ODE in (10) forward in time:
(76)

(

∂t ψ(t, x) = ν t, ψ(t, x)
ψ(s, x) = x

for x ∈ Ω and 0 ≤ s, t ≤ 1

with s a fixed time point. Integrating (76) w.r.t. time t gives
(77)

φνs,t (x)

Z
=x+

t


ν τ, φνs,τ (x) dτ

for 0 ≤ t ≤ 1 and x ∈ Ω.

s

A numerical implementation needs to discretize time, which can be done by sub-dividing
the time interval [0, 1] uniformly into N parts, i.e., ti = i/N for i = 0, 1, . . . , N . Let s = ti ,

26

C. CHEN AND O. ÖKTEM

t = ti+1 and ti−1 , the expressions for small deformations φνti ,ti+1 and φνti ,ti−1 would be derived
as the approximations
1
ν(ti , ·),
N
1
≈ Id − ν(ti , ·).
N

(78)

φνti ,ti+1 ≈ Id +

(79)

φνti ,ti−1

Furthermore, (11) implies that φνti ,0 = φνti−1 ,0 ◦ φνti ,ti−1 , which is combined with (79) yields
the following approximation:


1
φνti ,0 ≈ φνti−1 ,0 ◦ Id − ν(ti , · )
for i = 1, . . . , N .
N

(80)

Similarly, (11) implies φνti ,1 = φνti+1 ,1 ◦ φνti ,ti+1 , which is combined with (78) gives the following
approximation:


1
φνti ,1 ≈ φνti+1 ,1 ◦ Id + ν(ti , · )
for i = N − 1, . . . , 0.
N

(81)

Remark 9.1. The approximation used for deriving (80) and (81) from (79) and (78), respectively implies
(82)



Id +

−1
1
1
ν(ti , · )
≈ Id − ν(ti+1 , · )
N
N

for i = 0, . . . , N − 1.

Gradient descent algorithm. The approximations in the aforementioned paragraph play a
central role in the gradient descent algorithm that implements (75) for minimizing for finding
a local extrema to (24).
The numerical implementation is outlined in Algorithm 1 for the case with geometric group
action. The version for mass-preserving group action can be obtained from Algorithm 1 after
modifying two steps. The first is to replace step 8 with:
k
Update the Jacobian determinant Dφνti ,0 based on (80):




1
1
k
k
div ν k (ti , · ) Dφνti−1 ,0 ◦ Id − ν k (ti , · )
Dφνti ,0 ← 1 −
N
N
k

k

for i = 1, . . . , N , where Dφνt0 ,0 = Dφν0,0 = 1.
The second is to replace step 10 with:
Compute ∇E(ν k ) using (74) (use FFT based techniques for computing the
kernel):
∇E(ν k )(ti , , · ) ← 2γν k (ti , , · )
Z



k 
k 
k
+
D φνti ,0 (y) I ◦ φνti ,0 (y)K( · , y) · ∇ ∇L V(ν k , I) ◦ φνti ,1 (y) dy
Ω

for i = N, N − 1, . . . , 0.

INDIRECT IMAGE REGISTRATION WITH LARGE DIFFEOMORPHIC DEFORMATIONS

27

Remark 9.2. A similar gradient descent scheme is proposed in [11]. Here, the RKHS
structure that regulates the smoothness of the velocity fields is defined by a kernel function
instead of by a differential operator. This leads to a different implementation. Moreover, the
k
k
k
updating strategy in [11] that corresponds to in lines 7-9 is: compute φνti ,0 , φνti ,1 , and Dφνti ,1

k
k
k
first and then compute I ◦ φνti ,0 , ∇L V(ν k , I) ◦ φνti ,1 and Dφνti ,1 , respectively. Compared to
the updating strategy in Algorithm 1 (lines 7-9), the approach in [11] is more costly in terms
of computational effort and memory usage. The reason is that in each iteration Algorithm 1
just needs to update three scalar fields instead of two vector fields, and one tensor field, and
then three scalar fields on the considered domain.
Complexity analysis. The complexity analysis applies to Algorithm 1. Since the main part
(each iteration) of the algorithm is located to lines 7-11, the analysis focuses on those lines.
Suppose that Ω ⊂ R2 and the size of the image to be reconstructed is n × n. In line 7, one
k
updates I ◦ φνti ,0 for i = 1, . . . , N . Moreover, each of them needs to be used to compute the
gradient of the objective functional in line 10 at each time point. Hence they need to be stored
at hand. Through simple analysis, in this step, the computational cost is O(n2 N ) and the
k
space complexity is O(n2 N ) . Similarly, for lines 8 and 9, the Jacobian determinant Dφνti ,1

k
and ∇L V(ν k , I) ◦ φνti ,1 need to be updated and then stored for i = N − 1, . . . , 0. Therefore,
for each of the two steps, the computational cost is O(n2 N ) and the space complexity is
O(n2 N ). Evidently, lines 10-11 can be combined into one step. Line 10 uses FFT to compute
the gradient of the objective functional at each time point. Hence the computational cost for
this line is O(N n3 log n). In line 11, a vector field is updated at each time point and this
part requires twice the memory compared to a scalar field on 2D domain. In conclusion, the
computational cost is O(N n3 log n) and the memory requirements are O(n2 N ).
10. Application to 2D tomography. Indirect registration with geometric group action is
here applied to 2D parallel beam tomography with very sparse or highly noisy data. Although
this is not a full evaluation, it does still illustrate the performance of indirect registration.
T
10.1. The indirect registration problem. Let X = SBV (Ω, R) L ∞ (Ω, R) for a fixed
image bounded image domain Ω ⊂ R2 . Elements in X represent 2D images and this space
is equipped with an L 2 –inner product. Diffeomorphisms act on such elements through a
geometric group action and the goal is to register a given template I ∈ X against a target
that observed indirectly as in (1). Here, g ∈ Y = L 2 (M, R) with M denoting a fixed manifold
of parallel lines in R2 (parallel beam data). The template is assumed to be differentiable. The
forward operator T : X → Y is here the 2D ray/Radon transform, see [45] for further details.
The indirect registration scheme is given as the solution to (24) where L : X → R in (23)
is given as
2
L(f ) := T (f ) − g L 2 (M,R) for f ∈ X.
Hence, µ = 0 in (23), i.e., there is no additional regularization and (24) reduces to


2

2
ν
(83)
min
γkνkL 2 ([0,1],V ) + T I ◦ φ0,1 − g 2
.
ν∈L 2 ([0,1],V )

L (M,R)

In the above, γ > 0 is a fixed regularization parameter that weights the need for minimal
deformation against the need to register against the indirectly observed target.

28

C. CHEN AND O. ÖKTEM

Algorithm 1 Gradient descent scheme for minimizing E in (42) with geometric group action
1:
2:
3:
4:
5:
6:
7:

Initialize:
k ← 0.
ti ← i/N for i = 0, 1, . . . , N .
ν k (ti , · ) ← ν 0 (ti , · ) where ν 0 (ti , · ) is a given initial vector field.
Error tolerance  > 0, step size α > 0, and maximum iterations K > 0.
while ∇E(ν k ) >  and k < K do
k
Compute I ◦ φνti ,0 using (80) by

 
1
k
k
I ◦ φνti ,0 ← I ◦ φνti−1 ,0 ◦ Id − ν k (ti , · )
N
k

k

for i = 1, . . . , N , where I ◦ φνt0 ,0 = I ◦ φν0,0 = I.
8:

k

Update the Jacobian determinant Dφνti ,1 based on (81):




1
1
k
k
div ν k (ti , · ) Dφνti+1 ,1 ◦ Id + ν k (ti , · )
Dφνti ,1 ← 1 +
N
N
k

9:

k

for i = N − 1, . . . , 0, where DφtνN ,1 = Dφν1,1 = 1.

k
Compute ∇L V(ν k , I) ◦ φνti ,1 using (81) by

 



1
k
k
∇L V(ν k , I) ◦ φνti ,1 ← ∇L V(ν k , I) ◦ φνti+1 ,1 ◦ Id + ν k (ti , · )
N
for i = N − 1, . . . , 0, where



k
k
∇L V(ν k , I) ◦ φνtN ,1 = ∇L V(ν k , I) ◦ φν1,1 = ∇L V(ν k , I) .

10:

Compute ∇E(ν k ) using (73) (use FFT based techniques for computing the kernel):
∇E(ν k )(ti , , · ) ← 2γν k (ti , , · )
Z


k
k
k
−
∇L V(ν k , I) φνti ,1 (y) K( · , y) · ∇(I ◦ φνti ,0 )(y) D(φνti ,1 )(y) dy
Ω

11:

for i = N, N − 1, . . . , 0.
Update ν k :
ν k (ti , · ) ← ν k (ti , · ) − α∇E(ν k )(ti , · )

for i = N, N − 1, . . . , 0.
k ← k + 1.
13: end while
k
14: return I ◦ φν
ti ,0 for i = 1, . . . , N .
12:

INDIRECT IMAGE REGISTRATION WITH LARGE DIFFEOMORPHIC DEFORMATIONS

29

Next, consider a set V of vector fields that is a RKHS with a reproducing kernel represented
by symmetric and positive definite Gaussian function K : Ω×Ω → M2×2
+ . Then V is admissible
2
2
and is continuously embedded in L (Ω, R ). The Gaussian is parameterised by σ > 0 as
 1 0
 1
for x, y ∈ R2 .
(84)
K(x, y) := exp − 2 kx − yk2
0 1
2σ
The kernel width σ also acts as a regularization parameter.
10.2. Phantoms and data acquisition protocol. All images in the shape space X, such
as the target (phantom) and template, are supported in a fixed rectangular image domain
Ω = [−16, 16] × [−16, 16]. Data is obtained by first evaluating the 2D parallel beam ray
transform on M , which consists of parallel lines whose directions are equally distributed
within [0◦ , 180◦ ]. Finally, additive Gaussian white noise at varying noise levels is added to
data.
The noise level in data is quantified in terms of the signal-to-noise ratio (SNR) expressed
using the logarithmic decibel (dB) scale:
"
2 #
gideal − µideal Y
(85)
SNR(g) = 10 log10
for g = gideal + ∆g.
2
∆g − µnoise Y
In the above, gideal := T (ftrue ) is the noise-free component of data and ∆g ∈ Y is the noise
component. Furthermore, µideal is the average of g ideal and µnoise is the average of ∆g. Note
here that computing the data noise level requires access to the noise-free component of data,
which is the case for simulated data.
10.3. Reconstruction protocol. Test suites 1 and 2 compares Algorithm 1 against usual
reconstruction of the target from tomographic data. For the latter, filtered back projection
(FBP) and TV methods are used. The rationale behind comparing against direct reconstruction is to assess the influence of a priori information contained in (18), which states that
the unknown target can be obtained by a diffeomorphic deformation of a given template.
Hence, to some extent the indirect registration problem can be seen as a shape reconstruction
problem. Furthermore, the corresponding regularization parameters are set separately each
method and data noise level as to give the best visual agreement with the target.
The reconstruction methods FBP and TV are the ones available in Operator Discretization
Library (http://github.com/odlgroup/odl).
• FBP reconstruction is obtained using linear interpolation for the back projection and
the Hamming filter at varying frequency scaling.
• TV reconstruction is obtained by solving (3) with R as the total variation functional:
Z
R(f ) :=
∇f (x) dx.
Ω

The non-smooth minimization (3) is solved using the primal-dual method in [19].
10.4. Test suites and results. The test suites seek to assess robustness and performance
against different kind of objects, sensitivity of indirect registration against choice of regularization parameters γ, σ, and impact of using a template with different topology from target.
The targets include both single- and multi-object images.

30

C. CHEN AND O. ÖKTEM

Test suite 1: Single-object indirect registration. The aim here is to investigate the performance of indirect registration against highly noisy sparse tomographic data from a single
object target. The template and target are chosen similar to the direct registration test in [9,
section 4.2].
Images in shape space are digitized using 64 × 64 pixels. The parallel beam tomographic
data consists of 92 parallel line projections along ten uniformly distributed directions in
[0◦ , 180◦ ]. Hence, in a fully digitized setting, the inverse problem in (1) amounts to solving for
4 096 unknowns from 920 linear equations. Hence, this problem is highly under-determined
and data is also highly noisy with a noise level (SNR) of 4.87 dB.
When solving the indirect registration problem, the time interval [0, 1] for the flow of
diffeomorphisms is sampled uniformly at N = 20 time points. The gradient step size is set as
α = 0.02 and the regularization parameters are γ = 10−7 and σ = 6.0. The gradient descent
is stopped after 200 iterations. Finally, the comparison is against direct reconstructions from
FBP (Hamming filter with frequency scaling 0.4) and TV (1000 iterations, regularization
parameter µ = 3.0).
The results, which are shown in Figure 1, show that indirect registration performs fine
against a single-object target as long as the template has the topology (e.g., no holes).
Test suite 2: Multi-object indirect registration. The previous test used a target that only
contains one single object. This test considers a multi-object target, namely a target consisting
of six separately star-like objects with grey-values in [0, 1], which is digitized using 438 × 438
pixels. The choice of target is similar to the direct registration test in [13].
The data is highly sparse parallel beam tomographic projections taken along only six
directions that are uniformly distributed in [0, 180◦ ]. At each direction, one samples the ray
transform in a set of uniformly spaced 620 parallel lines. Hence, the inverse problem in (1)
amounts to reconstructing 191 844 unknowns from 3 720 linear equations. Besides being highly
under-determined, data is also highly noisy with noise level (SNR) of 4.75 dB.
As in test suite 1, the time interval [0, 1] for the flow of diffeomorphisms is sampled uniformly at N = 20 time points, the gradient step size is set to α = 0.04, and the regularization
parameter is γ = 10−7 . On the other hand, the other regularization parameter regulating the
width of the RKHS kernel is set as σ = 2.0 and the gradient descent was stopped after 200
iterations. Finally, comparison is against direct reconstructions from FBP (Hamming filter
with frequency scaling 0.4) and TV (1000 iterations, regularization parameter µ = 1.0).
The results, which are shown in Figure 2, show that indirect registration performs fine
also against a multi-object target as long as the template has the same number of objects.
Test suite 3: Sensitivity w.r.t. choice of regularization parameters. Choosing the regularization parameters is a well-known issue in most regularization schemes. Even though there is
some theory for how to do this, in practice it is often chosen using heuristic methods and
especially so when data is highly noisy and/or under-sampled. A natural question is therefore
to empirically investigate the sensitivity of the indirect registration against variations in the
regularization parameters γ and σ that enter in (83) and (84), respectively.
Here the target is the well-known Shepp-Logan phantom consisting of ten ellipsoids with
grey-values in [0, 1]. Images in shape space are digitized using 256 × 256 pixels. The parallel
beam tomographic data consists of 362 parallel line projections along ten uniformly distributed
directions in [0◦ , 180◦ ]. Hence, in a fully digitized setting, the inverse problem in (1) amounts

INDIRECT IMAGE REGISTRATION WITH LARGE DIFFEOMORPHIC DEFORMATIONS

γ
σ
1.0
2.0
2.5
3.0
4.0
8.0

10−7

10−5

10−3

10−1

10

0.8958
26.69
0.9039
26.96
0.9018
26.90
0.9007
26.90
0.8992
26.88
0.8305
24.15

0.8958
26.69
0.9039
26.96
0.9018
26.90
0.9007
26.90
0.8992
26.88
0.8305
24.15

0.8957
26.69
0.9038
26.96
0.9017
26.89
0.9006
26.89
0.8992
26.87
0.8307
24.16

0.8820
26.18
0.8958
26.72
0.8964
26.74
0.8968
26.77
0.8979
26.77
0.8438
24.84

0.6800
17.82
0.6967
18.45
0.7057
18.71
0.7164
19.03
0.7426
19.79
0.7960
21.27

31

Table 1: Test suite 3: SSIM and PSNR values of indirectly registered images as compared to
target for varying values of the regularization parameters γ and σ, see Figure 3 for selected
images. Each table entry has two values, the upper is the SSIM and the bottom is the PSNR.

to solving for 65 536 unknowns from 3 620 linear equations. Besides being highly underdetermined, data is also highly noisy with noise level (SNR) of 7.06 dB.
As in test suite 1, the time interval [0, 1] for the flow of diffeomorphism is sampled uniformly
at N = 20 time points, the gradient step size is set as α = 0.02, and the regularization
parameter is γ = 10−7 .
The regularization parameters (γ, σ) are varied and each choice results in a different indirectly registered image. Each of these are then matched against the target and the registration
is quantitatively assessed using structural similarity (SSIM) [62] and peak signal-to-noise ratio
(PSNR) [35] as figure of merits.
The SSIM and PSNR values are tabulated in Table 1, see also Figure 3 for some selected
reconstructions. As given in Table 1, the SSIM and PSNR values are quite similar for γ =
10−7 , 10−5 , 10−3 , 10−1 and σ = 1.0, 2.0, 2.5, 3.0, 4.0. The corresponding reconstructed results
in Figure 3 are almost the same and close to the target. But these values are obviously reduced
when γ = 10 or σ = 8.0. As shown in Figure 3, the resulting performance is not so good
compared with the target. Since the choosing intervals for γ and σ are relatively large, this
test shows that to some extent the proposed method is not sensitive to the precise choice of
the regularization parameters γ and σ under quantitative figures of merits (SSIM and PSNR).
However those values are chosen too big that would cause over-regularized results.
Test suite 4: Topology of the template. This test investigates the influence of a template
with a topology that differs from the target. The test suite involves two tests, one where the
the template lacks an object when compared to the target and the other where there is an
additional object in the template that does not exist in the target.
In both cases, the template and the target are modifications of the Shepp-Logan phantom

32

C. CHEN AND O. ÖKTEM

used in test suite 3 and tomographic data is generated following the same protocol as in test
suite 3, albeit with slightly differing noise levels (7.06 dB and 6.46 dB, respectively). Finally,
apart from the number of iterations which here is 1000, indirect registration in both cases was
performed using the same parameter setting as in test suite 3.
Results are shown in Figure 4. It is clear that the final indirectly registered template
retains the same topology, which is to be expected since diffeomorphic registration using
geometric group action can never introduce or remove an object that is not in the template.
This also points to the importance of having a template with the correct topology when using
diffeomorphic (indirect) registration with a geometric group action.
11. Discussion.
11.1. Topology and manifold structure on groups of diffeomorphisms. The LDDMM
framework can be seen as part of an even more abstract formulation as outlined in [16] that
involves replacing the group of diffeomorphisms in (12) with an abstract Lie group G of
transformations that acts on the shape space X, and the vector space V is the corresponding
Lie algebra of G. The situation becomes intricate when one seeks to introduce a useful
manifold structure on G in the infinite dimensional setting as outlined in [8].
As an example, Diff 1 (Ω) and Diff ∞ (Ω) are both Fréchet–Lie groups, but neither is a
Banach manifold. As nicely summarized in [15], working with such deformation groups allows
using methods from geometry, but less tools from analysis are applicable. The other choice is
to consider a group like Diff p (Ω). This group is a Banach manifold and a topological group
so there are powerful tools from analysis for prove existence results. On the other hand, it
is not a Lie group since the group operations are continuous, but not differentiable. Hence,
many of the tools from differential geometry are not available.
One may next consider groups generated as in (12) from an admissible Hilbert space of
vector fields. The class of admissible vector spaces is however very large, which in turn limits
how far one can develop a theory for image registration in this setting. As nicely outlined
in [17, p. 1538], GV in (12) does not need to have a differentiable structure in the general
setting. Furthermore, it does not have be a topological group under the topology induced by
the metric in (13). In fact, there is no natural way to define a topology on GV in the general
setting that makes it a topological group. Hence, one has studied specific classes of admissible
spaces. On the other hand, restricting attention to velocity fields with Sobolev regularity
generate classical groups of Sobolev diffeomorphisms. Such a group is a Hilbert manifold as
well as a topological group, but it is not necessarily a Lie group [17, p. 1512 and Theorem 8.3]
but the exponential map (Remark 5.4) is continuous [17, Theorem 4.4].
11.2. Alternative frameworks for large deformations. LDDMM is capable of generating
a vast range of diffeomorphic deformations, but it comes with some drawbacks coupled to its
non-parametric nature. There is currently no natural parametrisation of deformations using
variables that are easy to interpret from an image registration point-of-view. Furthermore,
representing these non-parametric diffeomorphic deformations in software in a computationally feasible way is challenging and requires introducing further structure. This need for
further structure on the group of diffeomorphic deformations has resulted in a number of
parametric frameworks, some of which are mentioned below. See the excellent introduction

INDIRECT IMAGE REGISTRATION WITH LARGE DIFFEOMORPHIC DEFORMATIONS

33

in [33, p. 30-35] for a more informative survey.
The poly-affine framework [6, 5] considers deformations obtained by integrating trajectories of vector fields that include local affine deformations. The weights for the each local affine
transformation and corresponding regions for its action remain fixed during the integration
of the flow. The approach is successfully used for registration of bones [56] and images with
cardiac motion [41], but the approach is not suitable for very large deformations. In the
GRID model [32, 48] the idea is to fuse local deformations by fusing corresponding velocity
fields. The local deformations can be more generic than in the poly-affine case, but each small
deformation is located around a point (seed) that defines a local region for its action. The
global deformation is a discrete temporal integration of a trajectory of corresponding vector
fields. The framework has been used to model biological growth [32]. The diffeons framework [66] considers deformations are defined as the solution of a flow equation that integrates
the identity map along a trajectory of vector fields. The difference to LDDMM is that one
only considers trajectories of vector fields that can be written as a finite linear combination
of a finite family of vector fields, called diffeons, which form a finite dimensional subspace.
A further development of diffeons is deformation modules [33]. Here, one introduces a new
deformation model that can generate diffeomorphisms that can be locally constrained to a
certain type of deformations with a natural interpretation. The diffeomorphisms are as in
LDDMM defined as final values of a flows of trajectories of vector fields and the constraints
corresponds to setting of a family generators of vector fields similar to the diffeons framework.
In these parametrised frameworks, the trajectory of the vector field is constrained to
be a combination of a few generators. Hence, by choosing these generators it should in
principle be possible to build trajectories of diffeomorphisms that correspond to specific a
priori deformations. It is however very difficult, if not impossible, to define generators that
correspond to more complex deformations that typically arise in many imaging applications.
Another issue that is important for image registration is to have a metric on the shape space.
Another aspect relevant for (indirect) image registration is to let the metric on the group of
deformations induce a metric on the shape space in a canonical manner. This is a key feature
of LDDMM and it also holds for the diffeons framework, whereas this is not possible within
the polyaffine and GRID frameworks.
A final recent development is to construct diffeomorphisms by using techniques from machine learning. Assume one has training data {(Ii , fi )}i ⊂ X × X where fi = W(Ii , φi ) for
some (unknown) deformation φi ∈ GV and W given by (15). One may then use principles from
machine learning to learn how to estimate the deformations from training data. The papers
[63, 64] implements such as scheme. More precisely, they use machine learning in order to
learn how to estimate the optimal diffeomorphic deformations. The idea is to learn the initial
momenta from the training data. Once these initial momenta are learned, one can use the
shooting method from LDDMM theory to generating the smallest diffeomorphic deformation
[65, Section 11.6.4]. The learning is done by using a (deep) convolution neural network that
learns an operator ΛΘ : X × X → GV whose parameters Θ are obtained by minimising the
mean absolute error. The trained parameter yields the initial momenta, so it is now possible
to generate elements in GV by means of the shooting method. The approach extends readily
to the indirect registration
setting. Training data is now {(Ii , gi )}i ⊂ X × Y where (18) holds,

i.e., gi = T W(Ii , φi ) for some (unknown) deformation φi ∈ GV . One can then perform

34

C. CHEN AND O. ÖKTEM

learning using a loss function similar to the direct matching case. This time the learning is
done by using a (deep) convolution neural network that learns an operator ΛΘ : X × Y → GV
whose parameters Θ are obtained by minimising a natural distance notion in data space. Such
an approach is possible especially when combined with [1] for reconstruction.
11.3. Spatiotemporal image reconstruction. The goal here is to recover time dependent
images ft ∈ X from time dependent data gt ∈ Y where gt = T (ft ) + ∆gt with ∆gt denoting

the noise term. A spatiotemporal analogue of (18) is to assume that ft = W φν0,t , I for
some unknown curve t 7→ φν0,t in GV . This allows one to formulate the following variational
approach to the spatiotemporal image reconstruction problem:
"Z  Z
 #
t
1

2
ν(s, ·) V ds + L ◦ W φν0,t , I dt .
γ
(86)
inf
ν∈L 2 ([0,1],V )
I∈X

0

0

In the above, I ∈ X is the starting image. The curve of diffeomorphism t 7→ φν0,t is generated
by the velocity field ν as a solution to (10), W : GV × X → X is given by the group action,
and L : X → R is the spatiotemporal analogue of (23), i.e.,

L(f ) := µR(f ) + D T (f ), gt
for f ∈ X.
The optimization problem in (86) is quite challenging and a natural approach is to consider
an intertwined scheme where a reconstruction (of the template) step is followed by an indirect
registration step, see [46, section 12.3] for further details. In such case, indirect registration
becomes a key element in spatiotemporal image reconstruction.
11.4. Shape based reconstruction. It is possible to view indirect registration as a variational reconstruction scheme that makes use of a priori shape information encoded by the
template. The tests in section 10 indicate that if the template has the correct topology, then
indirect registration performs fairly well. Furthermore, it is also robust against choice of
template and regularization parameters.
When using geometric group action, indirect registration is mostly useful as a reconstruction scheme in imaging problems where the aim is to recover the shape, so intensity variations
are of little, or no, importance. One example of such an imaging problem in applications is
when electron tomography (ET) is used to image the internal 3D structures at nano-scale of a
specimen [13]. This is however a rather limited category of inverse problems and many of the
medical imaging applications do not fall under this category. One approach to address this
is to consider mass-preserving group action instead of the geometric one. Another is to keep
the geometric group action, but let it act on the intensity map. This leads to metamorphosis
that is briefly described in subsection 11.5.
11.5. Metamorphosis. The LDDMM approach with geometric deformations only moves
intensities, it does not change them. Metamorphosis extends LDDMM by allowing the diffeomorphisms to also act on intensities. This is achieved by coupling the flow equation in (10)
to a similar flow equation for the intensities.
More precisely, in metamorphosis one needs to consider time dependent images. Hence,
introduce X as the set of time dependent images that are square integrable in time and

INDIRECT IMAGE REGISTRATION WITH LARGE DIFFEOMORPHIC DEFORMATIONS

35

contained in X at all time points. Hence, if ζ ∈ X then ζ(t, · ) ∈ X and
s
Z 1
2
kζkX :=
ζ(t, · ) X dt < ∞.
0

The Hilbert space structure of X induces a Hilbert space structures on X . Furthermore,
introduce I : [0, 1] × Ω → R as the solution to

 d I(t, x) = ζ t, φ(t, x) (t, x) ∈ [0, 1] × Ω,
dt
(87)

I(0, x) = I(x) x ∈ Ω.
In the above, φ is the diffeomorphism that solves (10), I ∈ X is the template, and ζ ∈ X is a
“source” term that is contained in the space X of time dependent images introduced above.
Bearing in mind the above, metamorphosis based image registration in [49] easily extends
to the indirect setting, which in turn extends the LDDMM based indirect registration in (24):
h

i
2
2

inf
γkνk
+
τ
kζk
+
L
◦
V
ν,
I(1,
·
)

X
L 2 ([0,1],V )
ν∈L 2 ([0,1],V )



ζ∈X

(88)

d


I(t, x) = ζ t, φν0,t (x) (t, x) ∈ [0, 1] × Ω,


dt



I(0, x) = I(x) x ∈ Ω.
Here, L : X → R is given by (23) and V : L 2 ([0, 1], V ) × X → X is given by (45) (geometric
group action). Note that the velocity field ν ∈ L 2 ([0, 1], V ) in the objective also generates,
through (10), the diffeomorphism φν0,t : Ω → Ω in the ODE constraint. Furthermore, I(1, · ) ∈
X in the objective functional depends on the optimization variable ζ ∈ X through the ODE
constraint. Finally, as in subsection 6.1.1 for LDDMM, one can re-phrased (88) as a PDE
constrained problem:

h
i
2
2

inf
γkνk
+
τ
kζk
+
L
f
(1,
·
)

X
L 2 ([0,1],V )

2 ([0,1],V )


ν∈Lζ∈X
D
E
(89)
∂
f
(t,
x)
+
∇
f
(t,
x),
ν(t,
x)
= ζ(t, x) (t, x) ∈ [0, 1] × Ω,

t
x


Rn


f (0, x) = I(x) x ∈ Ω.
11.6. Relation to optimal transport. Optimal transport is a well developed mathematical theory that defines a family of metrics between probability distributions. These metrics
measure the mass of an optimal displacement according to a base cost defined on the space
supporting the distributions. Wasserstein distance refers to using L p -norm as base costs. It
is a deep mathematical theory [61] with many applications. Applications to imaging is nicely
surveyed in [47], the focus here is on its usage in diffeomorphic image registration.
Optimal transport metrics are suitable for shape registration since they are sensitive to
spatial displacements of the shape. This framework is [20, 37] extended to handle measures

36

C. CHEN AND O. ÖKTEM

of different total mass (unbalanced optimal transport), which makes it applicable to a wider
range of image registration problems. Another aspect associated with using optimal transport
distances is computability. Without approximations, it is computationally feasible only for
1D images. This situation has however changed during the last few years, much due to using
entropic approximation schemes resulting in Sinkhorn-type iterations [23]. Such an approach
is developed in [36] for indirect image registration where the authors used Sinkhorn-type
iterations to compute the proximal operator of the transport problem for large problems,
such as indirect image registration in tomography. This scheme also applies to unbalanced
optimal transport in [29] to quantify the similarity in diffeomorphic registration of imaging
data. Finally, [40] combines metamorphosis with optimal transport for image registration.
The numerical method is based on adding a transport penalty to the objective in (89), which
is the PDE constrained formulation of metamorphosis.
Acknowledgments. The authors would like to thank Barbara Gris and Oliver Verdier
for valuable comments and remarks. We especially thank Barbara Gris for contributions in
stating and proving many of the results in section 7.
REFERENCES
[1] J. Adler and O. Öktem, Learned primal-dual reconstruction, E-print math.OC 1707.06474, ArXiv
e-prints, 2017. Accepted for publication in IEEE Transactions on Medical Imaging.
[2] L. Ambrosio, Transport equation and Cauchy problem for BV vector fields, Inventiones Mathematicae,
158 (2004), pp. 227—260.
[3] L. Ambrosio and E. De Giorgi, New functionals in the calculus of variations, Atti della Accademia
Nazionale dei Lincei: Rendiconti Lincei - Matematica e Applicazioni, 82:8 (1988), pp. 199—210.
[4] N. Aronszajn, Theory of reproducing kernels, Transactions of the American Mathematical Society, 68
(1950), pp. 337–404.
[5] V. Arsigny, O. Commowick, N. Ayache, and X. Pennec, A fast and log-Euclidean polyaffine framework for locally linear registration, Journal of Mathematical Imaging and Vision, 33 (2009), pp. 222—
238.
[6] V. Arsigny, X. Pennec, and N. Ayache, Polyrigid and polyaffine transformations: a novel geometrical
tool to deal with non-rigid deformations–application to the registration of histological slices, Medical
image analysis, 9 (2005), pp. 507—523.
[7] A. I. Awad and M. Hassaballah, eds., Image Feature Detectors and Descriptors: Foundations and
Applications, vol. 630 of Studies in Computational Intelligence, Springer-Verlag, 2016.
[8] A. Banyaga, The Structure of Classical Diffeomorphism Groups, Kluwer, 1997.
[9] M. Bauer, S. Joshi, and K. Modin, Diffeomorphic density matching by optimal information transport,
SIAM Journal on Imaging Sciences, 8 (2015), pp. 1718–1751.
[10] F. Becker, S. Petra, and C. Schnörr, Optical flow, in Handbook of Mathematical Methods in
Imaging, O. Scherzer, ed., Springer-Verlag, 2nd ed., 2015, pp. 1945–2004.
[11] F. M. Beg, M. I. Miller, A. Trouvé, and L. Younes, Computing large deformation metric mappings
via geodesic flow of diffeomorphisms, International Journal of Computer Vision, 61 (2005), pp. 139—
157.
[12] A. Berlinet and C. Thomas-Agnan, Reproducing Kernel Hilbert Spaces in Probability and Statisties,
Springer-Verlag, 2004.
[13] E. Bladt, D. M. Pelt, S. Bals, and K. J. Batenburg, Electron tomography based on highly limited
data using a neural network reconstruction technique, Ultramicroscopy, 158 (2015), pp. 81–88.
[14] F. Bonin-Font, A. Ortiz, and G. Oliver, Visual navigation for mobile robots: A survey, Journal of
Intelligent and Robotic Systems, 53 (2008), pp. 263–296.
[15] M. Bruveris, Notes on Riemannian geometry on manifolds of maps.
Lecture notes for

INDIRECT IMAGE REGISTRATION WITH LARGE DIFFEOMORPHIC DEFORMATIONS

[16]

[17]
[18]

[19]
[20]
[21]
[22]

[23]
[24]
[25]

[26]

[27]
[28]
[29]
[30]
[31]
[32]
[33]
[34]

[35]
[36]

37

IMS Summer School on Mathematics of Shapes, 4–15 July 2016. Available online from
http://www2.ims.nus.edu.sg/Programs/016shape/files/martins.pdf, 2016.
M. Bruveris and D. D. Holm, Geometry of image registration: The diffeomorphism group and momentum maps, in Geometry, Mechanics, and Dynamics: The Legacy of Jerry Marsden, C. D. E., D. D.
Holm, G. Patrick, and T. Ratiu, eds., vol. 73 of Fields Institute Communications, Springer-Verlag,
2015, pp. 19–56.
M. Bruveris and F.-X. Vialard, On completeness of groups of diffeomorphisms, Journal of the European Mathematical Society, 19 (2017), pp. 1507–1544.
M. Burger and S. Osher, A guide to the TV zoo, in Level Set and PDE Based Reconstruction Methods
in Imaging, M. Burger and S. Osher, eds., vol. 2090 of Lecture Notes in Mathematics, Springer-Verlag,
2013, pp. 1–70.
A. Chambolle and T. Pock, A first-order primal-dual algorithm for convex problems with applications
to imaging, Journal of Mathematical Imaging and Vision, 40 (2011), pp. 120–145.
L. Chizat, B. Schmitzer, G. Peyré, and F.-X. Vialard, An interpolating distance between optimal
transport and Fisher-Rao, Foundations of Computational Mathematics, (2016).
G. E. Christensen, R. D. Rabbitt, and M. I. Miller, Deformable template model using large deformation kinematics, IEEE Transactions on Image Processing, 5 (1996), pp. 1435–1447.
M. Colombo, Flows of Non-smooth Vector Fields and Degenerate Elliptic Equations with Applications
to the Vlasov-Poisson and Semigeostrophic Systems, Publications of the Scuola Normale Superiore,
Springer-Verlag, 2017.
M. Cuturi, Sinkhorn distances: Lightspeed computation of optimal transportation, in Proceedings of
NIPS, vol. 26, 2013, pp. 2292—2300.
T. D’Arcy, On Growth and Form, Cambridge University Press, New York, 1945.
S. Dawn, V. Saxena, and B. Sharma, Remote sensing image registration techniques: A survey, in
Image and Signal Processing. Proceedings of the 4th International Conference, ICISP 2010, TroisRivières, QC, Canada, June 30-July 2, 2010, A. Elmoataz, O. Lezoray, F. Nouboud, D. Mammass,
and J. Meunier, eds., vol. 6134 of Lecture Notes in Computer Science, 2010, pp. 103–112.
E. De Giorgi and L. Ambrosio, New functionals in calculus of variations, in Nonsmooth Optimization
and Related Topics, F. H. Clarke, V. F. Dem’yanov, and F. Giannessi, eds., vol. 43 of Part of the Ettore
Majorana International Science Series book series (EMISS), Springer-Verlag, 1989, ch. 4, pp. 49–59.
C. Demant, B. Streicher-Abel, and C. Garnica, Industrial Image Processing: Visual Quality Control
in Manufacturing, Springer-Verlag, 2013.
R. J. DiPerna and P.-L. Lions, Ordinary differential equations, transport theory and Sobolev spaces,
Inventiones Mathematicae, 98 (1989), pp. 511—547.
J. Feydy, B. Charlier, F.-X. Vialard, and G. Peyré, Optimal transport for diffeomorphic registration, E-print math.NA 1706.05218, ArXiv e-prints, 2017.
M. Grasmair, Generalized Bregman distances and convergence rates for non-convex regularization methods, Inverse Problems, 26 (2010), p. 115014.
U. Grenander and M. I. Miller, Pattern Theory. From Representation to Inference, Oxford University
Press, 2007.
U. Grenander, A. Srivastava, and A. Saini, A pattern-theoretic characterization of biological growth,
IEEE Transactions on Medical Imaging, 26 (2007), pp. 648—659.
B. Gris, Modular approach on shape spaces, sub-Riemannian geometry and computational anatomy, PhD
thesis, École Normale Supérieure de Paris-Saclay, 2016.
Y. Hong, S. Joshi, M. Sanchez, M. Styner, and M. Niethammer, Metamorphic geodesic regression,
in Medical Image Computing and Computer-Assisted Intervention – MICCAI 2012 15th International
Conference, Nice, France, October 1-5, 2012, Proceedings, Part II, N. Ayache, P. Delingette, H. Golland, and K. Mori, eds., vol. 7512 of Lecture Notes in Computer Science, Springer-Verlag, 2012,
pp. 197–205.
Q. Huynh-Thu and M. Ghanbari, Scope of validity of PSNR in image/video quality assessment, Electronics Letters, 44 (2008), pp. 800–801.
J. Karlsson and A. Ringh, Generalized sinkhorn iterations for regularizing inverse problems using
optimal mass transport, E-print math.OC 1612.02273, ArXiv e-prints, 2016. Accepted for publication
in SIAM Journal on Imaging Sciences.

38

C. CHEN AND O. ÖKTEM

[37] M. Liero, A. Mielke, and G. Savaré, Optimal entropy-transport problems and a new HellingerKantorovich distance between positive measures, E-print math.OC 1508.07941, ArXiv e-prints, 2016.
[38] T. Lindeberg, Generalized axiomatic scale-space theory, Advances in Imaging and Electron Physics, 178
(2013), pp. 1–96.
[39] T. Lindeberg, Image matching using generalized scale-space interest points, Journal of Mathematical
Imaging and Vision, 52 (2015), pp. 3–36.
[40] J. Maas, M. Rumpf, C. Schönlieb, and S. Simon, A generalized model for optimal transport of
images including dissipation and density modulation, ESAIM: Mathematical Modelling and Numerical
Analysis, 49 (2015), pp. 1745–1769.
[41] K. McLeod, M. Sermesant, P. Beerbaum, and X. Pennec, Spatio-temporal tensor decomposition
of a polyaffine motion model for a better analysis of pathological left ventricular dynamics, IEEE
Transactions on Medical Imaging, 34 (2015), pp. 1562—1575.
[42] M. Micheli and J. A. Glaunès, Matrix-valued kernels for shape deformation analysis, Geometry, Imaging, and Computing, 1 (2013), pp. 57–139.
[43] M. I. Miller, A. Trouvé, and L. Younes, Hamiltonian systems and optimal control in computational
anatomy: 100 years since D’Arcy Thompson, Annual Review of Biomedical Engineering, 17 (2015),
pp. 447–509.
[44] K. P. Murphy, Machine Learning: A Probabilistic Perspective, MIT Press, Cambridge, MA, USA, 2012.
[45] F. Natterer and F. Wübbeling, Mathematical Methods in Image Reconstruction, Mathematical Modeling and Computation, Society for Industrial and Applied Mathematics, 2001, http://epubs.siam.
org/doi/book/10.1137/1.9780898718324.
[46] O. Öktem, C. Chen, N. O. Domaniç, P. Ravikumar, and C. Bajaj, Shape based image reconstruction
using linearised deformations, Inverse Problems, 33 (2017), p. 035004.
[47] N. Papadakis, Optimal transport for image processing, hal id: tel-01246096, Université de Bordeaux,
2016. Avaialble from https://hal.archives-ouvertes.fr/tel-01246096v6.
[48] N. Portman and E. R. Vrscay, Existence and uniqueness of solutions to the grid macroscopic growth
equation, Applied Mathematics and Computation, 217 (2011), pp. 8318—8327.
[49] C. L. Richardson and L. Younes, Metamorphosis of images in reproducing kernel Hilbert spaces,
Advances in Computational Mathematics, 42 (2016), pp. 573—603.
[50] W. Ring and B. Wirth, Optimization methods on Riemannian manifolds and their application to shape
space, SIAM Journal of Optimization, 22 (2012), pp. 596–627.
[51] E. Sariyanidi, H. Güneş, and A. Cavallaro, Automatic analysis of facial affect: A survey of registration, representation, and recognition, IEEE Transactions on Pattern Analysis and Machine Intelligence, 37 (2015), pp. 1113–1133.
[52] O. Scherzer, M. Grasmair, H. Grossauer, M. Haltmeier, and F. Lenzen, Variational Methods in
Imaging, vol. 167 of Applied Mathematical Sciences, Springer-Verlag, New York, 2009.
[53] J. A. Schnabel, M. P. Heinrich, B. W. Papieżc, and M. J. Brady, Advances and challenges
in deformable image registration: From image fusion to complex motion modelling, Medical Image
Analysis, 33 (2016), pp. 145–148.
[54] B. Scholkopf and A. J. Smola, Learning with Kernels: Support Vector Machines, Regularization,
Optimization, and Beyond, MIT Press, Cambridge, MA, USA, 2001.
[55] L. Schwartz, Sous-espaces Hilbertiens d’espaces vectoriels topologiques et noyaux associés (noyaux reproduisants), Journal d’Analyse Mathématique, 13 (1964), pp. 115–256.
[56] C. Seiler, X. Pennec, and M. Reyes, Capturing the multi-scale anatomical shape variability with
polyaffine transformation trees, Medical image analysis, 16 (2012), pp. 1371—1384.
[57] A. Sotiras, C. Davatzikos, and N. Paragios, Deformable medical image registration: A survey, IEEE
Transactions on Medical Imaging, 32 (2013), pp. 1153–1190.
[58] A. Trouvé and L. Younes, Shape spaces, in Handbook of Mathematical Methods in Imaging, S. Otmar,
ed., Springer-Verlag, 2011, ch. 30, pp. 1309–1362.
[59] F.-X. Vialard, Hamiltonian Approach to Shape Spaces in a Diffeomorphic Framework: From the Discontinuous Image Matching Problem to a Stochastic Growth Model, PhD thesis, École normale supérieure
de Cachan, Centre de Mathématiques et de Leurs Applications, 61, avenue du Président Wilson,
Cachan, 2009.
[60] M. A. Viergever, J. B. Maintz, S. Klein, K. Murphy, M. Staring, and J. P. Pluim, A survey of

INDIRECT IMAGE REGISTRATION WITH LARGE DIFFEOMORPHIC DEFORMATIONS

39

medical image registration, Medical Image Analysis, 33 (2016), pp. 140–144.
[61] C. Villani, Optimal Transport: Old and New, Grundlehren der mathematischen Wissenschaften,
Springer-Verlag, 2009.
[62] Z. Wang, A. C. Bovik, H. R. Sheikh, and E. P. Simoncelli, Image quality assessment: from error
visibility to structural similarity, IEEE Transactions on Image Processing, 13 (2004), pp. 600–612.
[63] X. Yang, R. Kwitt, and M. Niethammer, Fast predictive image registration, in Deep Learning and
Data Labeling for Medical Applications: Proceedings of the First International Workshop (LABELS
2016) and Second International Workshop (DLMIA 2016) Held in Conjunction with MICCAI 2016,
G. Carneiro, D. Mateus, L. Peter, A. Bradley, J. M. R. S. Tavares, V. Belagiannis, J. P. Papa, J. C.
Nascimento, M. Loog, Z. Lu, J. S. Cardoso, and J. Cornebise, eds., vol. 10008 of Lecture Notes in
Computer Science, 2016, pp. 48–47.
[64] X. Yang, R. Kwitt, M. Styner, and M. Niethammer, Quicksilver: Fast predictive image registration
– a deep learning approach, NeuroImage, 158 (2017), pp. 378–396.
[65] L. Younes, Shapes and Diffeomorphisms, vol. 171 of Applied Mathematical Sciences, Springer-Verlag,
2010.
[66] L. Younes, Constrained diffeomorphic shape evolution, Foundations of Computational Mathematics, 12
(2012), pp. 295—325.
[67] B. Zitová and J. Flusser, Image registration methods: a survey, Image and Vision Computing, 21
(2003), pp. 977–1000.

40

C. CHEN AND O. ÖKTEM

Template

t = 0.25

t = 0.5

t = 0.75

t = 1.0

Target

Data at 0◦

FBP

TV

Figure 1: Test suite 1: Single-object indirect registration. The template (top left) is deformed
by a diffeomorphism to match the target (rightmost, second row). Images labeled with time
t show the final flow of diffeomorphic deformations obtained from indirect registration with
t = 1.0 denoting the final (indirectly) registered template. Data is highly noisy and undersampled, bottom leftmost image shows data at angle 0◦ (blue smooth curve is noise-free data,
red jagged curve is noisy data). The FBP and TV reconstructions are for comparison.

INDIRECT IMAGE REGISTRATION WITH LARGE DIFFEOMORPHIC DEFORMATIONS

Template

t = 0.25

t = 0.5

t = 0.75

t = 1.0

Target

Data at 0◦

FBP

TV

41

Figure 2: Test suite 2: Multi-object indirect registration. The template (top left) is a target
consisting of six separately triangle-like objects. Images labeled with time t show the final
flow of diffeomorphic deformations obtained from indirect registration with t = 1.0 denoting
the final (indirectly) registered template. Data is highly noisy and under-sampled, bottom
leftmost image shows data at angle 0◦ (blue smooth curve is noise-free data, red jagged curve
is noisy data). FBP and TV reconstructions are shown for comparison.

42

C. CHEN AND O. ÖKTEM

γ = 10−7 , σ = 2.0

γ = 10−5 , σ = 2.0

γ = 10−3 , σ = 2.0

γ = 10−7 , σ = 3.0

γ = 10−5 , σ = 3.0

γ = 10−3 , σ = 3.0

γ = 10−7 , σ = 8.0

γ = 10.0, σ = 3.0

γ = 10.0, σ = 8.0

Template

Target

Figure 3: Test suite 3: Sensitivity w.r.t. choice of regularization parameters γ and σ. Images
show final registered result (i.e., t = 1) for different values of these parameters. See Table 1
for a more quantitative comparison that includes a wider range of parameter values.

INDIRECT IMAGE REGISTRATION WITH LARGE DIFFEOMORPHIC DEFORMATIONS

Template

t = 0.25

t = 0.50

t = 0.75

t = 1.0

Target

Template

t = 0.25

t = 0.50

t = 0.75

t = 1.0

Target

43

Figure 4: Test suite 4: Topology of the template. First two rows is the case when the
template lacks one object as compared to target, the two following rows is the case when
template has one extra object. Images labeled with time t show the final flow of diffeomorphic
deformations obtained from indirect registration with t = 1.0 denoting the final (indirectly)
registered template.

