Next Article in Journal
Coupled SU(3)-Structures and Supersymmetry
Next Article in Special Issue
A Framework for Symmetric Part Detection in Cluttered Scenes
Previous Article in Journal
Group Theory of Wannier Functions Providing the Basis for a Deeper Understanding of Magnetism and Superconductivity
Previous Article in Special Issue
Unsupervised Object Modeling and Segmentation with Symmetry Detection for Human Activity Recognition
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Reduction by Lie Group Symmetries in Diffeomorphic Image Registration and Deformation Modelling

1
Department of Computer Science, University of Copenhagen, Universitetsparken 1, DK-2100 Copenhagen E, Denmark
2
Department of Mathematics, Imperial College, SW7 2AZ London, UK
*
Authors to whom correspondence should be addressed.
Symmetry 2015, 7(2), 599-624; https://doi.org/10.3390/sym7020599
Submission received: 30 November 2014 / Revised: 15 April 2015 / Accepted: 27 April 2015 / Published: 7 May 2015
(This article belongs to the Special Issue Symmetry: Theory and Applications in Vision)

Abstract

: We survey the role of reduction by symmetry in the large deformation diffeomorphic metric mapping framework for registration of a variety of data types (landmarks, curves, surfaces, images and higher-order derivative data). Particle relabelling symmetry allows the equations of motion to be reduced to the Lie algebra allowing the equations to be written purely in terms of the Eulerian velocity field. As a second use of symmetry, the infinite dimensional problem of finding correspondences between objects can be reduced for a range of concrete data types, resulting in compact representations of shape and spatial structure. Using reduction by symmetry, we describe these models in a common theoretical framework that draws on links between the registration problem and geometric mechanics. We outline these constructions and further cases where reduction by symmetry promises new approaches to the registration of complex data types.

1. Introduction

Registration, the task of establishing correspondences between multiple instances of objects, such as images, landmarks, curves and surfaces, plays a fundamental role in a range of computer vision applications, including shape modelling [1], motion compensation and optical flow [2], remote sensing [3] and medical imaging [4]. In the subfield of computational anatomy, establishing inter-subject correspondences between organs allows the statistical study of organ shape and shape variability [5]. Examples of the fundamental role of registration include quantifying developing Alzheimer’s disease by establishing correspondences between brain tissue at different stages of the disease [6]; measuring the effect of chronic obstructive pulmonary disease on lung tissue after removing variability caused by the respiratory process [7]; and correlating the shape of the hippocampus to schizophrenia after inter-subject registration [8].

In this paper, we survey the role of reduction by symmetry in diffeomorphic registration and deformation modelling, linking symmetry as seen from the field of geometric mechanics with the image registration problem. All of our calculations will be formal and void of functional analytic detail, although citations will be used when available. We focus on large deformations modelled in subgroups of the group of diffeomorphic mappings on the spatial domain in the context of large deformation diffeomorphic metric mapping (LDDMM) [1,911]. Connections with geometric mechanics [12,13] have highlighted the role of symmetry, and properties that were previously known to be connected with the specific data types have been described in a common theoretical framework [14]. We wish to describe these connections in a form that highlights the role of symmetry and points towards future applications of the ideas.

1.1. Symmetry and Information

One of the main reasons symmetry is useful in data analysis and numerics is its ability to reduce the complexity of information that represents data. Lower information complexity can lead to more stable statistical analysis and the reduced need of computational resources.

As a toy example, consider a spinning top. Upon choosing a reference configuration, the orientation of the top is given by a rotation matrix, i.e., an element R ∈ SO(3) (see Figure 1). To describe the direction of the tip of the top, it suffices to provide the orientation matrix R. However, R is contained in SO(3), a three-dimensional space, while the space of possible directions is the two-sphere, S2, which is only of dimension two. Therefore, providing the full matrix R is an over-representation of the tip direction. It suffices to solely provide the vector R · k ∈ S2 where k = (0, 0, 1) is the direction of the tip in a reference configuration. Note that if R ˜ k = k, then R k = R R ˜ k. Therefore, given only the direction k = R · k, we can only reconstruct R up to an element R ˜, which preserves k. The subgroup of rotations that preserve k can be identified with SO(2). Specifically, this identification comes from perceiving a rotation about k as a rotation of the plane, which is perpendicular to k. This insight allows us to express the space of directions S2 as a homogeneous space S2 SO(3)/SO(2). In terms of information, we can cartoonishly express this by:

orientation = direction of tip + orientation around the tip

This picture is typical for many group quotients. Generally speaking, if X is a manifold and G acts freely and properly on X, then:

dim ( X ) = dim ( X / G ) + dim ( G )

When X is infinite dimensional, this formula is less insightful. However, X/G is smaller than X in that there exists a surjective map from X to X/G with non-trivial level sets. Reduction by symmetry can be implemented when a problem posed on X has G symmetry and can be rewritten as a problem posed on X/G. As X/G is smaller than X, reduction by symmetry can yield more stable subsequent statistical analysis of observed data and more tractable algorithms, as will be shown later in the article. This reduction is particularly dramatic when dim(X) = and dim(X/G)<∞.

1.2. Symmetry in Registration

Registration of objects contained in a spatial domain, e.g., the volume to be imaged by a scanner, can be formulated as the search for a deformation that transforms both domain and objects to establish an inter-object match. The data available when solving a registration problem is generally insufficient for determining the displacement of every point of the domain. This is the case when images to be matched have areas of constant intensity and no derivative information can guide the registration. For example, the “best” deformation for matching the two discs in Figure 2 is ambiguous, except at the boundary of the discs, where the images are non-constant. Similarly, when 3D shapes are matched based on the similarity of their surfaces, the deformation of the interior cannot be derived from the available information. In these cases, the deformation model is over-complete, and a range of deformations can provide equally adequate matches for the data. The registration problem or the registration cost-function is thus symmetric with respect to the subset of transformations just described. When the deformation model is a Lie group, the deformations for which the registration is symmetric form a subgroup. The quotient by this subgroup of symmetries of the registration cost-function can provide vastly more compact representations. This situation arises in several cases with the LDDMM framework: when registering images, only displacements orthogonal to the level lines of an image are needed, and when registering shapes, the information left in the quotient is supported on the shape surface only.

1.3. Content and Outline

Although a degree of comfort with differential geometry will be assumed, it is the aim of this paper to make the role of symmetry in registration and deformation modelling clear to the non-expert in geometric mechanics and symmetry groups in image registration. We begin the paper by presenting the background for the registration problem and the large deformation approach before outlining some necessary notions from differential geometry. For more information on the Riemannian geometry behind the LDDMM approach to image registration, we refer the reader to [5]. We continue by describing how reduction by symmetry leads to an Eulerian formulation of the equations of motion when reducing to the Lie algebra. The symmetry of dissimilarity measures allows additional reductions, and we use isotropy subgroups to reduce the complexity of the registration problem further. Lastly, we survey the effect of symmetry in a range of concrete registration problems. The paper ends with concluding remarks.

2. Registration

The registration problem consists of finding correspondences between objects that are typically point sets (landmarks), curves, surfaces, images or more complicated spatially-dependent data, such as diffusion weighted images (DWI). The problem can be approached by letting M be a spatial domain containing the objects to be registered. M can be a compact finite dimensional differentiable manifold without boundary or ℝd itself with d = 2, 3. It is common to consider manifolds with boundaries, as well. In such cases, care must be taken with regards to boundary conditions. For example, vector-fields must be tangential to the boundary. Here, we consider only manifolds without a boundary.

A map φ : M → M can deform or warp the domain by mapping each xM to φ(x). The deformation encoded in the warp will apply to the objects in M, as well as the domain itself. For example, if the objects to be registered consist of point sets {x1,…, xN}, xiM, the set will be mapped to {φ(x1),…, φ(xN)}. For surfaces SM, φ similarly results in the warped surface φ(S). Because those operations are associative, the mapping φ acts on {xi} or S, and we write φ · {xi} and φ · S for the warped objects. An image is a function I : M → ℝ, and φ acts on I as well, in this case by composition with its inverse φ · I = I ○ φ1; see Figure 3. For this, φ must be is invertible, and commonly, we restrict to the set of invertible and differentiable mappings Diff(M). For various other types of data objects, the action of a warp on objects can be defined in a manner similar to that of point sets, surfaces and images. This fact relates a range registration problems to the common case of finding appropriate warps φ, which bring the objects into correspondence. Different shape instances can be realized by letting warps act on a base shape, and a class of shape models can thereby be obtained by using deformations as shape representations [1].

The search for appropriate warps can be formulated in a variational formulation with an energy:

E ( φ ) = R ( φ ) + F ( φ )
where F is a dissimilarity measure of the difference between the deformed objects and R is a regularization term that penalizes unwanted properties of φ, such as spatial irregularity. If two objects o1 and o2 are to be matched, F can take the form F (φ · o0, o1) using the action of φ on o0; for image matching, an often used dissimilarity measure is the L2-difference or sum of square differences (SSD) having the form F ( φ I 0 , I 1 ) = M | I 0 φ 1 ( x ) I 1 ( x ) | 2 d x.

The regularization term can take various forms often modelling physical properties, such as elasticity [15], and derivatives of φ are often penalized to ensure smoothness. For some choices of R, existence and analytical properties of minimizers of Equation (1) have been derived [16]; however, in general, it is difficult to ensure that solutions are diffeomorphic by penalizing φ in itself. The free-form-deformation (FFD; [17]) and related approaches model the deformation by a displacement vector field u on M = ℝd, so that φ(x) = x + u(x). Smoothness is here ensured by the choice of basis functions, e.g., B-splines, or by applying a regularization term on u. Smooth and invertible mappings can be obtained by integrating flows [9,11] to obtain one-parameter families or paths of mappings φt, t ∈ [0, 1]. The warp φ0 at t = 0 is here the identity mapping id ∈ Diff(M), and the dissimilarity is measured at the endpoint φ1. The time evolution of φt can be described by the differential equation:

d d t φ t ( x ) = u t ( φ t ( x ) )
with the flow field ut being a vector field on M. Numerically, the map φ can be represented by how it maps a finite set of points, and a numerical scheme might simply implement Euler-integration on each point, i.e., “xi+1xi + h · ut(xi)” with time-step size 0 < h ≪ 1. A relaxation of this idea is now a standard method in optical flows [18]. The space of flow fields is denoted by V. In the LDDMM [1] framework, the regularization is applied to the flow field and integrated over time giving the energy:
E ( φ ) = 0 1 | | u t | | V 2 d t + F ( φ 1 )

Here, the time-dependent diffeomorphism φt is related to ut through Equation (2). If the norm ‖·V that measures the irregularity of ut is sufficiently strong (e.g., Hk with k > d 2 + 1), then φt will be a diffeomorphism for all t [19]. This approach thus gives a direct way of enforcing properties of the generated warp: instead of regularizing φ directly, the analysis is lifted to a normed space V that is much easier control. The energy E in Equation (3) has the same minimizers as the geometric formulation of LDDMM used in the next section.

Direct approaches to solving the optimization problem in Equation (3) must handle the fact that the problem of finding a warp is now expanded to that of finding a time-dependent family of warps. This is a huge increase in dimensionality. This formulation of registration is thus very difficult to represent numerically and to optimize and analyse statistically. For several data types, it has been shown how optimal paths for Equation (3) have specific properties that reduce the dimensionality of the problem, making practical solutions feasible. In the next section, we outline the geometric framework that is needed when we, in the later sections, use reduction by symmetry to describe these data-dependent results in a common theoretical framework.

3. Notions from Differential Geometry

In this section, we will introduce a number of notions from differential geometry in a fairly informal manner. We will use conventions from [20,21] where a more rigorous understanding of differential geometry can be found.

We will assume that the reader has at least an intuitive picture of the notion of a smooth manifold M. For the purpose of this paper, M will either be assumed to be compact without a boundary or ℝn. The tangent bundle of M is the space of velocity vectors tangential to M. Notationally, the tangent bundle of M, denoted TM, is the set of pairs (x, v) where xM and v is a vector tangential to M at the point x (see Figure 4). A vector-field is a continuous map u : M → TM, such that u(x) ∈ TM is a vector above x for all xM. We will denote the space of vector-fields by X ( M ).

Given a vector-field u X ( M ), we may consider the initial value problem:

{ x ( 0 ) = x 0 d x d t = u ( x ( t ) )
for t ∈ [0, 1]. Given an initial condition x0, the point x1 = x(1) given by solving this initial value problem is uniquely determined if it exists. Under many circumstances (e.g., if M is compact or if M = ℝd and ‖u(x)‖ grows sub-linearly), an x1 exists for each x0, and there is a continuous invertible map Φ t u : x 0 M x 1 M, which we call the flow of u. Given a time-dependent vector-field, u t X ( M ) for t ∈ [0, 1], we can consider the initial value problem with d x d t = u t ( x ( t ) ). This will yield a flow map, Φ t 0 , t 1 u which is the flow from time t = t0 to t = t1. If ut is smooth, the flow map will be smooth, as well, in particular a diffeomorphism. We denote the set of diffeomorphisms by Diff(M).

Conversely, let φt ∈ Diff(M) be a time-dependent diffeomorphism. For any xM, we observe that φt(x) is a curve in M. If this curve is differentiable we may consider its time-derivative, d φ t d t ( x ) T M, which is a vector above the point φt(x) ∈ M. From these observations, it follows that d φ t d t [ φ t 1 ( x ) ] is a vector above x. Therefore, the map ut : M → TM, given by u t ( x ) : = ( d d t φ t ) [ φ t 1 ( x ) ], is a vector-field called the Eulerian velocity field of φt.

As will be described shortly, the Eulerian velocity field contains less data than d φ t d t. This reduction in data can be viewed from the perspective of symmetry. Given any ψ ∈ Diff(M), the curve φt can be transformed to the curve φt ○ ψ. We observe that:

u t : = d φ t d t φ t 1 = d φ t d t ψ ψ 1 φ t 1 = ( d d t ( φ t ψ ) ) ( φ t ψ ) 1

Thus, φt and φt ○ ψ both have the same Eulerian velocity fields. In other words, the Eulerian velocity field, ut, is invariant under particle relabellings. More precisely, we may view Diff(M) as a manifold in its own right, and view d φ t d t as a vector in the infinite-dimensional tangent bundle T Diff(M) above the “point” φt ∈ Diff(M). Thus, the vector d φ t d t contains both velocities and a base diffeomorphism φt. Given ut and φt, we can reconstruct d φ t d t via d φ t d t = u t φ t . As has been shown, we can also construct ut from d φ t d t by its own definition. However, we cannot reconstruct d φ t d t from ut, which is why ut contains less data.

Finally, we will denote some linear operators on the space of vector-fields. Let Φ ∈ Diff(M), and let u X ( M ). The push-forward of u by Φ is the vector-field given by:

Φ * u ( x ) : = D Φ | Φ 1 ( x ) u ( Φ 1 ( x ) )

In local coordinates (x1,…, xn), this looks like:

[ Φ * u ] i ( x ) = j = 1 n ( j | Φ 1 ( x ) Φ i ) u j ( Φ 1 ( x ) )

By inspection, we see that Φ is a linear operator on X ( M ). One can view Φu as “u in a new coordinate system”, because any geometric property of u is also inherited by Φu. For example, if S is invariant under u, then Φ(S) is invariant under Φu. We define the pull-back by Φu := (Φ1)u. Note that Φ(u)) = u.

As Φ is a linear operator, a well-defined operator exists, which is dual to Φ. Let X ( M ) * denote the dual space to X ( M ), i.e., the set of linear maps X ( M ) , which are continuous with respect to a chosen vector-space topology on X ( M ). Given m X ( M ) *, we define Φ m X ( M ) * by the identity:

Φ * m , Φ * u : = m , u
for all u X ( M ), where 〈m, u〉 denotes the evaluation of m on v. We can define Φm := (Φ1)m, which yields the identity:
Φ * m , u = m , Φ * u

In local coordinates, we may represent m as a one-form density, given by mi(x)dxi⊗(dx1∧…∧dxn) with components m1(x),…, mn(x). In this local coordinate description, the i-th component of the push-forward looks like:

[ Φ m ] i ( x ) : = det ( D Φ 1 ) | x i Φ j | Φ 1 ( x ) m j ( Φ 1 ( x ) )

Finally, we define the Lie derivative. Let w X ( M ). The Lie derivative operator, with respect to w, is the linear operator £ w : X ( M ) X ( M ) defined by:

£ w [ u ] = d d ϵ | ϵ = 0 ( Φ ϵ w ) * u

Again, as £w is a linear operator on X ( M ), we can define a dual operator on X ( M ) *. If m X ( M ) *, we can define £ w [ m ] X ( M ) £w[m] ∈ X(M) by the equation:

£ w [ m ] , u + m , £ w [ u ] = 0
for all u X ( M ).

We conclude the section with a table of notation for the reader’s convenience, see Table 1.

4. Reduction by Symmetry in LDDMM

In this section, we will present necessary conditions satisfied by local extremizers of the variational problem Equation (3). The resulting conditions will first involve the computation of a curve in Diff(M), as well as its time-derivative in T Diff(M). We then invoke a Diff(M) symmetry of the problem to reduce this computation to a computation on X ( M ) instead of T Diff(M) ≅ Diff(M) × X ( M ). Secondly, we describe how the symmetry of the dissimilarity measure allows further reductions.

4.1. Reduction to the Lie Algebra

The variational formulation Equation (3) of LDDMM is equivalent to minimizing the energy:

E = d ( i d , φ ) + F ( φ )
where d : Diff(M) × Diff(M) ℝ is a Riemannian distance metric on Diff(M) induced by the norm ‖·‖V, id is the identity diffeomorphism, and F : Diff(M) ℝ is a dissimilarity measure, i.e., a function measuring the disparity between the deformed template and the target object.

Example 1. Given images I0, I1L2(M), we consider the dissimilarity measure:

F ( φ ) = | | ( I 0 φ 1 ) I 1 | | L 2 ( M ) 2

In this article, we will consider the metric on connected components of Diff(M) given by:

d ( φ 0 , φ 1 ) = inf u C 0 ( [ 0 , 1 ] X ( M ) ) Φ 0 , 1 v φ 0 = φ 1 ( 0 1 | | u t | | d t )
where ut denotes a one-parameter family of vector fields and ‖·‖ is a norm on X ( M ), the Lie algebra of Diff(M). The norm is generally assumed to be admissible, i.e., embedded in C 0 1 ( M , n + k ) for sufficiently large k, so that a constant C exists satisfying ‖u‖ ≥ C‖u1,∞ for all u ∈ X(M) ([1], Chapter 9). In the case where M ≠ ℝn, we can define ‖u1,∞ chart-wise by a partition of unity (e.g., see the construction of Hk norms on M in [22]) or intrinsically in terms of the Riemannian gradient. Both choices would yield identical topologies, and so, this choice has no significance as far as the article is concerned. From now on, we will overload notation and let X ( M ) denote the set of Ck vector fields with finite norm. For finite k, this makes X ( M ) a Banach space, and it breaks the Lie algebra structure. The consequences of this breakage will not be explored here, and we will continue to treat X ( M ) formally as a Lie algebra. We will later be using the space of homeomorphisms generated by X ( M ), which is a subspace of Ck-diffeomorphisms. Again, we will overload the notation and call this space Diff(M), even though this is usually reserved for smooth diffeomorphisms. In the case where M = ℝd, we assume that our norm is such that decay conditions at infinity for u X ( M ) arise naturally as a result of requiring the norm to be finite. If ‖·‖ is induced by an inner-product, the inner-product is formally a Riemannian metric on Diff(M) given by:
d φ t d t ϕ t 2 : = d φ t d t φ t 1 2
and d is the Riemannian distance with respect to this metric. The norm is often defined in terms of an operator P : X ( M ) X ( M ) as ‖u2 = 〈P [u], u〉, and the assumed admissibility implies that ( X ( M ) , | | · | | )has a reproducible kernel Hilbert space structure (RKHS; [1], Chapter 8). For example, we could consider M = ℝ, P = d x ( 1 x 2 ), and the vector-field u(x) = exp(−|x|) is mapped to the one-form density dxδ, where δ is the Dirac delta distribution (see [23]. In particular, in the case that M = ℝn, a matrix-valued kernel function K : ℝn ×nn×n exists satisfying the reproducing property 〈P [K(·, x)a], u〉 = aT u(x) for all x ∈ ℝn and a ∈ Rn (see [24]). We will denote RKHSs by V and the norms by ‖·‖V.

Given P, minimizers of E Equation (6) must satisfy:

{ m t = ( Φ t , 1 u t ) * m 1 = P [ u t ] t Φ t , 1 u = u t Φ t , 1 u m 1 , w = d d ϵ | ϵ = 0 F ( Φ ϵ w Φ 0 , 1 u ) , w X ( M ) .

That Equation (7) is a necessary condition satisfied by the minimizers of Equation (6) (under certain analytic assumptions) is stated and proven in [1] (Proposition 11.6) for the case of M = ℝn. The proof is based on the definition of the curve energy Equation (6) as a sum a Riemannian distance d and the dissimilarity measure F. Minimization of the Riemannian distance d yields a search for a geodesic in Diff(M), because geodesics are distance minimizing. The corresponding geodesic equations are given by the first two lines of Equation (7). The term F only penalizes the end-point of the geodesic, and the minimization condition manifests as the third line of Equation (7).

Issues regarding the well-posedness of Equation (7) are non-trivial, because P is merely injective, but not bijective, and so, there is no guarantee that P can be inverted on a given m t X ( M ) at each time in order to obtain a vector-field u t X ( M ). Fortunately, safety guards for well-posedness are known (e.g., [19], Theorem 1, or [25]).

Using Equation (7) for computational purposes is difficult because Diff(M) is a non-linear infinite dimensional space. Moreover, the dissimilarity measure F only comes into play at time t = 1, and the distance function is an integral over the vector-space X ( M ). It would be beneficial if we could rewrite the extremizers in terms of the Eulerian velocity field u and the flow at t = 1. In fact, this is often the case. One (formally) must take the time-derivative of the term “ ( Φ t , 1 u ) m ( 1 )” and apply Equation (5). Explicitly, this computation is performed as follows. Let w X ( M ), and observe:

t m t , w = d d t ( Φ t , 1 u ) * m ( 1 ) , w = m , d d t ( Φ t , 1 u ) * w = m , £ u t [ w ] = £ u t [ m ] , w

As w is arbitrary, we find tmt + £ut[mt] = 0. This allows us to reformulate Equation (7) as:

{ t m t + £ u t [ m t ] = 0 , m t = P [ u t ] t [ 0 , 1 ] d d ϵ | ϵ = 0 [ F ( Φ ϵ w Φ 0 , 1 u ) ] + P [ u 1 ] , w = 0 , w X ( M )

The advantage of this formulation is that the bulk of the computation occurs on the vector-space X ( M )rather than on the space T Diff(M). Registration algorithms based on Equation (8) differ from the algorithm proposed by Beg et al. in [26]. In [26], a gradient descent on the time-dependent Eulerian vector field ut is used to minimize the curve energy. The algorithm is posed on the velocity field ut instead of the momentum field mt as Equation (8) suggests. The momentum at time t is retrieved by a transport equation similar to the first equation in Equation (7). The evolution Equation (8) effectively allows one to search over the space of initial conditions, X ( M ) *, rather than over the larger space of curves, C1([0, 1]; X ( M )).

This reduction of the problem from T Diff(M) given in Equation (7) to the problem over space of vector-fields, X ( M )given in Equation (8) is the first instance of reduction by symmetry. In particular, this corresponds to the fact that the space of vector-fields X ( M )is identifiable as a quotient space:

X ( M ) T Diff ( M ) / Diff ( M )

Additionally, the map ( φ t , d φ t d t ) T Diff ( M ) d φ t d t φ t 1 X ( M )is the quotient projection.

4.2. Isotropy Subgroups

The reduction to dynamics on T Diff(M) to dynamics on X ( M )occurs primarily because the distance function is Diff(M) invariant. However, one cannot completely abandon Diff(M), because the solution requires one to compute the Time 1 flow, Φ 0 , 1 u. Fortunately, there is a second reduction, which allows us to avoid computing Φ 0 , 1 uin its entirety. This second reduction corresponds to the invariance properties of the dissimilarity measure F. Let GF ⊂ Diff(M) denote the set of diffeomorphisms, which leave F invariant, i.e.:

G F : = { ψ Diff ( M ) | F ( φ ψ ) = F ( φ ) , φ Diff ( M ) }

One can readily verify that GF is a subgroup of Diff(M), and so, we call GF the isotropy subgroup of F.

Having defined GF, we can now consider the homogeneous space Q = Diff(M)/GF, which is the quotient space induced by the action of the right composition of GF on Diff(M). This quotient space is “smaller” in the sense of the amount of data required to describe an element of it. In terms of maps, this can be seen by defining the map φ ∈ Diff(M) ↦q = [φ]/GFQ, where [φ]/GF denotes the equivalence class of φ. We call this mapping the quotient projection, because it sends Diff(M) to Q surjectively. While these notions are theoretically quite complicated, often they manifest less so in practice.

Example 2. In this example, we consider a simple aspect of the landmark matching problem. Let M ⊂ ℝn be the closure of some open set. Let x1, x2, y1, y2M with x1x2, and consider the dissimilarity measure:

F ( φ ) = | | φ ( x 1 ) y 1 | | 2 + | | φ ( x 2 ) y 2 | | 2

We see that:

G F { ψ Diff ( M ) | ψ ( x 1 ) = x 1 , ψ ( x 2 ) = x 2 }
and:
Q = Diff ( M ) / G F { ( z 1 , z 2 ) M × M | z 1 z 2 } = M × M Δ M × M
where ΔM×M denotes the diagonal of M × M. The quotient projection is φ ∈ Diff(M) ↦ (φ(x1), φ(x2)) ∈ M × M − ΔM×M

Note that Diff(M) is infinite dimensional, while Q is of dimension 2 dim(M). Two examples of diffeomorphisms contained in GF are shown in Figure 5.

Example 3. In this example, we consider the matching problem for greyscale images. Let I0, I1Hk(M) be images. There is a natural action of Diff(M) on Hk given by sending each image IHk(M) to I ○ φHk(M). We could consider the matching function F : Diff(M) given by:

F ( φ ) = | | ( I 0 φ 1 ) I 1 | | H k

This function measures the difference between a deformed version of I0 and a fixed image, I1. The isotropy subgroup GF is the group of images that preserve 00. Such a diffeomorphism would preserve each of the level lets of I0, but could permute the points within a given level set.

If one is able to understand Q, then one can use this insight to reformulate the dissimilarity measure F as a function on Q, rather than Diff(M). In particular, there exists a unique function FQ : Q → ℝ defined by the property FQ([φ]/GF ) = F(φ). Again, this is useful in the sense of data, as illustrated in the following example.

Example 4. Consider the dissimilarity measure F of Example 2. The function, FQ : Q →is:

F Q ( q 1 , q 2 ) = | | q 1 y 1 | | 2 + | | q 2 y 2 | | 2

Finally, note that Diff(M) acts upon Q by the left action:

[ φ ] G F Diff ( M ) / G F ψ Diff ( M ) [ ψ φ ] / G F Diff ( M ) / G F

Usually, we will simply write ψ · q for the action of ψ ∈ Diff(M) on a given qQ. This means that X ( M )acts upon Q infinitesimally, as it is the Lie algebra of Diff(M).

Example 5. Consider the setup of Example 2. Here, Q = M × M − ΔM×M, and the left action of Diff(M) is given by:

ψ ( q 1 , q 2 ) = ( ψ ( q 1 ) , ψ ( q 2 ) )
for ψ ∈ Diff(M) and q = (q1, q2) ∈ Q. The infinitesimal action of u X ( M )on Q is:
u ( q 1 , q 2 ) = ( u ( q 1 ) , u ( q 2 ) ) T q Q

These constructions allow us to rephrase the initial optimization problem using a reduced curve energy. Minimization of E is equivalent to minimization of:

E Q = 0 1 | | u t | | d t + F Q ( q ( 1 ) )
where q(1) is obtained by integrating the ODE, d q ( t ) d t = u t q ( t )with the initial condition q ( 0 ) = [ i d ] / G F, where id ∈ Diff(M) is the identity transformation. We see that this curve energy only depends on the Eulerian velocity field and the equivalence class q(1). Minimizers of EQ must necessarily satisfy:
{ t m t + £ u t [ m t ] = 0 , m t = P [ u t ] m 1 , w = D F ( q ( 1 ) ) ( w q ( 1 ) ) , w X ( M )

A geometric derivation of this formula can be found in [13] (Lemma 2.8, (2.12), and (2.13)). Again, the solution only depends on the Eulerian velocity and q(1). For this reason, we see that the GF symmetry of F provides a second reduction in the data needed to solve our original problem.

4.3. Orthogonality

In addition to reducing the amount of data that we must keep track of, there is an additional consequence to the GF-symmetry of F. In particular, there is a potentially massive constraint satisfied by the Eulerian velocity u.

To describe this, we must introduce an isotropy algebra. Given q ( t ) = [ Φ 0 , t u ] / G F, we can define the (time-dependent) isotropy algebra:

g q ( t ) = { w X ( M ) | w q ( t ) = 0 }

This is nothing, but the “Lie algebra” associated with the isotropy group Gq(t) = {ψ ∈ Diff(M) | ψ·q(t) = q(t)}. The use of quotes here is deliberate. If we let X ( M )denote an RKHS obtained from the space of vector-fields, then some of these are permitted to be non-smooth, which means that the standard Lie bracket of vector-fields does not close.

It turns out that the velocity field ut that minimizes E (or EQ) is orthogonal to g q ( t )with respect to the chosen inner-product. Intuitively, this is quite sensible, because velocities that do not change q(t) do not alter the data and simply waste control effort. Equivalently, said from the perspective Lagrange-multipliers, we know that the Lagrange-multipliers used to enforce optimality should point in a direction orthogonal to those which leave the cost functional unaltered. These variations on the same statement are formalized below.

Proposition 1. Let u satisfy Equation (8) or Equation (9). Then, m = P[u] annihilates g q ( t ).

Proof. Let u be the solution to Equation (9). We will first prove that u1 is orthogonal to g q ( t ). Let w 1 g q ( t ). We observe:

P [ u 1 ] , w 1 = by ( 9 ) d d ϵ | ϵ = 0 F Q ( Φ ϵ w ( 1 ) q ( 1 ) )

However, w1 leaves q(1) fixed, so Φ ϵ w ( 1 ) · q ( 1 ) = 0. Therefore, 〈P [u1], w1〉 = 0. Let w t = [ Φ t , 1 u ] w 1. In coordinates, this means:

w t i ( x ) = j | [ Φ t , 1 u ] 1 ( x ) [ Φ t , 1 u ] i w t j ( 1 , [ Φ t , 1 u ] 1 ( x ) )

One can directly verify that w t g q ( t )for all t ∈ [0, 1]. Denoting mt = P [ut], as in Equation (9), we find:

d d t P [ u t ] , w t = d d t m t , w t = t m t , w t + m t , t w t = £ u t [ m t ] , w t + m t , £ u t [ w t ] = 0

The last equality follows from Equation (5). Thus, 〈P [ut], wt〉 is constant. We have already verified that at t = 1, this inner-product is zero, and therefore, 〈P [ut], wt〉 = 0 for all time. That w1 is an arbitrary element of g q ( 1 )makes wt an arbitrary element of g q ( 1 )at each time. Thus, ut is orthogonal to g q ( 1 )for all time. □

At this point, we should return to our example to illustrate this idea.

Example 6. Again, consider the setup of Example 2. In this case, q(t) = (q1(t), q2(t)) ∈ M × M − ΔM×M. The space g q ( 1 )is the space of vector-fields that vanish at q1(t) and q2(t). Therefore, ut is orthogonal to q(t) if and only if mt = P [ut] satisfies:

m t , v = p 1 v ( z 1 ( t ) ) + p 2 v ( z 2 ( t ) )
for some covectors p1, p 2 X ( M ) and for any v X ( M ). In other words:
m t = p 1 ( t ) δ q 1 ( t ) ( ) + p 2 δ q 2 ( t ) ( )
where δx(·) denotes the Dirac delta functional centred at x.

This orthogonality constraint allows one to reduce the evolution equation on X ( M )to an evolution equation on Q, which might be finite dimensional if GF is large enough. In particular, there is a horizontal lift, h : T Q X ( M ), uniquely defined by the conditions h ( q , q ˙ ) q = 0, and h ( q , q ˙ ) g qwith respect to the chosen inner-product on vector-fields.

Example 7. Consider the setup of Example 2 with M = ℝn. Then, Q = n × n Δ n × n. Let K : ℝn×nn×n be the reproducing kernel of P. Then, h : T Q X ( n )is given by:

h ( q , q ˙ ) ( x ) = K ( x q 1 ) p 1 + K ( x q 2 ) p 2
where p1, p2 ∈ℝn are such that p 1 + K ( q 1 q 2 ) p 2 = q ˙ 1and K ( q 2 q 1 ) p 1 + p 2 = q ˙ 2.

One can immediately observe that h is injective and linear in q ˙. In other words, h ( q , ) : T q Q X ( M )is an injective linear map for fixed qQ. Because the optimal ut is orthogonal to g q ( t ), we may invert h(q(t), ·) on ut. In particular, we may often write the equation of motion on TQ, rather than on X ( M ). This is a massive reduction if Q is finite dimensional. In particular, the inner-product structure on X ( M )induces a Riemannian metric on Q given by:

g q ( v 1 , v 2 ) = P [ h ( q , v 1 ) ] , h ( q , v 2 ) .

The equations of motion in Equations (8) and (9) map to the geodesic equations on Q.

Proposition 2. Let u extremize E or EQ. Then, there exists a unique trajectory q(t) ∈ Q, such that u = h ( d q ( t ) d t ). Moreover, q(t) is a geodesic with respect to the metric g.

Proof. Let u minimize E. Thus, u satisfies Equation (9). By the previous proposition, ut is orthogonal to g q ( t ). As h ( q ( t ) , ) : T q ( t ) Q X ( M )is injective on g q ( t ) , there exists a unique q ˙ ( t ), such that h ( q ( t ) , q ˙ ( t ) ) = u tNote that E can be written as:

E = | | h ( q ( t ) , q ˙ ( t ) ) | | v d t + F ( q ( 1 ) ) = g q ( q ˙ , q ˙ ) 1 / 2 d t + F ( q ( 1 ) )

Thus, minimizers of E correspond to geodesics in Q with respect to the metric g. □

If we let H: TQ → ℝ be the Hamiltonian induced by the metric on Q, we obtain the most data-efficient form or Equations (8) and (9). Minimizers of E (or EQ) are:

{ ( q , p ) ( t ) T * Q satisfies Hamilton’s equation p ( 1 ) = D F Q ( q ) q ( 0 ) = [ e ] / G F

We see that this is a boundary value problem posed entirely on Q. If Q is finite dimensional, this is a massive reduction in terms of data requirements.

Example 8. Consider the setup of Example 2 with M = ℝn. The metric on Q = M × M Δ M 2is most easily expressed on the cotangent bundle TQ. If K is the matrix valued kernel of P, the metric on TQ takes the form:

g q * ( p , p ) = i , j = 1 2 p i T K ( q i q j ) p j

4.4. Descending Group Action

A related approach to defining distances on a space of objects to be registered consists of defining an object space Oupon which Diff(M) acts transitively (this means that for any o1, o2 O, there exists a φ ∈ Diff(M) such that φ·o1 = o2) with distance:

d O ( o 1 , o 2 ) = inf φ Diff ( M ) { d ( i d , φ ) | φ o 1 = o 2 }

Here, the distance on Ois defined directly from the distance in the group acting on the objects; see for example [1,5]. With this approach, the Riemannian metric descends from Diff(M) to a Riemannian metric on O, and geodesics on Olift by horizontality to geodesics on Diff(M). The quotient spaces Q obtained by reduction by symmetry and their geometric structure correspond to the object spaces and geometries defined with this approach. Intuitively, reduction by symmetry can be considered a removal of redundant information to obtain compact representations, while letting the metric descend to the object space Oconstitutes an approach to defining a geometric structure on an already known space of objects. The resulting solutions are equivalent to the ones presented in this article, because O≅ Diff(M)/Go, where Go = {ψ ∈ Diff(M) | ψ(o) = o} for some fixed reference object o O.

5. Examples

Here, we present a number of concrete examples of how reduction by symmetry can reduce the infinite dimensional registration problem over Diff(M) to lower, in some cases finite, dimensional problems. In all examples, the symmetry of the dissimilarity measure with respect to a subgroup of Diff(M) gives a reduced space by quotienting out the isotropy subgroups.

5.1. Landmark Matching

The space Q used in the examples in Section 4 constitutes a special case of the landmark matching problem, where sets of landmarks Q = {(x1,…, xN)| xiM, xixjij} are placed into spatial correspondence trough the left action φ·(x1,…, xN) = (φ(x1),…, φ(xN)} of Diff(M) by minimizing the dissimilarity measure F ( φ ) = i = 1 N φ ( x i ) x i 2. The landmark space Q arises as a quotient of Diff(M) from the isotropy group GF, as in Example 2.

Reduction from Diff(M) to Q in the landmark case has been used in a series of papers starting with [27]. Hamilton’s equations (Equation (10)) take the form:

q ˙ i = j = 1 N K ( q i q j ) p j , p ˙ i = j = 1 N ( D K ( q i q j ) p j ) T p i
on TQ, where DK denotes the spatial derivative of the reproducing kernel K. Generalizing the situation in Example 6, the momentum field is a finite sum of Dirac measures j = 1 N p j δ q j, and the Eulerian velocity field is the corresponding finite linear combination of the kernel evaluated at q i : u ( ) = j = 1 N K ( q j ) p j. Registration of landmarks is often in practice done by optimizing over the initial value of the momentum p in the ODE to minimize E, a strategy called shooting [28]. Using symmetry, the optimization problem is thus reduced from an infinite dimensional time-dependent problem to an N dim(M) dimensional optimization problem involving integration of a 2N dim(M) dimensional ODE on TQ.

5.2. Curve and Surface Matching

The space of smooth non-intersecting closed parametrized curves in ℝn is also known as the space of embeddings, denoted Emb(S1, ℝn). The parametrization can be removed by considering the right action of Diff(S1) on Emb(S2, ℝn) given by:

c Emb ( S 1 , n ) ψ Diff ( S 1 ) c ψ Emb ( S 1 , n )

Then, the quotient space Gr(S1, ℝn):=Emb(S1, ℝn)/Diff(S1) is the space of unparameterized curves. The space Gr(S1, ℝn) is a special case of a non-linear Grassmannian [29], and it has a manifold structure under certain conditions on the space of embeddings and the space of diffeomorphisms [30].

When the parametrization is not removed, embedded curves and surfaces can be matched with the current dissimilarity measure [31,32]. If M is a volume manifold, then the objects are considered elements of the dual space of Ωk(M), the space of differential k-forms on M. In the surface case, a bounded submanifold SM can be seen as an element of [Ωk(M)] by its evaluation on a k-form, w ∈ Ω2(M), given by:

S ( w ) : = S w | S
where w|S ∈ Ωk(S) is the restriction of w to S. The dual space (Ωk(M)) can be equipped with a norm that enables surfaces to be quantitatively compared as elements of the vector-space (Ωk(M)). Note that the evaluation Equation (11) does not depend on the parametrization of S, as it is written in a coordinate-free form. Coordinate-based formulations of Equation (11) are available in [31,32]. This technique is computationally much more tractable than using the Hausdorff distance, which requires pairwise comparisons between all points between two surfaces.

The isotropy groups for curves and surfaces generalize the isotropy groups of landmarks by consisting of warps that keep the objects fixed, i.e.,

G F { ψ Diff ( M ) | ψ ( S ) = S } .

The momentum field will be supported on the transported curves/surfaces φ(t).S for optimal paths for E in Diff(M).

5.3. Image Matching

Images can be registered using either the L2-difference defined in Example 1 or with other dissimilarity measures, such as mutual information or the correlation ratio [33,34]. The similarity will be invariant to any infinitesimal deformation orthogonal to the gradient of dissimilarity measure. In the L2 case, this is equivalent to any infinitesimal deformation orthogonal to the level sets of the moving image [35]. The momentum field thus has the form p(t) = α(t)∇φ(t).I0 for a smooth function α(t) on M (see Figure 6), and the registration problem can be reduced to a search over the scalar field α(t) instead of vector field p(t).

Minimizers for E follow the PDE [5]:

u ˙ t M K ( y ) α ( y ) m t ( y ) d y , m ˙ t = m t T u t , α ˙ = ( α u t )
with mt representing the deformed image at time t.

In particular, the isotropy group of a source image f0C (M) is the subgroup of diffeomorphisms, which preserve the level sets of f0. The quotient space Diff(M)=Iso(f0) can be identified with the orbit of f0 by diffeomorphisms, i.e., Diff(M)=Iso(f0) ≅ Orb(f0):= {φ*f0 | φ ∈ Diff(M)}. This orbit is difficult to identify with a more concrete object, in contrast to, e.g., the case of landmark matching. However, it can be characterized by various properties. For example, for a function f ∈ Orb(f0) and any cf0 (M) ⊂ ℝ, the level sets f−1(c) and f0−1(c) have the same topology.

5.4. Jet Matching

In [14,36], an extension of the landmark case has been developed where higher-order spatial information is advected with the landmarks. The spaces of jet-particles arise as extensions of the reduced landmark space Q by quotienting out smaller isotropy subgroups known as jet-groups. A thorough account of jet-groups, including Lie group and algebra structures, can be found in [37] (Chapter 4). We provide a brief introduction here. Let G(0) be the isotropy subgroup for a single landmark:

G ( 0 ) : = { ψ G | ψ ( q ) = q }

Let now k be a positive integer. For any k-differentiable map f from a neighbourhood of q, the k-jet of f is denoted J q ( k ) ( f ). In coordinates, J q ( k ) ( f )consists of the coefficients of the k-th order Taylor expansions of f aboutx. The higher-order isotropy subgroups are then given by:

G ( k ) : = { ψ G ( 0 ) | J q ( k ) ψ = J q ( k ) i d }

That is, the elements of G(k) fix the Taylor expansion of the deformation φ up to order k. The definition naturally extends to finite numbers of landmarks, and the quotients Q(k) = G/G(k) can be identified as the sets:

Q ( 0 ) = { ( q 1 , , q N ) | q i M } Q ( 1 ) = { ( ( q i ( 0 ) , q i ( 1 ) ) , ) | ( q i ( 0 ) , q i ( 1 ) , ) M × GL ( d ) } Q ( 2 ) = { ( ( q i ( 0 ) , q i ( 1 ) , q i ( 2 ) ) , ) | q i ( 0 ) , q i ( 1 ) , q i ( 2 ) , ) M × GL ( d ) × S 2 1 }
with S 2 1being the space of rank (1, 2) tensors symmetric in the lower indices. Intuitively, the space Q(0) is the regular landmark space with information about the position of the points; the one-jet space Q(1) carries for each jet information about the position and the Jacobian matrix of the warp at the jet position; and the two-jet space Q(2) carries in addition the Hessian matrix of the warp at the jet position. The momentum for Q(0) in coordinates consists of N vectors representing the local displacement of the points. With the one-jet space Q(0), the momentum in addition contains d × d matrices that can be interpreted as locally linear deformations at the jet positions [36]. In combination with the displacement, the one-jet momenta can thus be regarded as locally affine transformations. The momentum fields for Q(2) add symmetric tensors encoding local second order deformation. The local effect of the jet-particles is sketched in Figure 7.

When the dissimilarity measure F is dependent not just on positions, but also on higher-order information around the points, reduction by symmetry implies that optimal solutions for E will be parametrized by k-jets in the same way as Q(0) parametrizesoptimal paths for E in the landmark case. The higher-order jets can thus be used for landmark matching when the dissimilarity measure is dependent on the local geometry around the landmarks. For example, matching of the first order structure, such as image gradients, leads to first-order jets, and matching of local curvature leads to second-order jets.

5.5. Discrete Image Matching

The image matching problem can be discretized by evaluating the L2-difference at a finite number of points. In practice, this always happens when the integral M | I 0 φ 1 ( x ) I 1 ( x ) | 2 d xis evaluated at finitely many pixels of the image. In [36,38], it is shown how this reduces the image matching PDE (12) to a finite dimensional system on Q when the integral is approximated by pointwise evaluation at a grid Λh:

F ( 0 ) ( φ ) x Λ h h d | I 0 ( φ 1 ( x ) I 1 ( x ) | 2
where h > 0 denotes the grid spacing. F(0) approximates F to order O(hd), d = dim(M). The reduced space Q encodes the position of the points φ1(x), x ∈ Λh, and the lifted Eulerian momentum field is a finite sum of point measures p = x Λ h a x δ φ 1 ( x ). For each grid point, the momentum encodes the local displacement of the point.

In [38], a discretization scheme with higher-order accuracy is in addition introduced with an O(hd+2) approximation F(2) of F. The increased accuracy results in the entire energy E being approximated to order O(hd+2). The solution space in these cases becomes the jet-space Q(2). For a given order of approximation, a corresponding reduction in the number of required discretization points is obtained. The reduction is countered by the increased information encoded in each two-jet. The momentum field thus encodes both local displacement, local linear deformation and second order deformation; see Figure 7. The discrete solutions will converge to solutions of the non-discretized problem as h → 0.

5.6. DWI/DTI Matching

Image matching is invariant with respect to variations parallel to the level lines of the images. With diffusion weighted images (DWI) and the variety of models for the diffusion information (e.g., diffusion tensor imaging (DTI) [39], Gaussian mixture fields [40]), first or higher-order information can be reintroduced into the matching problem. In essence, by letting the dissimilarity measure depend on the diffusion information, the isotropy subgroup of the matching problem becomes smaller.

The exact form of the of DWI matching problem depends on the diffusion model and how Diff(M) acts on the diffusion image. In [41], the diffusion is represented by the principal direction of the diffusion tensor, and the data objects to be matched are thus vector fields. The action by elements of Diff(M) is defined by:

φ I ( x ) = I φ 1 D φ I φ 1 D φ I φ 1 .

The action rotates the diffusion vector by the Jacobian of the warp, keeping its length fixed. Similar models can be applied to DTI with the preservation of principle direction scheme (PPD, [42,43]) and to GMF-based models [44]. The dependency on the Jacobian matrix implies that a reduced model must carry first order information in a similar fashion to the one-jet space Q(1); however, any irrotational part of the Jacobian can be removed by symmetry. The full effect of this has yet to be explored.

As in the case of image matching, the quotient can be identified with the orbit of the source data under diffeomorphisms.

5.7. Fluid Mechanics

Incidentally, the equation of motion:

t m t + L u [ m t ] = 0 , u t = K * m t
is an eccentric way of writing Euler’s equation for an inviscid incompressible fluid if we assume u t X ( n )is initially in the space of divergence free vector-fields and K∗ is the Riemannian flat map (which implies that mt and ut can be identified as functions on ℝn) [45]. This fact was exploited in [46] to create a sequence of regularized models to Euler’s equations by considering a sequence of kernels, such that the operator K∗ (viewed as a map to one-form densities) converges to a surjection onto the annihilator gradient vector-fields (this is written as a projection onto divergence free vector-fields in [46]). Moreover, if one replaces Diff(M) by the subgroup of volume preserving diffeomorphisms Diffvol(M), then (formally) one can produce incompressible particle methods using the same reduction arguments presented here. In fact, jet-particles were independently discovered in this context as a means of simulating fluids in [47]. It is notable that [47] is a mechanics paper, and the particle methods that were produced were approached from the perspective of reduction by symmetry without any knowledge of the related work being done in image registration.

In [48], one of the kernel parameters in [46], which controls the compressibility of the u, was taken to the incompressible limit. This allowed a realization of the particle methods described in [47]. The constructions of [48] are the same as presented in this article, but with Diff(M) replaced by the group of volume-preserving diffeomorphisms of ℝd. Velocity fields induced by first order jet-particles are visualized in Figure 8.

6. Discussion and Conclusions

The information available for solving the registration problem is in practice not sufficient for uniquely encoding the deformation between the objects to be registered. Symmetry arises in both particle relabelling symmetry that gives the Eulerian formulation of the equations of motion and in symmetry groups for specific dissimilarity measures.

For landmark matching, reduction by symmetry reduces the infinite dimensional registration problem to a finite dimensional problem on the reduced landmark space Q. For matching curves and surfaces, symmetry implies that the momentum stays concentrated at the curves and surfaces allowing a reduction by the isotropy groups of warps that leave the objects fixed. In image matching, symmetry allows reduction by the group of warps that do not change the level sets of the image. Jet-particles arise from smaller isotropy subgroups and, hence, larger reduced spaces Q(1) and Q(2) that encode locally affine and second order information.

Reduction by symmetry allows these cases to be handled in one theoretical framework. We have surveyed the mathematical construction behind the reduction approach and its relation to the above-mentioned examples. As data complexity rises both in terms of resolution and structure, symmetry will continue to be an important tool for removing redundant information and achieving compact data representations.

Acknowledgments

Henry O. Jacobs would like to thank Darryl Holm for providing him with a bridge from geometric mechanics to the field of image registration algorithms. Henry O. Jacobs is supported by the European Research Council Advanced Grant 267382 FCCA. Stefan Sommer is supported by the Danish Council for Independent Research with the project “Image-Based Quantification of Anatomical Change”.

Author Contributions

Henry O. Jacobs and Stefan Sommer wrote the paper.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Younes, L. Shapes and Diffeomorphisms; Springer-Berlin: Heidelberg, Germany, 2010. [Google Scholar]
  2. Brox, T.; Bruhn, A.; Papenberg, N.; Weickert, J. High Accuracy Optical Flow Estimation Based on a Theory for Warping. In Computer Vision-ECCV 2004; aPajdla, T., Matas, J., Eds.; Springer-Berlin: Heidelberg, Germany, 2004; pp. 25–36. [Google Scholar]
  3. Dawn, S.; Saxena, V.; Sharma, B. Remote Sensing Image Registration Techniques: A Survey. In Image and Signal Processing; Elmoataz, A., Lezoray, O., Nouboud, F., Mammass, D., Meunier,, J., Eds.; Springer-Berlin: Heidelberg, Germany, 2010; pp. 103–112. [Google Scholar]
  4. Sotiras, A.; Davatzikos, C.; Paragios, N. Deformable Medical Image Registration: A Survey. IEEE Trans. Med. Imag 2013, 32, 1153–1190. [Google Scholar]
  5. Younes, L.; Arrate, F.; Miller, M.I. Evolutions equations in computational anatomy. NeuroImage 2009, 45, S40–S50. [Google Scholar]
  6. Boyes, R.G.; Rueckert, D.; Aljabar, P.; Whitwell, J.; Schott, J.M.; Hill, D.L.G.; Fox, N.C. Cerebral atrophy measurements using Jacobian integration: Comparison with the boundary shift integral. NeuroImage 2006, 32, 159–169. [Google Scholar]
  7. Gorbunova, V.; Jacobs, S.S.A.M.; Lo, P.; Dirksen, A.; Nielsen, M.; Bab-Hadiashar, A.; de Bruijne, M. Early detection of emphysema progression. In Medical Image Computing and Computer-Assisted Intervention–MICCAI 2010; Springer-Berlin: Heidelberg, Germany, 2010; Volume 13, pp. 193–200. [Google Scholar]
  8. Joshi, S.C.; Miller, M.I.; Grenander, U. On the Geometry and Shape of Brain Sub-Manifolds. Int. J. Patt. Recogn. Artif. Intell 1997, 11, 1317–1343. [Google Scholar]
  9. Dupuis, P.; Grenander, U.; Miller, M.I. Variational Problems on Flows of Diffeomorphisms for Image Matching. Q. Appl. Math. 1998, 56, 587–600. [Google Scholar]
  10. Trouvé, A. An Infinite Dimensional Group Approach for Physics Based Models in Pattern Recognition, 1995. Available online: http://cis.jhu.edu/publications/papers_in_database/alain/trouve1995.pdf accessed on 28 April 2015.
  11. Christensen, G.; Rabbitt, R.; Miller, M. Deformable templates using large deformation kinematics. IEEE Trans. Image Process 1996, 5, 1435–1447. [Google Scholar]
  12. Holm, D.D.; Ratnanather, J.T.; Trouvé, A.; Younes, L. Soliton Dynamics in Computational Anatomy. NeuroImage 2004, 23, S170–S178. [Google Scholar]
  13. Bruveris, M.; Gay-Balmaz, F.; Holm, D.D.; Ratiu, T.S. The momentum map representation of images. J. Nonlinear Sci 2011, 21, 115–150. [Google Scholar]
  14. Jacobs, H. Symmetries in LDDMM with higher order momentum distributions, 2013, arXiv, p. 1306.3309. arXiv.org e-Print archive. Available online: http://arxiv.org/abs/1306.3309 accessed on 28 April 2015.
  15. Pennec, X.; Stefanescu, R.; Arsigny, V.; Fillard, P.; Ayache, N. Riemannian Elasticity: A Statistical Regularization Framework for Non-linear Registration. In MICCAI 2005; Springer-Berlin: Heidelberg, Germany, 2005; pp. 943–950. [Google Scholar]
  16. Derfoul, R.; Le Guyader, C. A. Relaxed Problem of Registration Based on the Saint Venant–Kirchhoff Material Stored Energy for the Mapping of Mouse Brain Gene Expression Data to a Neuroanatomical Mouse Atlas. SIAM J. Imag. Sci. 2014, 7, 2175–2195. [Google Scholar]
  17. Rueckert, D.; Sonoda, L.I.; Hayes, C.; Hill, D.L.; Leach, M.O.; Hawkes, D.J. Nonrigid registration using free-form deformations: Application to breast MR images. IEEE Trans. Med. Imag. 1999, 18, 712–721. [Google Scholar]
  18. Horn, B.K.; Schunck, B.G. Determining optical flow. Artif. Intell. 1981, 17, 185–203. [Google Scholar]
  19. Trouvé, A.; Younes, L. Local geometry of deformable templates. SIAM J. Math. Anal. 2005, 37, 17–59. [Google Scholar]
  20. Abraham, R.; Marsden, J.E. Foundations of Mechanics; Benjamin/Cummings Publishing Co. Inc.: Reading, MA, USA, 1978. [Google Scholar]
  21. Kobayashi, S.; Nomizu, K. Foundations of Differential Geometry, Volume I; Wiley Classics Library, John Wiley & Sons, Inc: New York, NY, USA, 1996. [Google Scholar]
  22. Taylor, M. Pseudo Differential Operators; Springer-Verlag: Berlin, Germany; New York, NY, USA, 1974. [Google Scholar]
  23. Holm, D.D.; Marsden, J.E. Momentum maps and measure-valued solutions (peakons, filaments, and sheets) for the EPDiff equation. In The Breadth of Symplectic and Poisson Geometry; Birkhäuser Boston: Boston, MA, USA, 2005; Volume 232, pp. 203–235. [Google Scholar]
  24. Micheli, M.; Glaunès, J.A. Matrix-valued Kernels for Shape Deformation Analysis. Geometry Imag. Comput 2014, 1, 57–39. [Google Scholar]
  25. Dupuis, P.; Grenander, U.; Miller, M.I. Variational problems on flows of diffeomorphisms for image matching. Quart. Appl. Math. 1998, 56, 587–600. [Google Scholar]
  26. Beg, M.F.; Miller, M.I.; Trouvé, A.; Younes, L. Computing Large Deformation Metric Mappings via Geodesic Flows of Diffeomorphisms. IJCV 2005, 61, 139–157. [Google Scholar]
  27. Joshi, S.; Miller, M. Landmark matching via large deformation diffeomorphisms. IEEE Trans. Image Process 2000, 9, 1357–1370. [Google Scholar]
  28. Vaillant, M.; Miller, M.; Younes, L.; Trouvé, A. Statistics on diffeomorphisms via tangent space representations. NeuroImage 2004, 23, S161–S169. [Google Scholar]
  29. Gay-Balmaz, F.; Vizman, C. Dual pairs in fluid dynamics. Ann. Glob. Anal. Geometry. 2012, 41, 1–24. [Google Scholar]
  30. Bauer, M.; Bruveris, M.; Michor, P.W. Overview of the Geometries of Shape Spaces and Diffeomorphism Groups. J. Math. Imag. Vis. 2014, 50, 60–97. [Google Scholar]
  31. Vaillant, M.; Glaunès, J. Surface Matching via Currents. In Information Processing in Medical Imaging; Springer-Berlin: Heidelberg, Germany, 2010; pp. 381–392. [Google Scholar]
  32. Glaunès, J. Transport Par Difféomorphismes de Points, de Mesures et de Courants Pour La Comparaison de Formes et L’anatomie Numérique. Ph.D. Thesis, Université Paris 13, Villetaneuse, France, 2005. [Google Scholar]
  33. Wells, W.M., III; Viola,, P.; Kikinis,, R. Multi-Modal Volume Registration by Maximization of Mutual Information. Med. Image Anal 1996, 1, 35–51. [Google Scholar]
  34. Roche, A.; Malandain, G.; Pennec, X.; Ayache, N. The Correlation Ratio as a New Similarity Measure for Multimodal Image Registration. In Proceedings of the First International Conference on Medical Image Computing and Computer-Assisted Intervention. In MICCAI 1998; Springer-Berlin: Heidelberg, Germany, 1998; pp. 1115–1124, ACM ID: 709612. [Google Scholar]
  35. Miller, M.I.; Trouvé, A.; Younes, L. Geodesic Shooting for Computational Anatomy. J. Math. Imag. Vis. 2006, 24, 209–228. [Google Scholar]
  36. Sommer, S.; Nielsen, M.; Darkner, S.; Pennec, X. Higher-Order Momentum Distributions and Locally Affine LDDMM Registration. SIAM J. Imag. Sci. 2013, 6, 341–367. [Google Scholar]
  37. Kolář, I.; Michor, P.W.; Slovák, J. Natural Operations in Differential Geometry; Springer Verlag: Berlin, Germany, 1999. [Google Scholar]
  38. Jacobs, H.O.; Sommer, S. Higher-order Spatial Accuracy in Diffeomorphic Image Registration. J. Geom. Imaging Comput, 2015. in press. Available online: http://arxiv.org/abs/1412.7504 accessed on 28 April 2015. [Google Scholar]
  39. Basser, P.J.; Mattiello, J.; LeBihan, D. Estimation of the effective self-diffusion tensor from the NMR spin echo. J. Magn. Resonance Ser. B 1994, 103, 247–254. [Google Scholar]
  40. Jian, B.; Vemuri, B.C.A. Unified Computational Framework for Deconvolution to Reconstruct Multiple Fibers From Diffusion Weighted MRI. IEEE Trans. Med. Imag 2007, 26, 1464–1471. [Google Scholar]
  41. Cao, Y.; Miller, M.I.; Winslow, R.L.; Younes, L. Large deformation diffeomorphic metric mapping of vector fields. IEEE Trans. Med. Imag. 2005, 24, 1216–1230. [Google Scholar]
  42. Alexander, D.C.; Gee, J.C.; Bajcsy, R. Strategies for Data Reorientation during Non-rigid Warps of Diffusion Tensor Images. In Medical Image Computing and Computer-Assisted Intervention—MICCAI’99; Taylor, C., Colchester, A., Eds.; Number 1679 in Lecture Notes in Computer Science; Springer-Verlag: Berlin, Germany, 1999; pp. 463–472. [Google Scholar]
  43. Alexander, D.C.; Pierpaoli, C.; Basser, P.J.; Gee, J.C. Spatial transformations of diffusion tensor magnetic resonance images. IEEE Trans. Med. Imag. 2001, 20, 1131–1139. [Google Scholar]
  44. Cheng, G.; Vemuri, B.C.; Carney, P.R.; Mareci, T.H. Non-rigid Registration of High Angular Resolution Diffusion Images Represented by Gaussian Mixture Fields. In Medical Image Computing and Computer-Assisted Intervention—MICCAI 2009; Springer-Verlag: Berlin, Germany, 2009; pp. 190–197. [Google Scholar]
  45. Arnold, V.I. Sur la géométrie différentielle des groupes de Lie de dimension infinie et ses applications à l’hydrodynamique des fluides parfaits. Ann. l’Instit. Fourier 1966, 16, 316–361. [Google Scholar]
  46. Mumford, D.; Michor, P.W. On Euler’s equation and “EPDiff”. J. Geometr. Mech 2013, 5, 319–344. [Google Scholar]
  47. Jacobs, H.O.; Ratiu, T.S.; Desbrun, M. On the coupling between an ideal fluid and immersed particles. Phys. D Nonlinear Phenome 2013, 265, 40–56. [Google Scholar]
  48. Cotter, C.J.; Holm, D.D.; Jacobs, H.O.; Meier, D.M. A jetlet hierarchy for ideal fluid dynamics. J. Phys. A Math. Theor. 2014, 47, 352001. [Google Scholar]
Figure 1. A diagram relating a top in the reference configuration to its current configuration via a rotation matrix R ∈ SO(3).
Figure 1. A diagram relating a top in the reference configuration to its current configuration via a rotation matrix R ∈ SO(3).
Symmetry 07 00599f1 1024
Figure 2. A registration of two discs of different sizes (a,b) with one example of a warp that brings (b) into correspondence with (a) visualized by its effect on an initially regular grid (c). Using symmetry, the dimensionality of the registration problem can be reduced from infinite to finite. In this case, six parameters of a one-jet particle (see Section 5.4) in the centre of the moving image encode the entire deformation. The six parameters can roughly be described as a position in ℝ2, a scaling, a stretch, a shear and a rotation. (a) Fixed image; (b) moving image; (c) warp.
Figure 2. A registration of two discs of different sizes (a,b) with one example of a warp that brings (b) into correspondence with (a) visualized by its effect on an initially regular grid (c). Using symmetry, the dimensionality of the registration problem can be reduced from infinite to finite. In this case, six parameters of a one-jet particle (see Section 5.4) in the centre of the moving image encode the entire deformation. The six parameters can roughly be described as a position in ℝ2, a scaling, a stretch, a shear and a rotation. (a) Fixed image; (b) moving image; (c) warp.
Symmetry 07 00599f2 1024
Figure 3. A warp φ ∈ Diff(M) acts on an image I: M → ℝ by composition with the inverse warp, φ·I = I ○φ−1. Given two images I0, I1: M → ℝ, image registration involves finding a warp φ, such that φ·I0 is close to I1 as measured by a dissimilarity measure F (φ·I0, I1).
Figure 3. A warp φ ∈ Diff(M) acts on an image I: M → ℝ by composition with the inverse warp, φ·I = I ○φ−1. Given two images I0, I1: M → ℝ, image registration involves finding a warp φ, such that φ·I0 is close to I1 as measured by a dissimilarity measure F (φ·I0, I1).
Symmetry 07 00599f3 1024
Figure 4. The tangent bundle T M of the manifold M consists of pairs (x, v) of points xM and tangent vectors vTxM. It is a fibre bundle over M with fibres TxM for each x ∈M.
Figure 4. The tangent bundle T M of the manifold M consists of pairs (x, v) of points xM and tangent vectors vTxM. It is a fibre bundle over M with fibres TxM for each x ∈M.
Symmetry 07 00599f4 1024
Figure 5. Examples of elements of the isotropy subgroup GF in Example 2 visualized by their effect on an initially square grid. The isotropy subgroup leaves the dissimilarity measure F invariant by not moving the points x1 and x2.
Figure 5. Examples of elements of the isotropy subgroup GF in Example 2 visualized by their effect on an initially square grid. The isotropy subgroup leaves the dissimilarity measure F invariant by not moving the points x1 and x2.
Symmetry 07 00599f5 1024
Figure 6. In image matching, the gradient of the L2-difference will be orthogonal to level lines of the image, and symmetry implies that the momentum field will be orthogonal to the level lines, so that p(t) = α(tφ(t):I0 for a time-dependent scalar field α.
Figure 6. In image matching, the gradient of the L2-difference will be orthogonal to level lines of the image, and symmetry implies that the momentum field will be orthogonal to the level lines, so that p(t) = α(tφ(t):I0 for a time-dependent scalar field α.
Symmetry 07 00599f6 1024
Figure 7. With discrete image matching, the image is sampled at a regular grid Λh, h > 0, and the image matching PDE (12) is reduced to an ODE on a finite dimensional reduced space Q. With the approximation F(0) (13), the momentum field will encode local displacement, as indicated by the horizontal arrows (top row). With a first order expansion, the solution space will be jet space Q(1), and locally affine motion is encoded around each grid point (middle row). The O(hd+2) approximation F(2) includes second order information, and the system reduces to the jet space Q(2) with second order motion encoded at each grid point (lower row).
Figure 7. With discrete image matching, the image is sampled at a regular grid Λh, h > 0, and the image matching PDE (12) is reduced to an ODE on a finite dimensional reduced space Q. With the approximation F(0) (13), the momentum field will encode local displacement, as indicated by the horizontal arrows (top row). With a first order expansion, the solution space will be jet space Q(1), and locally affine motion is encoded around each grid point (middle row). The O(hd+2) approximation F(2) includes second order information, and the system reduces to the jet space Q(2) with second order motion encoded at each grid point (lower row).
Symmetry 07 00599f7 1024
Figure 8. Velocity fields induced by first order incompressible jet-particles.
Figure 8. Velocity fields induced by first order incompressible jet-particles.
Symmetry 07 00599f8 1024
Table 1. Notation.
Table 1. Notation.
NotionNotation
manifoldM
tangent bundleTM
the group of diffeomorphismsDiff(M)
the space of vector-fields X ( M )
the push-forward of u u X(M) with respect to Φ ∈ Diff(M)Φ*u
the push-forward of m X(M)* with respect to Φ 2∈ Diff(M)Φ*m
The Lie derivative (w.r.t. u X(M)) L u
the dual space to X(M) X(M)*
Evaluation of m X(M)* on u X(M)m, u

Share and Cite

MDPI and ACS Style

Sommer, S.; Jacobs, H.O. Reduction by Lie Group Symmetries in Diffeomorphic Image Registration and Deformation Modelling. Symmetry 2015, 7, 599-624. https://doi.org/10.3390/sym7020599

AMA Style

Sommer S, Jacobs HO. Reduction by Lie Group Symmetries in Diffeomorphic Image Registration and Deformation Modelling. Symmetry. 2015; 7(2):599-624. https://doi.org/10.3390/sym7020599

Chicago/Turabian Style

Sommer, Stefan, and Henry O. Jacobs. 2015. "Reduction by Lie Group Symmetries in Diffeomorphic Image Registration and Deformation Modelling" Symmetry 7, no. 2: 599-624. https://doi.org/10.3390/sym7020599

APA Style

Sommer, S., & Jacobs, H. O. (2015). Reduction by Lie Group Symmetries in Diffeomorphic Image Registration and Deformation Modelling. Symmetry, 7(2), 599-624. https://doi.org/10.3390/sym7020599

Article Metrics

Back to TopTop